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

Prototyping latent heat, other misc fixes

Browse Source
This commit is contained in:
jtclemm
2023-12-22 10:46:44 -07:00
parent 7403426046
commit 08d2dd2699
16 changed files with 611 additions and 126 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
src/.gitignore vendored
View File

@ -197,40 +197,7 @@
/pair_tdpd.cpp
/pair_tdpd.h
/atom_vec_rheo.cpp
/atom_vec_rheo.h
/atom_vec_rheo_thermal.cpp
/atom_vec_rheo_thermal.h
/compute_rheo_grad.cpp
/compute_rheo_grad.h
/compute_rheo_interface.cpp
/compute_rheo_interface.h
/compute_rheo_kernel.cpp
/compute_rheo_kernel.h
/compute_rheo_property_atom.cpp
/compute_rheo_property_atom.h
/compute_rheo_rho_sum.cpp
/compute_rheo_rho_sum.h
/compute_rheo_surface.cpp
/compute_rheo_surface.h
/compute_rheo_vshift.cpp
/compute_rheo_vshift.h
/fix_rheo.cpp
/fix_rheo.h
/fix_rheo_pressure.cpp
/fix_rheo_pressure.h
/fix_rheo_stress.cpp
/fix_rheo_stress.h
/fix_rheo_tension.cpp
/fix_rheo_tension.h
/fix_rheo_thermal.cpp
/fix_rheo_thermal.h
/fix_rheo_viscosity.cpp
/fix_rheo_viscosity.h
/pair_rheo.cpp
/pair_rheo.h
/pair_rheo_react.cpp
/pair_rheo_react.h
/*rheo*
/compute_grid.cpp
/compute_grid.h

46
src/RHEO/atom_vec_rheo_thermal.cpp
View File

@ -35,6 +35,7 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
atom->status_flag = 1;
atom->conductivity_flag = 1;
atom->temperature_flag = 1;
atom->esph_flag = 1;
atom->heatflow_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
@ -45,17 +46,17 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"};
fields_copy = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"};
fields_comm = {"status", "rho", "temperature"};
fields_comm_vel = {"status", "rho", "temperature"};
fields_grow = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_copy = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_comm = {"status", "rho", "esph"};
fields_comm_vel = {"status", "rho", "esph"};
fields_reverse = {"drho", "heatflow"};
fields_border = {"status", "rho", "temperature"};
fields_border_vel = {"status", "rho", "temperature"};
fields_exchange = {"status", "rho", "temperature"};
fields_restart = {"status", "rho", "temperature"};
fields_create = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"};
fields_data_atom = {"id", "type", "status", "rho", "temperature", "x"};
fields_border = {"status", "rho", "esph"};
fields_border_vel = {"status", "rho", "esph"};
fields_exchange = {"status", "rho", "esph"};
fields_restart = {"status", "rho", "esph"};
fields_create = {"status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_data_atom = {"id", "type", "status", "rho", "esph", "x"};
fields_data_vel = {"id", "v"};
setup_fields();
@ -71,6 +72,7 @@ void AtomVecRHEOThermal::grow_pointers()
status = atom->status;
conductivity = atom->conductivity;
temperature = atom->temperature;
esph = atom->esph;
heatflow = atom->heatflow;
pressure = atom->pressure;
rho = atom->rho;
@ -98,6 +100,7 @@ void AtomVecRHEOThermal::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
heatflow[ilocal] = 0.0;
temperature[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
conductivity[ilocal] = 0.0;
@ -114,10 +117,11 @@ int AtomVecRHEOThermal::property_atom(const std::string &name)
if (name == "rho") return 1;
if (name == "drho") return 2;
if (name == "temperature") return 3;
if (name == "heatflow") return 4;
if (name == "conductivity") return 5;
if (name == "pressure") return 6;
if (name == "viscosity") return 7;
if (name == "esph") return 4;
if (name == "heatflow") return 5;
if (name == "conductivity") return 6;
if (name == "pressure") return 7;
if (name == "viscosity") return 8;
return -1;
}
@ -167,7 +171,7 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues,
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = heatflow[i];
buf[n] = esph[i];
else
buf[n] = 0.0;
n += nvalues;
@ -175,7 +179,7 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues,
} else if (index == 5) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = conductivity[i];
buf[n] = heatflow[i];
else
buf[n] = 0.0;
n += nvalues;
@ -183,12 +187,20 @@ void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues,
} else if (index == 6) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = pressure[i];
buf[n] = conductivity[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 7) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = pressure[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 8) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];

2
src/RHEO/atom_vec_rheo_thermal.h
View File

@ -36,7 +36,7 @@ class AtomVecRHEOThermal : virtual public AtomVec {
private:
int *status;
double *conductivity, *temperature, *heatflow;
double *conductivity, *temperature, *heatflow, *esph;
double *pressure, *rho, *drho, *viscosity;
};

63
src/RHEO/compute_rheo_grad.cpp
View File

@ -42,15 +42,15 @@ enum{COMMGRAD, COMMFIELD};
ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), rho0(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
gradv(nullptr), gradr(nullptr), gradt(nullptr), gradn(nullptr)
gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
velocity_flag = temperature_flag = rho_flag = eta_flag = 0;
velocity_flag = energy_flag = rho_flag = eta_flag = 0;
for (int iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg], "velocity") == 0) velocity_flag = 1;
else if (strcmp(arg[iarg], "rho") == 0) rho_flag = 1;
else if (strcmp(arg[iarg],"temperature") == 0) temperature_flag = 1;
else if (strcmp(arg[iarg], "energy") == 0) energy_flag = 1;
else if (strcmp(arg[iarg], "viscosity") == 0) eta_flag = 1;
else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]);
}
@ -72,7 +72,7 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
comm_reverse += dim;
}
if (temperature_flag) {
if (energy_flag) {
ncomm_grad += dim;
ncomm_field += 1;
comm_reverse += dim;
@ -96,7 +96,7 @@ ComputeRHEOGrad::~ComputeRHEOGrad()
{
memory->destroy(gradv);
memory->destroy(gradr);
memory->destroy(gradt);
memory->destroy(grade);
memory->destroy(gradn);
}
@ -130,7 +130,7 @@ void ComputeRHEOGrad::compute_peratom()
int i, j, k, ii, jj, jnum, itype, jtype, a, b, fluidi, fluidj;
double xtmp, ytmp, ztmp, delx, dely, delz;
double rsq, imass, jmass;
double rhoi, rhoj, Voli, Volj, drho, dT, deta;
double rhoi, rhoj, Voli, Volj, drho, de, deta;
double vi[3], vj[3], vij[3];
double wp, *dWij, *dWji;
@ -141,7 +141,7 @@ void ComputeRHEOGrad::compute_peratom()
double **x = atom->x;
double **v = atom->v;
double *rho = atom->rho;
double *temperature = atom->temperature;
double *energy = atom->esph;
double *viscosity = atom->viscosity;
int *status = atom->status;
int *type = atom->type;
@ -166,9 +166,9 @@ void ComputeRHEOGrad::compute_peratom()
for (k = 0; k < dim; k++)
gradr[i][k] = 0.0;
}
if (temperature_flag) {
if (energy_flag) {
for (k = 0; k < dim; k++)
gradt[i][k] = 0.0;
grade[i][k] = 0.0;
}
if (eta_flag) {
for (k = 0; k < dim; k++)
@ -234,7 +234,7 @@ void ComputeRHEOGrad::compute_peratom()
vij[2] = vi[2] - vj[2];
if (rho_flag) drho = rhoi - rhoj;
if (temperature_flag) dT = temperature[i] - temperature[j];
if (energy_flag) de = energy[i] - energy[j];
if (eta_flag) deta = viscosity[i] - viscosity[j];
wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));
@ -250,8 +250,8 @@ void ComputeRHEOGrad::compute_peratom()
if (rho_flag) // P,x P,y P,z
gradr[i][a] -= drho * Volj * dWij[a];
if (temperature_flag) // T,x T,y T,z
gradt[i][a] -= dT * Volj * dWij[a];
if (energy_flag) // e,x e,y e,z
grade[i][a] -= de * Volj * dWij[a];
if (eta_flag) // n,x n,y n,z
gradn[i][a] -= deta * Volj * dWij[a];
@ -267,8 +267,8 @@ void ComputeRHEOGrad::compute_peratom()
if (rho_flag) // P,x P,y P,z
gradr[j][a] += drho * Voli * dWji[a];
if (temperature_flag) // T,x T,y T,z
gradt[j][a] += dT * Voli * dWji[a];
if (energy_flag) // e,x e,y e,z
grade[j][a] += de * Voli * dWji[a];
if (eta_flag) // n,x n,y n,z
gradn[j][a] += deta * Voli * dWji[a];
@ -308,7 +308,7 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
int i,j,k,m;
int *mask = atom->mask;
double *rho = atom->rho;
double *temperature = atom->temperature;
double *energy = atom->esph;
double **v = atom->v;
int dim = domain->dimension;
double *h_rate = domain->h_rate;
@ -335,9 +335,9 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
for (k = 0; k < dim; k++)
buf[m++] = gradr[j][k];
if (temperature_flag)
if (energy_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradt[j][k];
buf[m++] = grade[j][k];
if (eta_flag)
for (k = 0; k < dim; k++)
@ -358,8 +358,8 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
if (rho_flag)
buf[m++] = rho[j];
if (temperature_flag)
buf[m++] = temperature[j];
if (energy_flag)
buf[m++] = energy[j];
}
}
return m;
@ -371,8 +371,7 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, m, last;
double *rho = atom->rho;
double *temperature = atom->temperature;
double **v = atom->v;
double *energy = atom->esph; double **v = atom->v;
int dim = domain->dimension;
m = 0;
@ -387,9 +386,9 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
for (k = 0; k < dim; k++)
gradr[i][k] = buf[m++];
if (temperature_flag)
if (energy_flag)
for (k = 0; k < dim; k++)
gradt[i][k] = buf[m++];
grade[i][k] = buf[m++];
if (eta_flag)
for (k = 0; k < dim; k++)
@ -403,8 +402,8 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
if (rho_flag)
rho[i] = buf[m++];
if (temperature_flag)
temperature[i] = buf[m++];
if (energy_flag)
energy[i] = buf[m++];
}
}
}
@ -427,9 +426,9 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf)
for (k = 0; k < dim; k++)
buf[m++] = gradr[i][k];
if (temperature_flag)
if (energy_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradt[i][k];
buf[m++] = grade[i][k];
if (eta_flag)
for (k = 0; k < dim; k++)
@ -456,9 +455,9 @@ void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf)
for (k = 0; k < dim; k++)
gradr[j][k] += buf[m++];
if (temperature_flag)
if (energy_flag)
for (k = 0; k < dim; k++)
gradt[j][k] += buf[m++];
grade[j][k] += buf[m++];
if (eta_flag)
for (k = 0; k < dim; k++)
@ -477,8 +476,8 @@ void ComputeRHEOGrad::grow_arrays(int nmax)
if (rho_flag)
memory->grow(gradr, nmax, dim, "rheo:grad_rho");
if (temperature_flag)
memory->grow(gradt, nmax, dim, "rheo:grad_temp");
if (energy_flag)
memory->grow(grade, nmax, dim, "rheo:grad_energy");
if (eta_flag)
memory->grow(gradn, nmax, dim, "rheo:grad_eta");
@ -498,7 +497,7 @@ double ComputeRHEOGrad::memory_usage()
if (rho_flag)
bytes = (size_t) nmax_store * dim * sizeof(double);
if (temperature_flag)
if (energy_flag)
bytes = (size_t) nmax_store * dim * sizeof(double);
if (eta_flag)

4
src/RHEO/compute_rheo_grad.h
View File

@ -40,7 +40,7 @@ class ComputeRHEOGrad : public Compute {
void forward_fields();
double **gradv;
double **gradr;
double **gradt;
double **grade;
double **gradn;
class FixRHEO *fix_rheo;
@ -48,7 +48,7 @@ class ComputeRHEOGrad : public Compute {
int comm_stage, ncomm_grad, ncomm_field, nmax_store;
double cut, cutsq, *rho0;
int velocity_flag, temperature_flag, rho_flag, eta_flag;
int velocity_flag, energy_flag, rho_flag, eta_flag;
int interface_flag, remap_v_flag;
class ComputeRHEOKernel *compute_kernel;

7
src/RHEO/compute_rheo_kernel.cpp
View File

@ -86,6 +86,8 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
}
comm_forward_save = comm_forward;
corrections_calculated = 0;
gsl_error_flag = 0;
}
/* ---------------------------------------------------------------------- */
@ -152,6 +154,10 @@ int ComputeRHEOKernel::check_corrections(int i)
if (coordination[i] < zmin)
return 0;
// Skip if corrections not yet calculated
if (!corrections_calculated)
return 0;
return 1;
}
@ -479,6 +485,7 @@ void ComputeRHEOKernel::compute_peratom()
gsl_error_tags.clear();
if (kernel_style == QUINTIC) return;
corrections_calculated = 1;
int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error;
double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj;

1
src/RHEO/compute_rheo_kernel.h
View File

@ -53,6 +53,7 @@ class ComputeRHEOKernel : public Compute {
int gsl_error_flag;
std::unordered_set<tagint> gsl_error_tags;
int corrections_calculated;
int kernel_style, zmin, dim, Mdim, ncor;
int nmax_store;
double h, hsq, hinv, hsqinv, pre_w, pre_wp;

8
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -104,6 +104,14 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
} else if (utils::strmatch(arg[iarg], "^grad/v")) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_gradv;
col_index[i] = get_tensor_index(arg[iarg]);
} else if (strcmp(arg[iarg], "energy") == 0) {
avec_index[i] = atom->avec->property_atom("esph");
if (avec_index[i] < 0)
error->all(FLERR,
"Invalid keyword {} for atom style {} in compute rheo/property/atom command ",
arg[iarg], atom->get_style());
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
thermal_flag = 1;
} else {
avec_index[i] = atom->avec->property_atom(arg[iarg]);
if (avec_index[i] < 0)

4
src/RHEO/compute_rheo_vshift.cpp
View File

@ -32,6 +32,8 @@
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace RHEO_NS;
@ -195,6 +197,7 @@ void ComputeRHEOVShift::compute_peratom()
if (mask[i] & groupbit) {
vmag = sqrt(vi[0] * vi[0] + vi[1] * vi[1] + vi[2] * vi[2]);
prefactor = vmag * volj * dr;
vshift[i][0] += prefactor * dx[0];
vshift[i][1] += prefactor * dx[1];
vshift[i][2] += prefactor * dx[2];
@ -204,6 +207,7 @@ void ComputeRHEOVShift::compute_peratom()
if (mask[j] & groupbit) {
vmag = sqrt(vj[0] * vj[0] + vj[1] * vj[1] + vj[2] * vj[2]);
prefactor = vmag * voli * dr;
vshift[j][0] -= prefactor * dx[0];
vshift[j][1] -= prefactor * dx[1];
vshift[j][2] -= prefactor * dx[2];

3
src/RHEO/fix_rheo.cpp
View File

@ -164,7 +164,7 @@ void FixRHEO::post_constructor()
compute_kernel->fix_rheo = this;
std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity";
if (thermal_flag) cmd += " temperature";
if (thermal_flag) cmd += " energy";
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(cmd));
compute_grad->fix_rheo = this;
@ -401,6 +401,7 @@ void FixRHEO::pre_force(int /*vflag*/)
compute_rhosum->compute_peratom();
compute_kernel->compute_peratom();
if (interface_flag) {
// Note on first setup, have no forces for pressure to reference
compute_interface->compute_peratom();

135
src/RHEO/fix_rheo_thermal.cpp
View File

@ -47,7 +47,8 @@ enum {NONE, CONSTANT};
FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr),
Tc(nullptr), kappa(nullptr), cv(nullptr), Tc_style(nullptr), kappa_style(nullptr), cv_style(nullptr),
Tc(nullptr), kappa(nullptr), cv(nullptr), L(nullptr),
Tc_style(nullptr), kappa_style(nullptr), cv_style(nullptr), L_style(nullptr),
fix_update_special_bonds(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");
@ -62,13 +63,22 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
int i, nlo, nhi;
int n = atom->ntypes;
memory->create(Tc_style, n + 1, "rheo:Tc_style");
memory->create(kappa_style, n + 1, "rheo:kappa_style");
memory->create(cv_style, n + 1, "rheo:cv_style");
memory->create(L_style, n + 1, "rheo:L_style");
memory->create(Tc, n + 1, "rheo:Tc");
memory->create(kappa, n + 1, "rheo:kappa");
memory->create(cv, n + 1, "rheo:cv");
memory->create(L, n + 1, "rheo:L");
for (i = 1; i <= n; i++) {
Tc_style[i] = NONE;
kappa_style[i] = NONE;
cv_style[i] = NONE;
L_style[i] = NONE;
}
int iarg = 3;
@ -81,9 +91,9 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg + 2], "constant") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error);
double kappa_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
double kappa_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (kappa_one < 0.0) error->all(FLERR, "The conductivity must be positive");
iarg += 1;
iarg += 2;
for (i = nlo; i <= nhi; i++) {
kappa_style[i] = CONSTANT;
@ -102,9 +112,9 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg + 2], "constant") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error);
double cv_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
double cv_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (cv_one < 0.0) error->all(FLERR, "The specific heat must be positive");
iarg += 1;
iarg += 2;
for (i = nlo; i <= nhi; i++) {
cv_style[i] = CONSTANT;
@ -124,8 +134,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg + 2], "constant") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error);
double Tc_one = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
iarg += 1;
double Tc_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
iarg += 2;
for (i = nlo; i <= nhi; i++) {
Tc_style[i] = CONSTANT;
@ -136,6 +146,28 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]);
}
iarg += 2;
} else if (strcmp(arg[iarg], "latent/heat") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat", error);
utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
// Cv arguments
if (strcmp(arg[iarg + 2], "constant") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error);
double L_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
if (L_one < 0.0) error->all(FLERR, "The latent heat must be positive");
iarg += 2;
for (i = nlo; i <= nhi; i++) {
L_style[i] = CONSTANT;
L[i] = L_one;
}
} else {
error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]);
}
iarg += 2;
} else if (strcmp(arg[iarg], "react") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal react", error);
@ -158,6 +190,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Must specify specific/heat for atom type {} in fix/rheo/thermal", i);
if (kappa_style[i] == NONE)
error->all(FLERR, "Must specify conductivity for atom type {} in fix/rheo/thermal", i);
if (Tc_style[i] == NONE && L_style[i] != NONE)
error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i);
}
}
@ -173,9 +207,11 @@ FixRHEOThermal::~FixRHEOThermal()
memory->destroy(cv_style);
memory->destroy(Tc_style);
memory->destroy(kappa_style);
memory->destroy(L_style);
memory->destroy(cv);
memory->destroy(Tc);
memory->destroy(kappa);
memory->destroy(L);
}
/* ---------------------------------------------------------------------- */
@ -210,6 +246,8 @@ void FixRHEOThermal::init()
dtf = 0.5 * update->dt * force->ftm2v;
if (atom->esph_flag != 1)
error->all(FLERR,"fix rheo/thermal command requires atom property esph");
if (atom->temperature_flag != 1)
error->all(FLERR,"fix rheo/thermal command requires atom property temperature");
if (atom->heatflow_flag != 1)
@ -267,9 +305,9 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
int i, a;
int *status = atom->status;
double *temperature = atom->temperature;
double **gradt = compute_grad->gradt;
double **vshift = compute_vshift->array_atom;
double *energy = atom->esph;
double **grade = compute_grad->grade;
double **vshift = compute_vshift->vshift;
int nlocal = atom->nlocal;
int dim = domain->dimension;
@ -279,9 +317,8 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
for (i = 0; i < nlocal; i++) {
if (status[i] & STATUS_NO_SHIFT) continue;
for (a = 0; a < dim; a++)
temperature[i] += dtv * vshift[i][a] * gradt[i][a];
energy[i] += dtv * vshift[i][a] * grade[i][a];
}
}
@ -290,29 +327,34 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
void FixRHEOThermal::post_integrate()
{
int i, itype;
double cvi, Tci, Ti;
double cvi, Tci, Ti, Li;
int *status = atom->status;
double *energy = atom->esph;
double *temperature = atom->temperature;
double *heatflow = atom->heatflow;
double *rho = atom->rho;
int *type = atom->type;
int n_melt = 0;
int n_freeze = 0;
//Integrate temperature and check status
//Integrate energy and check status
for (i = 0; i < atom->nlocal; i++) {
if (status[i] & STATUS_NO_INTEGRATION) continue;
itype = type[i];
cvi = calc_cv(i, type[i]);
temperature[i] += dtf * heatflow[i] / cvi;
cvi = calc_cv(i, itype);
energy[i] += dtf * heatflow[i];
temperature[i] = energy[i] / cvi;
if (Tc_style[itype] != NONE) {
Ti = temperature[i];
if (Tc_style[itype] == CONSTANT) {
Tci = Tc[itype];
Tci = calc_Tc(i, itype);
if (L_style[itype] != NONE) {
Li = calc_L(i, itype);
if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi);
temperature[i] = Ti;
}
if (Ti > Tci) {
@ -370,11 +412,39 @@ void FixRHEOThermal::post_neighbor()
}
/* ----------------------------------------------------------------------
Calculate temperature
In the future, update & forward evolving conductivity styles every timestep
------------------------------------------------------------------------- */
void FixRHEOThermal::pre_force(int /*vflag*/)
{
int i, itype;
double cvi, Tci, Ti, Li;
double *energy = atom->esph;
double *temperature = atom->temperature;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
// Calculate temperature
for (i = 0; i < nall; i++) {
itype = type[i];
cvi = calc_cv(i, itype);
temperature[i] = energy[i] / cvi;
if (Tc_style[itype] != NONE) {
Ti = temperature[i];
Tci = calc_Tc(i, itype);
if (L_style[itype] != NONE) {
Li = calc_L(i, itype);
if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi);
temperature[i] = Ti;
}
}
}
}
/* ---------------------------------------------------------------------- */
@ -382,17 +452,13 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
void FixRHEOThermal::final_integrate()
{
int *status = atom->status;
int *type = atom->type;
double *temperature = atom->temperature;
double *energy = atom->esph;
double *heatflow = atom->heatflow;
double cvi;
//Integrate temperature and check status
//Integrate energy
for (int i = 0; i < atom->nlocal; i++) {
if (status[i] & STATUS_NO_INTEGRATION) continue;
cvi = calc_cv(i, type[i]);
temperature[i] += dtf * heatflow[i] / cvi;
energy[i] += dtf * heatflow[i];
}
}
@ -535,6 +601,23 @@ double FixRHEOThermal::calc_cv(int i, int itype)
}
}
/* ---------------------------------------------------------------------- */
double FixRHEOThermal::calc_Tc(int i, int itype)
{
if (Tc_style[itype] == CONSTANT) {
return Tc[itype];
}
}
/* ---------------------------------------------------------------------- */
double FixRHEOThermal::calc_L(int i, int itype)
{
if (L_style[itype] == CONSTANT) {
return L[itype];
}
}
/* ---------------------------------------------------------------------- */

6
src/RHEO/fix_rheo_thermal.h
View File

@ -43,12 +43,14 @@ class FixRHEOThermal : public Fix {
void unpack_forward_comm(int, int, double *) override;
void reset_dt() override;
double calc_cv(int, int);
double calc_Tc(int, int);
double calc_L(int, int);
private:
double *cv, *Tc, *kappa;
double *cv, *Tc, *kappa, *L;
double dtf, dtv;
double cut_kernel, cut_bond, cutsq_bond;
int *cv_style, *Tc_style, *kappa_style;
int *cv_style, *Tc_style, *kappa_style, *L_style;
int btype;
class NeighList *list;

1
src/RHEO/pair_rheo.cpp
View File

@ -96,7 +96,6 @@ void PairRHEO::compute(int eflag, int vflag)
ev_init(eflag, vflag);
double **gradv = compute_grad->gradv;
double **gradt = compute_grad->gradt;
double **gradr = compute_grad->gradr;
double **v = atom->v;
double **x = atom->x;

351
src/RHEO/pair_rheo_solid.cpp Normal file
View File

@ -0,0 +1,351 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "pair_rheo_solid.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix_rheo.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */
PairRHEOSolid::PairRHEOSolid(LAMMPS *_lmp) : Pair(_lmp)
{
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairRHEOSolid::~PairRHEOSolid()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(k);
memory->destroy(cut);
memory->destroy(gamma);
}
}
/* ---------------------------------------------------------------------- */
void PairRHEOSolid::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double r, rsq, rinv, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double vxtmp, vytmp, vztmp, delvx, delvy, delvz, dot, smooth;
evdwl = 0.0;
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *status = atom->status;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double *special_lj = force->special_lj;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (!(status[i] & STATUS_SOLID)) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (!(status[j] & STATUS_SOLID)) continue;
factor_lj = special_lj[sbmask(j)];
if (factor_lj == 0) continue;
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
rinv = 1.0 / r;
fpair = k[itype][jtype] * (cut[itype][jtype] - r);
smooth = rsq / cutsq[itype][jtype];
smooth *= smooth;
smooth *= smooth;
smooth = 1.0 - smooth;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
dot = delx * delvx + dely * delvy + delz * delvz;
fpair -= gamma[itype][jtype] * dot * smooth * rinv;
fpair *= factor_lj * rinv;
if (eflag) evdwl = 0.0;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairRHEOSolid::allocate()
{
allocated = 1;
const int np1 = atom->ntypes + 1;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(k, np1, np1, "pair:k");
memory->create(cut, np1, np1, "pair:cut");
memory->create(gamma, np1, np1, "pair:gamma");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairRHEOSolid::settings(int narg, char ** /*arg*/)
{
if (narg != 0) error->all(FLERR, "Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairRHEOSolid::coeff(int narg, char **arg)
{
if (narg != 5) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double k_one = utils::numeric(FLERR, arg[2], false, lmp);
double cut_one = utils::numeric(FLERR, arg[3], false, lmp);
double gamma_one = utils::numeric(FLERR, arg[4], false, lmp);
if (cut_one <= 0.0) error->all(FLERR, "Incorrect args for pair coefficients");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
k[i][j] = k_one;
cut[i][j] = cut_one;
gamma[i][j] = gamma_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairRHEOSolid::init_style()
{
if (comm->ghost_velocity == 0)
error->all(FLERR,"Pair rheo/solid requires ghost atoms store velocity");
if (!atom->status_flag)
error->all(FLERR,"Pair rheo/solid requires atom_style rheo");
neighbor->add_request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairRHEOSolid::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]);
gamma[i][j] = mix_energy(gamma[i][i], gamma[j][j], cut[i][i], cut[j][j]);
}
cut[j][i] = cut[i][j];
k[j][i] = k[i][j];
gamma[j][i] = gamma[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairRHEOSolid::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&k[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
fwrite(&gamma[i][j], sizeof(double), 1, fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairRHEOSolid::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairRHEOSolid::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp, "%d %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairRHEOSolid::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp, "%d %d %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairRHEOSolid::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/,
double factor_lj, double &fforce)
{
double fpair, r, rinv;
double delx, dely, delz, delvx, delvy, delvz, dot, smooth;
if (rsq > cutsq[itype][jtype]) return 0.0;
double **x = atom->x;
double **v = atom->v;
r = sqrt(rsq);
rinv = 1.0 / r;
fpair = k[itype][jtype] * (cut[itype][jtype] - r);
smooth = rsq / cutsq[itype][jtype];
smooth *= smooth;
smooth = 1.0 - smooth;
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
delvx = v[i][0] - v[j][0];
delvy = v[i][1] - v[j][1];
delvz = v[i][2] - v[j][2];
dot = delx * delvx + dely * delvy + delz * delvz;
fpair -= gamma[itype][jtype] * dot * rinv * smooth;
fpair *= factor_lj;
fforce = fpair;
return 0.0;
}

51
src/RHEO/pair_rheo_solid.h Normal file
View File

@ -0,0 +1,51 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(rheo/solid,PairRHEOSolid);
// clang-format on
#else
#ifndef LMP_PAIR_RHEO_SOLID_H
#define LMP_PAIR_RHEO_SOLID_H
#include "pair.h"
namespace LAMMPS_NS {
class PairRHEOSolid : public Pair {
public:
PairRHEOSolid(class LAMMPS *);
~PairRHEOSolid() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
double single(int, int, int, int, double, double, double, double &) override;
protected:
double **k, **cut, **gamma;
void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif