tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Converting history indices to an enum

Browse Source
This commit is contained in:
jtclemm
2024-12-24 10:33:13 -07:00
parent aecbc21123
commit a9ce245527
3 changed files with 89 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
src/GRANULAR/fix_granular_mdr.cpp
View File

@ -40,6 +40,7 @@
using namespace LAMMPS_NS;
using namespace Granular_NS;
using namespace Granular_MDR_NS;
using namespace FixConst;
using MathConst::MY_PI;
@ -547,10 +548,10 @@ void FixGranularMDR::mean_surf_disp()
NeighList * list = pair->list;
const int size_history = pair->get_size_history();
int i,j,k,ii,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch;
double *history,*allhistory,**firsthistory;
int i, j, k, ii, jj, inum, jnum, itype, jtype;
int *ilist, *jlist, *numneigh, **firstneigh;
int *touch, **firsttouch;
double *history, *allhistory, **firsthistory;
bool touchflag = false;
class GranularModel* model;
@ -612,22 +613,14 @@ void FixGranularMDR::mean_surf_disp()
const double delta = model->radsum - sqrt(model->rsq);
const int delta_offset_0 = 0; // apparent overlap
const int delta_offset_1 = 1;
const int Ac_offset_0 = 18; // contact area
const int Ac_offset_1 = 19;
const int deltamax_offset_ = 23;
const int deltap_offset_0 = 24;
const int deltap_offset_1 = 25;
double deltamax = history[deltamax_offset_];
double deltap0 = history[deltap_offset_0];
double deltap1 = history[deltap_offset_1];
double deltamax = history[DELTA_MAX];
double deltap0 = history[DELTAP_0];
double deltap1 = history[DELTAP_1];
if (delta > deltamax) deltamax = delta;
double delta0old = history[delta_offset_0];
double delta1old = history[delta_offset_1];
double delta0old = history[DELTA_0];
double delta1old = history[DELTA_1];
int i0;
int i1;
@ -644,37 +637,37 @@ void FixGranularMDR::mean_surf_disp()
double delta_geo0;
double delta_geo1;
double deltaOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1));
double deltaOpt2 = deltamax*(deltamax - 2.0*R0)/(2.0*(deltamax - R0 - R1));
(R0 < R1) ? delta_geo0 = MAX(deltaOpt1,deltaOpt2) : delta_geo0 = MIN(deltaOpt1,deltaOpt2);
(R0 < R1) ? delta_geo1 = MIN(deltaOpt1,deltaOpt2) : delta_geo1 = MAX(deltaOpt1,deltaOpt2);
double deltaOpt1 = deltamax * (deltamax - 2.0 * R1) / (2.0 * (deltamax - R0 - R1));
double deltaOpt2 = deltamax * (deltamax - 2.0 * R0) / (2.0 * (deltamax - R0 - R1));
(R0 < R1) ? delta_geo0 = MAX(deltaOpt1, deltaOpt2) : delta_geo0 = MIN(deltaOpt1, deltaOpt2);
(R0 < R1) ? delta_geo1 = MIN(deltaOpt1, deltaOpt2) : delta_geo1 = MAX(deltaOpt1, deltaOpt2);
double overlap_limit = 0.75;
if (delta_geo0/R0 > overlap_limit) {
delta_geo0 = R0*overlap_limit;
if (delta_geo0 / R0 > overlap_limit) {
delta_geo0 = R0 * overlap_limit;
delta_geo1 = deltamax - delta_geo0;
} else if (delta_geo1/R1 > overlap_limit) {
delta_geo1 = R1*overlap_limit;
} else if (delta_geo1 / R1 > overlap_limit) {
delta_geo1 = R1 * overlap_limit;
delta_geo0 = deltamax - delta_geo1;
}
double deltap = deltap0 + deltap1;
double delta0 = delta_geo0 + (deltap0 - delta_geo0)/(deltap - deltamax)*(delta-deltamax);
double delta1 = delta_geo1 + (deltap1 - delta_geo1)/(deltap - deltamax)*(delta-deltamax);
double delta0 = delta_geo0 + (deltap0 - delta_geo0) / (deltap - deltamax) * (delta - deltamax);
double delta1 = delta_geo1 + (deltap1 - delta_geo1) / (deltap - deltamax) * (delta - deltamax);
double ddel0 = delta0 - delta0old;
double ddel1 = delta1 - delta1old;
if (Acon0[i0] != 0.0) {
const double Ac_offset0 = history[Ac_offset_0];
ddelta_bar[i0] += Ac_offset0/Acon0[i0]*ddel0;
const double Ac_offset0 = history[AC_0];
ddelta_bar[i0] += Ac_offset0 / Acon0[i0] * ddel0;
}
if (Acon0[i1] != 0.0) {
const double Ac_offset1 = history[Ac_offset_1];
ddelta_bar[i1] += Ac_offset1/Acon0[i1]*ddel1;
const double Ac_offset1 = history[AC_1];
ddelta_bar[i1] += Ac_offset1 / Acon0[i1] * ddel1;
}
}
}
@ -698,7 +691,7 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in)
model->history_update = history_update;
int regiondynamic = region->dynamic_check();
if (!regiondynamic) vwall[0] = vwall[1] = vwall[2] = 00;
if (!regiondynamic) vwall[0] = vwall[1] = vwall[2] = 0;
double **x = atom->x;
double *radius = atom->radius;
@ -762,7 +755,7 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in)
const double delta_offset0 = fix->history_many[i][fix->c2r[ic]][0];
const double ddelta = delta - delta_offset0;
const double Ac_offset0 = fix->history_many[i][fix->c2r[ic]][18];
ddelta_bar[i] += wij*Ac_offset0/Acon0[i]*ddelta;
ddelta_bar[i] += wij * Ac_offset0 / Acon0[i] * ddelta;
}
}
}

32
src/GRANULAR/fix_granular_mdr.h
View File

@ -23,6 +23,38 @@ FixStyle(GRANULAR/MDR,FixGranularMDR);
#include "fix.h"
namespace LAMMPS_NS {
namespace Granular_MDR_NS {
enum HistoryIndex {
DELTA_0 = 0, // apparent overlap
DELTA_1,
DELTAO_0, // displacement
DELTAO_1,
DELTA_MDR_0, // MDR apparent overlap
DELTA_MDR_1,
DELTA_BULK_0, // bulk displacement
DELTA_BULK_1,
DELTAMAX_MDR_0, // maximum MDR apparent overlap
DELTAMAX_MDR_1,
YFLAG_0, // yield flag
YFLAG_1,
DELTAY_0, // yield displacement
DELTAY_1,
CA_0, // contact area intercept
CA_1,
AADH_0, // adhesive contact radius
AADH_1,
AC_0, // contact area
AC_1,
EPS_BAR_0, // volume-averaged infinitesimal sor
EPS_BAR_1,
PENALTY, // contact penalty
DELTA_MAX,
DELTAP_0,
DELTAP_1
};
} // namespace Granular_MDR_NS
class FixGranularMDR : public Fix {
public:

88
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -16,6 +16,7 @@
#include "atom.h"
#include "error.h"
#include "citeme.h"
#include "fix_granular_mdr.h"
#include "granular_model.h"
#include "math_const.h"
#include "modify.h"
@ -494,6 +495,7 @@ void GranSubModNormalMDR::init()
double GranSubModNormalMDR::calculate_forces()
{
using namespace Granular_MDR_NS;
// To understand the structure of the overall code it is important to consider
// the following:
//
@ -545,41 +547,13 @@ double GranSubModNormalMDR::calculate_forces()
int i1 = 0;
double delta = gm->delta; // apparent overlap
// initialize indexing in history array of different constact history variables
const int delta_offset_0 = 0; // apparent overlap
const int delta_offset_1 = 1;
const int deltao_offset_0 = 2; // displacement
const int deltao_offset_1 = 3;
const int delta_MDR_offset_0 = 4; // MDR apparent overlap
const int delta_MDR_offset_1 = 5;
const int delta_BULK_offset_0 = 6; // bulk displacement
const int delta_BULK_offset_1 = 7;
const int deltamax_MDR_offset_0 = 8; // maximum MDR apparent overlap
const int deltamax_MDR_offset_1 = 9;
const int Yflag_offset_0 = 10; // yield flag
const int Yflag_offset_1 = 11;
const int deltaY_offset_0 = 12; // yield displacement
const int deltaY_offset_1 = 13;
const int cA_offset_0 = 14; // contact area intercept
const int cA_offset_1 = 15;
const int aAdh_offset_0 = 16; // adhesive contact radius
const int aAdh_offset_1 = 17;
const int Ac_offset_0 = 18; // contact area
const int Ac_offset_1 = 19;
const int eps_bar_offset_0 = 20; // volume-averaged infinitesimal strain tensor
const int eps_bar_offset_1 = 21;
const int penalty_offset_ = 22; // contact penalty
const int deltamax_offset_ = 23;
const int deltap_offset_0 = 24;
const int deltap_offset_1 = 25;
double * history = & gm->history[history_index]; // load in all history variables
// Rigid flat placement scheme
double * deltamax_offset = & history[deltamax_offset_];
double * deltamax_offset = & history[DELTA_MAX];
double deltamax = *deltamax_offset;
double * deltap_offset0 = & history[deltap_offset_0];
double * deltap_offset1 = & history[deltap_offset_1];
double * deltap_offset0 = & history[DELTAP_0];
double * deltap_offset1 = & history[DELTAP_1];
double deltap0 = *deltap_offset0;
double deltap1 = *deltap_offset1;
@ -627,18 +601,18 @@ double GranSubModNormalMDR::calculate_forces()
delta = delta_geo + (deltap0 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
}
delta_offset = & history[delta_offset_0];
deltao_offset = & history[deltao_offset_0];
delta_MDR_offset = & history[delta_MDR_offset_0];
delta_BULK_offset = & history[delta_BULK_offset_0];
deltamax_MDR_offset = & history[deltamax_MDR_offset_0];
Yflag_offset = & history[Yflag_offset_0];
deltaY_offset = & history[deltaY_offset_0];
cA_offset = & history[cA_offset_0];
aAdh_offset = & history[aAdh_offset_0];
Ac_offset = & history[Ac_offset_0];
eps_bar_offset = & history[eps_bar_offset_0];
deltap_offset = & history[deltap_offset_0];
delta_offset = & history[DELTA_0];
deltao_offset = & history[DELTAO_0];
delta_MDR_offset = & history[DELTA_MDR_0];
delta_BULK_offset = & history[DELTA_BULK_0];
deltamax_MDR_offset = & history[DELTAMAX_MDR_0];
Yflag_offset = & history[YFLAG_0];
deltaY_offset = & history[DELTAY_0];
cA_offset = & history[CA_0];
aAdh_offset = & history[AADH_0];
Ac_offset = & history[AC_0];
eps_bar_offset = & history[EPS_BAR_0];
deltap_offset = & history[DELTAP_0];
} else {
if (gm->contact_type != PAIR) break; // contact with particle-wall requires only one evaluation
if (itag_true < jtag_true) {
@ -668,18 +642,18 @@ double GranSubModNormalMDR::calculate_forces()
double deltap = deltap0 + deltap1;
delta = delta_geo + (deltap1 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
delta_offset = & history[delta_offset_1];
deltao_offset = & history[deltao_offset_1];
delta_MDR_offset = & history[delta_MDR_offset_1];
delta_BULK_offset = & history[delta_BULK_offset_1];
deltamax_MDR_offset = & history[deltamax_MDR_offset_1];
Yflag_offset = & history[Yflag_offset_1];
deltaY_offset = & history[deltaY_offset_1];
cA_offset = & history[cA_offset_1];
aAdh_offset = & history[aAdh_offset_1];
Ac_offset = & history[Ac_offset_1];
eps_bar_offset = & history[eps_bar_offset_1];
deltap_offset = & history[deltap_offset_1];
delta_offset = & history[DELTA_1];
deltao_offset = & history[DELTAO_1];
delta_MDR_offset = & history[DELTA_MDR_1];
delta_BULK_offset = & history[DELTA_BULK_1];
deltamax_MDR_offset = & history[DELTAMAX_MDR_1];
Yflag_offset = & history[YFLAG_1];
deltaY_offset = & history[DELTAY_1];
cA_offset = & history[CA_1];
aAdh_offset = & history[AADH_1];
Ac_offset = & history[AC_1];
eps_bar_offset = & history[EPS_BAR_1];
deltap_offset = & history[DELTAP_1];
}
// temporary i and j indices
@ -834,7 +808,7 @@ double GranSubModNormalMDR::calculate_forces()
adhesive_length[i] += aAdh;
// contact penalty scheme
penalty_offset = & history[penalty_offset_];
penalty_offset = & history[PENALTY];
double pij = *penalty_offset;
const double wij = std::max(1.0-pij,0.0);
@ -886,7 +860,7 @@ double GranSubModNormalMDR::calculate_forces()
gm->radi = radi_true;
gm->radj = radj_true;
double * penalty_offset = & history[penalty_offset_];
double * penalty_offset = & history[PENALTY];
const double pij = *penalty_offset;
const double wij = std::max(1.0-pij,0.0);
*penalty_offset = 0.0;