Dam example, generalizing non-quintic kernels

doc/src/Howto_rheo.rst

@@ -17,6 +17,10 @@ these calculations. Additional numerical details can be found in
 :ref:`(Palermo) <howto_rheo_palermo>` and
 :ref:`(Clemmer) <howto_rheo_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 <fix_rheo>` which integrates

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

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           cut equal 2.2
-variable           dx equal 4.0
+variable           dx equal 3.0
 
-region             box block -1 200 -1 80 -0.1 0.1 units box
+region             box block -1 150 -1 80 -0.1 0.1 units box
-#region             box block -1 100 -1 40 -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+1.0) $(xlo+40.0) $(ylo+v_dx+1.0) $(yhi-20.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             walls1 block $(xlo+v_dx) $(xhi-v_dx) $(ylo+v_dx) $(yhi-v_dx) EDGE EDGE side out units box
-region             walls block $(xlo+v_dx) $(xhi-v_dx) $(ylo+v_dx) $(yhi-v_dx) EDGE EDGE units box side out
+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           eta equal 0.1
-variable           Dr equal 0.05
+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
+                              rho/sum
-#fix                1 all rheo ${cut} RK1 10 surface/detection coordination 36 10 interface/reconstruct shift
 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                30000
-run   300000      #        400000

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

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

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

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
 

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 || kernel_style == WENDLANDC4) {
-  if (kernel_style == QUINTIC) {
     correction_order = -1;
   } else if (kernel_style == RK0) {
     correction_order = 0;
@@ -115,13 +116,23 @@ void ComputeRHEOKernel::init()
   hsqinv = hinv * hinv;
 
 
+  if (kernel_style != WENDLANDC4) {
     if (dim == 3) {
-    pre_w = 0.002652582384864922 * 27.0 * hsqinv * hinv;
+      pre_w = 1.0 / (120.0 * MY_PI) * 27.0 * hsqinv * hinv;
       pre_wp = pre_w * 3.0 * hinv;
     } else {
-    pre_w = 0.004661441847879780 * 9 * hsqinv;
+      pre_w = 7.0 / (478.0 * MY_PI)  * 9 * hsqinv;
       pre_wp = pre_w * 3.0 * hinv;
     }
+  } else {
+    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;
   memory->create(coordination, nmax_store, "rheo:coordination");
@@ -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;
+}
 
 /* ---------------------------------------------------------------------- */
 

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;

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]];
   }
 

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) {

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