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

Merge pull request #4091 from jtclemm/BPM

Browse Source 
Small fixes to GRANULAR/BPM packages
This commit is contained in:
Axel Kohlmeyer
2024-02-26 20:01:57 -05:00
committed by GitHub
parent e9deaa83d1 73ae3200d0
commit 4d89741d8c
13 changed files with 101 additions and 94 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
doc/src/Howto_granular.rst
View File

@ -45,10 +45,15 @@ atoms, and should be used for granular system instead of the fix style
To model heat conduction, one must add the temperature and heatflow
atom variables with:
* :doc:`fix property/atom <fix_property_atom>`
a temperature integration fix
* :doc:`fix heat/flow <fix_heat_flow>`
and a heat conduction option defined in both
* :doc:`pair_style granular <pair_granular>`
* :doc:`fix wall/gran <fix_wall_gran>`

8
doc/src/compute_fabric.rst
View File

@ -64,7 +64,7 @@ tangential force tensor. The contact tensor is calculated as
.. math::
C_{ab} = \frac{15}{2} (\phi_{ab} - \mathrm{Tr}(\phi) \delta_{ab})
C_{ab} = \frac{15}{2} (\phi_{ab} - \frac{1}{3} \mathrm{Tr}(\phi) \delta_{ab})
where :math:`a` and :math:`b` are the :math:`x`, :math:`y`, :math:`z`
directions, :math:`\delta_{ab}` is the Kronecker delta function, and
@ -83,7 +83,7 @@ The branch tensor is calculated as
.. math::
B_{ab} = \frac{15}{6 \mathrm{Tr}(D)} (D_{ab} - \mathrm{Tr}(D) \delta_{ab})
B_{ab} = \frac{15}{2\, \mathrm{Tr}(D)} (D_{ab} - \frac{1}{3} \mathrm{Tr}(D) \delta_{ab})
where the tensor :math:`D` is defined as
@ -101,7 +101,7 @@ The normal force fabric tensor is calculated as
.. math::
F^n_{ab} = \frac{15}{6\, \mathrm{Tr}(N)} (N_{ab} - \mathrm{Tr}(N) \delta_{ab})
F^n_{ab} = \frac{15}{2\, \mathrm{Tr}(N)} (N_{ab} - \frac{1}{3} \mathrm{Tr}(N) \delta_{ab})
where the tensor :math:`N` is defined as
@ -119,7 +119,7 @@ as
.. math::
F^t_{ab} = \frac{15}{9\, \mathrm{Tr}(N)} (T_{ab} - \mathrm{Tr}(T) \delta_{ab})
F^t_{ab} = \frac{5}{\mathrm{Tr}(N)} (T_{ab} - \frac{1}{3} \mathrm{Tr}(T) \delta_{ab})
where the tensor :math:`T` is defined as

3
doc/src/compute_property_atom.rst
View File

@ -23,6 +23,7 @@ Syntax
spx, spy, spz, sp, fmx, fmy, fmz,
nbonds,
radius, diameter, omegax, omegay, omegaz,
temperature, heatflow,
angmomx, angmomy, angmomz,
shapex, shapey, shapez,
quatw, quati, quatj, quatk, tqx, tqy, tqz,
@ -56,6 +57,8 @@ Syntax
*nbonds* = number of bonds assigned to an atom
*radius,diameter* = radius,diameter of spherical particle
*omegax,omegay,omegaz* = angular velocity of spherical particle
*temperature* = internal temperature of spherical particle
*heatflow* = internal heat flow of spherical particle
*angmomx,angmomy,angmomz* = angular momentum of aspherical particle
*shapex,shapey,shapez* = 3 diameters of aspherical particle
*quatw,quati,quatj,quatk* = quaternion components for aspherical or body particles

3
doc/src/dump.rst
View File

@ -104,7 +104,6 @@ Syntax
q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz,
heatflow, temperature,
c_ID, c_ID[I], f_ID, f_ID[I], v_name,
i_name, d_name, i2_name[I], d2_name[I]
@ -131,8 +130,6 @@ Syntax
omegax,omegay,omegaz = angular velocity of spherical particle
angmomx,angmomy,angmomz = angular momentum of aspherical particle
tqx,tqy,tqz = torque on finite-size particles
heatflow = rate of heat flow into particle
temperature = temperature of particle
c_ID = per-atom vector calculated by a compute with ID
c_ID[I] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below)
f_ID = per-atom vector calculated by a fix with ID

17
doc/src/fix_property_atom.rst
View File

@ -22,6 +22,8 @@ Syntax
*mol* = molecule IDs
*q* = charge
*rmass* = per-atom mass
*temperature* = internal temperature of atom
*heatflow* = internal heat flow of atom
i_name = new integer vector referenced by name
d_name = new floating-point vector referenced by name
i2_name = new integer array referenced by name
@ -59,14 +61,18 @@ these properties for each atom in the system when a data file is read.
This fix augments the set of per-atom properties with new custom
ones. This can be useful in several scenarios.
If the atom style does not define molecule IDs, per-atom charge, or
per-atom mass, they can be added using the *mol*\ , *q* or *rmass*
If the atom style does not define molecule IDs, per-atom charge,
per-atom mass, internal temperature, or internal heat flow, they can
be added using the *mol*\ , *q*, *rmass*, *temperature*, or *heatflow*
keywords. This could be useful to define "molecules" to use as rigid
bodies with the :doc:`fix rigid <fix_rigid>` command, or to carry
around an extra flag with atoms (stored as a molecule ID) that can be
used by various commands like :doc:`compute chunk/atom
<compute_chunk_atom>` to group atoms without having to use the group
command (which is limited to a total of 32 groups including *all*\ ).
For finite-size particles, an internal temperature and heat flow can
be used to model heat conduction as in the
:doc:`GRANULAR package <Howto_granular>`.
Another application is to use the *rmass* flag in order to have
per-atom masses instead of per-type masses. This could be used to
@ -85,9 +91,10 @@ properties that are not needed such as bond lists, which incurs some
overhead when there are no bonds.
In the future, we may add additional existing per-atom properties to
fix property/atom, similar to *mol*\ , *q* or *rmass*\ , which
"turn-on" specific properties defined by some atom styles, so they can
be easily used by atom styles that do not define them.
fix property/atom, similar to *mol*\ , *q*, *rmass*\ , *temperature*\ ,
or *heatflow* which "turn-on" specific properties defined by some atom
styles, so they can be easily used by atom styles that do not define
them.
More generally, the *i_name* and *d_name* options allow one or more
new custom per-atom vectors to be defined. Likewise the *i2_name* and

3
examples/granular/in.pour.heat
View File

@ -73,7 +73,8 @@ thermo 100
timestep 0.001
#dump 1 all custom 1000 ${name}.dump id type radius mass x y z temperature heatflow
compute 1 all property/atom temperature heatflow
#dump 1 all custom 1000 ${name}.dump id type radius mass x y z c_1[*]
run 100000

2
src/BPM/bond_bpm.cpp
View File

@ -357,7 +357,6 @@ void BondBPM::process_broken(int i, int j)
if (i < nlocal) {
for (m = 0; m < num_bond[i]; m++) {
if (bond_atom[i][m] == tag[j]) {
bond_type[i][m] = 0;
n = num_bond[i];
bond_type[i][m] = bond_type[i][n - 1];
bond_atom[i][m] = bond_atom[i][n - 1];
@ -372,7 +371,6 @@ void BondBPM::process_broken(int i, int j)
if (j < nlocal) {
for (m = 0; m < num_bond[j]; m++) {
if (bond_atom[j][m] == tag[i]) {
bond_type[j][m] = 0;
n = num_bond[j];
bond_type[j][m] = bond_type[j][n - 1];
bond_atom[j][m] = bond_atom[j][n - 1];

52
src/GRANULAR/compute_fabric.cpp
View File

@ -184,7 +184,7 @@ void ComputeFabric::compute_vector()
double nx, ny, nz;
double ncinv, denom, fn, ft, prefactor;
double br_tensor[6], ft_tensor[6], fn_tensor[6];
double trace_phi, trace_D, trace_Xfn, trace_Xft;
double trace_third_phi, trace_third_D, trace_third_Xfn, trace_third_Xft;
double phi_ij[6] = {0.0};
double Ac_ij[6] = {0.0};
double D_ij[6] = {0.0};
@ -295,11 +295,11 @@ void ComputeFabric::compute_vector()
MPI_Allreduce(phi_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
for (i = 0; i < 6; i++) phi_ij[i] = temp_dbl[i] * ncinv;
trace_phi = (1.0 / 3.0) * (phi_ij[0] + phi_ij[1] + phi_ij[2]);
trace_third_phi = (1.0 / 3.0) * (phi_ij[0] + phi_ij[1] + phi_ij[2]);
Ac_ij[0] = (15.0 / 2.0) * (phi_ij[0] - trace_phi);
Ac_ij[1] = (15.0 / 2.0) * (phi_ij[1] - trace_phi);
Ac_ij[2] = (15.0 / 2.0) * (phi_ij[2] - trace_phi);
Ac_ij[0] = (15.0 / 2.0) * (phi_ij[0] - trace_third_phi);
Ac_ij[1] = (15.0 / 2.0) * (phi_ij[1] - trace_third_phi);
Ac_ij[2] = (15.0 / 2.0) * (phi_ij[2] - trace_third_phi);
Ac_ij[3] = (15.0 / 2.0) * (phi_ij[3]);
Ac_ij[4] = (15.0 / 2.0) * (phi_ij[4]);
Ac_ij[5] = (15.0 / 2.0) * (phi_ij[5]);
@ -419,14 +419,14 @@ void ComputeFabric::compute_vector()
MPI_Allreduce(D_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
for (i = 0; i < 6; i++) D_ij[i] = temp_dbl[i];
trace_D = (1.0 / 3.0) * (D_ij[0] + D_ij[1] + D_ij[2]);
trace_third_D = (1.0 / 3.0) * (D_ij[0] + D_ij[1] + D_ij[2]);
br_tensor[0] = (15.0 / (6.0 * trace_D)) * (D_ij[0] - trace_D);
br_tensor[1] = (15.0 / (6.0 * trace_D)) * (D_ij[1] - trace_D);
br_tensor[2] = (15.0 / (6.0 * trace_D)) * (D_ij[2] - trace_D);
br_tensor[3] = (15.0 / (6.0 * trace_D)) * (D_ij[3]);
br_tensor[4] = (15.0 / (6.0 * trace_D)) * (D_ij[4]);
br_tensor[5] = (15.0 / (6.0 * trace_D)) * (D_ij[5]);
br_tensor[0] = (15.0 / (6.0 * trace_third_D)) * (D_ij[0] - trace_third_D);
br_tensor[1] = (15.0 / (6.0 * trace_third_D)) * (D_ij[1] - trace_third_D);
br_tensor[2] = (15.0 / (6.0 * trace_third_D)) * (D_ij[2] - trace_third_D);
br_tensor[3] = (15.0 / (6.0 * trace_third_D)) * (D_ij[3]);
br_tensor[4] = (15.0 / (6.0 * trace_third_D)) * (D_ij[4]);
br_tensor[5] = (15.0 / (6.0 * trace_third_D)) * (D_ij[5]);
for (i = 0; i < ntensors; i++) {
if (tensor_style[i] == BR) {
@ -439,17 +439,17 @@ void ComputeFabric::compute_vector()
MPI_Allreduce(Xfn_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
for (i = 0; i < 6; i++) Xfn_ij[i] = temp_dbl[i];
trace_Xfn = (1.0 / 3.0) * (Xfn_ij[0] + Xfn_ij[1] + Xfn_ij[2]);
trace_third_Xfn = (1.0 / 3.0) * (Xfn_ij[0] + Xfn_ij[1] + Xfn_ij[2]);
}
if (fn_flag) {
fn_tensor[0] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[0] - trace_Xfn);
fn_tensor[1] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[1] - trace_Xfn);
fn_tensor[2] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[2] - trace_Xfn);
fn_tensor[3] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[3]);
fn_tensor[4] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[4]);
fn_tensor[5] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[5]);
fn_tensor[0] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[0] - trace_third_Xfn);
fn_tensor[1] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[1] - trace_third_Xfn);
fn_tensor[2] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[2] - trace_third_Xfn);
fn_tensor[3] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[3]);
fn_tensor[4] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[4]);
fn_tensor[5] = (15.0 / (6.0 * trace_third_Xfn)) * (Xfn_ij[5]);
for (i = 0; i < ntensors; i++) {
if (tensor_style[i] == FN) {
@ -462,14 +462,14 @@ void ComputeFabric::compute_vector()
MPI_Allreduce(Xft_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
for (i = 0; i < 6; i++) Xft_ij[i] = temp_dbl[i];
trace_Xft = (1.0 / 3.0) * (Xft_ij[0] + Xft_ij[1] + Xft_ij[2]);
trace_third_Xft = (1.0 / 3.0) * (Xft_ij[0] + Xft_ij[1] + Xft_ij[2]);
ft_tensor[0] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[0] - trace_Xft);
ft_tensor[1] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[1] - trace_Xft);
ft_tensor[2] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[2] - trace_Xft);
ft_tensor[3] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[3]);
ft_tensor[4] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[4]);
ft_tensor[5] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[5]);
ft_tensor[0] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[0] - trace_third_Xft);
ft_tensor[1] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[1] - trace_third_Xft);
ft_tensor[2] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[2] - trace_third_Xft);
ft_tensor[3] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[3]);
ft_tensor[4] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[4]);
ft_tensor[5] = (15.0 / (9.0 * trace_third_Xfn)) * (Xft_ij[5]);
for (i = 0; i < ntensors; i++) {
if (tensor_style[i] == FT) {

8
src/GRANULAR/gran_sub_mod_damping.cpp
View File

@ -130,6 +130,12 @@ void GranSubModDampingTsuji::init()
double GranSubModDampingTsuji::calculate_forces()
{
damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta);
// in case argument <= 0 due to precision issues
double sqrt1;
if (gm->delta > 0.0)
sqrt1 = MAX(0.0, gm->meff * gm->Fnormal / gm->delta);
else
sqrt1 = 0.0;
damp_prefactor = damp * sqrt(sqrt1);
return -damp_prefactor * gm->vnnr;
}

38
src/compute_property_atom.cpp
View File

@ -205,6 +205,14 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (!atom->omega_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
pack_choice[i] = &ComputePropertyAtom::pack_omegaz;
} else if (strcmp(arg[iarg],"temperature") == 0) {
if (!atom->temperature_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
pack_choice[i] = &ComputePropertyAtom::pack_temperature;
} else if (strcmp(arg[iarg],"heatflow") == 0) {
if (!atom->heatflow_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
pack_choice[i] = &ComputePropertyAtom::pack_heatflow;
} else if (strcmp(arg[iarg],"angmomx") == 0) {
if (!atom->angmom_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
@ -1213,6 +1221,36 @@ void ComputePropertyAtom::pack_omegaz(int n)
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_temperature(int n)
{
double *temperature = atom->temperature;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = temperature[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_heatflow(int n)
{
double *heatflow = atom->heatflow;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = heatflow[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_angmomx(int n)
{
double **angmom = atom->angmom;

2
src/compute_property_atom.h
View File

@ -95,6 +95,8 @@ class ComputePropertyAtom : public Compute {
void pack_omegax(int);
void pack_omegay(int);
void pack_omegaz(int);
void pack_temperature(int);
void pack_heatflow(int);
void pack_angmomx(int);
void pack_angmomy(int);
void pack_angmomz(int);

50
src/dump_custom.cpp
View File

@ -41,7 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI,
IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,HEATFLOW,TEMPERATURE,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
TQX,TQY,TQZ,
COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY};
@ -929,18 +929,6 @@ int DumpCustom::count()
for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == HEATFLOW) {
if (!atom->heatflow_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = atom->heatflow;
nstride = 1;
} else if (thresh_array[ithresh] == TEMPERATURE) {
if (!atom->temperature_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = atom->temperature;
nstride = 1;
} else if (thresh_array[ithresh] == OMEGAX) {
if (!atom->omega_flag)
error->all(FLERR,
@ -1395,16 +1383,6 @@ int DumpCustom::parse_fields(int narg, char **arg)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[iarg] = &DumpCustom::pack_diameter;
vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"heatflow") == 0) {
if (!atom->heatflow_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[iarg] = &DumpCustom::pack_heatflow;
vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"temperature") == 0) {
if (!atom->temperature_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[iarg] = &DumpCustom::pack_temperature;
vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"omegax") == 0) {
if (!atom->omega_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
@ -1875,8 +1853,6 @@ int DumpCustom::modify_param(int narg, char **arg)
else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER;
else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = HEATFLOW;
else if (strcmp(arg[1],"temperature") == 0) thresh_array[nthresh] = TEMPERATURE;
else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;
else if (strcmp(arg[1],"omegaz") == 0) thresh_array[nthresh] = OMEGAZ;
@ -2791,30 +2767,6 @@ void DumpCustom::pack_diameter(int n)
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_heatflow(int n)
{
double *heatflow = atom->heatflow;
for (int i = 0; i < nchoose; i++) {
buf[n] = heatflow[clist[i]];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_temperature(int n)
{
double *temperature = atom->temperature;
for (int i = 0; i < nchoose; i++) {
buf[n] = temperature[clist[i]];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_omegax(int n)
{
double **omega = atom->omega;

2
src/dump_custom.h
View File

@ -188,8 +188,6 @@ class DumpCustom : public Dump {
void pack_mu(int);
void pack_radius(int);
void pack_diameter(int);
void pack_heatflow(int);
void pack_temperature(int);
void pack_omegax(int);
void pack_omegay(int);