From e3768f9ab05bcfd67244abde0c4fcc24098525fa Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 1 Oct 2024 10:01:12 -0600 Subject: [PATCH 01/27] Limiting coordination calculation in RHEO --- src/RHEO/compute_rheo_kernel.cpp | 3 +++ src/RHEO/compute_rheo_property_atom.cpp | 8 +++++++- src/RHEO/compute_rheo_property_atom.h | 2 +- src/RHEO/compute_rheo_surface.cpp | 2 ++ src/RHEO/fix_rheo.cpp | 22 ++++++++++++++-------- src/RHEO/fix_rheo.h | 1 + 6 files changed, 28 insertions(+), 10 deletions(-) diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 453f2b9028..658abd3189 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -137,6 +137,9 @@ void ComputeRHEOKernel::init() } } + if (correction_order != -1) + fix_rheo->coordination_flag = 1; + nmax_store = atom->nmax; memory->create(coordination, nmax_store, "rheo:coordination"); if (kernel_style == RK0) { diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index fb5848ac75..fd05c0373e 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -79,7 +79,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a size_peratom_cols = nvalues; pressure_flag = thermal_flag = interface_flag = 0; - surface_flag = shift_flag = shell_flag = 0; + surface_flag = shift_flag = shell_flag = coordination_flag = 0; // parse input values // customize a new keyword by adding to if statement @@ -109,6 +109,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a surface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; } else if (strcmp(arg[iarg], "coordination") == 0) { + coordination_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; } else if (strcmp(arg[iarg], "pressure") == 0) { pressure_flag = 1; @@ -223,6 +224,11 @@ void ComputeRHEOPropertyAtom::compute_peratom() { invoked_peratom = update->ntimestep; + // calculate optional values, if needed + + if (coordination_flag && !(fix_rheo->coordination_flag)) + compute_kernel->compute_coordination(); + // grow vector or array if necessary if (atom->nmax > nmax) { diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index 14f28075a3..745423b4ae 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -36,7 +36,7 @@ class ComputeRHEOPropertyAtom : public Compute { private: int nvalues, nmax; int pressure_flag, thermal_flag, interface_flag; - int surface_flag, shift_flag, shell_flag; + int surface_flag, shift_flag, shell_flag, coordination_flag; int *avec_index; int *col_index, *col_t_index; double *buf; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 4ff9587f6c..b4f7fe7c1c 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -81,6 +81,8 @@ void ComputeRHEOSurface::init() threshold_splash = fix_rheo->zmin_splash; interface_flag = fix_rheo->interface_flag; + fix_rheo->coordination_flag = 1; + cutsq = cut * cut; // need an occasional half neighbor list diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 4a25a46501..a978333aba 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -70,6 +70,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : surface_flag = 0; oxidation_flag = 0; self_mass_flag = 0; + coordination_flag = 0; int i; int n = atom->ntypes; @@ -398,34 +399,39 @@ void FixRHEO::initial_integrate(int /*vflag*/) void FixRHEO::pre_force(int /*vflag*/) { - compute_kernel->compute_coordination(); // Needed for rho sum + if (coordination_flag) + compute_kernel->compute_coordination(); - if (rhosum_flag) compute_rhosum->compute_peratom(); + if (rhosum_flag) + compute_rhosum->compute_peratom(); compute_kernel->compute_peratom(); - if (interface_flag) { - // Note on first setup, have no forces for pressure to reference + // Note on first setup, have no forces for pressure to reference + if (interface_flag) compute_interface->compute_peratom(); - } // No need to forward v, rho, or T for compute_grad since already done compute_grad->compute_peratom(); compute_grad->forward_gradients(); - if (shift_flag) compute_vshift->compute_peratom(); + // Depends on NO_SHIFT status + if (shift_flag) + compute_vshift->compute_peratom(); // Remove temporary options int *mask = atom->mask; int *status = atom->rheo_status; int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) - if (mask[i] & groupbit) status[i] &= OPTIONSMASK; + if (mask[i] & groupbit) + status[i] &= OPTIONSMASK; // Calculate surfaces, update status if (surface_flag) { compute_surface->compute_peratom(); - if (shift_flag) compute_vshift->correct_surfaces(); + if (shift_flag) + compute_vshift->correct_surfaces(); } } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 24e86b21d6..4b75299f27 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -53,6 +53,7 @@ class FixRHEO : public Fix { int interface_flag; int surface_flag; int oxidation_flag; + int coordination_flag; int viscosity_fix_defined; int pressure_fix_defined; From a3c23f2ce837e41811bacdf807233e03aa297fee Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 2 Oct 2024 17:35:27 -0600 Subject: [PATCH 02/27] Cleaning up rheo kernels --- src/RHEO/compute_rheo_kernel.cpp | 134 +++++++++++++++---------------- src/RHEO/compute_rheo_kernel.h | 9 ++- 2 files changed, 71 insertions(+), 72 deletions(-) diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 658abd3189..711995fc3e 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -189,27 +189,26 @@ double ComputeRHEOKernel::calc_w_self() double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r) { + if (kernel_style == WENDLANDC4) + return calc_w_wendlandc4(r); + if (kernel_style == QUINTIC) + return calc_w_quintic(r); + double w = 0.0; - int corrections_i, corrections_j, corrections; - - if (kernel_style == WENDLANDC4) return calc_w_wendlandc4(r); - - if (kernel_style != QUINTIC) { - corrections_i = check_corrections(i); - corrections_j = check_corrections(j); - corrections = corrections_i & corrections_j; - } else { - corrections = 0; - } + int corrections_i = check_corrections(i); + int corrections_j = check_corrections(j); + int corrections = corrections_i && corrections_j; if (!corrections) - w = calc_w_quintic(r); - else if (kernel_style == RK0) + return calc_w_quintic(r); + + double dx[3] = {delx, dely, delz}; + if (kernel_style == RK0) w = calc_w_rk0(i, j, r); else if (kernel_style == RK1) - w = calc_w_rk1(i, j, delx, dely, delz, r); + w = calc_w_rk1(i, j, dx, r); else if (kernel_style == RK2) - w = calc_w_rk2(i, j, delx, dely, delz, r); + w = calc_w_rk2(i, j, dx, r); return w; } @@ -218,26 +217,27 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double delz, double r) { + if (kernel_style == WENDLANDC4) + return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji); + if (kernel_style == QUINTIC) + return calc_dw_quintic(delx, dely, delz, r, dWij, dWji); + double wp; - int corrections_i, corrections_j; + int corrections_i = check_corrections(i); + int corrections_j = check_corrections(j); - if (kernel_style == WENDLANDC4) return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji); - - if (kernel_style != QUINTIC) { - corrections_i = check_corrections(i); - corrections_j = check_corrections(j); - } - - // Calc wp and default dW's, a bit inefficient but can redo later - wp = calc_dw_quintic(delx, dely, delz, r, dWij, dWji); + wp = calc_dw_scalar_quintic(delx, dely, delz, r); // Overwrite if there are corrections + double dxij[3] = {delx, dely, delz}; + double dxji[3] = {-delx, -dely, -delz}; + if (kernel_style == RK1) { - if (corrections_i) calc_dw_rk1(i, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk1(j, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk1(i, dxij, r, dWij); + if (corrections_j) calc_dw_rk1(j, dxji, r, dWji); } else if (kernel_style == RK2) { - if (corrections_i) calc_dw_rk2(i, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk2(j, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk2(i, dxij, r, dWij); + if (corrections_j) calc_dw_rk2(j, dxji, r, dWji); } return wp; @@ -278,10 +278,9 @@ double ComputeRHEOKernel::calc_w_quintic(double r) /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r, - double *dW1, double *dW2) +double ComputeRHEOKernel::calc_dw_scalar_quintic(double delx, double dely, double delz, double r) { - double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv; + double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s; s = r * 3.0 * cutinv; @@ -303,14 +302,23 @@ double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, } wp *= pre_wp; - wprinv = wp / r; + + return wp; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r, + double *dW1, double *dW2) +{ + double wp = calc_dw_scalar_quintic(delx, dely, delz, r); + double wprinv = wp / r; + dW1[0] = delx * wprinv; dW1[1] = dely * wprinv; dW1[2] = delz * wprinv; - dW2[0] = -delx * wprinv; - dW2[1] = -dely * wprinv; - dW2[2] = -delz * wprinv; + scale3(-1.0, dW1, dW2); return wp; } @@ -364,9 +372,7 @@ double ComputeRHEOKernel::calc_dw_wendlandc4(double delx, double dely, double de dW1[1] = dely * wprinv; dW1[2] = delz * wprinv; - dW2[0] = -delx * wprinv; - dW2[1] = -dely * wprinv; - dW2[2] = -delz * wprinv; + scale3(-1.0, dW1, dW2); return wp; } @@ -387,14 +393,11 @@ double ComputeRHEOKernel::calc_w_rk0(int i, int j, double r) /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, double delz, double r) +double ComputeRHEOKernel::calc_w_rk1(int i, int j, double *dx, double r) { int b; - double w, dx[3], H[MAX_MDIM]; + double w, H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; w = calc_w_quintic(r); if (dim == 2) { @@ -429,13 +432,11 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, double delz, double r) +double ComputeRHEOKernel::calc_w_rk2(int i, int j, double *dx, double r) { int b; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; + w = calc_w_quintic(r); if (dim == 2) { @@ -479,14 +480,10 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou /* ---------------------------------------------------------------------- */ -void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz, double r, - double *dW) +void ComputeRHEOKernel::calc_dw_rk1(int i, double *dx, double r, double *dW) { int a, b; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; w = calc_w_quintic(r); @@ -504,8 +501,8 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz // dWij[] = dWx dWy (dWz) //compute derivative operators + zero3(dW); for (a = 0; a < dim; a++) { - dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) @@ -516,14 +513,10 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz /* ---------------------------------------------------------------------- */ -void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz, double r, - double *dW) +void ComputeRHEOKernel::calc_dw_rk2(int i, double *dx, double r, double *dW) { int a, b; - double w, dx[3], H[MAX_MDIM]; - dx[0] = delx; - dx[1] = dely; - dx[2] = delz; + double w, H[MAX_MDIM]; w = calc_w_quintic(r); @@ -550,8 +543,8 @@ void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz // dWij[] = dWx dWy (dWz) //compute derivative operators + zero3(dW); for (a = 0; a < dim; a++) { - dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) @@ -567,7 +560,7 @@ void ComputeRHEOKernel::compute_peratom() lapack_error_flag = 0; lapack_error_tags.clear(); - if (kernel_style == QUINTIC) return; + if (correction_order == -1) return; corrections_calculated = 1; int i, j, ii, jj, inum, jnum, a, b, lapack_error; @@ -669,6 +662,7 @@ void ComputeRHEOKernel::compute_peratom() w = calc_w_quintic(r); rhoj = rho[j]; + if (interface_flag) if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j); @@ -856,7 +850,9 @@ void ComputeRHEOKernel::grow_arrays(int nmax) int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { + int a, b; int m = 0; + for (int i = 0; i < n; i++) { int j = list[i]; if (comm_stage == 0) { @@ -865,8 +861,8 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb if (kernel_style == RK0) { buf[m++] = C0[j]; } else { - for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) buf[m++] = C[j][a][b]; + for (a = 0; a < ncor; a++) + for (b = 0; b < Mdim; b++) buf[m++] = C[j][a][b]; } } } @@ -877,8 +873,10 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) { + int a, b; int m = 0; int last = first + n; + for (int i = first; i < last; i++) { if (comm_stage == 0) { coordination[i] = buf[m++]; @@ -886,8 +884,8 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) if (kernel_style == RK0) { C0[i] = buf[m++]; } else { - for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; + for (a = 0; a < ncor; a++) + for (b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; } } } diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index 8b70509e6a..2f458d14c9 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -40,6 +40,7 @@ class ComputeRHEOKernel : public Compute { double calc_w(int, int, double, double, double, double); double calc_dw(int, int, double, double, double, double); double calc_w_quintic(double); + double calc_dw_scalar_quintic(double, double, double, double); double calc_dw_quintic(double, double, double, double, double *, double *); double calc_w_wendlandc4(double); double calc_dw_wendlandc4(double, double, double, double, double *, double *); @@ -69,10 +70,10 @@ class ComputeRHEOKernel : public Compute { int check_corrections(int); double calc_w_rk0(int, int, double); - double calc_w_rk1(int, int, double, double, double, double); - double calc_w_rk2(int, int, double, double, double, double); - void calc_dw_rk1(int, double, double, double, double, double *); - void calc_dw_rk2(int, double, double, double, double, double *); + double calc_w_rk1(int, int, double*, double); + double calc_w_rk2(int, int, double*, double); + void calc_dw_rk1(int, double*, double, double*); + void calc_dw_rk2(int, double*, double, double*); }; } // namespace LAMMPS_NS #endif From 47e9a7fe79715ad43811257d7a121d382d43d2b4 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 29 Oct 2024 21:03:47 -0600 Subject: [PATCH 03/27] Typo + style changes in rheo doc --- doc/src/Howto_rheo.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index 90f9ff82d0..4d5734271a 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -70,7 +70,7 @@ particles to solid (e.g. with the :doc:`set ` command), (b) create bpm bonds between the particles (see the :doc:`bpm howto ` page for more details), and (c) use :doc:`pair rheo/solid ` to apply repulsive contact forces between distinct solid bodies. Akin to pair rheo, -pair rheo/solid considers a particles fluid/solid phase to determine whether to +pair rheo/solid considers a particle's fluid/solid phase to determine whether to apply forces. However, unlike pair rheo, pair rheo/solid does obey special bond settings such that contact forces do not have to be calculated between two bonded solid particles in the same elastic body. @@ -79,10 +79,10 @@ In systems with thermal evolution, fix rheo/thermal can optionally set a melting/solidification temperature allowing particles to dynamically swap their state between fluid and solid when the temperature exceeds or drops below the critical temperature, respectively. Using the *react* option, one can specify a maximum -bond length and a bond type. Then, when solidifying, particles will search their +bond length and a bond type. Then, when solidifying, particles search their local neighbors and automatically create bonds with any neighboring solid particles -in range. For BPM bond styles, bonds will then use the immediate position of the two -particles to calculate a reference state. When melting, particles will delete any +in range. For BPM bond styles, bonds then use the immediate position of the two +particles to calculate a reference state. When melting, particles delete any bonds of the specified type when reverting to a fluid state. Special bonds are updated as bonds are created/broken. From e1f14e1ee81b5cee442d0e2518b40e14b8607cc9 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 30 Oct 2024 08:02:35 -0600 Subject: [PATCH 04/27] unifying comm for/rev method args --- src/comm.cpp | 8 ++-- src/comm.h | 16 +++---- src/comm_brick.cpp | 108 ++++++++++++++++++++++++++++++-------------- src/comm_brick.h | 22 ++++----- src/comm_tiled.cpp | 109 ++++++++++++++++++++++++++++++--------------- src/comm_tiled.h | 22 ++++----- 6 files changed, 182 insertions(+), 103 deletions(-) diff --git a/src/comm.cpp b/src/comm.cpp index 02999fd541..8879352009 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -227,13 +227,13 @@ void Comm::init() } for (const auto &compute : modify->get_compute_list()) { - maxforward = MAX(maxforward,compute->comm_forward); - maxreverse = MAX(maxreverse,compute->comm_reverse); + maxforward = MAX(maxforward, compute->comm_forward); + maxreverse = MAX(maxreverse, compute->comm_reverse); } for (const auto &dump: output->get_dump_list()) { - maxforward = MAX(maxforward,dump->comm_forward); - maxreverse = MAX(maxreverse,dump->comm_reverse); + maxforward = MAX(maxforward, dump->comm_forward); + maxreverse = MAX(maxreverse, dump->comm_reverse); } if (force->newton == 0) maxreverse = 0; diff --git a/src/comm.h b/src/comm.h index fde4c3b81f..d9b4ac0dbf 100644 --- a/src/comm.h +++ b/src/comm.h @@ -87,17 +87,17 @@ class Comm : protected Pointers { // forward/reverse comm from a Pair, Bond, Fix, Compute, Dump - virtual void forward_comm(class Pair *) = 0; - virtual void reverse_comm(class Pair *) = 0; - virtual void forward_comm(class Bond *) = 0; - virtual void reverse_comm(class Bond *) = 0; + virtual void forward_comm(class Pair *, int size = 0) = 0; + virtual void reverse_comm(class Pair *, int size = 0) = 0; + virtual void forward_comm(class Bond *, int size = 0) = 0; + virtual void reverse_comm(class Bond *, int size = 0) = 0; virtual void forward_comm(class Fix *, int size = 0) = 0; virtual void reverse_comm(class Fix *, int size = 0) = 0; virtual void reverse_comm_variable(class Fix *) = 0; - virtual void forward_comm(class Compute *) = 0; - virtual void reverse_comm(class Compute *) = 0; - virtual void forward_comm(class Dump *) = 0; - virtual void reverse_comm(class Dump *) = 0; + virtual void forward_comm(class Compute *, int size = 0) = 0; + virtual void reverse_comm(class Compute *, int size = 0) = 0; + virtual void forward_comm(class Dump *, int size = 0) = 0; + virtual void reverse_comm(class Dump *, int size = 0) = 0; // forward comm of an array diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index cf38271029..9b8fc588ea 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -980,16 +980,21 @@ void CommBrick::borders() /* ---------------------------------------------------------------------- forward communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Pair + size > 0 -> Pair passes max size per atom + the latter is only useful if Pair does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommBrick::forward_comm(Pair *pair) +void CommBrick::forward_comm(Pair *pair, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = pair->comm_forward; + if (size) nsize = size; + else nsize = pair->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1017,16 +1022,21 @@ void CommBrick::forward_comm(Pair *pair) /* ---------------------------------------------------------------------- reverse communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Pair + size > 0 -> Pair passes max size per atom + the latter is only useful if Pair does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommBrick::reverse_comm(Pair *pair) +void CommBrick::reverse_comm(Pair *pair, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); + if (size) nsize = size; + else nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); for (iswap = nswap-1; iswap >= 0; iswap--) { @@ -1054,16 +1064,21 @@ void CommBrick::reverse_comm(Pair *pair) /* ---------------------------------------------------------------------- forward communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommBrick::forward_comm(Bond *bond) +void CommBrick::forward_comm(Bond *bond, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = bond->comm_forward; + if (size) nsize = size; + else nsize = bond->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1091,16 +1106,21 @@ void CommBrick::forward_comm(Bond *bond) /* ---------------------------------------------------------------------- reverse communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommBrick::reverse_comm(Bond *bond) +void CommBrick::reverse_comm(Bond *bond, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); + if (size) nsize = size; + else nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); for (iswap = nswap-1; iswap >= 0; iswap--) { @@ -1171,10 +1191,10 @@ void CommBrick::forward_comm(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix size/nsize used only to set recv buffer limit - size = 0 (default) -> use comm_forward from Fix + size = 0 (default) -> use comm_reverse from Fix size > 0 -> Fix passes max size per atom the latter is only useful if Fix does several comm modes, - some are smaller than max stored in its comm_forward + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ void CommBrick::reverse_comm(Fix *fix, int size) @@ -1255,16 +1275,21 @@ void CommBrick::reverse_comm_variable(Fix *fix) /* ---------------------------------------------------------------------- forward communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommBrick::forward_comm(Compute *compute) +void CommBrick::forward_comm(Compute *compute, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = compute->comm_forward; + if (size) nsize = size; + else nsize = compute->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1293,16 +1318,21 @@ void CommBrick::forward_comm(Compute *compute) /* ---------------------------------------------------------------------- reverse communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommBrick::reverse_comm(Compute *compute) +void CommBrick::reverse_comm(Compute *compute, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = compute->comm_reverse; + if (size) nsize = size; + else nsize = compute->comm_reverse; for (iswap = nswap-1; iswap >= 0; iswap--) { @@ -1330,16 +1360,21 @@ void CommBrick::reverse_comm(Compute *compute) /* ---------------------------------------------------------------------- forward communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommBrick::forward_comm(Dump *dump) +void CommBrick::forward_comm(Dump *dump, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = dump->comm_forward; + if (size) nsize = size; + else nsize = dump->comm_forward; for (iswap = 0; iswap < nswap; iswap++) { @@ -1368,16 +1403,21 @@ void CommBrick::forward_comm(Dump *dump) /* ---------------------------------------------------------------------- reverse communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommBrick::reverse_comm(Dump *dump) +void CommBrick::reverse_comm(Dump *dump, int size) { - int iswap,n; + int iswap,n,nsize; double *buf; MPI_Request request; - int nsize = dump->comm_reverse; + if (size) nsize = size; + else nsize = dump->comm_reverse; for (iswap = nswap-1; iswap >= 0; iswap--) { diff --git a/src/comm_brick.h b/src/comm_brick.h index 518e290fa1..a321daee2b 100644 --- a/src/comm_brick.h +++ b/src/comm_brick.h @@ -32,17 +32,17 @@ class CommBrick : public Comm { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix - void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix - void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond + void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond + void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix + void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix + void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix + void forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump + void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump void forward_comm_array(int, double **) override; // forward comm of array void *extract(const char *, int &) override; diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index b864e0523d..01a5bb1f69 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -1380,14 +1380,19 @@ void CommTiled::borders() /* ---------------------------------------------------------------------- forward communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Pair + size > 0 -> Pair passes max size per atom + the latter is only useful if Pair does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiled::forward_comm(Pair *pair) +void CommTiled::forward_comm(Pair *pair, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = pair->comm_forward; + if (size) nsize = size; + else nsize = pair->comm_forward; for (int iswap = 0; iswap < nswap; iswap++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1424,14 +1429,19 @@ void CommTiled::forward_comm(Pair *pair) /* ---------------------------------------------------------------------- reverse communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Fix + size > 0 -> Fix passes max size per atom + the latter is only useful if Fix does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiled::reverse_comm(Pair *pair) +void CommTiled::reverse_comm(Pair *pair, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); + if (size) nsize = MAX(pair->comm_reverse, pair->comm_reverse_off); + else nsize = pair->comm_reverse; for (int iswap = nswap-1; iswap >= 0; iswap--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1464,14 +1474,18 @@ void CommTiled::reverse_comm(Pair *pair) /* ---------------------------------------------------------------------- forward communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ - -void CommTiled::forward_comm(Bond *bond) +void CommTiled::forward_comm(Bond *bond, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = bond->comm_forward; + if (size) nsize = size; + else nsize = bond->comm_forward; for (int iswap = 0; iswap < nswap; iswap++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1508,14 +1522,19 @@ void CommTiled::forward_comm(Bond *bond) /* ---------------------------------------------------------------------- reverse communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiled::reverse_comm(Bond *bond) +void CommTiled::reverse_comm(Bond *bond, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); + if (size) nsize = size; + else nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); for (int iswap = nswap-1; iswap >= 0; iswap--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1596,10 +1615,10 @@ void CommTiled::forward_comm(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix size/nsize used only to set recv buffer limit - size = 0 (default) -> use comm_forward from Fix + size = 0 (default) -> use comm_reverse from Fix size > 0 -> Fix passes max size per atom the latter is only useful if Fix does several comm modes, - some are smaller than max stored in its comm_forward + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ void CommTiled::reverse_comm(Fix *fix, int size) @@ -1652,14 +1671,19 @@ void CommTiled::reverse_comm_variable(Fix * /*fix*/) /* ---------------------------------------------------------------------- forward communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiled::forward_comm(Compute *compute) +void CommTiled::forward_comm(Compute *compute, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = compute->comm_forward; + if (size) nsize = size; + else nsize = compute->comm_forward; for (int iswap = 0; iswap < nswap; iswap++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1695,14 +1719,19 @@ void CommTiled::forward_comm(Compute *compute) /* ---------------------------------------------------------------------- reverse communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiled::reverse_comm(Compute *compute) +void CommTiled::reverse_comm(Compute *compute, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = compute->comm_reverse; + if (size) nsize = size; + else nsize = compute->comm_reverse; for (int iswap = nswap-1; iswap >= 0; iswap--) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1736,14 +1765,19 @@ void CommTiled::reverse_comm(Compute *compute) /* ---------------------------------------------------------------------- forward communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiled::forward_comm(Dump *dump) +void CommTiled::forward_comm(Dump *dump, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = dump->comm_forward; + if (size) nsize = size; + else nsize = dump->comm_forward; for (int iswap = 0; iswap < nswap; iswap++) { nsend = nsendproc[iswap] - sendself[iswap]; @@ -1779,14 +1813,19 @@ void CommTiled::forward_comm(Dump *dump) /* ---------------------------------------------------------------------- reverse communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiled::reverse_comm(Dump *dump) +void CommTiled::reverse_comm(Dump *dump, int size) { - int i,irecv,n,nsend,nrecv; + int i,irecv,n,nsize,nsend,nrecv; - int nsize = dump->comm_reverse; + if (size) nsize = size; + else nsize = dump->comm_reverse; for (int iswap = nswap-1; iswap >= 0; iswap--) { nsend = nsendproc[iswap] - sendself[iswap]; diff --git a/src/comm_tiled.h b/src/comm_tiled.h index 857cddf033..74fcc0d1ff 100644 --- a/src/comm_tiled.h +++ b/src/comm_tiled.h @@ -32,17 +32,17 @@ class CommTiled : public Comm { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix - void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix - void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond + void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond + void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix + void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix + void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix + void forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump + void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump void forward_comm_array(int, double **) override; // forward comm of array From 3fd4f9b7f3b6b65a5434163e5b0c4626ac9f3d2f Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 30 Oct 2024 08:56:07 -0600 Subject: [PATCH 05/27] Minor clean ups to BPM bonds --- src/BPM/bond_bpm_rotational.cpp | 41 +++++++++++++++++---------------- src/BPM/bond_bpm_spring.cpp | 4 ++-- 2 files changed, 23 insertions(+), 22 deletions(-) diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index 2d5cb0652b..acc162ecbf 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -454,14 +454,12 @@ void BondBPMRotational::damping_forces(int i1, int i2, int type, double *rhat, d void BondBPMRotational::compute(int eflag, int vflag) { - if (!fix_bond_history->stored_flag) { fix_bond_history->stored_flag = true; store_data(); } - if (hybrid_flag) - fix_bond_history->compress_history(); + if (hybrid_flag) fix_bond_history->compress_history(); int i1, i2, itmp, n, type; double r[3], r0[3], rhat[3]; @@ -544,33 +542,36 @@ void BondBPMRotational::compute(int eflag, int vflag) // Apply forces and torques to particles // ------------------------------------------------------// - if (newton_bond || i1 < nlocal) { - f[i1][0] -= force1on2[0] * smooth; - f[i1][1] -= force1on2[1] * smooth; - f[i1][2] -= force1on2[2] * smooth; + MathExtra::scale3(smooth, force1on2); + MathExtra::scale3(smooth, torque2on1); - torque[i1][0] += torque2on1[0] * smooth; - torque[i1][1] += torque2on1[1] * smooth; - torque[i1][2] += torque2on1[2] * smooth; + if (newton_bond || i1 < nlocal) { + f[i1][0] -= force1on2[0]; + f[i1][1] -= force1on2[1]; + f[i1][2] -= force1on2[2]; + + torque[i1][0] += torque2on1[0]; + torque[i1][1] += torque2on1[1]; + torque[i1][2] += torque2on1[2]; } if (newton_bond || i2 < nlocal) { - f[i2][0] += force1on2[0] * smooth; - f[i2][1] += force1on2[1] * smooth; - f[i2][2] += force1on2[2] * smooth; + f[i2][0] += force1on2[0]; + f[i2][1] += force1on2[1]; + f[i2][2] += force1on2[2]; - torque[i2][0] += torque1on2[0] * smooth; - torque[i2][1] += torque1on2[1] * smooth; - torque[i2][2] += torque1on2[2] * smooth; + MathExtra::scale3(smooth, torque1on2); + torque[i2][0] += torque1on2[0]; + torque[i2][1] += torque1on2[1]; + torque[i2][2] += torque1on2[2]; } if (evflag) - ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0] * smooth, -force1on2[1] * smooth, - -force1on2[2] * smooth, r[0], r[1], r[2]); + ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0], -force1on2[1], + -force1on2[2], r[0], r[1], r[2]); } - if (hybrid_flag) - fix_bond_history->uncompress_history(); + if (hybrid_flag) fix_bond_history->uncompress_history(); } /* ---------------------------------------------------------------------- */ diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index 37437d5b60..0c61c20bf6 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -149,8 +149,6 @@ void BondBPMSpring::store_data() void BondBPMSpring::compute(int eflag, int vflag) { - if (hybrid_flag) fix_bond_history->compress_history(); - int i, bond_change_flag; double *vol0, *vol; @@ -195,6 +193,8 @@ void BondBPMSpring::compute(int eflag, int vflag) } } + if (hybrid_flag) fix_bond_history->compress_history(); + int i1, i2, itmp, n, type; double delx, dely, delz, delvx, delvy, delvz; double e, rsq, r, r0, rinv, smooth, fbond, dot; From df882a9552b96e0724d9024d8241f1703076e310 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 31 Oct 2024 15:27:22 -0600 Subject: [PATCH 06/27] Adding optional size arg for forward/reverse comm methods to Kokkos --- src/KOKKOS/comm_kokkos.cpp | 83 +++++++++++++++++++------------ src/KOKKOS/comm_kokkos.h | 26 +++++----- src/KOKKOS/comm_tiled_kokkos.cpp | 84 ++++++++++++++++++++++---------- src/KOKKOS/comm_tiled_kokkos.h | 22 ++++----- 4 files changed, 135 insertions(+), 80 deletions(-) diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp index 8f821c3036..eea79248fe 100644 --- a/src/KOKKOS/comm_kokkos.cpp +++ b/src/KOKKOS/comm_kokkos.cpp @@ -455,10 +455,10 @@ void CommKokkos::forward_comm_device(Fix *fix, int size) /* ---------------------------------------------------------------------- reverse communication invoked by a Fix size/nsize used only to set recv buffer limit - size = 0 (default) -> use comm_forward from Fix + size = 0 (default) -> use comm_reverse from Fix size > 0 -> Fix passes max size per atom the latter is only useful if Fix does several comm modes, - some are smaller than max stored in its comm_forward + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ void CommKokkos::reverse_comm(Fix *fix, int size) @@ -482,72 +482,94 @@ void CommKokkos::reverse_comm_variable(Fix *fix) /* ---------------------------------------------------------------------- forward communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommKokkos::forward_comm(Compute *compute) +void CommKokkos::forward_comm(Compute *compute, int size) { k_sendlist.sync(); - CommBrick::forward_comm(compute); + CommBrick::forward_comm(compute, size); } /* ---------------------------------------------------------------------- forward communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommKokkos::forward_comm(Bond *bond) +void CommKokkos::forward_comm(Bond *bond, int size) { - CommBrick::forward_comm(bond); + CommBrick::forward_comm(bond, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommKokkos::reverse_comm(Bond *bond) +void CommKokkos::reverse_comm(Bond *bond, int size) { - CommBrick::reverse_comm(bond); + CommBrick::reverse_comm(bond, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommKokkos::reverse_comm(Compute *compute) +void CommKokkos::reverse_comm(Compute *compute, int size) { k_sendlist.sync(); - CommBrick::reverse_comm(compute); + CommBrick::reverse_comm(compute, size); } /* ---------------------------------------------------------------------- forward communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommKokkos::forward_comm(Pair *pair) +void CommKokkos::forward_comm(Pair *pair, int size) { if (pair->execution_space == Host || forward_pair_comm_classic) { k_sendlist.sync(); - CommBrick::forward_comm(pair); + CommBrick::forward_comm(pair, size); } else { k_sendlist.sync(); - forward_comm_device(pair); + forward_comm_device(pair, size); } } /* ---------------------------------------------------------------------- */ template -void CommKokkos::forward_comm_device(Pair *pair) +void CommKokkos::forward_comm_device(Pair *pair, int size) { - int iswap,n; + int iswap,n,nsize; MPI_Request request; DAT::tdual_xfloat_1d k_buf_tmp; - int nsize = pair->comm_forward; + if (size) nsize = size; + else nsize = pair->comm_forward; + KokkosBase* pairKKBase = dynamic_cast(pair); int nmax = max_buf_pair; @@ -623,21 +645,21 @@ void CommKokkos::grow_buf_fix(int n) { /* ---------------------------------------------------------------------- */ -void CommKokkos::reverse_comm(Pair *pair) +void CommKokkos::reverse_comm(Pair *pair, int size) { if (pair->execution_space == Host || !pair->reverse_comm_device || reverse_pair_comm_classic) { k_sendlist.sync(); - CommBrick::reverse_comm(pair); + CommBrick::reverse_comm(pair, size); } else { k_sendlist.sync(); - reverse_comm_device(pair); + reverse_comm_device(pair, size); } } /* ---------------------------------------------------------------------- */ template -void CommKokkos::reverse_comm_device(Pair *pair) +void CommKokkos::reverse_comm_device(Pair *pair, int size) { int iswap,n; MPI_Request request; @@ -645,7 +667,8 @@ void CommKokkos::reverse_comm_device(Pair *pair) KokkosBase* pairKKBase = dynamic_cast(pair); - int nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); + if (size) nsize = size; + else nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); int nmax = max_buf_pair; for (iswap = 0; iswap < nswap; iswap++) { @@ -702,18 +725,18 @@ void CommKokkos::reverse_comm_device(Pair *pair) /* ---------------------------------------------------------------------- */ -void CommKokkos::forward_comm(Dump *dump) +void CommKokkos::forward_comm(Dump *dump, int size) { k_sendlist.sync(); - CommBrick::forward_comm(dump); + CommBrick::forward_comm(dump, size); } /* ---------------------------------------------------------------------- */ -void CommKokkos::reverse_comm(Dump *dump) +void CommKokkos::reverse_comm(Dump *dump, int size) { k_sendlist.sync(); - CommBrick::reverse_comm(dump); + CommBrick::reverse_comm(dump, size); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/comm_kokkos.h b/src/KOKKOS/comm_kokkos.h index 4fb4dfbe29..42941ff517 100644 --- a/src/KOKKOS/comm_kokkos.h +++ b/src/KOKKOS/comm_kokkos.h @@ -45,24 +45,24 @@ class CommKokkos : public CommBrick { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix - void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix - void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond + void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond + void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix + void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix + void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix + void forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump + void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump void forward_comm_array(int, double **) override; // forward comm of array template void forward_comm_device(); template void reverse_comm_device(); - template void forward_comm_device(Pair *pair); - template void reverse_comm_device(Pair *pair); + template void forward_comm_device(Pair *pair, int size=0); + template void reverse_comm_device(Pair *pair, int size=0); template void forward_comm_device(Fix *fix, int size=0); template void exchange_device(); template void borders_device(); diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp index 2e4ca30bed..69c5e8f847 100644 --- a/src/KOKKOS/comm_tiled_kokkos.cpp +++ b/src/KOKKOS/comm_tiled_kokkos.cpp @@ -417,7 +417,11 @@ void CommTiledKokkos::borders() /* ---------------------------------------------------------------------- forward communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Fix + size > 0 -> Fix passes max size per atom + the latter is only useful if Fix does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ void CommTiledKokkos::forward_comm(Pair *pair) @@ -427,32 +431,44 @@ void CommTiledKokkos::forward_comm(Pair *pair) /* ---------------------------------------------------------------------- reverse communication invoked by a Pair - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Pair + size > 0 -> Pair passes max size per atom + the latter is only useful if Pair does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiledKokkos::reverse_comm(Pair *pair) +void CommTiledKokkos::reverse_comm(Pair *pair, int size) { - CommTiled::reverse_comm(pair); + CommTiled::reverse_comm(pair, size); } /* ---------------------------------------------------------------------- forward communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiledKokkos::forward_comm(Bond *bond) +void CommTiledKokkos::forward_comm(Bond *bond, int size) { - CommTiled::forward_comm(bond); + CommTiled::forward_comm(bond, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Bond - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Bond + size > 0 -> Bond passes max size per atom + the latter is only useful if Bond does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiledKokkos::reverse_comm(Bond *bond) +void CommTiledKokkos::reverse_comm(Bond *bond, int size) { - CommTiled::reverse_comm(bond); + CommTiled::reverse_comm(bond, size); } /* ---------------------------------------------------------------------- @@ -466,21 +482,21 @@ void CommTiledKokkos::reverse_comm(Bond *bond) void CommTiledKokkos::forward_comm(Fix *fix, int size) { - CommTiled::forward_comm(fix,size); + CommTiled::forward_comm(fix, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Fix size/nsize used only to set recv buffer limit - size = 0 (default) -> use comm_forward from Fix + size = 0 (default) -> use comm_reverse from Fix size > 0 -> Fix passes max size per atom the latter is only useful if Fix does several comm modes, - some are smaller than max stored in its comm_forward + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ void CommTiledKokkos::reverse_comm(Fix *fix, int size) { - CommTiled::reverse_comm(fix,size); + CommTiled::reverse_comm(fix, size); } /* ---------------------------------------------------------------------- @@ -497,42 +513,58 @@ void CommTiledKokkos::reverse_comm_variable(Fix *fix) /* ---------------------------------------------------------------------- forward communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiledKokkos::forward_comm(Compute *compute) +void CommTiledKokkos::forward_comm(Compute *compute, int size) { - CommTiled::forward_comm(compute); + CommTiled::forward_comm(compute, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Compute - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Compute + size > 0 -> Compute passes max size per atom + the latter is only useful if Compute does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiledKokkos::reverse_comm(Compute *compute) +void CommTiledKokkos::reverse_comm(Compute *compute, int size) { - CommTiled::reverse_comm(compute); + CommTiled::reverse_comm(compute, size); } /* ---------------------------------------------------------------------- forward communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_forward from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiledKokkos::forward_comm(Dump *dump) +void CommTiledKokkos::forward_comm(Dump *dump, int size) { - CommTiled::forward_comm(dump); + CommTiled::forward_comm(dump, size); } /* ---------------------------------------------------------------------- reverse communication invoked by a Dump - nsize used only to set recv buffer limit + size/nsize used only to set recv buffer limit + size = 0 (default) -> use comm_reverse from Dump + size > 0 -> Dump passes max size per atom + the latter is only useful if Dump does several comm modes, + some are smaller than max stored in its comm_reverse ------------------------------------------------------------------------- */ -void CommTiledKokkos::reverse_comm(Dump *dump) +void CommTiledKokkos::reverse_comm(Dump *dump, int size) { - CommTiled::reverse_comm(dump); + CommTiled::reverse_comm(dump, size); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/comm_tiled_kokkos.h b/src/KOKKOS/comm_tiled_kokkos.h index 9033714796..ef226489c8 100644 --- a/src/KOKKOS/comm_tiled_kokkos.h +++ b/src/KOKKOS/comm_tiled_kokkos.h @@ -46,17 +46,17 @@ class CommTiledKokkos : public CommTiled { void exchange() override; // move atoms to new procs void borders() override; // setup list of atoms to comm - void forward_comm(class Pair *) override; // forward comm from a Pair - void reverse_comm(class Pair *) override; // reverse comm from a Pair - void forward_comm(class Bond *) override; // forward comm from a Bond - void reverse_comm(class Bond *) override; // reverse comm from a Bond - void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix - void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix - void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix - void forward_comm(class Compute *) override; // forward from a Compute - void reverse_comm(class Compute *) override; // reverse from a Compute - void forward_comm(class Dump *) override; // forward comm from a Dump - void reverse_comm(class Dump *) override; // reverse comm from a Dump + void forward_comm(class Pair *, int size = 0) override; // forward comm from a Pair + void reverse_comm(class Pair *, int size = 0) override; // reverse comm from a Pair + void forward_comm(class Bond *, int size = 0) override; // forward comm from a Bond + void reverse_comm(class Bond *, int size = 0) override; // reverse comm from a Bond + void forward_comm(class Fix *, int size = 0) override; // forward comm from a Fix + void reverse_comm(class Fix *, int size = 0) override; // reverse comm from a Fix + void reverse_comm_variable(class Fix *) override; // variable size reverse comm from a Fix + void forward_comm(class Compute *, int size = 0) override; // forward from a Compute + void reverse_comm(class Compute *, int size = 0) override; // reverse from a Compute + void forward_comm(class Dump *, int size = 0) override; // forward comm from a Dump + void reverse_comm(class Dump *, int size = 0) override; // reverse comm from a Dump void forward_comm_array(int, double **) override; // forward comm of array From a15a3153414094b1ddb9268ca371fcb1fbc6f50a Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 31 Oct 2024 15:38:54 -0600 Subject: [PATCH 07/27] Missing variable declaration --- src/KOKKOS/comm_kokkos.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp index eea79248fe..57c9d7c383 100644 --- a/src/KOKKOS/comm_kokkos.cpp +++ b/src/KOKKOS/comm_kokkos.cpp @@ -661,7 +661,7 @@ void CommKokkos::reverse_comm(Pair *pair, int size) template void CommKokkos::reverse_comm_device(Pair *pair, int size) { - int iswap,n; + int iswap,n,nsize; MPI_Request request; DAT::tdual_xfloat_1d k_buf_tmp; From 14041bb20c78309a35a85f4484f0bac709203acc Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 31 Oct 2024 15:56:45 -0600 Subject: [PATCH 08/27] Missing argument --- src/KOKKOS/comm_kokkos.cpp | 6 +++--- src/KOKKOS/comm_tiled_kokkos.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp index 57c9d7c383..e638a33306 100644 --- a/src/KOKKOS/comm_kokkos.cpp +++ b/src/KOKKOS/comm_kokkos.cpp @@ -378,10 +378,10 @@ void CommKokkos::forward_comm(Fix *fix, int size) { if (fix->execution_space == Host || !fix->forward_comm_device || forward_fix_comm_classic) { k_sendlist.sync(); - CommBrick::forward_comm(fix,size); + CommBrick::forward_comm(fix, size); } else { k_sendlist.sync(); - forward_comm_device(fix,size); + forward_comm_device(fix, size); } } @@ -668,7 +668,7 @@ void CommKokkos::reverse_comm_device(Pair *pair, int size) KokkosBase* pairKKBase = dynamic_cast(pair); if (size) nsize = size; - else nsize = MAX(pair->comm_reverse,pair->comm_reverse_off); + else nsize = MAX(pair->comm_reverse, pair->comm_reverse_off); int nmax = max_buf_pair; for (iswap = 0; iswap < nswap; iswap++) { diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp index 69c5e8f847..afddc079f4 100644 --- a/src/KOKKOS/comm_tiled_kokkos.cpp +++ b/src/KOKKOS/comm_tiled_kokkos.cpp @@ -424,9 +424,9 @@ void CommTiledKokkos::borders() some are smaller than max stored in its comm_forward ------------------------------------------------------------------------- */ -void CommTiledKokkos::forward_comm(Pair *pair) +void CommTiledKokkos::forward_comm(Pair *pair, int size) { - CommTiled::forward_comm(pair); + CommTiled::forward_comm(pair, size); } /* ---------------------------------------------------------------------- From f483a90d8a5a60a843e71585248c2c2fad4f67e2 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 1 Nov 2024 15:04:30 -0600 Subject: [PATCH 09/27] Simplifying RHEO compute comm calls --- src/RHEO/compute_rheo_grad.cpp | 6 ++---- src/RHEO/compute_rheo_interface.cpp | 6 ++---- src/RHEO/compute_rheo_kernel.cpp | 5 +---- src/RHEO/compute_rheo_kernel.h | 2 +- src/RHEO/compute_rheo_surface.cpp | 12 ++++-------- 5 files changed, 10 insertions(+), 21 deletions(-) diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index 63c77e3f3b..589f85e195 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -290,8 +290,7 @@ void ComputeRHEOGrad::compute_peratom() void ComputeRHEOGrad::forward_gradients() { comm_stage = COMMGRAD; - comm_forward = ncomm_grad; - comm->forward_comm(this); + comm->forward_comm(this, ncomm_grad); } /* ---------------------------------------------------------------------- */ @@ -299,8 +298,7 @@ void ComputeRHEOGrad::forward_gradients() void ComputeRHEOGrad::forward_fields() { comm_stage = COMMFIELD; - comm_forward = ncomm_field; - comm->forward_comm(this); + comm->forward_comm(this, ncomm_field); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index eec5659051..abfe6c2b0e 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -218,8 +218,7 @@ void ComputeRHEOInterface::compute_peratom() } comm_stage = 1; - comm_forward = 2; - comm->forward_comm(this); + comm->forward_comm(this, 2); } /* ---------------------------------------------------------------------- */ @@ -369,9 +368,8 @@ void ComputeRHEOInterface::store_forces() } // Forward comm forces - comm_forward = 3; comm_stage = 0; - comm->forward_comm(this); + comm->forward_comm(this, 3); } /* ---------------------------------------------------------------------- diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index d037894bd1..3a6629bfef 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -90,7 +90,6 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : comm_forward = ncor * Mdim; } - comm_forward_save = comm_forward; corrections_calculated = 0; lapack_error_flag = 0; } @@ -778,7 +777,6 @@ void ComputeRHEOKernel::compute_peratom() // communicate calculated quantities comm_stage = 1; - comm_forward = comm_forward_save; comm->forward_comm(this); } @@ -826,8 +824,7 @@ void ComputeRHEOKernel::compute_coordination() // communicate calculated quantities comm_stage = 0; - comm_forward = 1; - comm->forward_comm(this); + comm->forward_comm(this, 1); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index 2f458d14c9..037e9e3683 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -52,7 +52,7 @@ class ComputeRHEOKernel : public Compute { class FixRHEO *fix_rheo; private: - int comm_stage, comm_forward_save; + int comm_stage; int interface_flag; int lapack_error_flag; std::unordered_set lapack_error_tags; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index db78f0a4ff..74193f28f5 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -199,10 +199,8 @@ void ComputeRHEOSurface::compute_peratom() // reverse gradC and divr, forward divr comm_stage = 0; - comm_reverse = dim * dim + 1; - comm_forward = 1; - if (newton) comm->reverse_comm(this); - comm->forward_comm(this); + if (newton) comm->reverse_comm(this, dim * dim + 1); + comm->forward_comm(this, 1); // calculate nsurface for local atoms // Note, this isn't forwarded to ghosts @@ -299,10 +297,8 @@ void ComputeRHEOSurface::compute_peratom() // forward/reverse status and rsurface comm_stage = 1; - comm_reverse = 2; - comm_forward = 2; - if (newton) comm->reverse_comm(this); - comm->forward_comm(this); + if (newton) comm->reverse_comm(this, 2); + comm->forward_comm(this, 2); } /* ---------------------------------------------------------------------- */ From c1ac0cea1d06372510ca2ac7672386bbb7cdb485 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 1 Nov 2024 17:04:32 -0600 Subject: [PATCH 10/27] More rheo pressure EOS, allow variable cs --- doc/src/fix_rheo_pressure.rst | 25 +++++++++++- src/RHEO/compute_rheo_interface.cpp | 9 +++-- src/RHEO/compute_rheo_property_atom.cpp | 4 +- src/RHEO/fix_rheo_pressure.cpp | 52 ++++++++++++++++++++++--- src/RHEO/fix_rheo_pressure.h | 5 ++- src/RHEO/pair_rheo.cpp | 26 ++++++++++--- src/RHEO/pair_rheo.h | 1 + 7 files changed, 103 insertions(+), 19 deletions(-) diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index 2a714b298b..51f00b6684 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -19,8 +19,10 @@ Syntax .. parsed-literal:: *linear* args = none - *taitwater* args = none + *tait/water* args = none + *tait/general* args = exponent :math:`gamma` (unitless) and background pressure :math:`P[b]` (pressure) *cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2) + *ideal/gas* args = heat capacity ratio :math:`gamma` (unitless) Examples """""""" @@ -65,12 +67,31 @@ is a cubic equation of state which has an extra argument :math:`A_3`, P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . -Style *taitwater* is Tait's equation of state: +Style *tait/water* is Tait's equation of state: .. math:: P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr]. +Style *tait/general* generalizes this equation of state + +.. math:: + + P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr] + P[b]. + +where :math:`\gamma` is an exponent and :math:`P[b]` is a constant background pressure. + +Style *ideal/gas* is the ideal gas equation of state + +.. math:: + + P = (\gamma - 1) \rho e + +where :math:`\gamma` is the heat capacity ratio and :math:`e` is the internal energy of +a particle per unit mass. This style is only compatible with atom style rheo/thermal. +Note that when using this style, the speed of sound is no longer constant such that the +value of :math:`c` specified in :doc:`fix rheo ` is not used. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index abfe6c2b0e..b5a295aad6 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -97,6 +97,9 @@ void ComputeRHEOInterface::init() auto fixes = modify->get_fix_by_style("rheo/pressure"); fix_pressure = dynamic_cast(fixes[0]); + if (!fix_pressure->invertable_pressure) + error->all(FLERR, "RHEO interface reconstruction incompatible with pressure equation of state"); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); } @@ -178,7 +181,7 @@ void ComputeRHEOInterface::compute_peratom() 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); + rho[i] += w * (fix_pressure->calc_pressure(rho[j], j) - rho[j] * dot); normwf[i] += w; } } @@ -192,7 +195,7 @@ void ComputeRHEOInterface::compute_peratom() 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); + rho[j] += w * (fix_pressure->calc_pressure(rho[i], i) + rho[i] * dot); normwf[j] += w; } } @@ -210,7 +213,7 @@ void ComputeRHEOInterface::compute_peratom() if (status[i] & PHASECHECK) { if (normwf[i] != 0.0) { // Stores rho for solid particles 1+Pw in Adami Adams 2012 - rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], type[i])); + rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], i)); } else { rho[i] = rho0[itype]; } diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index fd05c0373e..4805c21903 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -440,7 +440,7 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); + buf[n] = fix_pressure->calc_pressure(rho[i], i); else buf[n] = 0.0; n += nvalues; @@ -490,7 +490,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (index == index_transpose) - p = fix_pressure->calc_pressure(rho[i], type[i]); + p = fix_pressure->calc_pressure(rho[i], i); else p = 0.0; buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p; diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 10bf6ff2c9..d4f66482b3 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL }; +enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL , IDEAL }; static constexpr double SEVENTH = 1.0 / 7.0; @@ -39,12 +39,14 @@ static constexpr double SEVENTH = 1.0 / 7.0; FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr), - rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr), - fix_rheo(nullptr) + rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), gamma(nullptr), + pressure_style(nullptr), fix_rheo(nullptr) { if (narg < 4) error->all(FLERR, "Illegal fix command"); comm_forward = 1; + variable_csq = 0; + invertable_pressure = 1; // Currently can only have one instance of fix rheo/pressure if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all"); @@ -55,6 +57,8 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : memory->create(c_cubic, n + 1, "rheo:c_cubic"); memory->create(tpower, n + 1, "rheo:tpower"); memory->create(pbackground, n + 1, "rheo:pbackground"); + memory->create(gamma, n + 1, "rheo:gamma"); + for (i = 1; i <= n; i++) pressure_style[i] = NONE; int iarg = 3; @@ -68,7 +72,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg + 1], "tait/water") == 0) { for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; } else if (strcmp(arg[iarg + 1], "tait/general") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait", error); + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait/general", error); double tpower_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); double pbackground_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); @@ -89,6 +93,22 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : pressure_style[i] = CUBIC; c_cubic[i] = c_cubic_one; } + + invertable_pressure = 0; + } else if (strcmp(arg[iarg + 1], "ideal/gas") == 0) { + if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure ideal/gas", error); + + double c_gamma_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 1; + + for (i = nlo; i <= nhi; i++) { + pressure_style[i] = IDEAL; + gamma[i] = c_gamma_one; + } + + variable_csq = 1; + if (atom->esph_flag != 1) + error->all(FLERR, "fix rheo/pressure ideal gas equation of state requires atom property esph"); } else { error->all(FLERR, "Illegal fix command, {}", arg[iarg]); } @@ -110,6 +130,7 @@ FixRHEOPressure::~FixRHEOPressure() memory->destroy(c_cubic); memory->destroy(tpower); memory->destroy(pbackground); + memory->destroy(gamma); } /* ---------------------------------------------------------------------- */ @@ -197,10 +218,11 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -double FixRHEOPressure::calc_pressure(double rho, int type) +double FixRHEOPressure::calc_pressure(double rho, int i) { double p = 0.0; double dr, rr3, rho_ratio; + int type = atom->type[i]; if (pressure_style[type] == LINEAR) { p = csq[type] * (rho - rho0[type]); @@ -215,15 +237,18 @@ double FixRHEOPressure::calc_pressure(double rho, int type) rho_ratio = rho * rho0inv[type]; p = csq[type] * rho0[type] * (pow(rho_ratio, tpower[type]) - 1.0) / tpower[type]; p += pbackground[type]; + } else if (pressure_style[type] == IDEAL) { + p = (gamma[type] - 1.0) * rho * atom->esph[i] / atom->mass[type]; } return p; } /* ---------------------------------------------------------------------- */ -double FixRHEOPressure::calc_rho(double p, int type) +double FixRHEOPressure::calc_rho(double p, int i) { double rho = 0.0; + int type = atom->type[i]; if (pressure_style[type] == LINEAR) { rho = csqinv[type] * p + rho0[type]; @@ -239,6 +264,21 @@ double FixRHEOPressure::calc_rho(double p, int type) rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]); rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]); rho *= pow(csq[type], -1.0 / tpower[type]); + } else if (pressure_style[type] == IDEAL) { + rho = p * atom->mass[type] / ((gamma[type] - 1.0) * atom->esph[i]); } return rho; } + +/* ---------------------------------------------------------------------- */ + +double FixRHEOPressure::calc_csq(double p, int i) +{ + int type = atom->type[i]; + double csq2 = csq[type]; + + if (pressure_style[type] == IDEAL) { + csq2 = (gamma[type] - 1.0) * atom->esph[i] / atom->mass[type]; + } + return csq2; +} diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index ca165b1ed5..ac60d1644d 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -36,9 +36,12 @@ class FixRHEOPressure : public Fix { void unpack_forward_comm(int, int, double *) override; double calc_pressure(double, int); double calc_rho(double, int); + double calc_csq(double, int); + int variable_csq; + int invertable_pressure; private: - double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground; + double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground, *gamma; int *pressure_style; class FixRHEO *fix_rheo; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index aa0c685aaa..5e06c45a3c 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -80,8 +80,8 @@ void PairRHEO::compute(int eflag, int vflag) int pair_force_flag, pair_rho_flag, pair_avisc_flag; int fluidi, fluidj; double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave; - double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, - kappa_ave, dT_prefactor; + double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj; + double eta_ave, kappa_ave, dT_prefactor; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; double *dWij, *dWji; double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3]; @@ -185,7 +185,13 @@ void PairRHEO::compute(int eflag, int vflag) kappaj = conductivity[j]; } - cs_ave = 0.5 * (cs[itype] + cs[jtype]); + if (!variable_csq) { + cs_ave = 0.5 * (cs[itype] + cs[jtype]); + } else { + csqi = fix_pressure->calc_csq(rhoi, i); + csqj = fix_pressure->calc_csq(rhoj, j); + cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj)); + } csq_ave = cs_ave * cs_ave; pair_rho_flag = 0; @@ -221,7 +227,7 @@ void PairRHEO::compute(int eflag, int vflag) if (fluidi && (!fluidj)) { compute_interface->correct_v(vj, vi, j, i); rhoj = compute_interface->correct_rho(j); - Pj = fix_pressure->calc_pressure(rhoj, jtype); + Pj = fix_pressure->calc_pressure(rhoj, j); if ((chi[j] > 0.9) && (r < (cutk * 0.5))) fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv; @@ -229,7 +235,7 @@ void PairRHEO::compute(int eflag, int vflag) } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vi, vj, i, j); rhoi = compute_interface->correct_rho(i); - Pi = fix_pressure->calc_pressure(rhoi, itype); + Pi = fix_pressure->calc_pressure(rhoi, i); if (chi[i] > 0.9 && r < (cutk * 0.5)) fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv; @@ -238,6 +244,14 @@ void PairRHEO::compute(int eflag, int vflag) rhoi = rho0i; rhoj = rho0j; } + + // recalculate speed of sound, if necessary + if (variable_csq && ((!fluidi) || (!fluidj))) { + csqi = fix_pressure->calc_csq(rhoi, i); + csqj = fix_pressure->calc_csq(rhoj, j); + cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj)); + csq_ave = cs_ave * cs_ave; + } } // Repel if close to inner solid particle @@ -480,6 +494,8 @@ void PairRHEO::setup() csq = fix_rheo->csq; rho0 = fix_rheo->rho0; + variable_csq = fix_pressure->variable_csq; + if (cutk != fix_rheo->cut) error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, fix_rheo->cut); diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index 444fcc2cb4..f1fcd10bf8 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -45,6 +45,7 @@ class PairRHEO : public Pair { int rho_damp_flag; int thermal_flag; int interface_flag; + int variable_csq; int harmonic_means_flag; From a93a930f2ff692f09679c42fc6ce7dc440a94d1a Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 5 Nov 2024 15:53:48 -0700 Subject: [PATCH 11/27] Updating developer comm doc page to include bond and optional args --- doc/src/Developer_comm_ops.rst | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/doc/src/Developer_comm_ops.rst b/doc/src/Developer_comm_ops.rst index 116b4d6ad8..66947de74f 100644 --- a/doc/src/Developer_comm_ops.rst +++ b/doc/src/Developer_comm_ops.rst @@ -79,19 +79,19 @@ containing ``double`` values. To correctly store integers that may be 64-bit (bigint, tagint, imageint) in the buffer, you need to use the :ref:`ubuf union ` construct. -The *Fix*, *Compute*, and *Dump* classes can also invoke the same kind -of forward and reverse communication operations using the same *Comm* -class methods. Likewise, the same pack/unpack methods and +The *Fix*, *Bond*, *Compute*, and *Dump* classes can also invoke the +same kind of forward and reverse communication operations using the +same *Comm* class methods. Likewise, the same pack/unpack methods and comm_forward/comm_reverse variables must be defined by the calling -*Fix*, *Compute*, or *Dump* class. +*Fix*, *Bond*, *Compute*, or *Dump* class. -For *Fix* classes, there is an optional second argument to the +For all of these classes, there is an optional second argument to the *forward_comm()* and *reverse_comm()* call which can be used when the -fix performs multiple modes of communication, with different numbers -of values per atom. The fix should set the *comm_forward* and +class performs multiple modes of communication, with different numbers +of values per atom. The class should set the *comm_forward* and *comm_reverse* variables to the maximum value, but can invoke the communication for a particular mode with a smaller value. For this -to work, the *pack_forward_comm()*, etc methods typically use a class +to work, the *pack_forward_comm()*, etc. methods typically use a class member variable to choose which values to pack/unpack into/from the buffer. From 009a976ae2bb309bb3ab0da658451ebe201392ab Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sat, 9 Nov 2024 21:47:38 -0700 Subject: [PATCH 12/27] Adding multitype correction to rheo vshfit --- src/RHEO/compute_rheo_vshift.cpp | 344 +++++++++++++++++++++++++++++-- src/RHEO/compute_rheo_vshift.h | 10 +- src/RHEO/fix_rheo.cpp | 44 ++-- src/RHEO/fix_rheo.h | 4 +- 4 files changed, 365 insertions(+), 37 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 8581d9ea8f..0047521e37 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -27,6 +27,7 @@ #include "error.h" #include "fix_rheo.h" #include "force.h" +#include "math_extra.h" #include "memory.h" #include "neigh_list.h" #include "neigh_request.h" @@ -36,20 +37,22 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; +using namespace MathExtra; /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr), - compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) + Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(nullptr), + fix_rheo(nullptr), rho0(nullptr), list(nullptr), compute_interface(nullptr), + compute_kernel(nullptr), compute_surface(nullptr) { if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command"); + comm_forward = 0; comm_reverse = 3; surface_flag = 0; - nmax_store = atom->nmax; - memory->create(vshift, nmax_store, 3, "rheo:vshift"); + nmax_store = 0; } /* ---------------------------------------------------------------------- */ @@ -73,9 +76,19 @@ void ComputeRHEOVShift::init() compute_surface = fix_rheo->compute_surface; rho0 = fix_rheo->rho0; + shift_type = fix_rheo->shift_type; cut = fix_rheo->cut; cutsq = cut * cut; cutthird = cut / 3.0; + + multiphase_flag = fix_rheo->shift_multiphase_flag; + if (multiphase_flag) { + scale = fix_rheo->shift_scale; + wmin = fix_rheo->shift_wmin; + cmin = fix_rheo->shift_cmin; + comm_forward = 1; + comm_reverse = 4; + } } /* ---------------------------------------------------------------------- */ @@ -120,6 +133,13 @@ void ComputeRHEOVShift::compute_peratom() if (nmax_store < atom->nmax) { memory->grow(vshift, atom->nmax, 3, "rheo:vshift"); + + if (multiphase_flag) { + memory->grow(ct, atom->nmax, "rheo:ct"); + memory->grow(cgradt, atom->nmax, 3, "rheo:cgradt"); + memory->grow(wsame, atom->nmax, "rheo:wsame"); + } + nmax_store = atom->nmax; } @@ -224,7 +244,17 @@ void ComputeRHEOVShift::compute_peratom() } } - if (newton_pair) comm->reverse_comm(this); + comm_stage = 0; + if (newton_pair) comm->reverse_comm(this, 3); + + // Zero any excluded types + + for (i = 0; i < nlocal; i++) + if (!shift_type[type[i]]) + for (a = 0; a < dim; a++) + vshift[i][a] = 0.0; + + if (multiphase_flag) correct_interfaces(); } /* ---------------------------------------------------------------------- */ @@ -246,7 +276,6 @@ void ComputeRHEOVShift::correct_surfaces() if (status[i] & PHASECHECK) continue; - //if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) { if (status[i] & STATUS_SURFACE) { nx = nsurface[i][0]; ny = nsurface[i][1]; @@ -283,16 +312,277 @@ void ComputeRHEOVShift::correct_surfaces() /* ---------------------------------------------------------------------- */ -int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) +void ComputeRHEOVShift::correct_interfaces() { - int m, last; + int i, j, a, ii, jj, jnum, itype, jtype; + int fluidi, fluidj; + double xtmp, ytmp, ztmp, rsq, r, rinv; + double w, wp, dr, w0, prefactor; + double imass, jmass, voli, volj, rhoi, rhoj; + double dx[3]; + int dim = domain->dimension; + + int *jlist; + int inum, *ilist, *numneigh, **firstneigh; + + int *type = atom->type; + int *status = atom->rheo_status; + int *mask = atom->mask; + double **x = atom->x; + double *rho = atom->rho; + double *mass = atom->mass; + double *rmass = atom->rmass; + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + size_t nbytes = nmax_store * sizeof(double); + memset(&ct[0], 0, nbytes); + memset(&wsame[0], 0, nbytes); + memset(&cgradt[0][0], 0, 3 * nbytes); + double ctmp, *dWij, *dWji; + + // Calculate color gradient + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + fluidi = !(status[i] & PHASECHECK); + jlist = firstneigh[i]; + jnum = numneigh[i]; + if (rmass) + imass = rmass[i]; + else + imass = mass[itype]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + dx[0] = xtmp - x[j][0]; + dx[1] = ytmp - x[j][1]; + dx[2] = ztmp - x[j][2]; + + rsq = lensq3(dx); + + if (rsq > cutsq) continue; + + fluidj = !(status[j] & PHASECHECK); + jtype = type[j]; + if (rmass) + jmass = rmass[j]; + else + jmass = mass[jtype]; + r = sqrt(rsq); + + rhoi = rho[i]; + rhoj = rho[j]; + + // Add corrections for walls + if (interface_flag) { + if (fluidi && (!fluidj)) { + rhoj = compute_interface->correct_rho(j); + } else if ((!fluidi) && fluidj) { + rhoi = compute_interface->correct_rho(i); + } else if ((!fluidi) && (!fluidj)) { + rhoi = rho0[itype]; + rhoj = rho0[jtype]; + } + } + + voli = imass / rhoi; + volj = jmass / rhoj; + + w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r); + + if (itype != jtype) ctmp = 1; + else ctmp = 0; + + ct[i] += ctmp * volj * w; + if (newton_pair || j < nlocal) + ct[j] += ctmp * voli * w; + } + } + + comm_stage = 1; + if (newton_pair) comm->reverse_comm(this, 1); + + // Calculate color gradient + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + fluidi = !(status[i] & PHASECHECK); + jlist = firstneigh[i]; + jnum = numneigh[i]; + imass = mass[itype]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + dx[0] = xtmp - x[j][0]; + dx[1] = ytmp - x[j][1]; + dx[2] = ztmp - x[j][2]; + rsq = lensq3(dx); + + if (rsq > cutsq) continue; + + fluidj = !(status[j] & PHASECHECK); + jtype = type[j]; + if (rmass) + jmass = rmass[j]; + else + jmass = mass[jtype]; + r = sqrt(rsq); + + rhoi = rho[i]; + rhoj = rho[j]; + + // Add corrections for walls + if (interface_flag) { + if (fluidi && (!fluidj)) { + rhoj = compute_interface->correct_rho(j); + } else if ((!fluidi) && fluidj) { + rhoi = compute_interface->correct_rho(i); + } else if ((!fluidi) && (!fluidj)) { + rhoi = rho0[itype]; + rhoj = rho0[jtype]; + } + } + + voli = imass / rhoi; + volj = jmass / rhoj; + + w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r); + dWij = compute_kernel->dWij; + dWji = compute_kernel->dWji; + + if (itype != jtype) ctmp = 1; + else ctmp = 0; + + for (a = 0; a < dim; a++) { + cgradt[i][a] -= ctmp * volj * dWij[a]; + if (newton_pair || j < nlocal) + cgradt[j][a] -= ctmp * voli * dWji[a]; + } + + if (itype == jtype) { + wsame[i] += w * r; + if (newton_pair || j < nlocal) + wsame[j] += w * r; + } + } + } + + comm_stage = 2; + if (newton_pair) comm->reverse_comm(this, 4); + comm->forward_comm(this, 1); + + // Correct shifting at fluid-fluid interface + // remove normal shifting component for interfacial particles + // Based on Yang, Rakhsha, Hu, & Negrut 2022 + + double *vs, ntmp[3]; + double minv, dot; + + for (i = 0; i < nlocal; i++) { + + // If isolated, just don't shift + if (wsame[i] < wmin) { + for (a = 0; a < dim; a++) + vshift[i][a] = 0.0; + continue; + } + + if (ct[i] < cmin) continue; + + minv = cgradt[i][0] * cgradt[i][0] + cgradt[i][1] * cgradt[i][1]; + if (dim == 3) minv += cgradt[i][2] * cgradt[i][2]; + + if (minv != 0) + minv = 1 / sqrt(minv); + + for (a = 0; a < dim; a++) + ntmp[a] = cgradt[i][a] * minv; + + vs = vshift[i]; + dot = ntmp[0] * vs[0] + ntmp[1] * vs[1]; + if (dim == 3) + dot += ntmp[2] * vs[2]; + + // To allowing shifting into the bulk + // if (dot > 0.0) continue; + + vshift[i][0] -= (1.0 - scale) * ntmp[0] * dot; + vshift[i][1] -= (1.0 - scale) * ntmp[1] * dot; + if (dim == 3) + vshift[i][2] -= (1.0 - scale) * ntmp[2] * dot; + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOVShift::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + m = 0; + + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = wsame[j]; + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOVShift::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; m = 0; last = first + n; - for (int i = first; i < last; i++) { - buf[m++] = vshift[i][0]; - buf[m++] = vshift[i][1]; - buf[m++] = vshift[i][2]; + for (i = first; i < last; i++) + wsame[i] = buf[m++]; +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, a, last; + + m = 0; + last = first + n; + if (comm_stage == 0) { + for (i = first; i < last; i++) { + buf[m++] = vshift[i][0]; + buf[m++] = vshift[i][1]; + buf[m++] = vshift[i][2]; + } + } else if (comm_stage == 1) { + for (i = first; i < last; i++) + buf[m++] = ct[i]; + } else { + for (i = first; i < last; i++) { + for (a = 0; a < 3; a++) + buf[m++] = cgradt[i][a]; + buf[m++] = wsame[i]; + } } return m; } @@ -301,14 +591,28 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) { - int i, j, m; + int i, j, a, m; m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - vshift[j][0] += buf[m++]; - vshift[j][1] += buf[m++]; - vshift[j][2] += buf[m++]; + if (comm_stage == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + vshift[j][0] += buf[m++]; + vshift[j][1] += buf[m++]; + vshift[j][2] += buf[m++]; + } + } else if (comm_stage == 1) { + for (i = 0; i < n; i++) { + j = list[i]; + ct[j] += buf[m++]; + } + } else { + for (i = 0; i < n; i++) { + j = list[i]; + for (a = 0; a < 3; a++) + cgradt[j][a] += buf[m++]; + wsame[j] += buf[m++]; + } } } @@ -319,5 +623,9 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) double ComputeRHEOVShift::memory_usage() { double bytes = 3 * nmax_store * sizeof(double); + + if (multiphase_flag) + bytes += 5 * nmax_store * sizeof(double); + return bytes; } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 485c6525f3..5d763ab008 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -31,19 +31,25 @@ class ComputeRHEOVShift : public Compute { void init() override; void init_list(int, class NeighList *) override; void compute_peratom() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; int pack_reverse_comm(int, int, double *) override; void unpack_reverse_comm(int, int *, double *) override; double memory_usage() override; void correct_surfaces(); + void correct_interfaces(); double **vshift; class FixRHEO *fix_rheo; private: - int nmax_store; + int nmax_store, comm_stage; double dtv, cut, cutsq, cutthird; - int surface_flag, interface_flag; + double scale, wmin, cmin; + int surface_flag, interface_flag, multiphase_flag; double *rho0; + double *wsame, *ct, **cgradt; + int *shift_type; class NeighList *list; class ComputeRHEOInterface *compute_interface; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index a978333aba..2241dbc585 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -70,15 +70,17 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : surface_flag = 0; oxidation_flag = 0; self_mass_flag = 0; - coordination_flag = 0; + shift_multiphase_flag = 0; int i; int n = atom->ntypes; memory->create(rho0, n + 1, "rheo:rho0"); memory->create(csq, n + 1, "rheo:csq"); + memory->create(shift_type, n + 1, "rheo:shift_type"); for (i = 1; i <= n; i++) { rho0[i] = 1.0; csq[i] = 1.0; + shift_type[i] = 1; } if (igroup != 0) error->all(FLERR, "fix rheo command requires group all"); @@ -112,6 +114,19 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg], "shift") == 0) { shift_flag = 1; + } else if (strcmp(arg[iarg], "shift/multiphase/scale") == 0) { + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/multiphase/scale", error); + shift_multiphase_flag = 1; + shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + shift_wmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg], "shift/type") == 0) { + if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/type", error); + for (i = 1; i <= n; i++) { + shift_type[i] = utils::logical(FLERR, arg[iarg + i], false, lmp); + } + iarg += n; } else if (strcmp(arg[iarg], "thermal") == 0) { thermal_flag = 1; } else if (strcmp(arg[iarg], "surface/detection") == 0) { @@ -128,7 +143,6 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : } else { error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]); } - iarg += 3; } else if (strcmp(arg[iarg], "interface/reconstruct") == 0) { interface_flag = 1; @@ -153,6 +167,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : iarg += 1; } + if ((!shift_flag) && shift_multiphase_flag) + error->all(FLERR, "Cannot use shift/multiphase/scale without shifting"); + if (self_mass_flag && (!rhosum_flag)) error->all(FLERR, "Cannot use self/mass setting without rho/sum"); @@ -174,6 +191,7 @@ FixRHEO::~FixRHEO() memory->destroy(csq); memory->destroy(rho0); + memory->destroy(shift_type); } /* ---------------------------------------------------------------------- @@ -376,7 +394,6 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Shifting atoms if (shift_flag) { for (i = 0; i < nlocal; i++) { - if (status[i] & STATUS_NO_SHIFT) continue; if (status[i] & PHASECHECK) continue; @@ -399,39 +416,34 @@ void FixRHEO::initial_integrate(int /*vflag*/) void FixRHEO::pre_force(int /*vflag*/) { - if (coordination_flag) - compute_kernel->compute_coordination(); + compute_kernel->compute_coordination(); // Needed for rho sum - if (rhosum_flag) - compute_rhosum->compute_peratom(); + if (rhosum_flag) compute_rhosum->compute_peratom(); compute_kernel->compute_peratom(); - // Note on first setup, have no forces for pressure to reference - if (interface_flag) + if (interface_flag) { + // Note on first setup, have no forces for pressure to reference compute_interface->compute_peratom(); + } // No need to forward v, rho, or T for compute_grad since already done compute_grad->compute_peratom(); compute_grad->forward_gradients(); - // Depends on NO_SHIFT status - if (shift_flag) - compute_vshift->compute_peratom(); + if (shift_flag) compute_vshift->compute_peratom(); // Remove temporary options int *mask = atom->mask; int *status = atom->rheo_status; int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) - if (mask[i] & groupbit) - status[i] &= OPTIONSMASK; + if (mask[i] & groupbit) status[i] &= OPTIONSMASK; // Calculate surfaces, update status if (surface_flag) { compute_surface->compute_peratom(); - if (shift_flag) - compute_vshift->correct_surfaces(); + if (shift_flag) compute_vshift->correct_surfaces(); } } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 4b75299f27..bb2d670f33 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -41,10 +41,12 @@ class FixRHEO : public Fix { // Model parameters double cut; double *rho0, *csq; + int *shift_type; int self_mass_flag; int zmin_kernel, zmin_surface, zmin_splash; int kernel_style, surface_style; double divr_surface; + double shift_scale, shift_wmin, shift_cmin; // Accessory fixes/computes int thermal_flag; @@ -53,7 +55,7 @@ class FixRHEO : public Fix { int interface_flag; int surface_flag; int oxidation_flag; - int coordination_flag; + int shift_multiphase_flag; int viscosity_fix_defined; int pressure_fix_defined; From bd9bcb3bf729da26644fb6085f2f02a356854ff9 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sun, 10 Nov 2024 14:27:22 -0700 Subject: [PATCH 13/27] adding missing rheo coordination flag --- src/RHEO/fix_rheo.cpp | 1 + src/RHEO/fix_rheo.h | 1 + 2 files changed, 2 insertions(+) diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 2241dbc585..806ef0873c 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -70,6 +70,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : surface_flag = 0; oxidation_flag = 0; self_mass_flag = 0; + coordination_flag = 0; shift_multiphase_flag = 0; int i; diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index bb2d670f33..acad9f0321 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -55,6 +55,7 @@ class FixRHEO : public Fix { int interface_flag; int surface_flag; int oxidation_flag; + int coordination_flag; int shift_multiphase_flag; int viscosity_fix_defined; From c0edafc0cfbe6cf50843128648d4d3f9cd47755f Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 11:31:16 -0700 Subject: [PATCH 14/27] Reorganizing optional subargs for fix rheo --- src/RHEO/compute_rheo_vshift.cpp | 10 ++--- src/RHEO/compute_rheo_vshift.h | 2 +- src/RHEO/fix_rheo.cpp | 63 +++++++++++++++++--------------- src/RHEO/fix_rheo.h | 24 ++++++------ 4 files changed, 53 insertions(+), 46 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 0047521e37..801b65fe3b 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -81,8 +81,8 @@ void ComputeRHEOVShift::init() cutsq = cut * cut; cutthird = cut / 3.0; - multiphase_flag = fix_rheo->shift_multiphase_flag; - if (multiphase_flag) { + cross_type_flag = fix_rheo->shift_cross_type_flag; + if (cross_type_flag) { scale = fix_rheo->shift_scale; wmin = fix_rheo->shift_wmin; cmin = fix_rheo->shift_cmin; @@ -134,7 +134,7 @@ void ComputeRHEOVShift::compute_peratom() if (nmax_store < atom->nmax) { memory->grow(vshift, atom->nmax, 3, "rheo:vshift"); - if (multiphase_flag) { + if (cross_type_flag) { memory->grow(ct, atom->nmax, "rheo:ct"); memory->grow(cgradt, atom->nmax, 3, "rheo:cgradt"); memory->grow(wsame, atom->nmax, "rheo:wsame"); @@ -254,7 +254,7 @@ void ComputeRHEOVShift::compute_peratom() for (a = 0; a < dim; a++) vshift[i][a] = 0.0; - if (multiphase_flag) correct_interfaces(); + if (cross_type_flag) correct_interfaces(); } /* ---------------------------------------------------------------------- */ @@ -624,7 +624,7 @@ double ComputeRHEOVShift::memory_usage() { double bytes = 3 * nmax_store * sizeof(double); - if (multiphase_flag) + if (cross_type_flag) bytes += 5 * nmax_store * sizeof(double); return bytes; diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 5d763ab008..04512c9740 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -46,7 +46,7 @@ class ComputeRHEOVShift : public Compute { int nmax_store, comm_stage; double dtv, cut, cutsq, cutthird; double scale, wmin, cmin; - int surface_flag, interface_flag, multiphase_flag; + int surface_flag, interface_flag, cross_type_flag; double *rho0; double *wsame, *ct, **cgradt; int *shift_type; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 806ef0873c..06ffa654af 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -52,9 +52,9 @@ static const char cite_rheo[] = /* ---------------------------------------------------------------------- */ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), - compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr), - compute_rhosum(nullptr), compute_vshift(nullptr) + Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), shift_type(nullptr), + compute_grad(nullptr), compute_kernel(nullptr), compute_interface(nullptr), + compute_surface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr) { time_integrate = 1; @@ -68,20 +68,18 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : shift_flag = 0; interface_flag = 0; surface_flag = 0; - oxidation_flag = 0; - self_mass_flag = 0; coordination_flag = 0; - shift_multiphase_flag = 0; + + rhosum_self_mass_flag = 0; + shift_cross_type_flag = 0; int i; int n = atom->ntypes; memory->create(rho0, n + 1, "rheo:rho0"); memory->create(csq, n + 1, "rheo:csq"); - memory->create(shift_type, n + 1, "rheo:shift_type"); for (i = 1; i <= n; i++) { rho0[i] = 1.0; csq[i] = 1.0; - shift_type[i] = 1; } if (igroup != 0) error->all(FLERR, "fix rheo command requires group all"); @@ -115,19 +113,26 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg], "shift") == 0) { shift_flag = 1; - } else if (strcmp(arg[iarg], "shift/multiphase/scale") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/multiphase/scale", error); - shift_multiphase_flag = 1; - shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - shift_wmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg], "shift/type") == 0) { - if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/type", error); - for (i = 1; i <= n; i++) { - shift_type[i] = utils::logical(FLERR, arg[iarg + i], false, lmp); + memory->create(shift_type, n + 1, "rheo:shift_type"); + for (i = 1; i <= n; i++) shift_type[i] = 1; + while (iarg < narg) { // optional sub-arguments + if (strcmp(arg[iarg], "scale/cross/type") == 0) { + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error); + shift_cross_type_flag = 1; + shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + shift_wmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg], "exclude/type") == 0) { + if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); + for (i = 1; i <= n; i++) + shift_type[i] = utils::logical(FLERR, arg[iarg + i], false, lmp); + iarg += n; + } else { + break; + } + iarg += 1; } - iarg += n; } else if (strcmp(arg[iarg], "thermal") == 0) { thermal_flag = 1; } else if (strcmp(arg[iarg], "surface/detection") == 0) { @@ -149,8 +154,14 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : interface_flag = 1; } else if (strcmp(arg[iarg], "rho/sum") == 0) { rhosum_flag = 1; - } else if (strcmp(arg[iarg], "self/mass") == 0) { - self_mass_flag = 1; + while (iarg < narg) { // optional sub-arguments + if (strcmp(arg[iarg], "self/mass") == 0) { + rhosum_self_mass_flag = 1; + } else { + break; + } + iarg += 1; + } } else if (strcmp(arg[iarg], "density") == 0) { if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error); for (i = 1; i <= n; i++) rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp); @@ -168,12 +179,6 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : iarg += 1; } - if ((!shift_flag) && shift_multiphase_flag) - error->all(FLERR, "Cannot use shift/multiphase/scale without shifting"); - - if (self_mass_flag && (!rhosum_flag)) - error->all(FLERR, "Cannot use self/mass setting without rho/sum"); - #if 0 if (lmp->citeme) lmp->citeme->add(cite_rheo); #endif @@ -212,7 +217,7 @@ void FixRHEO::post_constructor() if (rhosum_flag) { compute_rhosum = dynamic_cast( - modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag))); + modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", rhosum_self_mass_flag))); compute_rhosum->fix_rheo = this; } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index acad9f0321..aa80bc938e 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -41,23 +41,25 @@ class FixRHEO : public Fix { // Model parameters double cut; double *rho0, *csq; - int *shift_type; - int self_mass_flag; int zmin_kernel, zmin_surface, zmin_splash; int kernel_style, surface_style; double divr_surface; + + // Settings flags + int thermal_flag; + int rhosum_flag; + int interface_flag; + int surface_flag; + int shift_flag; + int coordination_flag; + + // Optional sub-flags/parameters + int rhosum_self_mass_flag; + int *shift_type; + int shift_cross_type_flag; double shift_scale, shift_wmin, shift_cmin; // Accessory fixes/computes - int thermal_flag; - int rhosum_flag; - int shift_flag; - int interface_flag; - int surface_flag; - int oxidation_flag; - int coordination_flag; - int shift_multiphase_flag; - int viscosity_fix_defined; int pressure_fix_defined; int thermal_fix_defined; From 59af94b8e36c1010516640a858e5ade51f2c3ca0 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 11:42:50 -0700 Subject: [PATCH 15/27] Updating rheo citation information --- doc/src/Howto_rheo.rst | 9 +++++++-- doc/src/fix_rheo.rst | 17 ++++++++++++++--- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index 4d5734271a..2407fd484d 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -15,8 +15,9 @@ details of the system, or develop new capabilities. For instance, the numerics associated with calculating gradients, reproducing kernels, etc. are separated into distinct classes to simplify the development of new integration schemes which can call these calculations. Additional numerical details can be found in -:ref:`(Clemmer) `. Example movies illustrating some of these -capabilities are found at https://www.lammps.org/movies.html#rheopackage. +:ref:`(Palermo) ` and :ref:`(Clemmer) `. +Example movies illustrating some of these capabilities are found at +https://www.lammps.org/movies.html#rheopackage. Note, if you simply want to run a traditional SPH simulation, the :ref:`SPH package ` package is likely better suited for your application. It has fewer advanced @@ -107,6 +108,10 @@ criteria for creating/deleting a bond or altering force calculations). ---------- +.. _howto_rheo_palermo: + +**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024). + .. _howto_rheo_clemmer: **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index eb88ef0536..617aa0eac5 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -46,8 +46,10 @@ Description .. versionadded:: 29Aug2024 Perform time integration for RHEO particles, updating positions, velocities, -and densities. For an overview of other features available in the RHEO package, -see :doc:`the RHEO howto `. +and densities. For a detailed breakdown of the integration timestep and +numerical details, see :ref:`(Palermo) `. For an overview +and list of other features available in the RHEO package, see +:doc:`the RHEO howto `. The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four kernels are currently available. The *quintic* kernel is a standard quintic @@ -70,7 +72,12 @@ and velocity of solid particles are alternatively reconstructed for every fluid-solid interaction to ensure no-slip and pressure-balanced boundaries. This is done by estimating the location of the fluid-solid interface and extrapolating fluid particle properties across the interface to calculate a -temporary apparent density and velocity for a solid particle. +temporary apparent density and velocity for a solid particle. The numerical +details are the same as those described in +:ref:`(Palermo) ` except there is an additional +restriction that the reconstructed solid density cannot be less than the +equilibrium density. This prevents fluid particles from sticking to solid +surfaces. A modified form of Fickian particle shifting can be enabled with the *shift* keyword. This effectively shifts particle positions to generate a @@ -163,6 +170,10 @@ Default ---------- +.. _rheo_palermo: + +**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024). + .. _fix_rheo_hu: **(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006). From 15767903b26ed37433c17f4b183a1c835e8645df Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 11:51:59 -0700 Subject: [PATCH 16/27] testing new argument description for fix rheo --- doc/src/fix_rheo.rst | 38 +++++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 617aa0eac5..4d3b30967a 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -20,17 +20,25 @@ Syntax .. parsed-literal:: - *thermal* values = none, turns on thermal evolution - *interface/reconstruct* values = none, reconstructs interfaces with solid particles - *surface/detection* values = *sdstyle* *limit* *limit/splash* - *sdstyle* = *coordination* or *divergence* - *limit* = threshold for surface particles - *limit/splash* = threshold for splash particles - *shift* values = none, turns on velocity shifting - *rho/sum* values = none, uses the kernel to compute the density of particles - *self/mass* values = none, a particle uses its own mass in a rho summation - *density* values = *rho01*, ... *rho0N* (density) - *speed/sound* values = *cs0*, ... *csN* (velocity) + *thermal* turns on thermal evolution + values = none + *interface/reconstruct* reconstructs interfaces with solid particles + values = none + *surface/detection* detects free-surfaces with an absence of particles + values = *sdstyle* *limit* *limit/splash* + *sdstyle* = *coordination* or *divergence* + *limit* = threshold for surface particles + *limit/splash* = threshold for splash particles + *shift* turns on velocity shifting + values = none + *rho/sum* density evolution performed by a kernel summation + values = none + optional args = *self/mass* + *self/mass* values = none, a particle uses its own mass in summation + *density* specify equilibrium densities for each atom type + values = *rho01*, ... *rho0N* (density) + *speed/sound* specify speeds of sound for each atom type + values = *cs0*, ... *csN* (velocity) Examples """""""" @@ -39,6 +47,7 @@ Examples fix 1 all rheo 3.0 quintic 0 thermal density 0.1 0.1 speed/sound 10.0 1.0 fix 1 all rheo 3.0 RK1 10 shift surface/detection coordination 40 + fix 1 all rheo 3.0 RK1 10 rhosum self/mass Description """"""""""" @@ -108,10 +117,9 @@ threshold for this classification is set by the numerical value of By default, RHEO integrates particles' densities using a mass diffusion equation. Alternatively, one can update densities every timestep by performing a kernel summation of the masses of neighboring particles by specifying the *rho/sum* -keyword. - -The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*. -Typically, the density :math:`\rho` of a particle is calculated as the sum over neighbors +keyword. The optional *self/mass* subargument then modifies the behavior of the +density summation. Typically, the density :math:`\rho` of a particle is calculated +as the sum over neighbors .. math:: \rho_i = \sum_{j} W_{ij} M_j From 0e90de7d95ec0517dd93900f8f247e0206e2052d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 11:52:21 -0700 Subject: [PATCH 17/27] Updaing rheo citeme --- src/RHEO/fix_rheo.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 06ffa654af..4c2ee735f0 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -38,16 +38,18 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; -#if 0 -// publication was removed from documentation static const char cite_rheo[] = - "@article{PalermoInPrep,\n" - " journal = {in prep},\n" - " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n" + "@article{Palermo2024,\n" + " journal = {Physics of Fluids},\n" + " title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\n" + " volume = {36},\n" + " number = {11},\n" + " pages = {113337},\n" " year = {2024},\n" - " author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n" + " issn = {1070-6631},\n" + " doi = {https://doi.org/10.1063/5.0228823},\n" + " author = {Palermo, Eric T. and Wolf, Ki T. and Clemmer, Joel T. and O'Connor, Thomas C.},\n" "}\n\n"; -#endif /* ---------------------------------------------------------------------- */ @@ -179,9 +181,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : iarg += 1; } -#if 0 if (lmp->citeme) lmp->citeme->add(cite_rheo); -#endif } /* ---------------------------------------------------------------------- */ From 8521ff05873e7db0bb485341415a9471d05ebc3b Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 15:20:54 -0700 Subject: [PATCH 18/27] Simplifying logic for shifting two types --- src/RHEO/compute_rheo_vshift.cpp | 31 ++++++++++++++++--------------- src/RHEO/compute_rheo_vshift.h | 2 +- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 801b65fe3b..248e857fb3 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -254,7 +254,7 @@ void ComputeRHEOVShift::compute_peratom() for (a = 0; a < dim; a++) vshift[i][a] = 0.0; - if (cross_type_flag) correct_interfaces(); + if (cross_type_flag) correct_type_interface(); } /* ---------------------------------------------------------------------- */ @@ -312,7 +312,7 @@ void ComputeRHEOVShift::correct_surfaces() /* ---------------------------------------------------------------------- */ -void ComputeRHEOVShift::correct_interfaces() +void ComputeRHEOVShift::correct_type_interface() { int i, j, a, ii, jj, jnum, itype, jtype; int fluidi, fluidj; @@ -417,6 +417,10 @@ void ComputeRHEOVShift::correct_interfaces() if (newton_pair) comm->reverse_comm(this, 1); // Calculate color gradient + // Note: in future might want to generalize this so color function can be used + // by other calculations (e.g. surface tension) + // maybe can create custom "calc_grad" method that takes an arbitrary field + // in ComputeRHEOGrad? for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -495,8 +499,7 @@ void ComputeRHEOVShift::correct_interfaces() // remove normal shifting component for interfacial particles // Based on Yang, Rakhsha, Hu, & Negrut 2022 - double *vs, ntmp[3]; - double minv, dot; + double ntmp[3], minv, dot; for (i = 0; i < nlocal; i++) { @@ -509,8 +512,9 @@ void ComputeRHEOVShift::correct_interfaces() if (ct[i] < cmin) continue; - minv = cgradt[i][0] * cgradt[i][0] + cgradt[i][1] * cgradt[i][1]; - if (dim == 3) minv += cgradt[i][2] * cgradt[i][2]; + minv = 0; + for (a = 0; a < dim; a++) + minv += cgradt[i][a] * cgradt[i][a]; if (minv != 0) minv = 1 / sqrt(minv); @@ -518,18 +522,15 @@ void ComputeRHEOVShift::correct_interfaces() for (a = 0; a < dim; a++) ntmp[a] = cgradt[i][a] * minv; - vs = vshift[i]; - dot = ntmp[0] * vs[0] + ntmp[1] * vs[1]; - if (dim == 3) - dot += ntmp[2] * vs[2]; + dot = 0.0; + for (a = 0; a < dim; a++) + dot += ntmp[a] * vshift[i][a]; - // To allowing shifting into the bulk + // To allowing shifting into the same phase bulk // if (dot > 0.0) continue; - vshift[i][0] -= (1.0 - scale) * ntmp[0] * dot; - vshift[i][1] -= (1.0 - scale) * ntmp[1] * dot; - if (dim == 3) - vshift[i][2] -= (1.0 - scale) * ntmp[2] * dot; + for (a = 0; a < dim; a++) + vshift[i][a] -= (1.0 - scale) * ntmp[a] * dot; } } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 04512c9740..b76d48d5e3 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -37,7 +37,7 @@ class ComputeRHEOVShift : public Compute { void unpack_reverse_comm(int, int *, double *) override; double memory_usage() override; void correct_surfaces(); - void correct_interfaces(); + void correct_type_interface(); double **vshift; class FixRHEO *fix_rheo; From 9e0047ff06ff282dc79c3484c1ba5e15361f0060 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 15:50:01 -0700 Subject: [PATCH 19/27] Documenting new multi-type shifting options --- doc/src/fix_rheo.rst | 61 ++++++++++++++++++++++++++++++++++++------- src/RHEO/fix_rheo.cpp | 14 +++++----- 2 files changed, 58 insertions(+), 17 deletions(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 4d3b30967a..089bfd29a0 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -28,9 +28,16 @@ Syntax values = *sdstyle* *limit* *limit/splash* *sdstyle* = *coordination* or *divergence* *limit* = threshold for surface particles - *limit/splash* = threshold for splash particles + *limit/splash* = threshold for splash particles (unitless) *shift* turns on velocity shifting values = none + optional args = *exclude/type* or *scale/cross/type* + *exclude/type* values = *types* + *types* = list of types + *scale/cross/type* values = *shiftscale* *cmin* *wmin* + *shiftscale* = fraction of shifting in normal direction to preserve (unitless) + *cmin* = minimum color function value required for scaling (unitless) + *wmin* = minimum local same-type support required for any shifting (unitless) *rho/sum* density evolution performed by a kernel summation values = none optional args = *self/mass* @@ -83,19 +90,49 @@ This is done by estimating the location of the fluid-solid interface and extrapolating fluid particle properties across the interface to calculate a temporary apparent density and velocity for a solid particle. The numerical details are the same as those described in -:ref:`(Palermo) ` except there is an additional +:ref:`(Palermo) ` except there is an additional restriction that the reconstructed solid density cannot be less than the equilibrium density. This prevents fluid particles from sticking to solid surfaces. A modified form of Fickian particle shifting can be enabled with the *shift* keyword. This effectively shifts particle positions to generate a -more uniform spatial distribution. Shifting currently does not consider the +more uniform spatial distribution. By default, shifting does not consider the type of a particle and therefore may be inappropriate in systems consisting -of multiple fluid phases. +of multiple atom types representing multiple fluid phases. However, two +optional subarguments can follow the *shift* keyword, *exclude/type* and +*scale/cross/type*. -In systems with free surfaces, the *surface/detection* keyword can be used -to classify the location of particles as being within the bulk fluid, on a +The *exclude/type* option lets the user specify a list of atom types which +are not shifted, *types*. A wild-card asterisk can be used in place +of or in conjunction with the *types* argument to set the coefficients for +multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with +no numeric values means all types from 1 to :math:`N`. A leading asterisk +means all types from 1 to n (inclusive). A trailing asterisk means all types +from m to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). + +The *scale/cross/type* option is designed to handle interfaces between fluids +made up of different atom types. Similar to the method by +:ref:`(Yang) `, a color function is calculated and used to +estimate a local interfacial normal vector. Shifting along this normal direction +is rescaled by a factor of *scaleshift*, such that a value of *scaleshift* of +zero implies there is no shifting in the normal direction and a value of +*scaleshift* of one implies no change in behavior. This scaling is only applied +to atoms with a color function value greater than *cmin*. To handle scenarios +of a small inclusion of one fluid type (e.g. a single atom) inside another, +the amount of same-type support is calculated + +.. math:: + W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij} + +where :math:`\delta_{ij}` is zero if atoms :math:`i` and :math:`j` have different +types but unity otherwise. If :math:`W_{i,\mathrm{same}}` is ever less than the +specified value of *wmin*, shifting is turned off for particle :math:`i` + +In systems with free surfaces (atom-vacuum), the *surface/detection* keyword +can classify the location of particles as being within the bulk fluid, on a free surface, or isolated from other particles in a splash or droplet. Shifting is then disabled in the normal direction away from the free surface to prevent particles from diffusing away. Surface detection can also be used @@ -117,9 +154,9 @@ threshold for this classification is set by the numerical value of By default, RHEO integrates particles' densities using a mass diffusion equation. Alternatively, one can update densities every timestep by performing a kernel summation of the masses of neighboring particles by specifying the *rho/sum* -keyword. The optional *self/mass* subargument then modifies the behavior of the -density summation. Typically, the density :math:`\rho` of a particle is calculated -as the sum over neighbors +keyword. Following this keyword, one may include the optional *self/mass* subargument +which modifies the behavior of the density summation. Typically, the density +:math:`\rho` of a particle is calculated as the sum over neighbors .. math:: \rho_i = \sum_{j} W_{ij} M_j @@ -182,6 +219,10 @@ Default **(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024). +.. _rheo_yang: + +**(Yang)** Yang, Rakhsha, Hu, Negrut, J. Comp. Physics, 458, 111079 (2022). + .. _fix_rheo_hu: -**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006). +**(Hu)** Hu, and Adams, J. Comp. Physics, 213, 844-861 (2006). diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 4c2ee735f0..56738c7a0c 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -75,7 +75,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : rhosum_self_mass_flag = 0; shift_cross_type_flag = 0; - int i; + int i, nlo, nhi; int n = atom->ntypes; memory->create(rho0, n + 1, "rheo:rho0"); memory->create(csq, n + 1, "rheo:csq"); @@ -122,14 +122,14 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error); shift_cross_type_flag = 1; shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - shift_wmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + shift_cmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + shift_wmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg += 3; } else if (strcmp(arg[iarg], "exclude/type") == 0) { - if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); - for (i = 1; i <= n; i++) - shift_type[i] = utils::logical(FLERR, arg[iarg + i], false, lmp); - iarg += n; + if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); + utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error); + for (i = nlo; i <= nhi; i++) shift_type[i] = 0; + iarg += 1; } else { break; } From 33128323d8aa656e1ea07a1866839dd798ed3ca2 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 16:31:09 -0700 Subject: [PATCH 20/27] Adding shift example to rheo doc, clarifying speed of sound exception --- doc/src/fix_rheo.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 089bfd29a0..6ec19ab6b5 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -54,6 +54,7 @@ Examples fix 1 all rheo 3.0 quintic 0 thermal density 0.1 0.1 speed/sound 10.0 1.0 fix 1 all rheo 3.0 RK1 10 shift surface/detection coordination 40 + fix 1 all rheo 3.0 RK1 10 shift exclude/type 2*4 scale/cross/type 0.05 0.02 0.5 fix 1 all rheo 3.0 RK1 10 rhosum self/mass Description @@ -172,7 +173,9 @@ equilibrium density *rho0*. The *speed/sound* keyword is used to specify the speed of sound of each of the N particle types. It must be followed by N numerical values specifying each type's -speed of sound *cs*. +speed of sound *cs*. These values may be ignored if the pressure equation of +state has a non-constant speed of sound, as discussed further in +:doc:`fix rheo/pressure `. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From cb5bad7ece106a99bf7d49a3750bed1ec5403502 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 Nov 2024 21:22:48 -0700 Subject: [PATCH 21/27] Dropping obsolete keyword --- doc/src/fix_rheo.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 6ec19ab6b5..24f7478566 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -16,7 +16,7 @@ Syntax * kstyle = *quintic* or *RK0* or *RK1* or *RK2* * zmin = minimal number of neighbors for reproducing kernels * zero or more keyword/value pairs may be appended to args -* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *self/mass* or *speed/sound* +* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *speed/sound* .. parsed-literal:: From 10c429fe217dbc83fc2cf05831006bc2e4fa4e9a Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 22 Nov 2024 14:37:47 -0700 Subject: [PATCH 22/27] Separating background pressure from EoS definition --- doc/src/fix_rheo_pressure.rst | 16 ++++++++++++---- src/RHEO/fix_rheo_pressure.cpp | 30 +++++++++++++++++++++++------- 2 files changed, 35 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index 51f00b6684..d06fadfa1c 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -14,15 +14,16 @@ Syntax * rheo/pressure = style name of this fix command * one or more types and pressure styles must be appended * types = lists of types (see below) -* pstyle = *linear* or *taitwater* or *cubic* +* pstyle = *linear* or *tait/water* or *tait/general* or *cubic* or *ideal/gas* or *background* .. parsed-literal:: *linear* args = none *tait/water* args = none - *tait/general* args = exponent :math:`gamma` (unitless) and background pressure :math:`P[b]` (pressure) + *tait/general* args = exponent :math:`gamma` (unitless) *cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2) *ideal/gas* args = heat capacity ratio :math:`gamma` (unitless) + *background* args = background pressure :math:`P[b]` (pressure) Examples """""""" @@ -31,6 +32,7 @@ Examples fix 1 all rheo/pressure * linear fix 1 all rheo/pressure 1 linear 2 cubic 10.0 + fix 1 all rheo/pressure * linear * background 0.1 Description """"""""""" @@ -77,9 +79,9 @@ Style *tait/general* generalizes this equation of state .. math:: - P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr] + P[b]. + P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr]. -where :math:`\gamma` is an exponent and :math:`P[b]` is a constant background pressure. +where :math:`\gamma` is an exponent. Style *ideal/gas* is the ideal gas equation of state @@ -92,6 +94,12 @@ a particle per unit mass. This style is only compatible with atom style rheo/the Note that when using this style, the speed of sound is no longer constant such that the value of :math:`c` specified in :doc:`fix rheo ` is not used. +The *background* style acts differently than the rest as it +only adds a constant background pressure shift :math:`P[b]` +to all atoms of the designated types. Therefore, this style +must be used in conjunction with another style that specifies +an equation of state. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index d4f66482b3..9364d08e47 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -47,6 +47,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : comm_forward = 1; variable_csq = 0; invertable_pressure = 1; + background_flag = 0; // Currently can only have one instance of fix rheo/pressure if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all"); @@ -59,7 +60,10 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : memory->create(pbackground, n + 1, "rheo:pbackground"); memory->create(gamma, n + 1, "rheo:gamma"); - for (i = 1; i <= n; i++) pressure_style[i] = NONE; + for (i = 1; i <= n; i++) { + pressure_style[i] = NONE; + pbackground[i] = 0.0; + } int iarg = 3; while (iarg < narg) { @@ -72,16 +76,14 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg + 1], "tait/water") == 0) { for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; } else if (strcmp(arg[iarg + 1], "tait/general") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait/general", error); + if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait/general", error); double tpower_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - double pbackground_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 2; + iarg += 1; for (i = nlo; i <= nhi; i++) { pressure_style[i] = TAITGENERAL; tpower[i] = tpower_one; - pbackground[i] = pbackground_one; } } else if (strcmp(arg[iarg + 1], "cubic") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure cubic", error); @@ -109,6 +111,15 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : variable_csq = 1; if (atom->esph_flag != 1) error->all(FLERR, "fix rheo/pressure ideal gas equation of state requires atom property esph"); + } else if (strcmp(arg[iarg + 1], "background") == 0) { + if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure background", error); + + double pbackground_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 1; + + for (i = nlo; i <= nhi; i++) + pbackground[i] = pbackground_one; + background_flag = 1; } else { error->all(FLERR, "Illegal fix command, {}", arg[iarg]); } @@ -236,10 +247,13 @@ double FixRHEOPressure::calc_pressure(double rho, int i) } else if (pressure_style[type] == TAITGENERAL) { rho_ratio = rho * rho0inv[type]; p = csq[type] * rho0[type] * (pow(rho_ratio, tpower[type]) - 1.0) / tpower[type]; - p += pbackground[type]; } else if (pressure_style[type] == IDEAL) { p = (gamma[type] - 1.0) * rho * atom->esph[i] / atom->mass[type]; } + + if (background_flag) + p += pbackground[type]; + return p; } @@ -250,6 +264,9 @@ double FixRHEOPressure::calc_rho(double p, int i) double rho = 0.0; int type = atom->type[i]; + if (background_flag) + p -= pbackground[type]; + if (pressure_style[type] == LINEAR) { rho = csqinv[type] * p + rho0[type]; } else if (pressure_style[type] == CUBIC) { @@ -260,7 +277,6 @@ double FixRHEOPressure::calc_rho(double p, int i) rho *= pow(rho0[type], 6.0 * SEVENTH); rho *= pow(csq[type], -SEVENTH); } else if (pressure_style[type] == TAITGENERAL) { - p -= pbackground[type]; rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]); rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]); rho *= pow(csq[type], -1.0 / tpower[type]); From e4683c3134041eb1ccc71a912e4705c9a80fd191 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 22 Nov 2024 14:40:45 -0700 Subject: [PATCH 23/27] typo in equation --- doc/src/fix_rheo_pressure.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index d06fadfa1c..5a732af647 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -58,7 +58,7 @@ equation of state with a particle pressure :math:`P` calculated as .. math:: - P = c (\rho - \rho_0) + P = c^2 (\rho - \rho_0) where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density, and :math:`\rho` is the current density of a particle. The numerical values of @@ -67,7 +67,7 @@ is a cubic equation of state which has an extra argument :math:`A_3`, .. math:: - P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . + P = c^2 ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) . Style *tait/water* is Tait's equation of state: From 8aba26bd71f09b4816ef5f3d86060e5f6951d24c Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 26 Nov 2024 15:55:08 -0700 Subject: [PATCH 24/27] doc fixes --- doc/src/fix_rheo.rst | 6 +++--- doc/src/fix_rheo_pressure.rst | 13 ++++++------- doc/src/fix_rheo_thermal.rst | 14 +++++++------- doc/src/fix_rheo_viscosity.rst | 13 ++++++------- 4 files changed, 22 insertions(+), 24 deletions(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 24f7478566..6341df1a1f 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -102,12 +102,12 @@ more uniform spatial distribution. By default, shifting does not consider the type of a particle and therefore may be inappropriate in systems consisting of multiple atom types representing multiple fluid phases. However, two optional subarguments can follow the *shift* keyword, *exclude/type* and -*scale/cross/type*. +*scale/cross/type* to adjust shifting at fluid interfaces. The *exclude/type* option lets the user specify a list of atom types which are not shifted, *types*. A wild-card asterisk can be used in place -of or in conjunction with the *types* argument to set the coefficients for -multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" +of or in conjunction with the *types* argument to toggle shifting for +multiple atom types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom types, then an asterisk with no numeric values means all types from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive). A trailing asterisk means all types diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index 5a732af647..7d5d08ddd1 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -44,13 +44,12 @@ define different equations of state for different atom types. An equation must be specified for every atom type. One first defines the atom *types*. A wild-card asterisk can be used in place -of or in conjunction with the *types* argument to set the coefficients for -multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" -or "m\*n". If :math:`N` is the number of atom types, then an asterisk with -no numeric values means all types from 1 to :math:`N`. A leading asterisk -means all types from 1 to n (inclusive). A trailing asterisk means all types -from m to :math:`N` (inclusive). A middle asterisk means all types from m to n -(inclusive). +of or in conjunction with the *types* argument to set values for multiple atom +types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is +the number of atom types, then an asterisk with no numeric values means all types +from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive). +A trailing asterisk means all types from m to :math:`N` (inclusive). A middle +asterisk means all types from m to n (inclusive). The *types* definition is followed by the pressure style, *pstyle*. Current options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear diff --git a/doc/src/fix_rheo_thermal.rst b/doc/src/fix_rheo_thermal.rst index 214bc1db86..bc0e8dbfd0 100644 --- a/doc/src/fix_rheo_thermal.rst +++ b/doc/src/fix_rheo_thermal.rst @@ -70,13 +70,13 @@ of the energy is used to shift energies. This may be inappropriate in systems with multiple atom types with different specific heats. For each property, one must first define a list of atom types. A wild-card -asterisk can be used in place of or in conjunction with the *types* argument -to set the coefficients for multiple pairs of atom types. This takes the -form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom -types, then an asterisk with no numeric values means all types from 1 to -:math:`N`. A leading asterisk means all types from 1 to n (inclusive). -A trailing asterisk means all types from m to :math:`N` (inclusive). A -middle asterisk means all types from m to n (inclusive). +asterisk can be used in place of or in conjunction with the *types* argument to +set values for multiple atom types. This takes the form "\*" or "\*n" or "m\*" +or "m\*n". If :math:`N` is the number of atom types, then an asterisk with no +numeric values means all types from 1 to :math:`N`. A leading asterisk means +all types from 1 to n (inclusive). A trailing asterisk means all types from m +to :math:`N` (inclusive). A middle asterisk means all types from m to n +(inclusive). The *types* definition for each property is followed by the style. Currently, the only option is *constant*. Style *constant* simply applies a constant value diff --git a/doc/src/fix_rheo_viscosity.rst b/doc/src/fix_rheo_viscosity.rst index 804059e6f8..4d72d9b209 100644 --- a/doc/src/fix_rheo_viscosity.rst +++ b/doc/src/fix_rheo_viscosity.rst @@ -45,13 +45,12 @@ viscosities for different atom types, but a viscosity must be specified for every atom type. One first defines the atom *types*. A wild-card asterisk can be used in place -of or in conjunction with the *types* argument to set the coefficients for -multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*" -or "m\*n". If :math:`N` is the number of atom types, then an asterisk with -no numeric values means all types from 1 to :math:`N`. A leading asterisk -means all types from 1 to n (inclusive). A trailing asterisk means all types -from m to :math:`N` (inclusive). A middle asterisk means all types from m to n -(inclusive). +of or in conjunction with the *types* argument to set values for multiple atom +types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is +the number of atom types, then an asterisk with no numeric values means all types +from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive). +A trailing asterisk means all types from m to :math:`N` (inclusive). A middle +asterisk means all types from m to n (inclusive). The *types* definition is followed by the viscosity style, *vstyle*. Two options are available, *constant* and *power*. Style *constant* simply From 09cd209c62ed6843f8dd41c96cacad19af2a3500 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 26 Nov 2024 15:56:03 -0700 Subject: [PATCH 25/27] typo --- src/RHEO/compute_rheo_interface.cpp | 2 +- src/RHEO/fix_rheo_pressure.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index b5a295aad6..bbc00bd8d8 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -97,7 +97,7 @@ void ComputeRHEOInterface::init() auto fixes = modify->get_fix_by_style("rheo/pressure"); fix_pressure = dynamic_cast(fixes[0]); - if (!fix_pressure->invertable_pressure) + if (!fix_pressure->invertible_pressure) error->all(FLERR, "RHEO interface reconstruction incompatible with pressure equation of state"); neighbor->add_request(this, NeighConst::REQ_DEFAULT); diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 9364d08e47..d972f19ffa 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -46,7 +46,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : comm_forward = 1; variable_csq = 0; - invertable_pressure = 1; + invertible_pressure = 1; background_flag = 0; // Currently can only have one instance of fix rheo/pressure @@ -96,7 +96,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : c_cubic[i] = c_cubic_one; } - invertable_pressure = 0; + invertible_pressure = 0; } else if (strcmp(arg[iarg + 1], "ideal/gas") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure ideal/gas", error); From e9073e4e1d20f3c73df9fe2f0abfea30cb24fe9d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 26 Nov 2024 16:19:19 -0700 Subject: [PATCH 26/27] Missed changes --- doc/src/fix_rheo.rst | 2 +- src/RHEO/fix_rheo_pressure.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 6341df1a1f..0214289dee 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -123,7 +123,7 @@ zero implies there is no shifting in the normal direction and a value of *scaleshift* of one implies no change in behavior. This scaling is only applied to atoms with a color function value greater than *cmin*. To handle scenarios of a small inclusion of one fluid type (e.g. a single atom) inside another, -the amount of same-type support is calculated +the degree of same-type support is calculated .. math:: W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij} diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index ac60d1644d..ea9aca3169 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -38,7 +38,7 @@ class FixRHEOPressure : public Fix { double calc_rho(double, int); double calc_csq(double, int); int variable_csq; - int invertable_pressure; + int invertible_pressure; private: double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground, *gamma; From 296941ede79a45569c74504e62eed45645d1157f Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 2 Dec 2024 10:05:12 -0700 Subject: [PATCH 27/27] Missed file, whitespace --- src/BPM/bond_bpm_rotational.cpp | 2 +- src/RHEO/fix_rheo_pressure.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index cb5518a506..efb84f1709 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -543,7 +543,7 @@ void BondBPMRotational::compute(int eflag, int vflag) // ------------------------------------------------------// MathExtra::scale3(smooth, force1on2); - + if (newton_bond || i1 < nlocal) { f[i1][0] -= force1on2[0]; f[i1][1] -= force1on2[1]; diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index ea9aca3169..c15e9b4ddc 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -42,7 +42,7 @@ class FixRHEOPressure : public Fix { private: double *c_cubic, *csq, *csqinv, *rho0, *rho0inv, *tpower, *pbackground, *gamma; - int *pressure_style; + int *pressure_style, background_flag; class FixRHEO *fix_rheo; };