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

cleaned up comments

Browse Source
This commit is contained in:
William Zunker
2024-12-17 11:18:26 -05:00
parent 97f19d9d54
commit a7ba185a4c
16 changed files with 56 additions and 552 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
examples/granular/in.mpfem.triaxial12particles
View File

@ -12,7 +12,7 @@ timestep 1e-6
######################### SIMULATION BOUNDING BOX, INTEGRATION, AND, GRAVITY ###########################
boundary f f f
read_data spheres.data
read_data spheres12.data
fix integr all nve/sphere
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################

0
examples/granular/spheres.data → examples/granular/spheres12.data
View File

21
examples/granular/spheresSTLgeneration.data
View File

@ -1,21 +0,0 @@
#LAMMPS data file created by matlab.
10 atoms
1 atom types
-0.005000 0.005000 xlo xhi
-0.005000 0.005000 ylo yhi
-0.001000 0.020000 zlo zhi
Atoms
1 1 0.000249 1560.000000 0.000469 -0.000120 0.008380
2 1 0.000260 1560.000000 -0.001811 0.000425 0.004090
3 1 0.000225 1560.000000 -0.001371 0.000426 0.002952
4 1 0.000252 1560.000000 0.000928 -0.000089 0.004118
5 1 0.000210 1560.000000 0.000208 0.000384 0.000842
6 1 0.000218 1560.000000 -0.002456 -0.000587 0.003834
7 1 0.000201 1560.000000 -0.000278 0.003459 0.003445
8 1 0.000231 1560.000000 -0.002123 0.003188 0.004799
9 1 0.000268 1560.000000 -0.000338 -0.001557 0.007382
10 1 0.000244 1560.000000 0.002537 0.002151 0.005999

1
src/GRANULAR/fix_wall_gran.cpp
View File

@ -496,7 +496,6 @@ void FixWallGran::post_force(int /*vflag*/)
model->dx[2] = dz;
model->radi = radius[i];
model->radj = rwall;
if (model->beyond_contact) model->touch = history_one[i][0];
touchflag = model->check_contact();

4
src/GRANULAR/fix_wall_gran.h
View File

@ -51,8 +51,8 @@ class FixWallGran : public Fix {
void reset_dt() override;
// for granular model choices
class Granular_NS::GranularModel *model; // MOVED HERE FROM PROTECTED FOR MDR MODEL
void clear_stored_contacts(); // MOVED HERE FROM PROTECTED FOR MDR MODEL
class Granular_NS::GranularModel *model;
void clear_stored_contacts();
protected:
int wallstyle, wiggle, wshear, axis;

7
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -29,7 +29,6 @@
#include "region.h"
#include "update.h"
#include "variable.h"
#include <iostream>
using namespace LAMMPS_NS;
using namespace Granular_NS;
@ -130,8 +129,6 @@ void FixWallGranRegion::post_force(int /*vflag*/)
if (update->setupflag) history_update = 0;
model->history_update = history_update;
//std::cout << (uint64_t)this << ", " << (uint64_t)model << std::endl;
// if just reneighbored:
// update rigid body masses for owned atoms if using FixRigid
// body[i] = which body atom I is in, -1 if none
@ -231,8 +228,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
model->radi = radius[i];
model->radj = region->contact[ic].radius;
model->r = region->contact[ic].r;
model->i = i; // Added for MDR
model->j = ic; // Added for MDR
model->i = i;
model->j = ic;
if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0];

16
src/GRANULAR/fix_wall_gran_region.h
View File

@ -44,13 +44,13 @@ class FixWallGranRegion : public FixWallGran {
int size_restart(int) override;
int maxsize_restart() override;
class Region *region; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
void update_contacts(int, int); // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
int *ncontact; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
int **walls; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
int *c2r; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
double ***history_many; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
int tmax; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
class Region *region;
void update_contacts(int, int);
int tmax; // max # of region walls one particle can touch
int *ncontact; // # of shear contacts per particle
int **walls; // which wall each contact is with
double ***history_many; // history per particle per contact
int *c2r; // contact to region mapping
private:
@ -59,8 +59,6 @@ class FixWallGranRegion : public FixWallGran {
// shear history for multiple contacts per particle
// c2r[i] = index of Ith contact in
// region-contact[] list of contacts
int motion_resetflag; // used by restart to indicate that region

282
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -17,11 +17,10 @@
#include "granular_model.h"
#include "math_const.h"
#include "atom.h"
#include "csv_writer.h"
#include "update.h"
#include "citeme.h"
#include <cmath>
#include <iostream>
#include <iomanip>
#include <sstream>
@ -38,6 +37,17 @@ static constexpr double INVROOT6 = 0.40824829046386307274; // 1/sqrt(6)
static constexpr double FOURTHIRDS = (4.0 / 3.0); // 4/3
static constexpr double JKRPREFIX = 1.2277228507842888; // cbrt(3*PI**2/16)
static const char cite_mdr[] =
"MDR contact model command: \n\n"
"@Article\n"
" author = \n"
" title = \n"
" journal = ,\n"
" year = ,\n"
" volume = ,\n"
" pages = \n"
"}\n\n";
/* ----------------------------------------------------------------------
Default normal model
------------------------------------------------------------------------- */
@ -399,6 +409,8 @@ void GranSubModNormalJKR::set_fncrit()
GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
GranSubModNormal(gm, lmp)
{
if (lmp->citeme) lmp->citeme->add(cite_mdr);
num_coeffs = 6; // Young's Modulus, Poisson's ratio, yield stress, effective surface energy, psi_b, coefficent of restitution
contact_radius_flag = 1;
size_history = 26;
@ -441,14 +453,15 @@ double GranSubModNormalMDR::calculate_forces()
// The MDR contact model was developed by imagining individual particles being
// squished between a number of rigid flats (references below). To allow
// for many interacting particles, we extend the idea of isolated particles surrounded
// by rigid flats. In particular, we imagine placing rigid flats at the overlap
// midpoints between particles. The force is calculated seperately on both sides
// by rigid flats. In particular, we imagine placing rigid flats at the overlaps
// between particles. The force is calculated seperately on both sides
// of the contact assuming interaction with a rigid flat. The two forces are then
// averaged on either side of the contact to determine the final force. If the
// contact is between a particle and wall then only one force evaluation is required.
//
// Zunker and Kamrin, 2024, Part I: https://doi.org/10.1016/j.jmps.2023.105492
// Zunker and Kamrin, 2024, Part II: https://doi.org/10.1016/j.jmps.2023.105493
// Zunker, Dunatunga, Thakur, Tang, and Kamrin, 2025:
const int itag_true = atom->tag[gm->i]; // true i particle tag
const int jtag_true = atom->tag[gm->j]; // true j particle tag
@ -466,10 +479,6 @@ double GranSubModNormalMDR::calculate_forces()
int i1 = 0;
double delta = gm->delta; // apparent overlap
//if (gm->contact_type == PAIR) delta = gm->delta/2.0; // half displacement to imagine interaction with rigid flat
//std::cout << "Normal force is called for: " << i_true << ", " << j_true << std::endl;
//std::cout << "Contact model has been entered " << gm->contact_type << ", " << PAIR << ", " << WALL << ", " << WALLREGION << ", " << gm->itype << ", " << gm->jtype << ", " << gm->delta << std::endl;
// initialize indexing in history array of different constact history variables
const int delta_offset_0 = 0; // apparent overlap
const int delta_offset_1 = 1;
@ -538,7 +547,6 @@ double GranSubModNormalMDR::calculate_forces()
double * contacts = atom->dvector[index_contacts];
double * adhesive_length = atom->dvector[index_adhesive_length];
double * history = & gm->history[history_index]; // load in all history variables
// Rigid flat placement scheme
@ -604,7 +612,6 @@ double GranSubModNormalMDR::calculate_forces()
double deltap = deltap0 + deltap1;
delta = delta_geo + (deltap0 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
//std::cout << "CS 0: " << gm->radi << ", " << gm->radj << ", gm->delta " << gm->delta << ", delta " << delta << ", deltamax " << deltamax << ", delta_geo " << delta_geo << ", delta_geo_alt " << delta_geo_alt << ", delta_geo/Ri " << delta_geo/gm->radi << ", delta_geo_alt/Rj " << delta_geo_alt/gm->radj << ", delta_geoOpt1 " << delta_geoOpt1 << ", delta_geoOpt2 " << delta_geoOpt2 << ", deltamax-delta " << (deltamax-gm->delta) << std::endl;
}
delta_offset = & history[delta_offset_0];
deltao_offset = & history[deltao_offset_0];
@ -648,8 +655,6 @@ double GranSubModNormalMDR::calculate_forces()
double deltap = deltap0 + deltap1;
delta = delta_geo + (deltap1 - delta_geo)/(deltap - deltamax)*(gm->delta-deltamax);
//std::cout << "CS 1: " << gm->radi << ", " << gm->radj << ", gm->delta " << gm->delta << ", delta " << delta << ", deltamax " << deltamax << ", delta_geo " << delta_geo << ", delta_geo_alt " << delta_geo_alt << ", delta_geo/Ri " << delta_geo/gm->radi << ", delta_geo_alt/Rj " << delta_geo_alt/gm->radj << ", delta_geoOpt1 " << delta_geoOpt1 << ", delta_geoOpt2 " << delta_geoOpt2 << ", deltamax-delta " << (deltamax-gm->delta) << std::endl;
delta_offset = & history[delta_offset_1];
deltao_offset = & history[deltao_offset_1];
delta_MDR_offset = & history[delta_MDR_offset_1];
@ -668,8 +673,6 @@ double GranSubModNormalMDR::calculate_forces()
const int i = gm->i;
const int j = gm->j;
//std::cout << lmp->update->ntimestep << std::endl;
// material and geometric property definitions
// E, nu, Y gamma , psi_b, and CoR are already defined.
const double G = E/(2.0*(1.0+nu)); // shear modulus
@ -714,12 +717,6 @@ double GranSubModNormalMDR::calculate_forces()
}
}
//if (i_true == 167 && j_true == 204) {
//std::cout << "i " << i << " | j " << j << " | delta_BULK: " << delta_BULK << " | delta_MDR " << delta_MDR << " | ddelta_BULK " << ddelta_BULK << " | ddelta_MDR " << ddelta_MDR << std::endl;
//}
//std::cout << "Yield Flag: " << *Yflag_offset << ", " << R << std::endl;
// MDR force calculation
double F_MDR;
double A; // height of elliptical indenter
@ -746,32 +743,21 @@ double GranSubModNormalMDR::calculate_forces()
// added for rigid flat placement
*deltap_offset = deltamax_MDR - (deltae1Dmax + deltaR);
//std::cout << *deltap_offset << std::endl;
}
//std::cout << psi_b << ", " << psi[i] << ", " << A << ", " << B << ", " << pY << ", " << amax << " || " << deltao << ", " << delta << ", " << ddelta << ", " << *delta_offset << ", " << ddelta_bar[i] << " || " << delta_MDR << ", " << ddelta_MDR << ", " << *delta_MDR_offset << ", " << deltamax_MDR << " || " << delta_BULK << ", " << ddelta_BULK << ", " << *delta_BULK_offset << " || " << R << std::endl;
//std::cout << i << ", " << j << ", " << A << ", " << B << " || " << deltao << ", " << delta << ", " << ddelta << ", " << R << ", " << MY_PI*pow(amax,2.0) << std::endl;
double a_na;
double a_fac = 0.99;
(deltae1D >= 0.0) ? a_na = B*sqrt(A - deltae1D)*sqrt(deltae1D)/A : a_na = 0.0;
double aAdh = *aAdh_offset;
if (aAdh > a_fac*amax) aAdh = a_fac*amax;
//if (i_true == 4 && j_true == 52){
//std::cout << "CS: " << contactSide << ", aAdh: " << aAdh << ", deltae1D: " << deltae1D << ", A: " << A << ", B:" << B << ", amax: " << amax << ", deltae1D: " << deltae1D << ", R: " << R << std::endl;
//}
if ( gamma > 0.0 ) { // adhesive contact
double g_aAdh;
if (delta_MDR == deltamax_MDR || a_na >= aAdh ) { // case 1: no tensile springs, purely compressive contact
(deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0)));
if ( std::isnan(F_MDR) ) {
std::cout << "F_MDR is NaN, case 1: no tensile springs" << std::endl;
//std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl;
std::cout << "i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", contact type: " << gm->contact_type << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << std::endl;
std::exit(1);
error->one(FLERR, "F_MDR is NaN, case 1: no tensile springs");
}
*aAdh_offset = a_fac*a_na;
} else {
@ -784,27 +770,19 @@ double GranSubModNormalMDR::calculate_forces()
pow(Eeff,2),0.3333333333333333))/pow(A,2))/6;
if ( (deltae1D + lmax - g_aAdh) >= 0.0) { // case 2: tensile springs, but not exceeding critical length --> deltae + lmax - g(aAdhes) >= 0
//if (contactSide == 0) {
//std::cout << "Case 2 tensile springs not exceeding critical length, R " << R << "deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh" << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl;
//}
const double deltaeAdh = g_aAdh;
const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0)));
const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh;
F_MDR = F_na + F_Adhes;
if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 2: tensile springs, but not exceeding critical length" << std::endl;
if ( std::isnan(F_MDR) ) error->one(FLERR, "F_MDR is NaN, case 2: tensile springs, but not exceeding critical length");
} else { // case 3: tensile springs exceed critical length --> deltae + lmax - g(aAdhes) = 0
//if (contactSide == 0) {
//std::cout << "Case 3 tensile springs exceed critical length, R " << R << " deltae1D " << deltae1D << " , lmax " << lmax << " , g_adh " << g_aAdh << " sum " << (deltae1D + lmax - g_aAdh) << std::endl;
//}
if ( aAdh < acrit ) {
aAdh = 0.0;
F_MDR = 0.0;
} else {
// newton-raphson to find aAdh
const double maxIterations = 100;
const double error = 1e-10;
const double error1 = 1e-10;
const double error2 = 1e-16;
double aAdh_tmp = aAdh;
double fa;
@ -812,7 +790,7 @@ double GranSubModNormalMDR::calculate_forces()
double dfda;
for (int lv1 = 0; lv1 < maxIterations; ++lv1) {
fa = deltae1D + sqrt(2.0*MY_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) );
if (abs(fa) < error) {
if (abs(fa) < error1) {
break;
}
dfda = -((aAdh_tmp*A)/(B*sqrt(-pow(aAdh_tmp,2.0) + pow(B,2.0)/4.0))) + (gamma*sqrt(MY_PI/2.0))/(Eeff*sqrt((aAdh_tmp*gamma)/Eeff));
@ -832,7 +810,7 @@ double GranSubModNormalMDR::calculate_forces()
const double F_na = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltaeAdh/A) - (1.0 - 2.0*deltaeAdh/A)*sqrt(4.0*deltaeAdh/A - 4.0*pow(deltaeAdh,2.0)/pow(A,2.0)));
const double F_Adhes = 2.0*Eeff*(deltae1D - deltaeAdh)*aAdh;
F_MDR = F_na + F_Adhes;
if ( std::isnan(F_MDR) ) std::cout << "F_MDR is NaN, case 3: tensile springs exceed critical length" << std::endl;
if ( std::isnan(F_MDR) ) error->one(FLERR, "F_MDR is NaN, case 3: tensile springs exceed critical length");
}
*aAdh_offset = aAdh;
}
@ -841,15 +819,10 @@ double GranSubModNormalMDR::calculate_forces()
*aAdh_offset = a_na;
(deltae1D <= 0.0) ? F_MDR = 0.0 : F_MDR = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1D/A) - (1.0 - 2.0*deltae1D/A)*sqrt(4.0*deltae1D/A - 4.0*pow(deltae1D,2.0)/pow(A,2.0)));
if ( std::isnan(F_MDR) ) {
std::cout << "F_MDR is NaN, non-adhesive case" << std::endl;
//std::cout << "Normal model: " << gm->delta << ", " << ddelta << ", " << gm->radi << ", " << gm->radj << " | delta: " << delta0 << ", " << delta1 << " | delta2_offset: " << *delta2_offset0 << ", " << *delta2_offset1 << "| dde: " << dde0 << ", " << dde1 << "| Fold: " << F0old << ", " << F1old << " | a: " << a0 << ", " << a1 << " | k_BULK: " << k_BULK0 << ", " << k_BULK1 << " | h_BULK: " << h_BULK0 << ", " << h_BULK1 << std::endl;
std::cout << "i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << std::endl;
std::exit(1);
error->one(FLERR, "F_MDR is NaN, non-adhesive case");
}
}
//std::cout << gm->i << ", " << gm->j << ", aAdh_offset: " << *aAdh_offset << ", aAdh: " << aAdh << ", a_na: " << a_na << std::endl;
contacts[i] += 1;
adhesive_length[i] += aAdh;
@ -858,8 +831,6 @@ double GranSubModNormalMDR::calculate_forces()
double pij = *penalty_offset;
const double wij = std::max(1.0-pij,0.0);
//std::cout << gm->i << ", " << gm->j << ", " << gm->contact_type << ", " << *gm->xi[1] << ", " << *gm->xj[1] << std::endl;
// area related calculations
double Ac;
(*Yflag_offset == 0.0) ? Ac = MY_PI*delta*R : Ac = MY_PI*((2.0*delta*R - pow(delta,2.0)) + cA/MY_PI);
@ -871,80 +842,9 @@ double GranSubModNormalMDR::calculate_forces()
double F_BULK;
(delta_BULK <= 0.0) ? F_BULK = 0.0 : F_BULK = (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac;
//if (atom->tag[i_true] == 4135 && lmp->update->ntimestep > 10935000 && gm->contact_type == 2) {
// std::cout << "CS: " << contactSide << ", i_true: " << i_true << ", j_true: " << j_true << ", i_tag: " << atom->tag[i_true] << ", j_tag: " << atom->tag[j_true] << ", deltae1D: " << deltae1D << ", A: " << A << ", B: " << B << ", amax: " << amax << ", deltamax_MDR: " << deltamax_MDR << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
//}
//if (i == 35 && lmp->update->ntimestep % 1000 == 0) {
//double **x = atom->x;
//const double xi = x[i_true][0];
// std::cout << "i: " << i << ", j: " << j << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", j_true: " << j_true << ", j_tag " << atom->tag[j_true] << ", radi_true: " << radi_true << ", radj_true " << radj_true << ", gm->radi: " << gm->radi << ", gm->radj: " << gm->radj << ", R: " << R << ", Ro: " << Ro << std::endl;
//}
//if (i_true == 35 && lmp->update->ntimestep >= 736002 && lmp->update->ntimestep <= 736005) {
// //double **x = atom->x;
// //const double xi = x[i_true][0];
// std::cout << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", R: " << R << ", Ro: " << Ro << std::endl;
//}
//if ( ((atom->tag[i_true] == 4 && atom->tag[j_true] == 11) || (atom->tag[i_true] == 11 && atom->tag[j_true] == 4)) && lmp->update->ntimestep == 45000) {
// std::cout << "CS: " << contactSide << ", contact_type: " << gm->contact_type << ", pair: " << PAIR << ", wall: " << WALL << ", WALLREGION " << WALLREGION << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
//}
//if ( ((atom->tag[i_true] == 4301) || (atom->tag[j_true] == 4076)) && lmp->update->ntimestep > 1016700) {
// std::cout << "i: " << i << ", j: " << j << "i_true: " << i_true << ", i_tag: " << atom->tag[i_true] << ", j_true: " << j_true << ", j_tag " << atom->tag[j_true] << ", radi_true: " << radi_true << ", radj_true " << radj_true << ", gm->radi: " << gm->radi << ", gm->radj: " << gm->radj << ", R: " << R << ", Ro: " << Ro << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << std::endl;
//}
//if ( ((atom->tag[i_true] == 4) || (atom->tag[j_true] == 4)) && lmp->update->ntimestep == 45000) {
// int rank = 0;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
// double **x = atom->x;
// const double xi = x[i][0];
// const double xj = x[j][0];
// const double yi = x[i][1];
// const double yj = x[j][1];
// const double zi = x[i][2];
// const double zj = x[j][2];
// std::cout << "CS: " << contactSide << ", rank, " << rank << ", CT: " << gm->contact_type << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", i: " << i << ", j: " << j << ", nlocal: " << atom->nlocal << ", nghost: " << atom->nghost << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << ", R: " << R << ", xi: " << xi << ", xj: " << xj << ", yi: " << yi << ", yj: " << yj << ", zi: " << zi << ", zj: " << zj << std::endl;
//}
//int rank = 0;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//double **x = atom->x;
//const double xi = x[i][0];
//const double xj = x[j][0];
//const double yi = x[i][1];
//const double yj = x[j][1];
//const double zi = x[i][2];
//const double zj = x[j][2];
//const double dis = sqrt(pow((xi-xj),2.0) + pow((yi-yj),2.0) + pow((zi-zj),2.0));
//const double delta_test = gm->radi + gm->radj - dis;
//if (delta_test < 0.0 && gm->contact_type != 2) {
// std::cout << "Particles are not touching but a force is evaluated, CS: " << contactSide << ", rank, " << rank << ", contact_type: " << gm->contact_type << ", pair: " << PAIR << ", wall: " << WALL << ", WALLREGION " << WALLREGION << ", itag_true: " << atom->tag[i_true] << ", jtag_true: " << atom->tag[j_true] << ", i: " << i << ", j: " << j << ", nlocal: " << atom->nlocal << ", nghost: " << atom->nghost << ", wij: " << wij << ", gm->delta: " << gm->delta << ", delta: " << delta << ", delmax: " << deltamax << ", deltap: " << *deltap_offset << ", R: " << R << ", xi: " << xi << ", xj: " << xj << ", yi: " << yi << ", yj: " << yj << ", zi: " << zi << ", zj: " << zj << std::endl;
// std::exit(1);
//}
//if (F_BULK > 0.0) {
// std::cout << "F_BULK is: " << F_BULK << std::endl;
//}
//std::cout << delta_BULK << ", " << F_BULK << ", " << (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac << std::endl;
//int rank = 0;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//std::cout << "Step: " << lmp->update->ntimestep << ", CS: " << contactSide << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", rank: " << rank << ", gm->delta: " << gm->delta << ", delta: " << delta << ", ddelta: " << ddelta << ", delta_bar: " << ddelta_bar[i] << ", R: " << R << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << std::endl;
//std::cout << gm->i << ", " << gm->j << ", " << Vgeo[i] << ", " << Acon0[i] << ", " << Acon1[i] << ", " << Ac << ", " << kappa << " || " << psi[i] << ", " << ddelta_bar[i] << ", " << ddelta << ", " << ddelta_MDR << ", " << ddelta_BULK << ", " << delta << ", " << delta_MDR << ", " << delta_BULK << ", " << F_MDR << ", " << F_BULK << ", " << R << " || " << deltae1D << ", " << A << ", " << B << std::endl;
//std::cout << gm->i << ", " << gm->j << ", " << (1.0/Vgeo[i])*Acon0[i]*delta_BULK*kappa*Ac << std::endl;
// total force calculation
(contactSide == 0) ? F0 = F_MDR + F_BULK : F1 = F_MDR + F_BULK;
//std::cout << gm->i << ", " << gm->j << " | " << deltao << ", " << ddelta_bar[i] << ", " << R << ", " << psi[i] << ", " << psi_b << ", " << Ac << " | " << pij << ", " << wij << std::endl;
// mean surface dipslacement calculation
*Ac_offset = wij*Ac;
@ -961,64 +861,17 @@ double GranSubModNormalMDR::calculate_forces()
const double eps_bar_contact = (1.0/(3*kappa*Velas[i]))*(fx*bx + fy*by + fz*bz);
eps_bar[i] += eps_bar_contact;
//double **x = atom->x;
//const double xi = x[i_true][0];
//const double xj = x[j_true][0];
//std::cout << i_true << ", " << j_true << ", " << xi << ", " << xj << ", " << gm->nx[0] << ", " << gm->nx[1] << ", " << gm->nx[2] << std::endl;
//if ( (i == 0 && j == 2 && gm->contact_type == 0) || (i == 2 && j == 0 && gm->contact_type == 0)) {
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << delta << ", " << *delta_offset << ", " << (uintptr_t)(delta_offset) << " || " << deltao << ", " << *deltao_offset << ", " << (uintptr_t)(deltao_offset) << " || " << delta_MDR << ", " << *delta_MDR_offset << ", " << (uintptr_t)(delta_MDR_offset) << " || " << *Yflag_offset << ", " << (uintptr_t)(Yflag_offset) << " || " << R << std::endl;
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << fx << ", " << fy << ", " << fz << " || " << bx << ", " << by << ", " << bz << ", " << Velas[i] << std::endl;
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << eps_bar_contact << ", " << *eps_bar_offset << ", " << (uintptr_t)(eps_bar_offset) << " || " << wij << ", " << ddeltao << ", " << deltao << ", " << delta << ", " << *delta_offset << " || " << Ro << ", " << R << std::endl;
//}
//if () {
// std::cout << j << ", " << -Vo*(eps_bar_contact - *eps_bar_offset) - wij*MY_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) ) << ", " << -Vo*(eps_bar_contact - *eps_bar_offset) << ", " << wij*MY_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) ) << std::endl;
//std::cout << i << ", " << j << ", " << gm->contact_type << " || " << eps_bar_contact << ", " << *eps_bar_offset << ", " << (uintptr_t)(eps_bar_offset) << " || " << wij << ", " << ddeltao << ", " << deltao << " || " << Ro << ", " << R << std::endl;
//}
double desp_bar_contact = eps_bar_contact - *eps_bar_offset; // && desp_bar_contact < 0.0
double desp_bar_contact = eps_bar_contact - *eps_bar_offset;
if(delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){
const double Vo = (4.0/3.0)*MY_PI*pow(Ro,3.0);
dRnumerator[i] += -Vo*(eps_bar_contact - *eps_bar_offset) - wij*MY_PI*ddeltao*( 2.0*deltao*Ro - pow(deltao,2.0) + pow(R,2.0) - pow(Ro,2.0) );
dRdenominator[i] += wij*2.0*MY_PI*R*(deltao + R - Ro);
//if ( (atom->tag[i] == 9) ) {
// std::cout << "CT: " << gm->contact_type << ", " << PAIR << "i_tag: " << atom->tag[i] << ", j_tag: " << atom->tag[j] << ", deltae1D: " << deltae1D << ", R: " << R << ", Ro: " << Ro << ", F_MDR: " << F_MDR << ", F_BULK: " << F_BULK << ", wij: " << wij << ", deltao: " << deltao << ", ddeltao: " << ddeltao << ", desp_bar: " << eps_bar_contact - *eps_bar_offset << std::endl;
//}
}
*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);
//std::cout << psi_b << ", " << psi[i] << ", " << A << ", " << B << ", " << pY << ", " << amax << " || " << deltao << ", " << delta << ", " << ddelta << ", " << *delta_offset << ", " << ddelta_bar[i] << " || " << delta_MDR << ", " << ddelta_MDR << ", " << *delta_MDR_offset << ", " << deltamax_MDR << " || " << delta_BULK << ", " << ddelta_BULK << ", " << *delta_BULK_offset << " || " << R << " || " << Ac << ", " << *Ac_offset << ", " << Acon0[i] << ", " << Acon1[i] << " || " << F_MDR << ", " << F_BULK << ", " << Vgeo[i] << std::endl;
//std::cout << gm->i << ", " << gm->j << ", " << gm->radi << ", " << gm->radj << " | " << delta << ", " << F_MDR << " | " << deltae1D << ", " << A << ", " << B << std::endl;
//if (atom->tag[i] == 1 && lmp->update->ntimestep == 45000) {
// double nx;
// double ny;
// double nz;
// if (i == j_true) {
// nx = gm->nx[0];
// ny = gm->nx[1];
// nz = gm->nx[2];
// } else {
// nx = -gm->nx[0];
// ny = -gm->nx[1];
// nz = -gm->nx[2];
// }
// double deltae = deltamax_MDR - *deltap_offset;
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/MPFEM/deformed_particle_shape_data.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision
// rowDataStream << nx << ", " << ny << ", " << nz << ", " << Ro << ", " << R << ", " << delta << ", " << deltae << ", " << A << ", " << B << ", " << a_na;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
}
gm->i = i_true;
@ -1031,24 +884,13 @@ double GranSubModNormalMDR::calculate_forces()
const double wij = std::max(1.0-pij,0.0);
*penalty_offset = 0.0;
//std::cout << gm->i << ", " << gm->j << ", " << xi << ", " << xj << std::endl;
// force magnifiers to prevent over penetration
double * deltao_offset = & history[deltao_offset_0];
//const double wallForceMagnifer = std::exp(10.0*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0;
const double wallForceMagnifer = 1.0;
// assign final force
//(gm->contact_type != PAIR) ? F = wij*F0*wallForceMagnifer : F = wij*(F0 + F1)/2; // F = 2*wij*pow(1/F0 + 1/F1,-1);
if (gm->contact_type != PAIR) {
F = wij*F0*wallForceMagnifer;
F = wij*F0;
} else {
F = wij*(F0 + F1)/2.0;
}
//std::cout << F << ", " << F0 << ", " << F1 << " | " << R0 << ", " << R1 << std::endl;
// calculate damping force
if (F > 0.0) {
double Eeff;
@ -1065,53 +907,9 @@ double GranSubModNormalMDR::calculate_forces()
const double damp_prefactor = beta*sqrt(gm->meff*kn);
const double F_DAMP = -damp_prefactor*(gm->vnnr);
//std:: cout << gm->contact_type << ", " << Eeff << " , " << Reff << ", " << gm->radi << ", " << gm->radj << " || " << kn << ", " << beta << ", " << gm->meff << " || " << F_DAMP << ", " << F << std::endl;
F += wij*F_DAMP;
}
//double **x = atom->x;
//const double xi = x[gm->i][0];
//const double xj = x[gm->j][0];
//const double del = 20.0 - abs(xi-xj);
//if (i_true == 146 && j_true == 152) {
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/rigid_flat_output.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << dde0 << ", " << dde1 << ", " << (dde0 + dde1) << ", " << F0old << ", " << F1old << ", " << a0 << ", " << a1 << ", " << h0 << ", " << h1 << ", " << k_BULK0 << ", " << k_BULK1 << ", " << h_BULK0 << ", " << h_BULK1;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//if (i_true == 0 && j_true == 1) {
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/twoParticleDifferingRadii/rigid_flat_output.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << dde0 << ", " << dde1 << ", " << (dde0 + dde1) << ", " << F0old << ", " << F1old << ", " << a0 << ", " << a1 << ", " << h0 << ", " << h1 << ", " << k_BULK0 << ", " << k_BULK1 << ", " << h_BULK0 << ", " << h_BULK1 << ", " << R0 << ", " << R1;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//if (i_true == 0 && j_true == 1) {
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/twoParticleDifferingRadii/rigid_flat_output_fit.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(6); // Set the format and precision
// rowDataStream << gm->delta << ", " << delta0 << ", " << delta1 << ", " << delta01;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//if (gm->i == 0 && gm->j == 2) {
// CSVWriter csvWriter("/Users/willzunker/lammps/sims/compressionSleeve/pairContactsBotCen.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
// rowDataStream << del << ", " << F;
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//std:: cout << "The force F is: " << F << std::endl;
return F;
}
@ -1121,33 +919,3 @@ void GranSubModNormalMDR::set_fncrit()
{
Fncrit = fabs(F);
}
//std::cout << sidata.i << ", " << sidata.j << ", " << R << ", " << deltan << ", " << deltao << ", " << dRsums_i[0] << ", " << dRsums_i[1] << ", " << numQuant << std::endl;
//std::cout << gm->i << ", " << gm->j << " | " << gm->nx[0] << ", " << gm->nx[1] << ", " << gm->nx[2] << std::endl;
//std::cout << "F is: " << F << std::endl;
//std::cout << "Ro from particle history is: " << Ro[gm->i] << std::endl;
//std::cout << "MDR contact model has been entered." << std::endl;
// std::cout << "F is: " << F << std::endl;
//std::cout << "gamma > 0.0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
//std::cout << "deltae1D <= 0.0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
//std::cout << "F_MDR should be > 0: " << F_MDR << ", " << gm->i << ", " << gm->j << std::endl;
//std::cout << gm->i << ", " << gm->j << " || " << delta << ", " << delta_MDR << ", " << deltamax_MDR << ", " << deltae1D << " || " << A << ", " << B << std::endl;
//std::cout << i_true << ", " << j_true << std::endl;
//std::cout << history_index << ", " << history[0] << ", " << history[1] << ", " << history[2] << std::endl;
// initialize all history variables
//double delta_offset;
//double deltao_offset;
//double delta_MDR_offset;
//double delta_BULK_offset;
//double deltamax_MDR_offset;
//double Yflag;
//double deltaY_offset;
//double Ac_offset;
//double aAdh_offset;
//double deltap_offset;
//double cA_offset;
//double eps_bar_offset;
//double wall_contact_flag_offset;

3
src/GRANULAR/gran_sub_mod_rolling.cpp
View File

@ -17,7 +17,6 @@
#include "gran_sub_mod_normal.h"
#include "granular_model.h"
#include "math_extra.h"
#include <iostream>
#include <cmath>
@ -82,8 +81,6 @@ void GranSubModRollingSDS::calculate_forces()
hist_temp[1] = gm->history[rhist1];
hist_temp[2] = gm->history[rhist2];
//std::cout << "Frcrit rolling is: " << Frcrit << std::endl;
if (gm->history_update) {
rolldotn = dot3(hist_temp, gm->nx);

6
src/GRANULAR/granular_model.cpp
View File

@ -243,7 +243,11 @@ 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 (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); TURNED OFF FOR MDR
if (normal_model->name == "mdr") {
if (damping_model->name != "none") error->all(FLERR, "Damping model must be set to 'none' for the MDR contact model. Specify a coefficient of restitution CoR < 1 if normal damping is desired.");
} else {
if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model");
}
if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model");
// Twisting, rolling, and heat are optional

28
src/GRANULAR/pair_granular.cpp
View File

@ -39,7 +39,6 @@
#include "gran_sub_mod_normal.h"
#include <cstring>
#include <iostream>
using namespace LAMMPS_NS;
using namespace Granular_NS;
@ -490,28 +489,7 @@ void PairGranular::init_style()
if (!fix_history) error->all(FLERR,"Could not find pair fix neigh history ID");
}
/*
// Example
//Store persistent per atom quantities
if (! fix_flag) {
int tmp1, tmp2;
id_fix = "BOND_BPM_PLASTIC_FIX_PROP_ATOM";
modify->add_fix(
fmt::format("{} all property/atom d_plastic d_plastic_temp d_viscous_temp d_plastic_heat d_viscous_heat ghost yes", id_fix));
index_plastic = atom->find_custom("plastic",tmp1,tmp2);
index_pq = atom->find_custom("plastic_heat",tmp1,tmp2);
index_pt = atom->find_custom("plastic_temp",tmp1,tmp2);
index_vq = atom->find_custom("viscous_heat",tmp1,tmp2);
index_vt = atom->find_custom("viscous_temp",tmp1,tmp2);
fix_flag = 1;
}
*/
if (model->normal_model->name == "mdr") {
std::cout << "MDR history variables have been initialized" << std::endl;
// FOR MDR CONTACT MODEL
//Store persistent per atom quantities
if (! fix_flag) {
int tmp1, tmp2;
@ -537,11 +515,7 @@ void PairGranular::init_style()
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_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
std::cout << "MDR history variables have been initialized 2" << ", " << index_Ro << std::endl;
//index_volSums = atom->find_custom("volSums",tmp1,tmp2);
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");

1
src/GRANULAR/pair_granular.h
View File

@ -49,7 +49,6 @@ class PairGranular : public Pair {
size_t get_size_history() const;
// granular models
// MOVED HERE FROM PRIVATE FOR MDR MODEL
class Granular_NS::GranularModel** models_list;
int **types_indices;
int nmodels, maxmodels;

18
src/SRD/fix_srd.cpp
View File

@ -40,7 +40,6 @@
#include <cmath>
#include <cstring>
#include <iostream>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -55,7 +54,7 @@ enum { SHIFT_NO, SHIFT_YES, SHIFT_POSSIBLE };
static constexpr double EINERTIA = 0.2; // moment of inertia prefactor for ellipsoid
static constexpr int ATOMPERBIN = 100;
static constexpr int ATOMPERBIN = 30;
static constexpr double BIG = 1.0e20;
static constexpr int VBINSIZE = 5;
static constexpr double TOLERANCE = 0.00001;
@ -1322,7 +1321,6 @@ void FixSRD::collisions_single()
#endif
if (t_remain > dt) {
print_collision(i, j, ibounce, t_remain, dt, xscoll, xbcoll, norm, type);
ninside++;
if (insideflag == INSIDE_ERROR || insideflag == INSIDE_WARN) {
std::string mesg;
@ -2605,17 +2603,14 @@ void FixSRD::parameterize()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) {
if (mass_srd == 0.0) {
if (mass_srd == 0.0)
mass_srd = rmass[i];
std::cout << "mass if #1: " << rmass[i] << std::endl;
std::cout << "srd radius is: " << radius[i] << std::endl;
} else if (rmass[i] != mass_srd)
else if (rmass[i] != mass_srd)
flag = 1;
} else {
if (mass_srd == 0.0) {
if (mass_srd == 0.0)
mass_srd = mass[type[i]];
std::cout << "mass if #2: " << mass[type[i]] << std::endl;
} else if (mass[type[i]] != mass_srd)
else if (mass[type[i]] != mass_srd)
flag = 1;
}
}
@ -2717,7 +2712,6 @@ void FixSRD::parameterize()
if (dimension == 3) {
volsrd = (domain->xprd * domain->yprd * domain->zprd) - volbig;
density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd * domain->zprd - volbig);
std::cout << nsrd << ", " << mass_srd << ", " << domain->xprd << ", " << domain->yprd << ", " << domain->zprd << ", " << volbig << std::endl;
} else {
volsrd = (domain->xprd * domain->yprd) - volbig;
density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd - volbig);
@ -2779,7 +2773,7 @@ void FixSRD::parameterize()
gridsrd, binsize3x, binsize3y, binsize3z);
mesg += fmt::format(" SRD per actual grid cell = {:.8}\n", srd_per_cell);
mesg += fmt::format(" SRD viscosity = {:.8}\n", viscosity);
mesg += fmt::format(" big/SRD mass density ratio = {:.8} | big density = {:.8}, small density = {:.8}\n", mdratio,density_big,density_srd);
mesg += fmt::format(" big/SRD mass density ratio = {:.8}\n", mdratio);
utils::logmesg(lmp, mesg);
}

25
src/csv_writer.h
View File

@ -1,25 +0,0 @@
#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_;
};

81
src/fix_mdr_mean_surf_disp.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -37,7 +36,6 @@
#include "update.h"
#include "fix_wall_gran_region.h"
#include "comm.h"
#include <iostream>
#include <algorithm>
using namespace LAMMPS_NS;
@ -53,8 +51,6 @@ FixMDRmeanSurfDisp::FixMDRmeanSurfDisp(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 1; // value needs to match number of values you communicate
}
// FOR MDR
int FixMDRmeanSurfDisp::setmask()
{
int mask = 0;
@ -79,7 +75,6 @@ int FixMDRmeanSurfDisp::pack_forward_comm(int n, int *list, double *buf, int /*p
for (int i = 0; i < n; i++) {
int j = list[i];
buf[m++] = ddelta_bar[j];
//buf[m++] = Acon0[j];
}
return m;
}
@ -90,7 +85,6 @@ void FixMDRmeanSurfDisp::unpack_forward_comm(int n, int first, double *buf)
int last = first + n;
for (int i = first; i < last; i++) {
ddelta_bar[i] = buf[m++];
//Acon0[i] = buf[m++];
}
}
@ -102,9 +96,6 @@ void FixMDRmeanSurfDisp::pre_force(int)
const int size_history = pair->get_size_history();
//std::cout << " " << std::endl;
//std::cout << "New Step" << std::endl;
{
int i,j,k,lv1,ii,jj,inum,jnum,itype,jtype,ktype;
@ -114,10 +105,6 @@ void FixMDRmeanSurfDisp::pre_force(int)
bool touchflag = false;
//class GranularModel* model;
//class GranularModel** models_list = pair->models_list;
//int ** types_indices = pair->types_indices;
double **x = atom->x;
int *type = atom->type;
double *radius = atom->radius;
@ -186,8 +173,8 @@ void FixMDRmeanSurfDisp::pre_force(int)
history_ik = &allhistory[size_history * kk];
double * pik = &history_ik[22]; // penalty for contact i and k
// we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa, so we need to manually search to figure out the owner
// check if k is in the neighbor list of j
// we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa,
// so we need to manually search to figure out the owner check if k is in the neighbor list of j
double * pjk = NULL;
int * const jklist = firstneigh[j];
const int jknum = numneigh[j];
@ -197,10 +184,6 @@ void FixMDRmeanSurfDisp::pre_force(int)
allhistory_j = firsthistory[j];
history_jk = &allhistory_j[size_history * jk];
pjk = &history_jk[22]; // penalty for contact j and k
//int rank = 0;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//std::cout << "Print 183 pjk: " << rank << ", " << pjk << std::endl;
//std::cout << "Print 183 pjk[0]: " << rank << ", " << pjk[0] << std::endl;
break;
}
}
@ -215,16 +198,11 @@ void FixMDRmeanSurfDisp::pre_force(int)
allhistory_k = firsthistory[k];
history_kj = &allhistory_k[size_history * kj];
pjk = &history_kj[22]; // penalty for contact j and k
//int rank = 0;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//std::cout << "Print 198 pjk: " << rank << ", " << pjk << std::endl;
//std::cout << "Print 198 pjk[0]: " << rank << ", " << pjk[0] << std::endl;
break;
}
}
}
//std::cout << "Print: " << __LINE__ << std::endl;
std::vector<double> distances = {r_ij,r_ik,r_jk};
auto maxElement = std::max_element(distances.begin(), distances.end());
double maxValue = *maxElement;
@ -256,9 +234,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double eny_ik = dely_ik * rinv_ik;
const double enz_ik = delz_ik * rinv_ik;
const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik);
//std::cout << "Print: " << __LINE__ << std::endl;
pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) );
//std::cout << "Print: " << __LINE__ << std::endl;
}
}
}
@ -334,8 +310,6 @@ void FixMDRmeanSurfDisp::pre_force(int)
history = &allhistory[size_history * jj];
model->history = history;
//const double pij = history[22];
//const double wij = std::max(1.0-pij,0.0);
const double delta = model->radsum - sqrt(model->rsq);
const int delta_offset_0 = 0; // apparent overlap
@ -365,18 +339,6 @@ void FixMDRmeanSurfDisp::pre_force(int)
i1 = i;
}
int rank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int i_ghost;
int j_ghost;
(i >= atom->nlocal) ? i_ghost = 1 : i_ghost = 0;
(j >= atom->nlocal) ? j_ghost = 1 : j_ghost = 0;
if (i_ghost == 1) {
std::cout << "rank: " << rank << ", i: " << i << ", j: " << j << ", i is ghost? " << i_ghost << ", j is ghost? " << j_ghost << ", nlocal: " << atom->nlocal << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << std::endl;
}
//int i0 = std::max(i,j);
//int i1 = std::min(i,j);
double R0 = radius[i0];
double R1 = radius[i1];
@ -407,32 +369,13 @@ void FixMDRmeanSurfDisp::pre_force(int)
if (Acon0[i0] != 0.0) {
const double Ac_offset0 = history[Ac_offset_0];
ddelta_bar[i0] += Ac_offset0/Acon0[i0]*ddel0; // Multiply by 0.5 since displacement is shared equally between deformable particles.
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;
}
//int rank = 0;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//std::cout << "delta_bar calc: Step: " << lmp->update->ntimestep << ", itag: " << atom->tag[i] << ", jtag: " << atom->tag[j] << ", i: " << i << ", j: " << j << ", rank: " << rank << ", Ac_offset0: " << history[Ac_offset_0] << ", Ac_offset1: " << history[Ac_offset_1] << ", Acon[i0]: " << Acon0[i0] << ", Acon[i1]: " << Acon0[i1] << ", ddel0: " << ddel0 << ", ddel1: " << ddel1 << ", ddelta_bar[i0]: " << ddelta_bar[i0] << ", ddelta_bar[i1]: " << ddelta_bar[i1] << std::endl;
//if (Acon0[j] != 0.0) {
// const double delta_offset0 = history[0];
// const double ddelta = delta/2.0 - delta_offset0; // Divide by 2.0 since we are storing 1/2 deltan in main MDR script
// const double Ac_offset0 = history[18];
// ddelta_bar[j] += Ac_offset0/Acon0[j]*ddelta; // Multiply by 0.5 since displacement is shared equally between deformable particles.
//}
//
//if (Acon0[i] != 0.0) {
// const double delta_offset1 = history[1];
// const double ddelta = delta/2.0 - delta_offset1; // Divide by 2.0 since we are storing 1/2 deltan in main MDR script
// const double Ac_offset1 = history[19];
// ddelta_bar[i] += Ac_offset1/Acon0[i]*ddelta;
//}
}
}
}
@ -515,27 +458,11 @@ void FixMDRmeanSurfDisp::pre_force(int)
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; // Multiply by 0.5 since displacement is shared equally between deformable particles.
//std::cout << delta << ", " << delta_offset0 << " || " << Ac_offset0 << ", " << Acon0[i] << ", " << ddelta << ", " << ddelta_bar[i] << std::endl;
ddelta_bar[i] += wij*Ac_offset0/Acon0[i]*ddelta;
}
}
}
}
}
comm->forward_comm(this);
//and the function delcarations in the header:
//int pack_forward_comm(int, int *, double *, int, int *) override;
//void unpack_forward_comm(int, int, double *) override;
//Then the methods would look like::
//where comm_stage is a public flag to control hich quantity is being communicated
}
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*MY_PI*pow(R,2.0) << std::endl;
//std::cout << "Fix radius update setup has been entered !!!" << std::endl;
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;

107
src/fix_mdr_radius_update.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -28,14 +27,12 @@
#include "memory.h"
#include "modify.h"
#include "variable.h"
#include "csv_writer.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 <iostream>
#include <iomanip>
#include <sstream>
@ -52,8 +49,6 @@ FixMDRradiusUpdate::FixMDRradiusUpdate(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 20; // value needs to match number of values you communicate
}
// FOR MDR
int FixMDRradiusUpdate::setmask()
{
int mask = 0;
@ -167,43 +162,10 @@ void FixMDRradiusUpdate::pre_force(int)
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");
//model = models_list[0];
//class GranSubModNormalMDR* norm_model = dynamic_cast<GranSubModNormalMDR *>(model->normal_model);
//std::cout << "Preforce was called radius update" << std::endl;
// assign correct value to initially non-zero MDR particle history variables
//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_Atot = atom->find_custom("Atot",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);
//double * Ro = atom->dvector[index_Ro];
//double * Vgeo = atom->dvector[index_Vgeo];
//double * Velas = atom->dvector[index_Velas];
//double * Atot = atom->dvector[index_Atot];
//double * psi = atom->dvector[index_psi];
//double * psi_b = atom->dvector[index_psi_b];
//double * sigmaxx = atom->dvector[index_sigmaxx];
//double * sigmayy = atom->dvector[index_sigmayy];
//double * sigmazz = atom->dvector[index_sigmazz];
//double * history_setup_flag = atom->dvector[index_history_setup_flag];
//double * contacts = atom->dvector[index_contacts];
//double * adhesive_length = atom->dvector[index_adhesive_length];
double *radius = atom->radius;
int nlocal = atom->nlocal;
//std::cout << "ntotal " << ntotal << ", nlocal " << atom->nlocal << ", nghost " << atom->nghost << std::endl;
for (int i = 0; i < nlocal; i++) {
if (history_setup_flag[i] < 1e-16) {
Ro[i] = radius[i];
@ -220,9 +182,7 @@ void FixMDRradiusUpdate::pre_force(int)
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*/)
@ -284,21 +244,9 @@ void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf)
void FixMDRradiusUpdate::end_of_step()
{
//comm->forward_comm(this);
//std::cout << "end_of_step() was called radius update" << std::endl;
// update the apparent radius of every particle
double *radius = atom->radius;
int nlocal = atom->nlocal;
double sigmaxx_sum = 0.0;
double sigmayy_sum = 0.0;
double sigmazz_sum = 0.0;
double Vparticles = 0.0;
//std::cout << "New step" << std::endl;
for (int i = 0; i < nlocal; i++) {
@ -308,32 +256,15 @@ void FixMDRradiusUpdate::end_of_step()
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);
//(Vgeoi < Vo) ? Vgeo[i] = Vgeoi : Vgeo[i] = Vo;
const double Afree = Atot[i] - Acon1[i];
psi[i] = Afree/Atot[i];
//if (psi[i] < 0.08) {
// std::cout << "psi is: " << psi[i] << std::endl;
//}
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;
//radius[i] += dR;
}
//if (dR > 1e-7) {
// std::cout << "big dR change" << std::endl;
//}
//if (atom->tag[i] == 9){
// std::cout << i << ", radius: " << radius[i] << ", dR: " << dR << ", dRnum: " << dRnumerator[i] << ", dRdenom " << dRdenominator[i] << ", dRdem_full " << dRdenominator[i] - 4.0*MY_PI*pow(R,2.0) << std::endl;
//}
//std::cout << i << ", " << radius[i] << " | " << psi_b[i] << ", " << psi[i] << " | " << Acon1[i] << ", " << Atot[i] << std::endl;
Velas[i] = Vo*(1.0 + eps_bar[i]);
Vcaps[i] = 0.0;
eps_bar[i] = 0.0;
@ -341,46 +272,8 @@ void FixMDRradiusUpdate::end_of_step()
dRdenominator[i] = 0.0;
Acon0[i] = Acon1[i];
Acon1[i] = 0.0;
//std::cout << "Acon reset: " << Acon0[i] << ", " << Acon1[i] << std::endl;
Atot_sum[i] = 0.0;
ddelta_bar[i] = 0.0;
//adhesive_length[i] = adhesive_length[i]/contacts[i]; // convert adhesive length to average aAdh for each particle
//if (lmp->update->ntimestep == 487500) {
// double **x = atom->x;
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/atom_stresses_at_mid_compression.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
// rowDataStream << sigmaxx[i] << ", " << sigmayy[i] << ", " << sigmazz[i] << ", " << Velas[i] << ", " << x[i][2];
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//
//if (lmp->update->ntimestep == 595750) {
// double **x = atom->x;
// CSVWriter csvWriter("/Users/willzunker/lammps_mdr_develop/sims/avicelTableting/atom_stresses_at_max_compression.csv");
// std::stringstream rowDataStream;
// rowDataStream << std::scientific << std::setprecision(4); // Set the format and precision
// rowDataStream << sigmaxx[i] << ", " << sigmayy[i] << ", " << sigmazz[i] << ", " << Velas[i] << ", " << x[i][2];
// std::string rowData = rowDataStream.str();
// csvWriter.writeRow(rowData);
//}
//sigmaxx_sum += sigmaxx[i];
//sigmayy_sum += sigmayy[i];
//sigmazz_sum += sigmazz[i];
//Vparticles += Velas[i];
}
//double sigmaxx_avg = sigmaxx_sum/nlocal;
//double sigmayy_avg = sigmayy_sum/nlocal;
//double sigmazz_avg = sigmazz_sum/nlocal;
comm->forward_comm(this);
}
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*MY_PI*pow(R,2.0) << std::endl;
//std::cout << "Fix radius update setup has been entered !!!" << std::endl;
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;