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

Combining fixes

Browse Source
This commit is contained in:
jtclemm
2024-12-18 12:58:37 -07:00
parent 61734a414c
commit ab6e3f6015
6 changed files with 261 additions and 360 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/.gitignore vendored
View File

@ -845,6 +845,8 @@
/fix_ffl.h
/fix_filter_corotate.cpp
/fix_filter_corotate.h
/fix_granular_mdr.cpp
/fix_granular_mdr.h
/fix_viscosity.cpp
/fix_viscosity.h
/fix_ehex.cpp

278
src/fix_mdr_mean_surf_disp.cpp → src/GRANULAR/fix_granular_mdr.cpp
View File

@ -17,78 +17,189 @@
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "fix_mdr_mean_surf_disp.h"
#include "fix_granular_mdr.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix_wall_gran_region.h"
#include "fix_neigh_history.h"
#include "force.h"
#include "granular_model.h"
#include "gran_sub_mod_normal.h"
#include "input.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "variable.h"
#include "fix_neigh_history.h"
#include "neigh_list.h"
#include "pair.h"
#include "pair_granular.h"
#include "granular_model.h"
#include "neigh_list.h"
#include "region.h"
#include "update.h"
#include "fix_wall_gran_region.h"
#include "comm.h"
#include "variable.h"
#include <algorithm>
#include <iomanip>
#include <sstream>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace Granular_NS;
using namespace FixConst;
using MathConst::MY_PI;
enum {MEAN_SURF_DISP, RADIUS_UPDATE};
/* ---------------------------------------------------------------------- */
FixMDRmeanSurfDisp::FixMDRmeanSurfDisp(LAMMPS *lmp, int narg, char **arg) :
FixGranularMDR::FixGranularMDR(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
comm_forward = 1; // value needs to match number of values you communicate
comm_forward = 20; // value needs to match number of values you communicate
}
int FixMDRmeanSurfDisp::setmask()
int FixGranularMDR::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
mask |= PRE_FORCE | END_OF_STEP;
return mask;
}
void FixMDRmeanSurfDisp::setup_pre_force(int /*vflag*/)
void FixGranularMDR::setup_pre_force(int /*vflag*/)
{
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
psi_b = atom->dvector[index_psi_b];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
pre_force(0);
}
int FixMDRmeanSurfDisp::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
void FixGranularMDR::setup(int /*vflag*/)
{
int m = 0;
for (int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = ddelta_bar[j];
}
return m;
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
psi_b = atom->dvector[index_psi_b];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
end_of_step();
}
void FixMDRmeanSurfDisp::unpack_forward_comm(int n, int first, double *buf)
void FixGranularMDR::pre_force(int)
{
int m = 0;
int last = first + n;
for (int i = first; i < last; i++) {
ddelta_bar[i] = buf[m++];
}
radius_update();
mean_surf_disp();
}
void FixMDRmeanSurfDisp::pre_force(int)
void FixGranularMDR::radius_update()
{
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
class GranularModel* model;
class GranularModel** models_list = pair->models_list;
class GranSubModNormalMDR* norm_model = nullptr;
for (int i = 0; i < pair->nmodels; i++) {
model = models_list[i];
if (model->normal_model->name == "mdr") norm_model = dynamic_cast<GranSubModNormalMDR *>(model->normal_model);
}
if (norm_model == nullptr) error->all(FLERR, "Did not find mdr model");
double *radius = atom->radius;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (history_setup_flag[i] < 1e-16) {
Ro[i] = radius[i];
Vgeo[i] = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
Velas[i] = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
Atot[i] = 4.0*MY_PI*pow(Ro[i],2.0);
psi[i] = 1.0;
psi_b[i] = norm_model->psi_b;
history_setup_flag[i] = 1.0;
}
sigmaxx[i] = 0.0;
sigmayy[i] = 0.0;
sigmazz[i] = 0.0;
contacts[i] = 0.0;
adhesive_length[i] = 0.0;
}
comm_stage = RADIUS_UPDATE;
comm->forward_comm(this, 20);
}
void FixGranularMDR::mean_surf_disp()
{
FixNeighHistory * fix_history = dynamic_cast<FixNeighHistory *>(modify->get_fix_by_id("NEIGH_HISTORY_GRANULAR"));
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
@ -464,5 +575,116 @@ void FixMDRmeanSurfDisp::pre_force(int)
}
}
}
comm->forward_comm(this);
comm_stage = MEAN_SURF_DISP;
comm->forward_comm(this, 1);
}
int FixGranularMDR::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
{
int m = 0;
if (comm_stage == RADIUS_UPDATE) {
for (int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = Ro[j]; // 1
buf[m++] = Vgeo[j]; // 2
buf[m++] = Velas[j]; // 3
buf[m++] = Vcaps[j]; // 4
buf[m++] = eps_bar[j]; // 5
buf[m++] = dRnumerator[j]; // 6
buf[m++] = dRdenominator[j]; // 7
buf[m++] = Acon0[j]; // 8
buf[m++] = Acon1[j]; // 9
buf[m++] = Atot[j]; // 10
buf[m++] = Atot_sum[j]; // 11
buf[m++] = ddelta_bar[j]; // 12
buf[m++] = psi[j]; // 13
buf[m++] = psi_b[j]; // 14
buf[m++] = sigmaxx[j]; // 15
buf[m++] = sigmayy[j]; // 16
buf[m++] = sigmazz[j]; // 17
buf[m++] = history_setup_flag[j]; // 18
buf[m++] = contacts[j]; // 19
buf[m++] = adhesive_length[j]; // 20
}
} else {
for (int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = ddelta_bar[j];
}
}
return m;
}
void FixGranularMDR::unpack_forward_comm(int n, int first, double *buf)
{
int m = 0;
int last = first + n;
if (comm_stage == RADIUS_UPDATE) {
for (int i = first; i < last; i++) {
Ro[i] = buf[m++]; // 1
Vgeo[i] = buf[m++]; // 2
Velas[i] = buf[m++]; // 3
Vcaps[i] = buf[m++]; // 4
eps_bar[i] = buf[m++]; // 5
dRnumerator[i] = buf[m++]; // 6
dRdenominator[i] = buf[m++]; // 7
Acon0[i] = buf[m++]; // 8
Acon1[i] = buf[m++]; // 9
Atot[i] = buf[m++]; // 10
Atot_sum[i] = buf[m++]; // 11
ddelta_bar[i] = buf[m++]; // 12
psi[i] = buf[m++]; // 13
psi_b[i] = buf[m++]; // 14
sigmaxx[i] = buf[m++]; // 15
sigmayy[i] = buf[m++]; // 16
sigmazz[i] = buf[m++]; // 17
history_setup_flag[i] = buf[m++]; // 18
contacts[i] = buf[m++]; // 19
adhesive_length[i] = buf[m++]; // 20
}
} else {
for (int i = first; i < last; i++) {
ddelta_bar[i] = buf[m++];
}
}
}
void FixGranularMDR::end_of_step()
{
// update the apparent radius of every particle
double *radius = atom->radius;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
const double R = radius[i];
Atot[i] = 4.0*MY_PI*pow(R,2.0) + Atot_sum[i];
const double Vo = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
const double Vgeoi = 4.0/3.0*MY_PI*pow(R,3.0) - Vcaps[i];
Vgeo[i] = std::min(Vgeoi,Vo);
const double Afree = Atot[i] - Acon1[i];
psi[i] = Afree/Atot[i];
const double dR = std::max(dRnumerator[i]/(dRdenominator[i] - 4.0*MY_PI*pow(R,2.0)),0.0);
if (psi_b[i] < psi[i]) {
if ((radius[i] + dR) < (1.5*Ro[i])) radius[i] += dR;
}
Velas[i] = Vo*(1.0 + eps_bar[i]);
Vcaps[i] = 0.0;
eps_bar[i] = 0.0;
dRnumerator[i] = 0.0;
dRdenominator[i] = 0.0;
Acon0[i] = Acon1[i];
Acon1[i] = 0.0;
Atot_sum[i] = 0.0;
ddelta_bar[i] = 0.0;
}
comm_stage = RADIUS_UPDATE;
comm->forward_comm(this, 20);
}

13
src/fix_mdr_radius_update.h → src/GRANULAR/fix_granular_mdr.h
View File

@ -13,18 +13,18 @@
#ifdef FIX_CLASS
// clang-format off
FixStyle(mdr/radius/update,FixMDRradiusUpdate);
FixStyle(GRANULAR/MDR,FixGranularMDR);
// clang-format on
#else
#ifndef LMP_FIX_MDR_RADIUS_UPDATE_H
#define LMP_FIX_MDR_RADIUS_UPDATE_H
#ifndef LMP_FIX_GRANULAR_MDR_H
#define LMP_FIX_GRANULAR_MDR_H
#include "fix.h"
namespace LAMMPS_NS {
class FixMDRradiusUpdate : public Fix {
class FixGranularMDR : public Fix {
public:
double * Ro;
double * Vgeo;
@ -47,7 +47,7 @@ class FixMDRradiusUpdate : public Fix {
double * contacts;
double * adhesive_length;
FixMDRradiusUpdate(class LAMMPS *, int, char **);
FixGranularMDR(class LAMMPS *, int, char **);
int setmask() override;
void setup(int) override;
void setup_pre_force(int) override;
@ -57,7 +57,10 @@ class FixMDRradiusUpdate : public Fix {
void unpack_forward_comm(int, int, double *) override;
private:
int comm_stage;
void radius_update();
void mean_surf_disp();
};
} // namespace LAMMPS_NS

3
src/GRANULAR/pair_granular.cpp
View File

@ -504,8 +504,7 @@ void PairGranular::init_style()
index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total length of adhesive contact on a particle
// Initiate MDR radius update fix
modify->add_fix("fix_mdr_radius_update all mdr/radius/update");
modify->add_fix("fix_mdr_mean_surf_disp all mdr/mean/surf/disp");
modify->add_fix("fix_mdr all GRANULAR/MDR");
fix_flag = 1;
}

46
src/fix_mdr_mean_surf_disp.h
View File

@ -1,46 +0,0 @@
/* -*- 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 FIX_CLASS
// clang-format off
FixStyle(mdr/mean/surf/disp,FixMDRmeanSurfDisp);
// clang-format on
#else
#ifndef LMP_FIX_MDR_MEAN_SURF_DISP_H
#define LMP_FIX_MDR_MEAN_SURF_DISP_H
#include "fix.h"
namespace LAMMPS_NS {
class FixMDRmeanSurfDisp : public Fix {
public:
double * Acon0;
double * ddelta_bar;
FixMDRmeanSurfDisp(class LAMMPS *, int, char **);
int setmask() override;
void setup_pre_force(int) override;
void pre_force(int) override; // FOR MDR
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
private:
};
} // namespace LAMMPS_NS
#endif
#endif

279
src/fix_mdr_radius_update.cpp
View File

@ -1,279 +0,0 @@
/* ----------------------------------------------------------------------
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:
William Zunker (MIT), Sachith Dunatunga (MIT),
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "fix_mdr_radius_update.h"
#include "atom.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "variable.h"
#include "granular_model.h"
#include "pair_granular.h"
#include "pair.h"
#include "gran_sub_mod_normal.h"
#include "update.h"
#include "comm.h"
#include <iomanip>
#include <sstream>
using namespace LAMMPS_NS;
using namespace Granular_NS;
using namespace FixConst;
using MathConst::MY_PI;
/* ---------------------------------------------------------------------- */
FixMDRradiusUpdate::FixMDRradiusUpdate(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
comm_forward = 20; // value needs to match number of values you communicate
}
int FixMDRradiusUpdate::setmask()
{
int mask = 0;
mask |= PRE_FORCE | END_OF_STEP;
return mask;
}
void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/)
{
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
psi_b = atom->dvector[index_psi_b];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
pre_force(0);
}
void FixMDRradiusUpdate::setup(int /*vflag*/)
{
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
psi_b = atom->dvector[index_psi_b];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
end_of_step();
}
void FixMDRradiusUpdate::pre_force(int)
{
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
class GranularModel* model;
class GranularModel** models_list = pair->models_list;
class GranSubModNormalMDR* norm_model = nullptr;
for (int i = 0; i < pair->nmodels; i++) {
model = models_list[i];
if (model->normal_model->name == "mdr") norm_model = dynamic_cast<GranSubModNormalMDR *>(model->normal_model);
}
if (norm_model == nullptr) error->all(FLERR, "Did not find mdr model");
double *radius = atom->radius;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (history_setup_flag[i] < 1e-16) {
Ro[i] = radius[i];
Vgeo[i] = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
Velas[i] = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
Atot[i] = 4.0*MY_PI*pow(Ro[i],2.0);
psi[i] = 1.0;
psi_b[i] = norm_model->psi_b;
history_setup_flag[i] = 1.0;
}
sigmaxx[i] = 0.0;
sigmayy[i] = 0.0;
sigmazz[i] = 0.0;
contacts[i] = 0.0;
adhesive_length[i] = 0.0;
}
comm->forward_comm(this);
}
int FixMDRradiusUpdate::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
{
int m = 0;
for (int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = Ro[j]; // 1
buf[m++] = Vgeo[j]; // 2
buf[m++] = Velas[j]; // 3
buf[m++] = Vcaps[j]; // 4
buf[m++] = eps_bar[j]; // 5
buf[m++] = dRnumerator[j]; // 6
buf[m++] = dRdenominator[j]; // 7
buf[m++] = Acon0[j]; // 8
buf[m++] = Acon1[j]; // 9
buf[m++] = Atot[j]; // 10
buf[m++] = Atot_sum[j]; // 11
buf[m++] = ddelta_bar[j]; // 12
buf[m++] = psi[j]; // 13
buf[m++] = psi_b[j]; // 14
buf[m++] = sigmaxx[j]; // 15
buf[m++] = sigmayy[j]; // 16
buf[m++] = sigmazz[j]; // 17
buf[m++] = history_setup_flag[j]; // 18
buf[m++] = contacts[j]; // 19
buf[m++] = adhesive_length[j]; // 20
}
return m;
}
void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf)
{
int m = 0;
int last = first + n;
for (int i = first; i < last; i++) {
Ro[i] = buf[m++]; // 1
Vgeo[i] = buf[m++]; // 2
Velas[i] = buf[m++]; // 3
Vcaps[i] = buf[m++]; // 4
eps_bar[i] = buf[m++]; // 5
dRnumerator[i] = buf[m++]; // 6
dRdenominator[i] = buf[m++]; // 7
Acon0[i] = buf[m++]; // 8
Acon1[i] = buf[m++]; // 9
Atot[i] = buf[m++]; // 10
Atot_sum[i] = buf[m++]; // 11
ddelta_bar[i] = buf[m++]; // 12
psi[i] = buf[m++]; // 13
psi_b[i] = buf[m++]; // 14
sigmaxx[i] = buf[m++]; // 15
sigmayy[i] = buf[m++]; // 16
sigmazz[i] = buf[m++]; // 17
history_setup_flag[i] = buf[m++]; // 18
contacts[i] = buf[m++]; // 19
adhesive_length[i] = buf[m++]; // 20
}
}
void FixMDRradiusUpdate::end_of_step()
{
// update the apparent radius of every particle
double *radius = atom->radius;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
const double R = radius[i];
Atot[i] = 4.0*MY_PI*pow(R,2.0) + Atot_sum[i];
const double Vo = 4.0/3.0*MY_PI*pow(Ro[i],3.0);
const double Vgeoi = 4.0/3.0*MY_PI*pow(R,3.0) - Vcaps[i];
Vgeo[i] = std::min(Vgeoi,Vo);
const double Afree = Atot[i] - Acon1[i];
psi[i] = Afree/Atot[i];
const double dR = std::max(dRnumerator[i]/(dRdenominator[i] - 4.0*MY_PI*pow(R,2.0)),0.0);
if (psi_b[i] < psi[i]) {
if ((radius[i] + dR) < (1.5*Ro[i])) radius[i] += dR;
}
Velas[i] = Vo*(1.0 + eps_bar[i]);
Vcaps[i] = 0.0;
eps_bar[i] = 0.0;
dRnumerator[i] = 0.0;
dRdenominator[i] = 0.0;
Acon0[i] = Acon1[i];
Acon1[i] = 0.0;
Atot_sum[i] = 0.0;
ddelta_bar[i] = 0.0;
}
comm->forward_comm(this);
}