Patching vremap in SPH

This commit is contained in:
jtclemm
2024-07-25 17:10:23 -06:00
parent e02ad1a3b2
commit 6abbbdba6c
2 changed files with 41 additions and 3 deletions

View File

@ -61,6 +61,10 @@ void FixSPH::init()
void FixSPH::setup_pre_force(int /*vflag*/) void FixSPH::setup_pre_force(int /*vflag*/)
{ {
remap_v_flag = domain->deform_vremap;
if (remap_v_flag && (!comm->ghost_velocity))
error->all(FLERR, "Fix sph requires ghost atoms store velocity when deforming with remap v");
// set vest equal to v // set vest equal to v
double **v = atom->v; double **v = atom->v;
double **vest = atom->vest; double **vest = atom->vest;
@ -119,9 +123,16 @@ void FixSPH::initial_integrate(int /*vflag*/)
rho[i] += dtf * drho[i]; // ... and density rho[i] += dtf * drho[i]; // ... and density
// extrapolate velocity for use with velocity-dependent potentials, e.g. SPH // extrapolate velocity for use with velocity-dependent potentials, e.g. SPH
vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; // if velocities are remapped, perform this extrapolation after communication
vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; if (remap_v_flag) {
vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; vest[i][0] = dtfm * f[i][0];
vest[i][1] = dtfm * f[i][1];
vest[i][2] = dtfm * f[i][2];
} else {
vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0];
vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1];
vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2];
}
v[i][0] += dtfm * f[i][0]; v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1]; v[i][1] += dtfm * f[i][1];
@ -136,6 +147,31 @@ void FixSPH::initial_integrate(int /*vflag*/)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixSPH::pre_force(int /*vflag*/)
{
// if velocities are remapped, calculate estimates here
// note that vest currently stores dtfm * force
if (!remap_v_flag) return;
double **v = atom->v;
double **vest = atom->vest;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) {
if (mask[i] & groupbit) {
vest[i][0] += v[i][0];
vest[i][1] += v[i][1];
vest[i][2] += v[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixSPH::final_integrate() void FixSPH::final_integrate()
{ {
// update v, rho, and e of atoms in group // update v, rho, and e of atoms in group

View File

@ -30,6 +30,7 @@ class FixSPH : public Fix {
int setmask() override; int setmask() override;
void init() override; void init() override;
void setup_pre_force(int) override; void setup_pre_force(int) override;
void pre_force(int) override;
void initial_integrate(int) override; void initial_integrate(int) override;
void final_integrate() override; void final_integrate() override;
void reset_dt() override; void reset_dt() override;
@ -41,6 +42,7 @@ class FixSPH : public Fix {
double dtv, dtf; double dtv, dtf;
double *step_respa; double *step_respa;
int mass_require; int mass_require;
int remap_v_flag;
class Pair *pair; class Pair *pair;
}; };