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

Adding output option for reorganized peratom values

Browse Source
This commit is contained in:
jtclemm
2023-04-25 16:29:29 -06:00
parent 7fc916a1d4
commit be568d257d
27 changed files with 1443 additions and 449 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
src/RHEO/atom_vec_rheo.cpp
View File

@ -33,15 +33,17 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp)
forceclearflag = 1;
atom->status_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
atom->viscosity_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// 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"};
fields_copy = {"status", "rho", "drho"};
fields_grow = {"status", "rho", "drho", "pressure", "viscosity"};
fields_copy = {"status", "rho", "drho", "pressure", "viscosity"};
fields_comm = {"status", "rho"};
fields_comm_vel = {"status", "rho"};
fields_reverse = {"drho"};
@ -49,7 +51,7 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp)
fields_border_vel = {"status", "rho"};
fields_exchange = {"status", "rho"};
fields_restart = {"status", "rho"};
fields_create = {"status", "rho", "drho"};
fields_create = {"status", "rho", "drho", "pressure", "viscosity"};
fields_data_atom = {"id", "type", "status", "rho", "x"};
fields_data_vel = {"id", "v"};
@ -64,8 +66,10 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp)
void AtomVecRHEO::grow_pointers()
{
status = atom->status;
pressure = atom->pressure;
rho = atom->rho;
drho = atom->drho;
viscosity = atom->viscosity;
}
/* ----------------------------------------------------------------------
@ -86,6 +90,8 @@ void AtomVecRHEO::force_clear(int n, size_t nbytes)
void AtomVecRHEO::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
@ -96,8 +102,10 @@ void AtomVecRHEO::data_atom_post(int ilocal)
int AtomVecRHEO::property_atom(const std::string &name)
{
if (name == "status") return 0;
if (name == "rho") return 1;
if (name == "drho") return 2;
if (name == "pressure") return 1;
if (name == "rho") return 2;
if (name == "drho") return 3;
if (name == "viscosity") return 4;
return -1;
}
@ -123,12 +131,20 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
buf[n] = pressure[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
@ -136,5 +152,13 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

2
src/RHEO/atom_vec_rheo.h
View File

@ -36,7 +36,7 @@ class AtomVecRHEO : virtual public AtomVec {
private:
int *status;
double *rho, *drho;
double *pressure, *rho, *drho, *viscosity;
};
} // namespace LAMMPS_NS

200
src/RHEO/atom_vec_rheo_thermal.cpp Normal file
View File

@ -0,0 +1,200 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
----------------------------------------------------------------------- */
#include "atom_vec_rheo_thermal.h"
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::ATOMIC;
mass_type = PER_TYPE;
forceclearflag = 1;
atom->status_flag = 1;
atom->conductivity_flag = 1;
atom->temperature_flag = 1;
atom->heatflow_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
atom->viscosity_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// 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_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_data_vel = {"id", "v"};
setup_fields();
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::grow_pointers()
{
status = atom->status;
conductivity = atom->conductivity;
temperature = atom->temperature;
heatflow = atom->heatflow;
pressure = atom->pressure;
rho = atom->rho;
drho = atom->drho;
viscosity = atom->viscosity;
}
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::force_clear(int n, size_t nbytes)
{
memset(&drho[n], 0, nbytes);
memset(&heatflow[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
heatflow[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
conductivity[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
assign an index to named atom property and return index
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecRHEOThermal::property_atom(const std::string &name)
{
if (name == "status") return 0;
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;
return -1;
}
/* ----------------------------------------------------------------------
pack per-atom data into buf for ComputePropertyAtom
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = status[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = temperature[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = heatflow[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 5) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = conductivity[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 6) {
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 == 7) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

46
src/RHEO/atom_vec_rheo_thermal.h Normal file
View File

@ -0,0 +1,46 @@
/* -*- 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 ATOM_CLASS
// clang-format off
AtomStyle(rheo/thermal,AtomVecRHEOThermal);
// clang-format on
#else
#ifndef LMP_ATOM_VEC_RHEO_THERMAL_H
#define LMP_ATOM_VEC_RHEO_THERMAL_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecRHEOThermal : virtual public AtomVec {
public:
AtomVecRHEOThermal(class LAMMPS *);
void grow_pointers() override;
void force_clear(int, size_t) override;
void data_atom_post(int) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:
int *status;
double *conductivity, *temperature, *heatflow;
double *pressure, *rho, *drho, *viscosity;
};
} // namespace LAMMPS_NS
#endif
#endif

85
src/RHEO/compute_rheo_grad.cpp
View File

@ -84,69 +84,34 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
}
comm_forward = ncomm_grad;
nmax_store = 0;
grow_arrays(atom->nmax);
}
/* ---------------------------------------------------------------------- */
ComputeRHEOGrad::~ComputeRHEOGrad()
{
int dim = domain->dimension;
int tmp1, tmp2, index;
index = atom->find_custom("rheo_grad_v", tmp1, tmp2);
if (index != 1) atom->remove_custom(index, 1, dim * dim);
index = atom->find_custom("rheo_grad_rho", tmp1, tmp2);
if (index != 1) atom->remove_custom(index, 1, dim);
index = atom->find_custom("rheo_grad_t", tmp1, tmp2);
if (index != 1) atom->remove_custom(index, 1, dim);
index = atom->find_custom("rheo_grad_eta", tmp1, tmp2);
if (index != 1) atom->remove_custom(index, 1, dim);
memory->destroy(gradv);
memory->destroy(gradr);
memory->destroy(gradt);
memory->destroy(gradn);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOGrad::init()
{
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
cut = fix_rheo->cut;
cutsq = cut * cut;
rho0 = fix_rheo->rho0;
interface_flag = fix_rheo->interface_flag;
compute_kernel = fix_rheo->compute_kernel;
compute_interface = fix_rheo->compute_interface;
int tmp1, tmp2;
index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
// Create coordination array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int index;
int dim = domain->dimension;
if (velocity_flag) {
index = atom->add_custom("rheo_grad_v", 1, dim * dim);
gradv = atom->darray[index];
}
if (rho_flag) {
index = atom->add_custom("rheo_grad_rho", 1, dim);
gradr = atom->darray[index];
}
if (temperature_flag) {
index= atom->add_custom("rheo_grad_temp", 1, dim);
gradt = atom->darray[index];
}
if (eta_flag) {
index = atom->add_custom("rheo_grad_eta", 1, dim);
gradn = atom->darray[index];
}
nmax_store = 0;
grow_arrays(atom->nmax);
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
}
/* ---------------------------------------------------------------------- */
@ -175,7 +140,7 @@ void ComputeRHEOGrad::compute_peratom()
double **v = atom->v;
double *rho = atom->rho;
double *temperature = atom->temperature;
double *viscosity = atom->dvector[index_visc];
double *viscosity = atom->viscosity;
int *status = atom->status;
int *type = atom->type;
double *mass = atom->mass;
@ -240,15 +205,17 @@ void ComputeRHEOGrad::compute_peratom()
vj[2] = v[j][2];
// Add corrections for walls
if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j, i);
} else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) {
compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i, j);
} else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
rhoi = rho0;
rhoj = rho0;
if (interface_flag) {
if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j, i);
} else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) {
compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i, j);
} else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
rhoi = rho0;
rhoj = rho0;
}
}
Voli = mass[itype] / rhoi;
@ -481,16 +448,16 @@ void ComputeRHEOGrad::grow_arrays(int nmax)
{
int dim = domain->dimension;
if (velocity_flag)
memory->grow(gradv, nmax, dim * dim, "atom:rheo_grad_v");
memory->grow(gradv, nmax, dim * dim, "rheo:grad_v");
if (rho_flag)
memory->grow(gradr, nmax, dim, "atom:rheo_grad_rho");
memory->grow(gradr, nmax, dim, "rheo:grad_rho");
if (temperature_flag)
memory->grow(gradt, nmax, dim, "atom:rheo_grad_temp");
memory->grow(gradt, nmax, dim, "rheo:grad_temp");
if (eta_flag)
memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta");
memory->grow(gradn, nmax, dim, "rheo:grad_eta");
nmax_store = nmax;
}

8
src/RHEO/compute_rheo_grad.h
View File

@ -46,14 +46,14 @@ class ComputeRHEOGrad : public Compute {
private:
int comm_stage, ncomm_grad, ncomm_field, nmax_store;
int index_visc;
double cut, cutsq, rho0;
class NeighList *list;
int velocity_flag, temperature_flag, rho_flag, eta_flag;
int interface_flag;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOInterface *compute_interface;
int velocity_flag, temperature_flag, rho_flag, eta_flag;
class NeighList *list;
void grow_arrays(int);
};

66
src/RHEO/compute_rheo_interface.cpp
View File

@ -46,24 +46,35 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
nmax_store = 0;
comm_forward = 3;
comm_reverse = 4;
nmax_store = atom->nmax;
memory->create(chi, nmax_store, "rheo:chi");
memory->create(norm, nmax_store, "rheo/interface:norm");
memory->create(normwf, nmax_store, "rheo/interface:normwf");
// For fp_store, create an instance of fix property atom
// Need restarts + exchanging with neighbors since it needs to persist
// between timesteps (fix property atom will handle callbacks)
int tmp1, tmp2;
int index = atom->find_custom("fp_store", tmp1, tmp2);
if (index == -1) {
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
index = atom->find_custom("fp_store", tmp1, tmp2);
}
fp_store = atom->darray[index];
}
/* ---------------------------------------------------------------------- */
ComputeRHEOInterface::~ComputeRHEOInterface()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_chi", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 0);
if (id_fix_pa && modify->nfix) modify->delete_fix(id_fix_pa);
delete[] id_fix_pa;
memory->destroy(chi);
memory->destroy(norm);
memory->destroy(normwf);
}
@ -80,37 +91,6 @@ void ComputeRHEOInterface::init()
cutsq = cut * cut;
wall_max = sqrt(3.0) / 12.0 * cut;
// Create chi array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int tmp1, tmp2;
int nmax = atom->nmax;
int index = atom->find_custom("rheo_chi", tmp1, tmp2);
if (index == -1) {
index = atom->add_custom("rheo_chi", 1, 0);
memory->destroy(norm);
memory->destroy(normwf);
memory->create(norm, nmax, "rheo/interface:norm");
memory->create(normwf, nmax, "rheo/interface:normwf");
nmax_store = nmax;
}
chi = atom->dvector[index];
// For fp_store, go ahead and create an instance of fix property atom
// Need restarts + exchanging with neighbors since it needs to persist
// between timesteps (fix property atom will handle callbacks)
index = atom->find_custom("fp_store", tmp1, tmp2);
if (index == -1) {
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
index = atom->find_custom("fp_store", tmp1, tmp2);
}
fp_store = atom->darray[index];
// need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
}
@ -145,11 +125,9 @@ void ComputeRHEOInterface::compute_peratom()
if (atom->nmax > nmax_store) {
nmax_store = atom->nmax;
memory->destroy(norm);
memory->destroy(normwf);
memory->create(norm, nmax_store, "rheo/interface:norm");
memory->create(normwf, nmax_store, "rheo/interface:normwf");
memory->grow(chi, nmax_store, "rheo/interface:chi");
memory->grow(norm, nmax_store, "rheo/interface:norm");
memory->grow(normwf, nmax_store, "rheo/interface:normwf");
memory->grow(chi, nmax_store, "rheo:chi");
}
for (i = 0; i < nall; i++) {

53
src/RHEO/compute_rheo_kernel.cpp
View File

@ -68,7 +68,6 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
correction_order = 2;
}
solid_flag = 0;
dim = domain->dimension;
comm_forward = 1;
@ -93,11 +92,7 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
ComputeRHEOKernel::~ComputeRHEOKernel()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_coordination", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 0, 0);
memory->destroy(coordination);
memory->destroy(C);
memory->destroy(C0);
}
@ -112,11 +107,8 @@ void ComputeRHEOKernel::init()
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/kernel");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
int icompute = modify->find_compute("rheo_interface");
if (icompute != -1) {
compute_interface = ((ComputeRHEOInterface *) modify->compute[icompute]);
solid_flag = 1;
}
interface_flag = fix_rheo->interface_flag;
compute_interface = fix_rheo->compute_interface;
zmin = fix_rheo->zmin_kernel;
h = fix_rheo->h;
@ -133,22 +125,8 @@ void ComputeRHEOKernel::init()
pre_wp = pre_w * 3.0 * hinv;
}
// Create coordination array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int tmp1, tmp2;
int nmax = atom->nmax;
int index = atom->find_custom("rheo_coordination", tmp1, tmp2);
if (index == -1) {
index = atom->add_custom("rheo_coordination", 0, 0);
nmax_store = nmax;
}
coordination = atom->ivector[index];
// Create local arrays for kernel arrays, I can't foresee a reason to print
nmax_store = atom->nmax;
memory->create(coordination, nmax_store, "rheo:coordination");
if (kernel_style == CRK0) {
memory->create(C0, nmax_store, "rheo/kernel:C0");
} else if (kernel_style == CRK1) {
@ -499,7 +477,7 @@ void ComputeRHEOKernel::compute_peratom()
if (kernel_style == QUINTIC) return;
int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error;
double xtmp, ytmp, ztmp, r, rsq, w, vj;
double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj;
double dx[3];
gsl_matrix_view gM;
@ -549,10 +527,12 @@ void ComputeRHEOKernel::compute_peratom()
if (rsq < hsq) {
r = sqrt(rsq);
w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
if (!(status[j] & STATUS_FLUID) && solid_flag) {
vj = mass[type[j]] / compute_interface->correct_rho(j,i);
} else vj = mass[type[j]] / rho[j];
rhoj = rho[j];
if (interface_flag)
if (!(status[j] & STATUS_FLUID))
rhoj = compute_interface->correct_rho(j,i);
vj = mass[type[j]] / rhoj;
M += w * vj;
}
}
@ -596,11 +576,12 @@ void ComputeRHEOKernel::compute_peratom()
r = sqrt(rsq);
w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
if (solid_flag)
rhoj = rho[j];
if (interface_flag)
if (!(status[j] & STATUS_FLUID))
vj = mass[type[j]]/compute_interface->correct_rho(j,i);
else
vj = mass[type[j]]/rho[j];
rhoj = compute_interface->correct_rho(j,i);
vj = mass[type[j]] / rhoj;
//Populate the H-vector of polynomials (2D)
if (dim == 2) {
@ -759,7 +740,7 @@ void ComputeRHEOKernel::compute_coordination()
void ComputeRHEOKernel::grow_arrays(int nmax)
{
memory->grow(coordination, nmax, "atom:rheo_coordination");
memory->grow(coordination, nmax, "rheo:coordination");
if (kernel_style == CRK0) {
memory->grow(C0, nmax, "rheo/kernel:C0");

2
src/RHEO/compute_rheo_kernel.h
View File

@ -49,7 +49,7 @@ class ComputeRHEOKernel : public Compute {
private:
int comm_stage, comm_forward_save;
int solid_flag;
int interface_flag;
int gsl_error_flag;
std::unordered_set<tagint> gsl_error_tags;

293
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -21,20 +21,29 @@
#include "atom.h"
#include "atom_vec.h"
#include "compute_rheo_interface.h"
#include "compute_rheo_kernel.h"
#include "compute_rheo_surface.h"
#include "compute_rheo_vshift.h"
#include "error.h"
#include "fix_rheo.h"
#include "fix_rheo_thermal.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */
ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
index(nullptr), colindex(nullptr), pack_choice(nullptr)
Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr),
compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr),
index(nullptr), pack_choice(nullptr)
{
if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error);
@ -43,31 +52,66 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
thermal_flag, interface_flag, surface_flag, shift_flag = 0;
// parse input values
// customize a new keyword by adding to if statement
pack_choice = new FnPtrPack[nvalues];
index = new int[nvalues];
int i;
for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3;
if (strcmp(arg[iarg],"id") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_id;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (!atom->molecule_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_molecule;
} else if (strcmp(arg[iarg],"proc") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_proc;
} else if (strcmp(arg[iarg],"type") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_type;
} else if (strcmp(arg[iarg],"mass") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_mass;
if (strcmp(arg[iarg],"phase") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase;
} else if (strcmp(arg[iarg],"chi") == 0) {
interface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi;
} else if (strcmp(arg[iarg],"surface") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface;
} else if (strcmp(arg[iarg],"surface/r") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r;
} else if (strcmp(arg[iarg],"surface/divr") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
} else if (strcmp(arg[iarg],"surface/nx") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nx;
} else if (strcmp(arg[iarg],"surface/ny") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_ny;
} else if (strcmp(arg[iarg],"surface/nz") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nz;
} else if (strcmp(arg[iarg],"coordination") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
} else if (strcmp(arg[iarg],"cv") == 0) {
thermal_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv;
} else if (strcmp(arg[iarg],"shift/vx") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vx;
} else if (strcmp(arg[iarg],"shift/vy") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vy;
} else if (strcmp(arg[iarg],"shift/vz") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vz;
} else {
error->all(FLERR,"Invalid keyword {} for compute rheo/property/atom command ", arg[iarg]);
index[i] = atom->avec->property_atom(arg[iarg]);
if (index[i] < 0)
error->all(FLERR,
"Invalid keyword {} for atom style {} in compute rheo/property/atom command ",
atom->get_style(), arg[iarg]);
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"heatflow") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"conductivity") == 0) thermal_flag = 1;
}
}
@ -79,12 +123,41 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
{
delete[] pack_choice;
delete[] index;
memory->destroy(vector_atom);
memory->destroy(array_atom);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::init()
{
auto fixes = modify->get_fix_by_style("rheo");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
if (interface_flag && !(fix_rheo->interface_flag))
error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo");
if (surface_flag && !(fix_rheo->surface_flag))
error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo");
if (shift_flag && !(fix_rheo->shift_flag))
error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo");
if (thermal_flag && !(fix_rheo->thermal_flag))
error->all(FLERR, "Cannot request thermal property without fix rheo/thermal");
compute_interface = fix_rheo->compute_interface;
compute_kernel = fix_rheo->compute_kernel;
compute_surface = fix_rheo->compute_surface;
compute_vshift = fix_rheo->compute_vshift;
if (thermal_flag) {
fixes = modify->get_fix_by_style("rheo/thermal");
fix_thermal = dynamic_cast<FixRHEOThermal *>(fixes[0]);
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
@ -133,14 +206,16 @@ double ComputeRHEOPropertyAtom::memory_usage()
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_rho(int n)
void ComputeRHEOPropertyAtom::pack_phase(int n)
{
double *rho = atom->rho;
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int inverse_mask = ~PHASEMASK;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = rho[i];
if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask);
else buf[n] = 0.0;
n += nvalues;
}
@ -148,15 +223,189 @@ void ComputeRHEOPropertyAtom::pack_rho(int n)
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_status(int n)
void ComputeRHEOPropertyAtom::pack_chi(int n)
{
double *chi = compute_interface->chi;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = chi[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface(int n)
{
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int inverse_mask = ~SURFACEMASK;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = status[i];
if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_r(int n)
{
double *rsurface = compute_surface->rsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = rsurface[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_divr(int n)
{
double *divr = compute_surface->divr;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = divr[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nx(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][0];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_ny(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nz(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_coordination(int n)
{
int *coordination = compute_kernel->coordination;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = coordination[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_cv(int n)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vx(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][0];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vy(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vz(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_atom_style(int n)
{
atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit);
}

39
src/RHEO/compute_rheo_property_atom.h
View File

@ -28,37 +28,40 @@ class ComputeRHEOPropertyAtom : public Compute {
public:
ComputeRHEOPropertyAtom(class LAMMPS *, int, char **);
~ComputeRHEOPropertyAtom() override;
void init() override;
void compute_peratom() override;
double memory_usage() override;
private:
int nvalues;
int nmax;
int nvalues, nmax;
int thermal_flag, interface_flag, surface_flag, shift_flag;
int *index;
double *buf;
typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_rho(int);
void pack_drho(int);
void pack_temperature(int);
void pack_heatflow(int);
void pack_status(int);
void pack_phase(int);
void pack_chi(int);
void pack_surface(int);
void pack_r_surface(int);
void pack_divr_surface(int);
void pack_nx_surface(int);
void pack_ny_surface(int);
void pack_nz_surface(int);
void pack_surface_r(int);
void pack_surface_divr(int);
void pack_surface_nx(int);
void pack_surface_ny(int);
void pack_surface_nz(int);
void pack_coordination(int);
void pack_viscosity(int);
void pack_pressure(int);
void pack_conductivity(int);
void pack_cv(int);
void pack_vx_shift(int);
void pack_vy_shift(int);
void pack_vz_shift(int);
void pack_shift_vx(int);
void pack_shift_vy(int);
void pack_shift_vz(int);
void pack_atom_style(int);
class FixRHEO *fix_rheo;
class FixRHEOThermal *fix_thermal;
class ComputeRHEOInterface *compute_interface;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSurface *compute_surface;
class ComputeRHEOVShift *compute_vshift;
};

4
src/RHEO/compute_rheo_surface.h
View File

@ -36,13 +36,13 @@ class ComputeRHEOSurface : public Compute {
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
double **nsurface, *rsurface;
double **nsurface, *rsurface, *divr;
class FixRHEO *fix_rheo;
private:
double cut, cutsq, rho0, threshold_divr;
int surface_style, nmax_store, threshold_z;
double **B, **gradC, *divr;
double **B, **gradC;
int threshold_style, comm_stage;
class NeighList *list;

69
src/RHEO/compute_rheo_vshift.cpp
View File

@ -22,6 +22,7 @@
#include "comm.h"
#include "compute_rheo_interface.h"
#include "compute_rheo_kernel.h"
#include "compute_rheo_surface.h"
#include "domain.h"
#include "error.h"
#include "fix_rheo.h"
@ -38,36 +39,22 @@ using namespace RHEO_NS;
ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), vshift(nullptr), fix_rheo(nullptr),
compute_kernel(nullptr), compute_interface(nullptr)
compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr)
{
if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");
comm_reverse = 3;
surface_flag = 0;
// Create vshift array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int tmp1, tmp2;
int index = atom->find_custom("rheo_vshift", tmp1, tmp2);
if (index == -1) {
index = atom->add_custom("rheo_vshift", 1, 3);
nmax_store = atom->nmax;
}
vshift = atom->darray[index];
nmax_store = atom->nmax;
memory->create(vshift, nmax_store, 3, "rheo:vshift");
}
/* ---------------------------------------------------------------------- */
ComputeRHEOVShift::~ComputeRHEOVShift()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_vshift", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 3);
memory->destroy(vshift);
}
/* ---------------------------------------------------------------------- */
@ -76,12 +63,12 @@ void ComputeRHEOVShift::init()
{
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
surface_flag = 0;
if (fix_rheo->surface_flag)
surface_flag = 1;
surface_flag = fix_rheo->surface_flag;
interface_flag = fix_rheo->interface_flag;
compute_kernel = fix_rheo->compute_kernel;
compute_interface = fix_rheo->compute_interface;
compute_surface = fix_rheo->compute_surface;
cut = fix_rheo->cut;
cutsq = cut * cut;
@ -128,7 +115,7 @@ void ComputeRHEOVShift::compute_peratom()
firstneigh = list->firstneigh;
if (nmax_store < atom->nmax) {
memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift");
memory->grow(vshift, atom->nmax, 3, "rheo:vshift");
nmax_store = atom->nmax;
}
@ -176,15 +163,17 @@ void ComputeRHEOVShift::compute_peratom()
rhoj = rho[j];
// Add corrections for walls
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j,i);
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i,j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;
rhoj = 1.0;
if (interface_flag) {
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j,i);
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i,j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;
rhoj = 1.0;
}
}
voli = imass / rhoi;
@ -227,30 +216,28 @@ void ComputeRHEOVShift::correct_surfaces()
{
if (!surface_flag) return;
int i, a, b;
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int i, a, b;
int dim = domain->dimension;
double **nsurface = compute_surface->nsurface;
int tmp1, tmp2;
int index_nsurf = atom->find_custom("rheo_nsurf", tmp1, tmp2);
if (index_nsurf == -1) error->all(FLERR, "Cannot find rheo nsurf");
double **nsurf = atom->darray[index_nsurf];
int nlocal = atom->nlocal;
int dim = domain->dimension;
double nx,ny,nz,vx,vy,vz;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
nx = nsurf[i][0];
ny = nsurf[i][1];
nx = nsurface[i][0];
ny = nsurface[i][1];
vx = vshift[i][0];
vy = vshift[i][1];
vz = vshift[i][2];
vshift[i][0] = (1 - nx * nx) * vx - nx * ny * vy;
vshift[i][1] = (1 - ny * ny) * vy - nx * ny * vx;
if (dim > 2) {
nz = nsurf[i][2];
nz = nsurface[i][2];
vshift[i][0] -= nx * nz * vz;
vshift[i][1] -= ny * nz * vz;
vshift[i][2] = (1 - nz * nz) * vz - nz * ny * vy - nx * nz * vx;

5
src/RHEO/compute_rheo_vshift.h
View File

@ -42,11 +42,12 @@ class ComputeRHEOVShift : public Compute {
private:
int nmax_store;
double dtv, cut, cutsq, cutthird;
int surface_flag;
int surface_flag, interface_flag;
class NeighList *list;
class ComputeRHEOInterface *compute_interface ;
class ComputeRHEOInterface *compute_interface;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSurface *compute_surface;
};
} // namespace LAMMPS_NS

4
src/RHEO/fix_rheo.cpp
View File

@ -60,8 +60,12 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
if (igroup != 0)
error->all(FLERR,"fix rheo command requires group all");
if (atom->pressure_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with pressure");
if (atom->rho_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with density");
if (atom->viscosity_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with viscosity");
if (atom->status_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with status");

2
src/RHEO/fix_rheo_pressure.cpp
View File

@ -38,7 +38,7 @@ static constexpr double SEVENTH = 1.0 / 7.0;
/* ---------------------------------------------------------------------- */
FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fix_rheo(nullptr), pressure(nullptr)
Fix(lmp, narg, arg), fix_rheo(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");

141
src/RHEO/fix_rheo_thermal.cpp
View File

@ -40,7 +40,7 @@ enum {NONE, CONSTANT, TYPE};
FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr),
Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr), conductivity(nullptr)
Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");
@ -48,9 +48,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
cv_style = NONE;
conductivity_style = NONE;
comm_forward = 1;
nmax_store = 0;
int ntypes = atom->ntypes;
int iarg = 3;
while (iarg < narg) {
@ -170,9 +167,11 @@ void FixRHEOThermal::init()
dtf = 0.5 * update->dt * force->ftm2v;
if (atom->temperature_flag != 1)
error->all(FLERR,"fix rheo/thermal command requires atoms store temperature property");
error->all(FLERR,"fix rheo/thermal command requires atom property temperature");
if (atom->heatflow_flag != 1)
error->all(FLERR,"fix rheo/thermal command requires atoms store heatflow property");
error->all(FLERR,"fix rheo/thermal command requires atom property heatflow");
if (atom->conductivity_flag != 1)
error->all(FLERR,"fix rheo/thermal command requires atom property conductivity");
}
/* ---------------------------------------------------------------------- */
@ -181,34 +180,6 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
{
fix_rheo->thermal_fix_defined = 1;
// Identify whether this is the first/last instance of fix thermal
// First will grow arrays, last will communicate
first_flag = 0;
last_flag = 0;
int i = 0;
auto fixlist = modify->get_fix_by_style("rheo/thermal");
for (const auto &fix : fixlist) {
if (strcmp(fix->id, id) == 0) break;
i++;
}
if (i == 0) first_flag = 1;
if ((i + 1) == fixlist.size()) last_flag = 1;
// Create conductivity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int tmp1, tmp2;
int index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index== -1) {
index = atom->add_custom("rheo_conductivity", 1, 0);
nmax_store = atom->nmax;
}
conductivity = atom->dvector[index];
post_neighbor();
pre_force(0);
}
@ -294,14 +265,10 @@ void FixRHEOThermal::post_neighbor()
int i;
int *type = atom->type;
int *mask = atom->mask;
double *conductivity = atom->conductivity;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
nmax_store = atom->nmax;
}
if (conductivity_style == CONSTANT) {
for (i = 0; i < nall; i++)
if (mask[i] & groupbit) conductivity[i] = kappa;
@ -312,39 +279,11 @@ void FixRHEOThermal::post_neighbor()
}
/* ----------------------------------------------------------------------
Update (and forward) evolving conductivity styles every timestep
Zero heat flow
In the future, update & forward evolving conductivity styles every timestep
------------------------------------------------------------------------- */
void FixRHEOThermal::pre_force(int /*vflag*/)
{
// send updated temperatures to ghosts if first instance of fix
// then clear heatflow for next force calculation
double *heatflow = atom->heatflow;
if (first_flag) {
comm->forward_comm(this);
for (int i = 0; i < atom->nmax; i++) heatflow[i] = 0.0;
}
// Not needed yet, when needed add stage check for (un)pack_forward_comm() methods
//int i;
//double *conductivity = atom->dvector[index_cond];
//int *mask = atom->mask;
//int nlocal = atom->nlocal;
//if (first_flag && (nmax_store < atom->nmax)) {
// memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
// nmax_store = atom->nmax;
//}
//if (conductivity_style == TBD) {
// for (i = 0; i < nlocal; i++) {
// if (mask[i] & groupbit) {
// }
// }
//}
//if (last_flag && comm_forward) comm->forward_comm(this);
}
/* ---------------------------------------------------------------------- */
@ -387,69 +326,3 @@ double FixRHEOThermal::calc_cv(int i)
return(cv_type[atom->type[i]]);
}
}
/* ---------------------------------------------------------------------- */
int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i, j, m;
double *temperature = atom->temperature;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = temperature[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRHEOThermal::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
double *temperature = atom->temperature;
for (i = first; i < last; i++) temperature[i] = buf[m++];
}
/* ---------------------------------------------------------------------- */
int FixRHEOThermal::pack_reverse_comm(int n, int first, double *buf)
{
int m = 0;
int last = first + n;
double *heatflow = atom->heatflow;
for (int i = first; i < last; i++) {
buf[m++] = heatflow[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRHEOThermal::unpack_reverse_comm(int n, int *list, double *buf)
{
int m = 0;
double *heatflow = atom->heatflow;
for (int i = 0; i < n; i++)
heatflow[list[i]] += buf[m++];
}
/* ---------------------------------------------------------------------- */
double FixRHEOThermal::memory_usage()
{
double bytes = 0.0;
bytes += (size_t) nmax_store * sizeof(double);
return bytes;
}

11
src/RHEO/fix_rheo_thermal.h
View File

@ -37,29 +37,22 @@ class FixRHEOThermal : public Fix {
void pre_force(int) override;
void final_integrate() override;
void reset_dt() override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
double memory_usage() override;
double calc_cv(int);
private:
double *cv_type, cv;
double *Tc_type, Tc;
double *kappa_type, kappa;
double dtf, dtv;
double *conductivity;
int Tc_style;
int cv_style;
int conductivity_style;
int first_flag, last_flag;
int nmax_store;
class FixRHEO *fix_rheo;
class ComputeRHEOGrad *compute_grad;
class ComputeRHEOVShift *compute_vshift;
double calc_cv(int);
void grow_array(int);
};
} // namespace LAMMPS_NS

54
src/RHEO/fix_rheo_viscosity.cpp
View File

@ -37,14 +37,13 @@ enum {NONE, CONSTANT, TYPE, POWER};
/* ---------------------------------------------------------------------- */
FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fix_rheo(nullptr), eta_type(nullptr), viscosity(nullptr), compute_grad(nullptr)
Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), eta_type(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");
viscosity_style = NONE;
comm_forward = 0;
nmax_store = 0;
int ntypes = atom->ntypes;
int iarg = 3;
@ -86,11 +85,6 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
FixRHEOViscosity::~FixRHEOViscosity()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 0);
memory->destroy(eta_type);
}
@ -121,9 +115,7 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
{
fix_rheo->viscosity_fix_defined = 1;
// Identify whether this is the first/last instance of fix viscosity
// First will grow arrays, last will communicate
first_flag = 0;
// Identify whether this is the last instance of fix viscosity
last_flag = 0;
int i = 0;
@ -133,22 +125,8 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
i++;
}
if (i == 0) first_flag = 1;
if ((i + 1) == fixlist.size()) last_flag = 1;
// Create viscosity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_store exceeded
int tmp1, tmp2;
int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index == -1) {
index = atom->add_custom("rheo_viscosity", 1, 0);
nmax_store = atom->nmax;
}
viscosity = atom->dvector[index];
post_neighbor();
pre_force(0);
}
@ -163,14 +141,9 @@ void FixRHEOViscosity::post_neighbor()
int *type = atom->type;
int *mask = atom->mask;
double *viscosity = atom->viscosity;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_store = atom->nmax;
}
int nall = atom->nlocal + atom->nghost;
if (viscosity_style == CONSTANT) {
for (i = 0; i < nall; i++)
@ -191,16 +164,12 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
double tmp, gdot;
int *mask = atom->mask;
double *viscosity = atom->viscosity;
double **gradv = compute_grad->gradv;
int nlocal = atom->nlocal;
int dim = domain->dimension;
if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_store = atom->nmax;
}
if (viscosity_style == POWER) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -231,7 +200,8 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
int i, j, k, m;
double *viscosity = atom->viscosity;
m = 0;
for (i = 0; i < n; i++) {
@ -246,6 +216,7 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, m, last;
double *viscosity = atom->viscosity;
m = 0;
last = first + n;
@ -253,12 +224,3 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
viscosity[i] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
double FixRHEOViscosity::memory_usage()
{
double bytes = 0.0;
bytes += (size_t) nmax_store * sizeof(double);
return bytes;
}

6
src/RHEO/fix_rheo_viscosity.h
View File

@ -35,15 +35,11 @@ class FixRHEOViscosity : public Fix {
void pre_force(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
double memory_usage() override;
private:
double *eta_type, eta;
double npow, K, gd0, tau0;
double *viscosity;
int viscosity_style;
int first_flag, last_flag;
int nmax_store;
int viscosity_style, last_flag;
class FixRHEO *fix_rheo;
class ComputeRHEOGrad *compute_grad;

19
src/RHEO/pair_rheo.cpp
View File

@ -99,6 +99,9 @@ void PairRHEO::compute(int eflag, int vflag)
double *rho = atom->rho;
double *mass = atom->mass;
double *drho = atom->drho;
double *pressure = atom->pressure;
double *viscosity = atom->viscosity;
double *conductivity = atom->conductivity;
double *temperature = atom->temperature;
double *heatflow = atom->heatflow;
double *special_lj = force->special_lj;
@ -112,22 +115,6 @@ void PairRHEO::compute(int eflag, int vflag)
chi = compute_interface->chi;
}
int tmp1, tmp2;
int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index == -1) error->all(FLERR, "Cannot find rheo viscosity");
double *viscosity = atom->dvector[index];
index = atom->find_custom("rheo_pressure", tmp1, tmp2);
if (index == -1) error->all(FLERR, "Cannot find rheo pressure");
double *pressure = atom->dvector[index];
double *conductivity;
if (thermal_flag) {
index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index == -1) error->all(FLERR, "Cannot find rheo conductivity");
conductivity = atom->dvector[index];
}
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;

14
src/atom.cpp
View File

@ -200,6 +200,9 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp)
// RHEO package
status = nullptr;
conductivity = nullptr;
pressure = nullptr;
viscosity = nullptr;
// SPH package
@ -533,6 +536,9 @@ void Atom::peratom_create()
// RHEO package
add_peratom("status",&status,INT,0);
add_peratom("conductivity",&conductivity,DOUBLE,0);
add_peratom("pressure",&pressure,DOUBLE,0);
add_peratom("viscosity",&viscosity,DOUBLE,0);
// SPH package
@ -639,7 +645,7 @@ void Atom::set_atomflag_defaults()
temperature_flag = heatflow_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
status_flag = 0;
status_flag = conductivity_flag = pressure_flag = viscosity_flag = 0;
rho_flag = esph_flag = cv_flag = vest_flag = 0;
dpd_flag = edpd_flag = tdpd_flag = 0;
sp_flag = 0;
@ -2943,6 +2949,9 @@ void *Atom::extract(const char *name)
// RHEO package
if (strcmp(name,"status") == 0) return (void *) status;
if (strcmp(name,"conductivity") == 0) return (void *) conductivity;
if (strcmp(name,"pressure") == 0) return (void *) pressure;
if (strcmp(name,"viscosity") == 0) return (void *) viscosity;
// SPH package
@ -3069,6 +3078,9 @@ int Atom::extract_datatype(const char *name)
// RHEO package
if (strcmp(name,"status") == 0) return LAMMPS_INT;
if (strcmp(name,"conductivity") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"pressure") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"viscosity") == 0) return LAMMPS_DOUBLE;
// SPH package

5
src/atom.h
View File

@ -158,6 +158,9 @@ class Atom : protected Pointers {
// RHEO package
int *status;
double *conductivity;
double *pressure;
double *viscosity;
// SPH package
@ -194,7 +197,7 @@ class Atom : protected Pointers {
int temperature_flag, heatflow_flag;
int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag;
int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag;
int status_flag;
int status_flag, conductivity_flag, pressure_flag, viscosity_flag;
int rho_flag, esph_flag, cv_flag, vest_flag;
int dpd_flag, edpd_flag, tdpd_flag;
int mesont_flag;

200
src/atom_vec_rheo_thermal.cpp Normal file
View File

@ -0,0 +1,200 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
----------------------------------------------------------------------- */
#include "atom_vec_rheo_thermal.h"
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::ATOMIC;
mass_type = PER_TYPE;
forceclearflag = 1;
atom->status_flag = 1;
atom->conductivity_flag = 1;
atom->temperature_flag = 1;
atom->heatflow_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
atom->viscosity_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// 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_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_data_vel = {"id", "v"};
setup_fields();
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::grow_pointers()
{
status = atom->status;
conductivity = atom->conductivity;
temperature = atom->temperature;
heatflow = atom->heatflow;
pressure = atom->pressure;
rho = atom->rho;
drho = atom->drho;
viscosity = atom->viscosity;
}
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::force_clear(int n, size_t nbytes)
{
memset(&drho[n], 0, nbytes);
memset(&heatflow[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
heatflow[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
conductivity[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
assign an index to named atom property and return index
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecRHEOThermal::property_atom(const std::string &name)
{
if (name == "status") return 0;
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;
return -1;
}
/* ----------------------------------------------------------------------
pack per-atom data into buf for ComputePropertyAtom
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = status[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = temperature[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = heatflow[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 5) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = conductivity[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 6) {
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 == 7) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

46
src/atom_vec_rheo_thermal.h Normal file
View File

@ -0,0 +1,46 @@
/* -*- 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 ATOM_CLASS
// clang-format off
AtomStyle(rheo/thermal,AtomVecRHEOThermal);
// clang-format on
#else
#ifndef LMP_ATOM_VEC_RHEO_THERMAL_H
#define LMP_ATOM_VEC_RHEO_THERMAL_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecRHEOThermal : virtual public AtomVec {
public:
AtomVecRHEOThermal(class LAMMPS *);
void grow_pointers() override;
void force_clear(int, size_t) override;
void data_atom_post(int) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:
int *status;
double *conductivity, *temperature, *heatflow;
double *pressure, *rho, *drho, *viscosity;
};
} // namespace LAMMPS_NS
#endif
#endif

411
src/compute_rheo_property_atom.cpp Normal file
View File

@ -0,0 +1,411 @@
// clang-format off
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
----------------------------------------------------------------------- */
#include "compute_rheo_property_atom.h"
#include "atom.h"
#include "atom_vec.h"
#include "compute_rheo_interface.h"
#include "compute_rheo_kernel.h"
#include "compute_rheo_surface.h"
#include "compute_rheo_vshift.h"
#include "error.h"
#include "fix_rheo.h"
#include "fix_rheo_thermal.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */
ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr),
compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr),
index(nullptr), pack_choice(nullptr)
{
if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error);
peratom_flag = 1;
nvalues = narg - 3;
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
thermal_flag, interface_flag, surface_flag, shift_flag = 0;
// parse input values
// customize a new keyword by adding to if statement
pack_choice = new FnPtrPack[nvalues];
index = new int[nvalues];
int i;
for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3;
if (strcmp(arg[iarg],"phase") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase;
} else if (strcmp(arg[iarg],"chi") == 0) {
interface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi;
} else if (strcmp(arg[iarg],"surface") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface;
} else if (strcmp(arg[iarg],"surface/r") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r;
} else if (strcmp(arg[iarg],"surface/divr") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
} else if (strcmp(arg[iarg],"surface/nx") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nx;
} else if (strcmp(arg[iarg],"surface/ny") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_ny;
} else if (strcmp(arg[iarg],"surface/nz") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nz;
} else if (strcmp(arg[iarg],"coordination") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
} else if (strcmp(arg[iarg],"cv") == 0) {
thermal_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv;
} else if (strcmp(arg[iarg],"shift/vx") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vx;
} else if (strcmp(arg[iarg],"shift/vy") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vy;
} else if (strcmp(arg[iarg],"shift/vz") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vz;
} else {
index[i] = atom->avec->property_atom(arg[iarg]);
if (index[i] < 0)
error->all(FLERR,
"Invalid keyword {} for atom style {} in compute rheo/property/atom command ",
atom->get_style(), arg[iarg]);
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"heatflow") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"conductivity") == 0) thermal_flag = 1;
}
}
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
{
delete[] pack_choice;
delete[] index;
memory->destroy(vector_atom);
memory->destroy(array_atom);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::init()
{
auto fixes = modify->get_fix_by_style("rheo");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
if (interface_flag && !(fix_rheo->interface_flag))
error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo");
if (surface_flag && !(fix_rheo->surface_flag))
error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo");
if (shift_flag && !(fix_rheo->shift_flag))
error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo");
if (thermal_flag && !(fix_rheo->thermal_flag))
error->all(FLERR, "Cannot request thermal property without fix rheo/thermal");
compute_interface = fix_rheo->compute_interface;
compute_kernel = fix_rheo->compute_kernel;
compute_surface = fix_rheo->compute_surface;
compute_vshift = fix_rheo->compute_vshift;
if (thermal_flag) {
fixes = modify->get_fix_by_style("rheo/thermal");
fix_thermal = dynamic_cast<FixRHEOThermal *>(fixes[0]);
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow vector or array if necessary
if (atom->nmax > nmax) {
nmax = atom->nmax;
if (nvalues == 1) {
memory->destroy(vector_atom);
memory->create(vector_atom,nmax,"rheo/property/atom:vector");
} else {
memory->destroy(array_atom);
memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array");
}
}
// fill vector or array with per-atom values
if (nvalues == 1) {
buf = vector_atom;
(this->*pack_choice[0])(0);
} else {
if (nmax) buf = &array_atom[0][0];
else buf = nullptr;
for (int n = 0; n < nvalues; n++)
(this->*pack_choice[n])(n);
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeRHEOPropertyAtom::memory_usage()
{
double bytes = (double)nmax * nvalues * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword compute rheo/property/atom can output
the atom property is packed into buf starting at n with stride nvalues
customize a new keyword by adding a method
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_phase(int n)
{
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int inverse_mask = ~PHASEMASK;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_chi(int n)
{
double *chi = compute_interface->chi;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = chi[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface(int n)
{
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int inverse_mask = ~SURFACEMASK;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_r(int n)
{
double *rsurface = compute_surface->rsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = rsurface[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_divr(int n)
{
double *divr = compute_surface->divr;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = divr[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nx(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][0];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_ny(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nz(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_coordination(int n)
{
int *coordination = compute_kernel->coordination;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = coordination[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_cv(int n)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vx(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][0];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vy(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vz(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_atom_style(int n)
{
atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit);
}

71
src/compute_rheo_property_atom.h Normal file
View File

@ -0,0 +1,71 @@
/* -*- 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 COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/property/atom,ComputeRHEOPropertyAtom);
// clang-format on
#else
#ifndef LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
#define LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeRHEOPropertyAtom : public Compute {
public:
ComputeRHEOPropertyAtom(class LAMMPS *, int, char **);
~ComputeRHEOPropertyAtom() override;
void init() override;
void compute_peratom() override;
double memory_usage() override;
private:
int nvalues, nmax;
int thermal_flag, interface_flag, surface_flag, shift_flag;
int *index;
double *buf;
typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_phase(int);
void pack_chi(int);
void pack_surface(int);
void pack_surface_r(int);
void pack_surface_divr(int);
void pack_surface_nx(int);
void pack_surface_ny(int);
void pack_surface_nz(int);
void pack_coordination(int);
void pack_cv(int);
void pack_shift_vx(int);
void pack_shift_vy(int);
void pack_shift_vz(int);
void pack_atom_style(int);
class FixRHEO *fix_rheo;
class FixRHEOThermal *fix_thermal;
class ComputeRHEOInterface *compute_interface;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSurface *compute_surface;
class ComputeRHEOVShift *compute_vshift;
};
} // namespace LAMMPS_NS
#endif
#endif