More rheo pressure EOS, allow variable cs

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 <fix_rheo>` is not used.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

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<FixRHEOPressure *>(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];
       }

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;

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),
+    rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), gamma(nullptr),
-    fix_rheo(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;
+}

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;

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,
+  double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj;
-      kappa_ave, dT_prefactor;
+  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];
         }
 
+        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);

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;