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

Adding setup check to MDR submodel

Browse Source
This commit is contained in:
jtclemm
2025-01-06 13:58:49 -07:00
parent 53084018ae
commit 09735b29e7
1 changed files with 15 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -545,7 +545,8 @@ double GranSubModNormalMDR::calculate_forces()
double R1 = 0.0;
double delta = gm->delta; // apparent overlap
double * history = & gm->history[history_index]; // load in all history variables
double *history = & gm->history[history_index]; // load in all history variables
int update = gm->history_update;
// Rigid flat placement scheme
double * deltamax_offset = & history[DELTA_MAX];
@ -816,9 +817,11 @@ double GranSubModNormalMDR::calculate_forces()
// area related calculations
double Ac;
(*Yflag_offset == 0.0) ? Ac = MY_PI * delta * R : Ac = MY_PI * ((2.0 * delta * R - pow(delta, 2)) + cA / MY_PI);
if (Ac < 0.0 ) Ac = 0.0;
Atot_sum[i] += wij * (Ac - 2.0 * MY_PI * R * (deltamax_MDR + delta_BULK));
Acon1[i] += wij * Ac;
if (Ac < 0.0) Ac = 0.0;
if (update) {
Atot_sum[i] += wij * (Ac - 2.0 * MY_PI * R * (deltamax_MDR + delta_BULK));
Acon1[i] += wij * Ac;
}
// bulk force calculation
double F_BULK;
@ -831,7 +834,7 @@ double GranSubModNormalMDR::calculate_forces()
*Ac_offset = wij * Ac;
// radius update scheme quantity calculation
Vcaps[i] += (MY_PI * THIRD) * pow(delta, 2) * (3.0 * R - delta);
if (update) Vcaps[i] += (MY_PI * THIRD) * pow(delta, 2) * (3.0 * R - delta);
const double Fntmp = wij * (F_MDR + F_BULK);
const double fx = Fntmp * gm->nx[0];
@ -841,19 +844,21 @@ double GranSubModNormalMDR::calculate_forces()
const double by = -(Ro - deltao) * gm->nx[1];
const double bz = -(Ro - deltao) * gm->nx[2];
const double eps_bar_contact = (1.0 / (3 * kappa * Velas[i])) * (fx * bx + fy * by + fz * bz);
eps_bar[i] += eps_bar_contact;
if (update) eps_bar[i] += eps_bar_contact;
double desp_bar_contact = eps_bar_contact - *eps_bar_offset;
if(delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){
if (update && elta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){
const double Vo = (4.0 * THIRD) * MY_PI * pow(Ro, 3);
dRnumerator[i] += -Vo * (eps_bar_contact - *eps_bar_offset) - wij * MY_PI * ddeltao * (2.0 * deltao * Ro - pow(deltao, 2) + pow(R, 2) - pow(Ro, 2));
dRdenominator[i] += wij * 2.0 * MY_PI * R * (deltao + R - Ro);
}
*eps_bar_offset = eps_bar_contact;
sigmaxx[i] += (1.0 / Velas[i]) * (fx * bx);
sigmayy[i] += (1.0 / Velas[i]) * (fy * by);
sigmazz[i] += (1.0 / Velas[i]) * (fz * bz);
if (update) {
sigmaxx[i] += (1.0 / Velas[i]) * (fx * bx);
sigmayy[i] += (1.0 / Velas[i]) * (fy * by);
sigmazz[i] += (1.0 / Velas[i]) * (fz * bz);
}
}
gm->i = i_true;