From b3de75da971382b79d8c5b585bb58ace365d48ee Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 23 Feb 2024 13:26:28 -0700 Subject: [PATCH] Cleaning up math, fixing tension bug, patching bond creation --- src/RHEO/compute_rheo_interface.cpp | 54 ++++++++++--------- src/RHEO/fix_rheo.cpp | 6 --- src/RHEO/fix_rheo_tension.cpp | 2 +- src/RHEO/fix_rheo_thermal.cpp | 80 +++++++++++++++++++++++------ 4 files changed, 90 insertions(+), 52 deletions(-) diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 001f15a472..3cb2fcf058 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -26,6 +26,7 @@ #include "force.h" #include "fix_rheo.h" #include "fix_rheo_pressure.h" +#include "math_extra.h" #include "memory.h" #include "modify.h" #include "neighbor.h" @@ -36,6 +37,7 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; +using namespace MathExtra; static constexpr double EPSILON = 1e-1; @@ -107,8 +109,8 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr) void ComputeRHEOInterface::compute_peratom() { - int i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match; - double xtmp, ytmp, ztmp, delx, dely, delz, rsq, w, dot; + int a, i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match; + double xtmp, ytmp, ztmp, rsq, w, dot, dx[3]; int inum, *ilist, *jlist, *numneigh, **firstneigh; int nlocal = atom->nlocal; @@ -153,15 +155,15 @@ void ComputeRHEOInterface::compute_peratom() j = jlist[jj]; j &= NEIGHMASK; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx * delx + dely * dely + delz * delz; + dx[0] = xtmp - x[j][0]; + dx[1] = ytmp - x[j][1]; + dx[2] = ztmp - x[j][2]; + rsq = lensq3(dx); if (rsq < cutsq) { jtype = type[j]; fluidj = !(status[j] & PHASECHECK); - w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq)); + w = compute_kernel->calc_w_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq)); status_match = 0; norm[i] += w; @@ -172,9 +174,9 @@ void ComputeRHEOInterface::compute_peratom() chi[i] += w; } else { if (!fluidi) { - dot = (-fp_store[j][0] + fp_store[i][0]) * delx; - dot += (-fp_store[j][1] + fp_store[i][1]) * dely; - dot += (-fp_store[j][2] + fp_store[i][2]) * delz; + dot = 0; + for (a = 0; a < 3; a++) + dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a]; rho[i] += w * (fix_pressure->calc_pressure(rho[j], jtype) - rho[j] * dot); normwf[i] += w; @@ -187,9 +189,9 @@ void ComputeRHEOInterface::compute_peratom() chi[j] += w; } else { if (!fluidj) { - dot = (-fp_store[i][0] + fp_store[j][0]) * delx; - dot += (-fp_store[i][1] + fp_store[j][1]) * dely; - dot += (-fp_store[i][2] + fp_store[j][2]) * delz; + dot = 0; + for (a = 0; a < 3; a++) + dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a]; rho[j] += w * (fix_pressure->calc_pressure(rho[i], itype) + rho[i] * dot); normwf[j] += w; @@ -225,7 +227,7 @@ void ComputeRHEOInterface::compute_peratom() int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - int i,j,k,m; + int i, j, k, m; m = 0; double *rho = atom->rho; @@ -267,7 +269,7 @@ void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf) int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf) { - int i,k,m,last; + int i, k, m, last; double *rho = atom->rho; m = 0; @@ -285,7 +287,7 @@ int ComputeRHEOInterface::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf) { - int i,k,j,m; + int i, k, j, m; double *rho = atom->rho; int *status = atom->status; m = 0; @@ -350,23 +352,19 @@ void ComputeRHEOInterface::store_forces() for (const auto &fix : fixlist) { for (int i = 0; i < atom->nlocal; i++) { minv = 1.0 / mass[type[i]]; - if (mask[i] & fix->groupbit) { - fp_store[i][0] = f[i][0] * minv; - fp_store[i][1] = f[i][1] * minv; - fp_store[i][2] = f[i][2] * minv; - } else { - fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv; - fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv; - fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv; - } + if (mask[i] & fix->groupbit) + for (int a = 0; a < 3; a++) + fp_store[i][a] = f[i][a] * minv; + else + for (int a = 0; a < 3; a++) + fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } } } else { for (int i = 0; i < atom->nlocal; i++) { minv = 1.0 / mass[type[i]]; - fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv; - fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv; - fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv; + for (int a = 0; a < 3; a++) + fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } } diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 3a999e12dd..beba940174 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -422,12 +422,6 @@ void FixRHEO::pre_force(int /*vflag*/) if (mask[i] & groupbit) status[i] &= OPTIONSMASK; - // Reinstate temporary options - for (int i = 0; i < nall; i++) - if (mask[i] & groupbit) - if (status[i] & STATUS_SOLID) - status[i] |= STATUS_NO_SHIFT; - // Calculate surfaces, update status if (surface_flag) { compute_surface->compute_peratom(); diff --git a/src/RHEO/fix_rheo_tension.cpp b/src/RHEO/fix_rheo_tension.cpp index 8d00c8b988..388b574365 100644 --- a/src/RHEO/fix_rheo_tension.cpp +++ b/src/RHEO/fix_rheo_tension.cpp @@ -467,7 +467,7 @@ void FixRHEOTension::pre_force(int vflag) if (wsame[i] < wmin) continue; - weight = MAX(1.0, wsame[i] * wmin_inv); + weight = MIN(1.0, wsame[i] * wmin_inv); //MAX -> MIN 2/14/24 itype = type[i]; if (norm[i] != 0) diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 635cf78c85..0640dd6827 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -423,6 +423,7 @@ void FixRHEOThermal::pre_force(int /*vflag*/) double *energy = atom->esph; double *temperature = atom->temperature; int *type = atom->type; + int *status = atom->status; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; @@ -444,6 +445,11 @@ void FixRHEOThermal::pre_force(int /*vflag*/) } } } + + // Add temporary options, wiped by preceding fix rheo preforce + for (int i = 0; i < nall; i++) + if (status[i] & STATUS_SOLID) + status[i] |= STATUS_NO_SHIFT; } /* ---------------------------------------------------------------------- */ @@ -481,39 +487,79 @@ void FixRHEOThermal::break_bonds() tagint **bond_atom = atom->bond_atom; int *num_bond = atom->num_bond; - int nlocal = atom->nlocal; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + // Rapidly delete all bonds for local atoms that melt (no shifting) for (int i = 0; i < nlocal; i++) { + if (!(status[i] & STATUS_MELTING)) continue; for (m = 0; m < num_bond[i]; m++) { j = atom->map(bond_atom[i][m]); - if (!(status[i] & STATUS_MELTING) && !(status[j] & STATUS_MELTING)) continue; + bond_type[i][m] = 0; if (n_histories > 0) for (auto &ihistory: histories) - dynamic_cast(ihistory)->delete_history(i, num_bond[i] - 1); + dynamic_cast(ihistory)->delete_history(i, m); - if (fix_update_special_bonds) fix_update_special_bonds->add_broken_bond(i, j); + if (fix_update_special_bonds) + fix_update_special_bonds->add_broken_bond(i, j); + } + num_bond[i] = 0; + } - // For non-melting neighbors, selectively delete bond if necessary - if (j >= nlocal || (status[j] & STATUS_MELTING)) continue; - for (n = 0; n < num_bond[j]; n++) { - if (bond_atom[j][n] == tag[i]) { - bond_type[j][n] = 0; - nmax = num_bond[j] - 1; - bond_type[j][n] = bond_type[j][nmax]; - bond_atom[j][n] = bond_atom[j][nmax]; - if (n_histories > 0) { + // Update bond list and break solid-melted bonds + for (n = 0; n < nbondlist; n++) { + + // skip bond if already broken + if (bondlist[n][2] <= 0) continue; + i = bondlist[n][0]; + j = bondlist[n][1]; + + if (!(status[i] & STATUS_MELTING) && !(status[j] & STATUS_MELTING)) continue; + + bondlist[n][2] = 0; + + // Delete bonds for non-melted local atoms (shifting) + if (i < nlocal) { + for (m = 0; m < num_bond[i]; m++) { + if (bond_atom[i][m] == tag[j]) { + nmax = num_bond[i] - 1; + bond_type[i][m] = bond_type[i][nmax]; + bond_atom[i][m] = bond_atom[i][nmax]; + if (n_histories > 0) for (auto &ihistory: histories) { - dynamic_cast(ihistory)->shift_history(j, n, nmax); - dynamic_cast(ihistory)->delete_history(j, nmax); + auto fix_bond_history = dynamic_cast (ihistory); + fix_bond_history->shift_history(i, m, nmax); + fix_bond_history->delete_history(i, nmax); } - } + bond_type[i][nmax] = 0; + num_bond[i]--; + break; + } + } + } + + if (j < nlocal) { + for (m = 0; m < num_bond[j]; m++) { + if (bond_atom[j][m] == tag[i]) { + nmax = num_bond[j] - 1; + bond_type[j][m] = bond_type[j][nmax]; + bond_atom[j][m] = bond_atom[j][nmax]; + if (n_histories > 0) + for (auto &ihistory: histories) { + auto fix_bond_history = dynamic_cast (ihistory); + fix_bond_history->shift_history(j, m, nmax); + fix_bond_history->delete_history(j, nmax); + } + bond_type[j][nmax] = 0; num_bond[j]--; break; } } } - num_bond[i] = 0; } }