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

Stablized radius update, added WFM, added outputs

Browse Source
This commit is contained in:
William Zunker
2025-03-26 10:43:51 -04:00
parent 285baf27b5
commit 9d01ac2caf
7 changed files with 236 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
src/GRANULAR/fix_granular_mdr.cpp
View File

@ -38,6 +38,10 @@
#include "update.h"
#include "variable.h"
// wzunker
#include "csv_writer.h"
#include <sstream>
using namespace LAMMPS_NS;
using namespace Granular_NS;
using namespace Granular_MDR_NS;
@ -45,7 +49,7 @@ using namespace FixConst;
using MathConst::MY_PI;
static constexpr double EPSILON = 1e-16;
static constexpr double OVERLAP_LIMIT = 0.75;
static constexpr double OVERLAP_LIMIT = 0.95;
enum { COMM_1, COMM_2 };
@ -85,7 +89,7 @@ void FixGranularMDR::post_constructor()
modify->add_fix(
fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator "
"d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi "
"d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz ghost yes",
"d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_dRavg ghost yes",
id_fix));
index_Ro = atom->find_custom("Ro", tmp1, tmp2);
@ -105,6 +109,7 @@ void FixGranularMDR::post_constructor()
index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2);
index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2);
index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2);
index_dRavg = atom->find_custom("dRavg", tmp1, tmp2);
}
/* ---------------------------------------------------------------------- */
@ -208,6 +213,7 @@ void FixGranularMDR::pre_force(int)
double *sigmayy = atom->dvector[index_sigmayy];
double *sigmazz = atom->dvector[index_sigmazz];
double *history_setup_flag = atom->dvector[index_history_setup_flag];
double *dRavg = atom->dvector[index_dRavg];
int new_atom;
int nlocal = atom->nlocal;
@ -250,10 +256,35 @@ void FixGranularMDR::pre_force(int)
psi[i] = (Atot[i] - Acon1[i]) / Atot[i];
if (psi_b_coeff < psi[i]) {
const double dR = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)), 0.0);
if ((radius[i] + dR) < (1.5 * Ro[i])) radius[i] += dR;
double w_confinement;
( psi[i] > 0.1 ) ? w_confinement = 1.0/(1.0 + exp(-75.0*(psi[i]-0.2))) : w_confinement = 0.0;
const double dR = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0))*w_confinement, 0.0);
const double N_window = 10.0;
if (dR > 0.0) dRavg[i] += (dR - dRavg[i]) / N_window;
if (((radius[i] + dR) < (1.5 * Ro[i])) && (dR > 0.0)) radius[i] += dRavg[i];
//(dR + dRavg[i])/2.0;
//dRavg[i] = (dR + dRavg[i])/2.0;
//// wzunker
//const double dRdenominatorTrue = (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0));
//if (i == 0) {
// CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRavg[i] << ", " << R << ", " << dRdenominatorTrue << ", " << lmp->update->ntimestep;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
}
Acon0[i] = Acon1[i];
// wzunker
//const double dR_tmp = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)), 0.0);
//printf("i = %d, psi_b = %f, psi = %f, Atot = %f, Acon = %f, dRnumerator = %f, dRdenominator = %f, R = %f, dR = %f \n", i, psi_b_coeff, psi[i], Atot[i], Acon1[i], dRnumerator[i], dRdenominator[i],radius[i],dR_tmp);
}
comm_stage = COMM_1;
@ -510,6 +541,9 @@ void FixGranularMDR::calculate_contact_penalty()
"(e.g. increase the skin distance).");
pjk[0] += 1.0 / (1.0 + std::exp(-50.0 * (alpha / MY_PI - 0.5)));
// wzunker
//printf("i = %d, j = %d, k = %d, pij = %f, pik = %f, pjk = %f \n", i, j, k, pij[0], pik[0], pjk[0]);
}
}
}

3
src/GRANULAR/fix_granular_mdr.h
View File

@ -97,8 +97,7 @@ class FixGranularMDR : public Fix {
int index_sigmayy; // yy-component of the stress tensor, not necessary forforce calculation
int index_sigmazz; // zz-component of the stress tensor, not necessary forforce calculation
int index_history_setup_flag; // flag to check if history variables have beeninitialized
int index_contacts; // total contacts on particle
int index_adhesive_length; // total length of adhesive contact on a particle
int index_dRavg; // total contacts on particle
};
} // namespace LAMMPS_NS

6
src/GRANULAR/gran_sub_mod_damping.cpp
View File

@ -18,6 +18,8 @@
#include "math_special.h"
#include "math_const.h"
#include "style_gran_sub_mod.h" // IWYU pragma: keep
#include <cmath>
using namespace LAMMPS_NS;
@ -74,9 +76,13 @@ GranSubModDampingVelocity::GranSubModDampingVelocity(GranularModel *gm, LAMMPS *
double GranSubModDampingVelocity::calculate_forces()
{
if (gm->normal_model->name == "mdr") {
return 0.0;
} else {
damp_prefactor = damp;
return -damp_prefactor * gm->vnnr;
}
}
/* ----------------------------------------------------------------------
Mass velocity damping

171
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -26,6 +26,9 @@
#include <iomanip>
#include <sstream>
// wzunker
#include "csv_writer.h"
using namespace LAMMPS_NS;
using namespace Granular_NS;
using namespace MathConst;
@ -47,7 +50,7 @@ static constexpr int MDR_MAX_IT = 100; // Newton-Raphs
static constexpr double MDR_EPSILON1 = 1e-10; // Newton-Raphson for MDR
static constexpr double MDR_EPSILON2 = 1e-16; // Newton-Raphson for MDR
static constexpr double MDR_EPSILON3 = 1e-20; // For precision checks
static constexpr double MDR_OVERLAP_LIMIT = 0.75; // Maximum contact overlap for MDR
static constexpr double MDR_OVERLAP_LIMIT = 0.95; // Maximum contact overlap for MDR
static const char cite_mdr[] =
"MDR contact model command: (i) https://doi.org/10.1016/j.jmps.2023.105492 || (ii) https://doi.org/10.1016/j.jmps.2023.105493 || (iii) https://doi.org/10.31224/4289\n\n"
@ -472,14 +475,16 @@ void GranSubModNormalMDR::coeffs_to_local()
Y = coeffs[2]; // yield stress
gamma = coeffs[3]; // effective surface energy
psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total}
CoR = coeffs[5]; // coefficent of restitution
//CoR = coeffs[5]; // coefficent of restitution
damp = coeffs[5]; // coefficent of restitution
if (E <= 0.0) error->all(FLERR, "Illegal MDR normal model, Young's modulus must be greater than 0");
if (nu < 0.0 || nu > 0.5) error->all(FLERR, "Illegal MDR normal model, Poisson's ratio must be between 0 and 0.5");
if (Y < 0.0) error->all(FLERR, "Illegal MDR normal model, yield stress must be greater than or equal to 0");
if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0");
if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0");
if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0");
//if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0");
if (damp < 0.0) error->all(FLERR, "Illegal MDR normal model, damping coefficent must be greater than 1");
G = E / (2.0 * (1.0 + nu)); // shear modulus
kappa = E / (3.0 * (1.0 - 2.0 * nu)); // bulk modulus
@ -490,6 +495,7 @@ void GranSubModNormalMDR::coeffs_to_local()
Eeffinv = 1.0 / Eeff;
Eeffsq = Eeff * Eeff;
Eeffsqinv = Eeffinv * Eeffinv;
Eeff2particle = 0.5 * Eeff;
gammasq = gamma * gamma;
gamma3 = gammasq * gamma;
@ -528,6 +534,7 @@ void GranSubModNormalMDR::init()
index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); // xx-component of the stress tensor, not necessary for force calculation
index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); // yy-component of the stress tensor, not necessary for force calculation
index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); // zz-component of the stress tensor, not necessary for force calculation
index_dRavg = atom->find_custom("dRavg", tmp1, tmp2); // radius update increment
}
/* ---------------------------------------------------------------------- */
@ -566,6 +573,7 @@ double GranSubModNormalMDR::calculate_forces()
double *sigmaxx = atom->dvector[index_sigmaxx];
double *sigmayy = atom->dvector[index_sigmayy];
double *sigmazz = atom->dvector[index_sigmazz];
double *dRavg = atom->dvector[index_dRavg];
const int itag_true = atom->tag[gm->i]; // true i particle tag
const int jtag_true = atom->tag[gm->j]; // true j particle tag
@ -579,6 +587,7 @@ double GranSubModNormalMDR::calculate_forces()
double F1 = 0.0; // force on contact side 1
double delta = gm->delta; // apparent overlap
double Ac_avg = 0.0; // average contact area across both sides
double a_damp = 0.0; // damping contact radius
double *history = & gm->history[history_index]; // load in all history variables
int history_update = gm->history_update;
@ -596,6 +605,11 @@ double GranSubModNormalMDR::calculate_forces()
if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta;
double deltamax = *deltamax_offset;
//wzunker
double F_MDR0;
double F_BULK0;
double F_MDR1;
double F_BULK1;
for (int contactSide = 0; contactSide < 2; contactSide++) {
@ -804,8 +818,10 @@ double GranSubModNormalMDR::calculate_forces()
F_MDR = calculate_nonadhesive_mdr_force(deltae1D, Ainv, Eeff, A, B);
}
if (std::isnan(F_MDR))
if (std::isnan(F_MDR)) {
printf("itag = %d, jtag = %d \n", itag_true, jtag_true);
error->one(FLERR, "F_MDR is NaN, case 1: no tensile springs");
}
if (history_update) *aAdh_offset = a_fac * a_na;
} else {
@ -861,6 +877,10 @@ double GranSubModNormalMDR::calculate_forces()
}
aAdh = aAdh_tmp;
if (aAdh < acrit) {
aAdh = 0.0;
F_MDR = 0.0;
} else {
g_aAdh = A * 0.5 - A * Binv * sqrt(Bsq * 0.25 - pow(aAdh, 2));
g_aAdh = round_up_negative_epsilon(g_aAdh);
@ -871,6 +891,7 @@ double GranSubModNormalMDR::calculate_forces()
if (std::isnan(F_MDR))
error->one(FLERR, "F_MDR is NaN, case 3: tensile springs exceed critical length");
}
}
if (history_update) *aAdh_offset = aAdh;
}
}
@ -891,13 +912,22 @@ double GranSubModNormalMDR::calculate_forces()
}
Ac_avg += wij * Ac;
// wzunker damping contact area related things
(gamma > 0.0) ? a_damp += aAdh : a_damp += a_na;
// bulk force calculation
double F_BULK;
(delta_BULK <= 0.0) ? F_BULK = 0.0 : F_BULK = (1.0 / Vgeo[i]) * Acon0[i] * delta_BULK * kappa * Ac;
// force-magnifier
// total force calculation
(contactSide == 0) ? F0 = F_MDR + F_BULK : F1 = F_MDR + F_BULK;
// wzunker
(contactSide == 0) ? F_MDR0 = F_MDR : F_MDR1 = F_MDR;
(contactSide == 0) ? F_BULK0 = F_BULK : F_BULK1 = F_BULK;
if (history_update) {
// mean surface displacement calculation
*Ac_offset = wij * Ac;
@ -921,6 +951,32 @@ double GranSubModNormalMDR::calculate_forces()
dRnumerator[i] -= Vo * (eps_bar_contact - *eps_bar_offset);
dRnumerator[i] -= wij * MY_PI * ddeltao * (2 * deltao * Ro - pow(deltao, 2) + pow(R, 2) - pow(Ro, 2));
dRdenominator[i] += wij * 2.0 * MY_PI * R * (deltao + R - Ro);
// wzunker
//if (gm->contact_type == PAIR && i == 0) {
// CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters_0.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << wij << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << eps_bar_contact << ", " << *eps_bar_offset << ", " << ddeltao << ", " << deltao << ", " << R << ", " << Fntmp << ", " << lmp->update->ntimestep;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//
//if (gm->contact_type == PAIR && i == 1) {
// CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters_1.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << wij << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << eps_bar_contact << ", " << *eps_bar_offset << ", " << ddeltao << ", " << deltao << ", " << R << ", " << Fntmp << ", " << lmp->update->ntimestep;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
// wzunker
//if (gm->contact_type == PAIR) {
// printf("i = %d, j = %d, contact_type = %d, dRnumerator = %f, dRdenominator = %f, eps_bar_contact = %f, eps_bar_offset = %f, wij = %f, ddelto = %f, deltao = %f, R = %f, Fntmp = %f \n", i, gm->j, gm->contact_type, dRnumerator[i], dRdenominator[i], eps_bar_contact, *eps_bar_offset, wij, ddeltao, deltao, R, Fntmp);
//}
}
if (history_update) {
@ -944,31 +1000,106 @@ double GranSubModNormalMDR::calculate_forces()
const double wij = MAX(1.0 - pij, 0.0);
// assign final force
double damp_prefactor;
if (gm->contact_type != PAIR) {
F = wij * F0;
a_damp = a_damp/2.0;
damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff2particle * a_damp);
double *deltao_offset = &history[DELTAO_0];
const double wfm = std::exp(10.7*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0; // wall force magnifier
//const double wfm = 1.0;
F = wij * F0 * wfm;
} else {
damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff * a_damp);
F = wij * (F0 + F1) * 0.5;
}
// calculate damping force
if (F > 0.0) {
double Eeff2;
double Reff2;
if (gm->contact_type == PAIR) {
Eeff2 = E / (2.0 * (1.0 - pow(nu, 2)));
Reff2 = 1.0 / ((1.0 / gm->radi + 1.0 / gm->radj));
} else {
Eeff2 = E / (1.0 - pow(nu, 2));
Reff2 = gm->radi;
}
const double kn = Eeff2 * Reff2;
const double beta = -log(CoR) / sqrt(pow(log(CoR), 2) + PISQ);
const double damp_prefactor = beta * sqrt(gm->meff * kn);
const double F_DAMP = -damp_prefactor * gm->vnnr;
//if (F > 0.0) {
// double Eeff2;
// double Reff2;
// if (gm->contact_type == PAIR) {
// Eeff2 = E / (2.0 * (1.0 - pow(nu, 2)));
// Reff2 = 1.0 / ((1.0 / gm->radi + 1.0 / gm->radj));
// } else {
// Eeff2 = E / (1.0 - pow(nu, 2));
// Reff2 = gm->radi;
// }
// const double kn = Eeff2 * Reff2;
// const double beta = -log(CoR) / sqrt(pow(log(CoR), 2) + PISQ);
// const double damp_prefactor = beta * sqrt(gm->meff * kn);
// const double F_DAMP = -damp_prefactor * gm->vnnr;
//
// F += wij * F_DAMP;
//}
F += wij * F_DAMP;
// calculate damping force
if (F > 0.0) F += -wij * damp * gm->vnnr;
//F += -wij * damp_prefactor * gm->vnnr;
//printf("i = %d, j = %d, wij = %f, damp = %f, vnnr = %e, F_DAMP = %e \n", gm->i, gm->j, wij, damp, gm->vnnr, -wij * damp * gm->vnnr);
// wzunker
double *delta_offset_0 = &history[DELTA_0];
double *delta_offset_1 = &history[DELTA_0 + 1];
const double delta0 = *delta_offset_0;
const double delta1 = *delta_offset_1;
double *delta_BULK_offset_0 = &history[DELTA_BULK_0];
double *delta_BULK_offset_1 = &history[DELTA_BULK_0 + 1];
const double deltaBULK0 = *delta_BULK_offset_0;
const double deltaBULK1 = *delta_BULK_offset_1;
double *aAdh_offset_0 = &history[AADH_0];
double *aAdh_offset_1 = &history[AADH_0 + 1];
const double aAdh0 = *aAdh_offset_0;
const double aAdh1 = *aAdh_offset_1;
double *Ac_offset_0 = &history[AC_0];
double *Ac_offset_1 = &history[AC_0 + 1];
const double Ac0 = *Ac_offset_0;
const double Ac1 = *Ac_offset_1;
double F_DAMP = -wij * damp * gm->vnnr;
const double i_0 = 0;
const double i_1 = 1;
double psi_0;
double psi_1;
double Acon_0;
double Acon_1;
double Vgeo_0;
double Vgeo_1;
if (itag_true == i_0) {
psi_0 = psi[gm->i];
psi_1 = psi[gm->j];
Acon_0 = Acon0[gm->i];
Acon_1 = Acon0[gm->j];
Vgeo_0 = Vgeo[gm->i];
Vgeo_1 = Vgeo[gm->j];
} else {
psi_0 = psi[gm->j];
psi_1 = psi[gm->i];
Acon_0 = Acon0[gm->j];
Acon_1 = Acon0[gm->i];
Vgeo_0 = Vgeo[gm->j];
Vgeo_1 = Vgeo[gm->i];
}
//if ( (itag_true == i_0 && jtag_true == i_1) || (itag_true == i_1 && jtag_true == i_0) ) {
// CSVWriter csvWriter("/Users/willzunker/simulations/lammps/avicel_tableting/avicel_tableting_pair_info.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << itag_true << ", " << jtag_true << ", " << delta0 << ", " << delta1 << ", " << F0 << ", " << F1 << ", " << F << ", " << gm->radi << ", " << gm->radj << ", " << wij << ", " << psi_0 << ", " << psi_1 << ", " << aAdh0 << ", " << aAdh1 << ", " << Ac0 << ", " << Ac1 << ", " << F_MDR0 << ", " << F_MDR1 << ", " << F_BULK0 << ", " << F_BULK1 << ", " << deltaBULK0 << ", " << deltaBULK1 << ", " << Acon_0 << ", " << Acon_1 << ", " << Vgeo_0 << ", " << Vgeo_1 << ", " << kappa << ", " << lmp->update->ntimestep;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//if ( (gm->contact_type == PAIR) ) {
// CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/sticky/pair_info.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << itag_true << ", " << jtag_true << ", " << delta0 << ", " << delta1 << ", " << F0 << ", " << F1 << ", " << F << ", " << gm->radi << ", " << gm->radj << ", " << wij << ", " << psi_0 << ", " << psi_1 << ", " << aAdh0 << ", " << aAdh1 << ", " << Ac0 << ", " << Ac1 << ", " << F_MDR0 << ", " << F_MDR1 << ", " << F_BULK0 << ", " << F_BULK1 << ", " << deltaBULK0 << ", " << deltaBULK1 << ", " << Acon_0 << ", " << Acon_1 << ", " << Vgeo_0 << ", " << Vgeo_1 << ", " << kappa << ", " << F_DAMP << ", " << gm->vnnr << ", " << delta << ", " << lmp->update->ntimestep;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
return F;
}

3
src/GRANULAR/gran_sub_mod_normal.h
View File

@ -147,7 +147,7 @@ namespace Granular_NS {
protected:
double G, kappa, Eeff; // derived coeffs
double Eeffsq, Eeffinv, Eeffsqinv;
double Eeffsq, Eeffinv, Eeffsqinv, Eeff2particle;
double gammasq, gamma3, gamma4;
int warn_flag;
@ -155,6 +155,7 @@ namespace Granular_NS {
int index_Ro, index_Vgeo, index_Velas, index_Vcaps, index_eps_bar, index_dRnumerator;
int index_dRdenominator, index_Acon0, index_Acon1, index_Atot, index_Atot_sum, index_ddelta_bar;
int index_psi, index_sigmaxx, index_sigmayy, index_sigmazz, index_contacts, index_adhesive_length;
int index_dRavg;
int fix_mdr_flag;
char *id_fix;

4
src/GRANULAR/granular_model.cpp
View File

@ -251,8 +251,8 @@ void GranularModel::init()
// Must have valid normal, damping, and tangential models
if (normal_model->name == "none") error->all(FLERR, "Must specify normal granular model");
if (normal_model->name == "mdr") {
if (damping_model->name != "none")
error->all(FLERR, "MDR require 'none' damping model. To damp, specify a coefficient of restitution < 1.");
//if (damping_model->name != "none")
//error->all(FLERR, "MDR require 'none' damping model. To damp, specify a coefficient of restitution < 1.");
} else {
if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model");
}

25
src/csv_writer.h Normal file
View File

@ -0,0 +1,25 @@
#include <iostream>
#include <fstream>
#include <string>
class CSVWriter {
public:
CSVWriter(const std::string& filename) : filename_(filename) {}
void writeRow(const std::string& data) {
std::ofstream file;
// Use the append mode to add data to the end of the file if it exists
file.open(filename_, std::ios::out | std::ios::app);
if (!file.is_open()) {
std::cerr << "Failed to open file: " << filename_ << std::endl;
return;
}
file << data << std::endl;
file.close();
}
private:
std::string filename_;
};