ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Separating background pressure from EoS definition

Browse Source
This commit is contained in:
jtclemm
2024-11-22 14:37:47 -07:00
parent cb5bad7ece
commit 10c429fe21
2 changed files with 35 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
doc/src/fix_rheo_pressure.rst
View File

@ -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 <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
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

30
src/RHEO/fix_rheo_pressure.cpp
View File

@ -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]);