diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index b09925d6ac..939754aa95 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -17,6 +17,10 @@ these calculations. Additional numerical details can be found in :ref:`(Palermo) ` and :ref:`(Clemmer) `. +Note, if you simply want to run a traditional SPH simulation, the SPH package +is likely better suited for your application. It has fewer advanced features +and therefore benefits from improved performance. + ---------- At the core of the package is :doc:`fix rheo ` which integrates diff --git a/examples/rheo/balloon/in.rheo.balloon b/examples/rheo/balloon/in.rheo.balloon index d24573b0fb..c82fa39df3 100644 --- a/examples/rheo/balloon/in.rheo.balloon +++ b/examples/rheo/balloon/in.rheo.balloon @@ -60,6 +60,7 @@ fix 1 all rheo ${cut} quintic 0 & fix 2 all rheo/viscosity * constant ${eta} fix 3 all rheo/pressure * linear fix 4 all wall/harmonic ylo EDGE 2.0 1.0 1.0 +fix 5 all enforce2d compute rho all rheo/property/atom rho compute phase all rheo/property/atom phase diff --git a/examples/rheo/dam-break/in.rheo.dam.break b/examples/rheo/dam-break/in.rheo.dam.break index b59b710451..870b3529b6 100644 --- a/examples/rheo/dam-break/in.rheo.dam.break +++ b/examples/rheo/dam-break/in.rheo.dam.break @@ -3,24 +3,24 @@ dimension 2 units lj atom_style rheo -boundary f f p +boundary f s p comm_modify vel yes newton off # ------ Create simulation box ------ # variable n equal 1.0 -variable cut equal 3.0 -variable dx equal 4.0 +variable cut equal 2.2 +variable dx equal 3.0 -region box block -1 200 -1 80 -0.1 0.1 units box -#region box block -1 100 -1 40 -0.1 0.1 units box +region box block -1 150 -1 80 -0.1 0.1 units box create_box 2 box lattice hex ${n} -region fluid block $(xlo+v_dx) $(xlo+44.0) $(ylo+v_dx) $(yhi-16.0) EDGE EDGE units box -#region fluid block $(xlo+v_dx) $(xlo+20.0) $(ylo+v_dx) $(yhi-10.0) EDGE EDGE units box -region walls block $(xlo+v_dx) $(xhi-v_dx) $(ylo+v_dx) $(yhi-v_dx) EDGE EDGE units box side out +region fluid block $(xlo+v_dx+1.0) $(xlo+40.0) $(ylo+v_dx+1.0) $(yhi-20.0) EDGE EDGE units box +region walls1 block $(xlo+v_dx) $(xhi-v_dx) $(ylo+v_dx) $(yhi-v_dx) EDGE EDGE side out units box +region walls2 block EDGE EDGE EDGE $(yhi-v_dx) EDGE EDGE side in units box +region walls intersect 2 walls1 walls2 create_atoms 1 region fluid create_atoms 2 region walls @@ -33,12 +33,13 @@ group rig type 2 variable rho0 equal 1.0 variable mp equal ${rho0}/${n} variable cs equal 1.0 -variable zeta equal 0.05 +variable zeta equal 0.1 variable dt_max equal 0.1*${cut}/${cs}/3 -variable eta equal 0.05 -variable Dr equal 0.05 +variable eta equal 0.1 +variable Dr equal 0.1 -mass * ${mp} +mass 1 ${mp} +mass 2 $(2*v_mp) set group all rheo/rho ${rho0} set group all rheo/status 0 @@ -46,36 +47,30 @@ set group rig rheo/status 1 timestep ${dt_max} -pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr} +pair_style rheo ${cut} artificial/visc ${zeta} #rho/damp ${Dr} pair_coeff * * - # ------ Fixes & computes ------ # fix 1 all rheo ${cut} quintic 10 & - shift & surface/detection coordination 22 8 & - interface/reconstruct -#fix 1 all rheo ${cut} RK1 10 surface/detection coordination 36 10 interface/reconstruct shift + rho/sum fix 2 all rheo/viscosity * constant ${eta} fix 3 all rheo/pressure * linear -fix 4 all gravity 5e-4 vector 0 -1 0 +fix 4 all gravity 1e-3 vector 0 -1 0 fix 5 rig setforce 0.0 0.0 0.0 +fix 6 all enforce2d compute rho all rheo/property/atom rho -compute phase all rheo/property/atom phase compute p all rheo/property/atom pressure compute surf all rheo/property/atom surface compute sn all rheo/property/atom surface/n/x surface/n/y -compute vshift all rheo/property/atom shift/v/x shift/v/y - # ------ Output & Run ------ # -thermo 100 +thermo 20 thermo_style custom step time ke press -dump 1 all custom 20 atomDump id type x y vx vy fx fy c_rho c_phase c_surf c_p c_vshift[*] c_sn[*] +#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho c_surf c_p c_sn[*] -#thermo 1 -run 300000 # 400000 +run 30000 diff --git a/examples/rheo/ice-cubes/in.rheo.ice.cubes b/examples/rheo/ice-cubes/in.rheo.ice.cubes index 12eb54c32e..b023a8bac9 100644 --- a/examples/rheo/ice-cubes/in.rheo.ice.cubes +++ b/examples/rheo/ice-cubes/in.rheo.ice.cubes @@ -64,6 +64,7 @@ fix 4 all rheo/thermal conductivity * constant ${kappa} & fix 5 all wall/region wall harmonic 1.0 1.0 1.0 fix 6 all gravity 5e-4 vector 0 -1 0 fix 7 all deposit 8 0 1000 37241459 mol my_mol region drop near 2.0 vy -0.02 -0.02 +fix 8 all enforce2d compute rho all rheo/property/atom rho compute phase all rheo/property/atom phase diff --git a/examples/rheo/oxidation/in.rheo.oxidation b/examples/rheo/oxidation/in.rheo.oxidation index 5c1456ad5a..d26b379adb 100644 --- a/examples/rheo/oxidation/in.rheo.oxidation +++ b/examples/rheo/oxidation/in.rheo.oxidation @@ -10,7 +10,6 @@ newton off region box block -60 60 0 80 -0.01 0.01 units box create_box 3 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50 - region lbar block -15 0 3 80 EDGE EDGE units box region rbar block 0 15 3 80 EDGE EDGE units box region bar union 2 lbar rbar @@ -83,6 +82,7 @@ fix 7 all gravity 5e-5 vector 0 -1 0 fix 8 floor setforce 0.0 0.0 0.0 fix 9 surf add/heat linear 1.1 0.05 fix 10 floor add/heat constant 0 overwrite yes # fix the temperature of the floor +fix 11 all enforce2d compute surf all rheo/property/atom surface compute rho all rheo/property/atom rho diff --git a/examples/rheo/poiseuille/in.rheo.poiseuille b/examples/rheo/poiseuille/in.rheo.poiseuille index ee8ad5523c..7c38c8a6e0 100644 --- a/examples/rheo/poiseuille/in.rheo.poiseuille +++ b/examples/rheo/poiseuille/in.rheo.poiseuille @@ -6,7 +6,6 @@ atom_style rheo boundary p p p comm_modify vel yes - # ------ Create simulation box ------ # variable n equal 1.0 @@ -53,7 +52,6 @@ timestep ${dt_max} pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr} pair_coeff * * - # ------ Fixes & computes ------ # fix 1 all rheo ${cut} quintic 0 shift @@ -62,10 +60,10 @@ fix 2 all rheo/viscosity * constant ${eta} fix 3 all rheo/pressure * linear fix 4 rig setforce 0.0 0.0 0.0 fix 5 fluid addforce ${fext} 0.0 0.0 +fix 6 all enforce2d compute rho all rheo/property/atom rho - # ------ Output & Run ------ # thermo 200 diff --git a/examples/rheo/taylor-green/in.rheo.taylor.green b/examples/rheo/taylor-green/in.rheo.taylor.green index e386d8262c..4485387440 100644 --- a/examples/rheo/taylor-green/in.rheo.taylor.green +++ b/examples/rheo/taylor-green/in.rheo.taylor.green @@ -7,7 +7,6 @@ boundary p p p comm_modify vel yes newton off - # ------ Create simulation box ------ # variable n equal 1.0 @@ -47,12 +46,12 @@ timestep ${dt_max} pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr} pair_coeff * * - # ------ Fixes & computes ------ # fix 1 all rheo ${cut} RK1 8 shift fix 2 all rheo/viscosity * constant ${eta} fix 3 all rheo/pressure * linear +fix 4 all enforce2d compute rho all rheo/property/atom rho diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index bd16a89b6a..d0f24850bd 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -25,6 +25,7 @@ #include "error.h" #include "fix_rheo.h" #include "force.h" +#include "math_const.h" #include "math_extra.h" #include "memory.h" #include "modify.h" @@ -43,6 +44,7 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; +using namespace MathConst; using namespace MathExtra; static constexpr int DELTA = 2000; @@ -57,8 +59,7 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : kernel_style = utils::inumeric(FLERR,arg[3],false,lmp); - - if (kernel_style == QUINTIC) { + if (kernel_style == QUINTIC || kernel_style == WENDLANDC4) { correction_order = -1; } else if (kernel_style == RK0) { correction_order = 0; @@ -115,12 +116,22 @@ void ComputeRHEOKernel::init() hsqinv = hinv * hinv; - if (dim == 3) { - pre_w = 0.002652582384864922 * 27.0 * hsqinv * hinv; - pre_wp = pre_w * 3.0 * hinv; + if (kernel_style != WENDLANDC4) { + if (dim == 3) { + pre_w = 1.0 / (120.0 * MY_PI) * 27.0 * hsqinv * hinv; + pre_wp = pre_w * 3.0 * hinv; + } else { + pre_w = 7.0 / (478.0 * MY_PI) * 9 * hsqinv; + pre_wp = pre_w * 3.0 * hinv; + } } else { - pre_w = 0.004661441847879780 * 9 * hsqinv; - pre_wp = pre_w * 3.0 * hinv; + if (dim == 3) { + pre_w = 495.0 / (32.0 * MY_PI * hsq * h); + pre_wp = pre_w * hinv; + } else { + pre_w = 9.0 / (MY_PI * hsq); + pre_wp = pre_w * hinv; + } } nmax_store = atom->nmax; @@ -163,11 +174,27 @@ int ComputeRHEOKernel::check_corrections(int i) /* ---------------------------------------------------------------------- */ +double ComputeRHEOKernel::calc_w_self(int i, int j) +{ + double w; + if (kernel_style == WENDLANDC4) + w = calc_w_wendlandc4(i, j, 0.0, 0.0, 0.0, 0.0); + else + w = calc_w_quintic(i, j, 0.0, 0.0, 0.0, 0.0); + + return w; +} + +/* ---------------------------------------------------------------------- */ + double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double delz, double r) { double w; int corrections_i, corrections_j, corrections; + if (kernel_style == WENDLANDC4) + return calc_w_wendlandc4(i,j,delx,dely,delz,r); + if (kernel_style != QUINTIC) { corrections_i = check_corrections(i); corrections_j = check_corrections(j); @@ -191,6 +218,9 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double double wp; int corrections_i, corrections_j; + if (kernel_style == WENDLANDC4) + return calc_dw_wendlandc4(i,j,delx,dely,delz,r,dWij,dWji); + if (kernel_style != QUINTIC) { corrections_i = check_corrections(i); corrections_j = check_corrections(j); @@ -288,6 +318,62 @@ double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely return wp; } +/* ---------------------------------------------------------------------- */ + +double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double dely, double delz, double r) +{ + double w, tmp6, s; + s = r * hinv; + + if (s > 1.0) { + w = 0.0; + } else { + tmp6 = (1.0 - s) * (1.0 - s); + tmp6 *= tmp6 * tmp6; + w = tmp6 * (1.0 + 6.0 * s + 35.0 * THIRD * s * s); + } + + w *= pre_w; + + Wij = w; + Wji = w; + + return w; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2) +{ + double wp, tmp1, tmp5, tmp6, s, wprinv; + double *mass = atom->mass; + int *type = atom->type; + + s = r * hinv; + + if (s > 1.0) { + wp = 0.0; + } else { + tmp1 = 1.0 - s; + tmp5 = tmp1 * tmp1; + tmp5 = tmp5 * tmp5 * tmp1; + tmp6 = tmp5 * tmp1; + wp = tmp6 * (6.0 + 70.0 * THIRD * s); + wp -= 6 * tmp5 * (1.0 + 6.0 * s + 35.0 * THIRD * s * s); + } + + wp *= pre_wp; + 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; + + return wp; +} /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index ed190c19ce..543e5d88c0 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -36,10 +36,13 @@ class ComputeRHEOKernel : public Compute { void unpack_forward_comm(int, int, double *) override; double memory_usage() override; void compute_coordination(); + double calc_w_self(int,int); double calc_w(int,int,double,double,double,double); double calc_dw(int,int,double,double,double,double); double calc_w_quintic(int,int,double,double,double,double); double calc_dw_quintic(int,int,double,double,double,double,double *,double *); + double calc_w_wendlandc4(int,int,double,double,double,double); + double calc_dw_wendlandc4(int,int,double,double,double,double,double *,double *); void grow_arrays(int); double dWij[3], dWji[3], Wij, Wji; diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index c7270897a4..d7f432a03d 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -94,7 +94,7 @@ void ComputeRHEORhoSum::compute_peratom() // initialize arrays, local with quintic self-contribution, ghosts are zeroed for (i = 0; i < nlocal; i++) { - w = compute_kernel->calc_w_quintic(i, i, 0.0, 0.0, 0.0, 0.0); + w = compute_kernel->calc_w_self(i, i); rho[i] = w * mass[type[i]]; } diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index cb2abfb938..c704677bec 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -90,6 +90,8 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : cut = h; if (strcmp(arg[4], "quintic") == 0) { kernel_style = QUINTIC; + } else if (strcmp(arg[4], "wendland/c4") == 0) { + kernel_style = WENDLANDC4; } else if (strcmp(arg[4], "RK0") == 0) { kernel_style = RK0; } else if (strcmp(arg[4], "RK1") == 0) { diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 0af9fa8d01..251f82a99a 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -72,7 +72,7 @@ class FixRHEO : public Fix { namespace RHEO_NS { - enum {QUINTIC, RK0, RK1, RK2}; + enum {QUINTIC, WENDLANDC4, RK0, RK1, RK2}; enum {COORDINATION, DIVR}; // Status variables