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

result of "make fix-whitespace"

Browse Source
This commit is contained in:
Sachith Dunatunga
2024-12-14 13:20:16 -08:00
parent 00ebe9a3e8
commit e48f288e2b
12 changed files with 198 additions and 198 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/GRANULAR/fix_wall_gran.h
View File

@ -87,7 +87,7 @@ class FixWallGran : public Fix {
int store;
};
} // namespace LAMMPS_NS

2
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -231,7 +231,7 @@ 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->i = i; // Added for MDR
model->j = ic; // Added for MDR
if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0];

10
src/GRANULAR/fix_wall_gran_region.h
View File

@ -53,20 +53,20 @@ class FixWallGranRegion : public FixWallGran {
int tmax; // MOVED FROM PRIVATE TO PUBLIC FOR MDR MODEL
private:
int nregion;
// 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
// vel info is to be reset
};
} // namespace LAMMPS_NS

204
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -22,7 +22,7 @@
#include <cmath>
#include <iostream>
#include <iomanip>
#include <iomanip>
#include <sstream>
using namespace LAMMPS_NS;
@ -401,12 +401,12 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
{
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;
size_history = 26;
nondefault_history_transfer = 1;
transfer_history_factor = new double[size_history];
//transfer_history_factor[0] = +1;
for (int i = 0; i < size_history; i++) {
for (int i = 0; i < size_history; i++) {
transfer_history_factor[i] = +1;
}
}
@ -435,18 +435,18 @@ void GranSubModNormalMDR::coeffs_to_local()
double GranSubModNormalMDR::calculate_forces()
{
// To understand the structure of the overall code it is important to consider
// To understand the structure of the overall code it is important to consider
// the following:
//
// 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
// 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
// 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
// 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
@ -456,8 +456,8 @@ double GranSubModNormalMDR::calculate_forces()
const int j_true = gm->j; // true j particle index
const double radi_true = gm->radi; // true i particle initial radius
const double radj_true = gm->radj; // true j particle initial radius
F = 0.0; // average force
F = 0.0; // average force
double F0 = 0.0; // force on contact side 0
double F1 = 0.0; // force on contact side 1
double R0 = 0.0;
@ -465,21 +465,21 @@ double GranSubModNormalMDR::calculate_forces()
int i0 = 0;
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;
const int deltao_offset_0 = 2; // displacement
const int deltao_offset_1 = 3;
//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;
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_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;
@ -493,32 +493,32 @@ double GranSubModNormalMDR::calculate_forces()
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 penalty_offset_ = 22; // contact penalty
const int deltamax_offset_ = 23;
const int deltap_offset_0 = 24;
const int deltap_offset_1 = 25;
// initialize particle history variables
// initialize particle history variables
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
int index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
int index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
int index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
int index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
int index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
int index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
int index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
double * Rinitial = atom->dvector[index_Ro];
double * Vgeo = atom->dvector[index_Vgeo];
double * Velas = atom->dvector[index_Velas];
@ -527,7 +527,7 @@ double GranSubModNormalMDR::calculate_forces()
double * dRnumerator = atom->dvector[index_dRnumerator];
double * dRdenominator = atom->dvector[index_dRdenominator];
double * Acon0 = atom->dvector[index_Acon0];
double * Acon1 = atom->dvector[index_Acon1];
double * Acon1 = atom->dvector[index_Acon1];
double * Atot = atom->dvector[index_Atot];
double * Atot_sum = atom->dvector[index_Atot_sum];
double * ddelta_bar = atom->dvector[index_ddelta_bar];
@ -537,9 +537,9 @@ double GranSubModNormalMDR::calculate_forces()
double * sigmazz = atom->dvector[index_sigmazz];
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
double * history = & gm->history[history_index]; // load in all history variables
// Rigid flat placement scheme
double * deltamax_offset = & history[deltamax_offset_];
@ -552,19 +552,19 @@ double GranSubModNormalMDR::calculate_forces()
if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta;
deltamax = *deltamax_offset;
for (int contactSide = 0; contactSide < 2; contactSide++) {
for (int contactSide = 0; contactSide < 2; contactSide++) {
double * delta_offset;
double * delta_offset;
double * deltao_offset;
double * delta_MDR_offset;
double * delta_BULK_offset;
double * deltamax_MDR_offset;
double * Yflag_offset;
double * deltaY_offset;
double * delta_MDR_offset;
double * delta_BULK_offset;
double * deltamax_MDR_offset;
double * Yflag_offset;
double * deltaY_offset;
double * cA_offset;
double * aAdh_offset;
double * Ac_offset;
double * eps_bar_offset;
double * aAdh_offset;
double * Ac_offset;
double * eps_bar_offset;
double * penalty_offset;
double * deltap_offset;
@ -631,7 +631,7 @@ double GranSubModNormalMDR::calculate_forces()
gm->radi = radj_true;
gm->radj = radi_true;
}
double delta_geo;
double delta_geo_alt;
double delta_geoOpt1 = deltamax*(deltamax - 2.0*R1)/(2.0*(deltamax - R0 - R1));
@ -649,7 +649,7 @@ double GranSubModNormalMDR::calculate_forces()
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];
@ -678,8 +678,8 @@ double GranSubModNormalMDR::calculate_forces()
const double Ro = Rinitial[i]; // initial radius
const double R = gm->radi; // apparent radius
// kinematics
// kinematics
const double ddelta = delta - *delta_offset;
*delta_offset = delta;
@ -689,7 +689,7 @@ double GranSubModNormalMDR::calculate_forces()
double ddelta_MDR;
double ddelta_BULK;
if ( psi[i] < psi_b ) { // if true, bulk response has triggered, split displacement increment between the MDR and BULK components
if ( psi[i] < psi_b ) { // if true, bulk response has triggered, split displacement increment between the MDR and BULK components
ddelta_MDR = std::min(ddelta-ddelta_bar[i], delta-*delta_MDR_offset);
ddelta_BULK = ddelta_bar[i];
} else { // if false, no bulk response, full displacement increment goes to the MDR component
@ -712,7 +712,7 @@ double GranSubModNormalMDR::calculate_forces()
*deltaY_offset = delta_MDR;
*cA_offset = M_PI*(pow(*deltaY_offset,2.0) - *deltaY_offset*R);
}
}
}
//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;
@ -725,25 +725,25 @@ double GranSubModNormalMDR::calculate_forces()
double A; // height of elliptical indenter
double B; // width of elliptical indenter
double deltae1D; // transformed elastic displacement
double deltaR; // displacement correction
double deltaR; // displacement correction
double amax; // maximum experienced contact radius
const double cA = *cA_offset; // contact area intercept
if ( *Yflag_offset == 0.0 ) { // elastic contact
A = 4.0*R;
B = 2.0*R;
deltae1D = delta_MDR;
(deltae1D > 0) ? amax = sqrt(deltae1D*R) : amax = 0.0;
A = 4.0*R;
B = 2.0*R;
deltae1D = delta_MDR;
(deltae1D > 0) ? amax = sqrt(deltae1D*R) : amax = 0.0;
} else { // plastic contact
amax = sqrt((2.0*deltamax_MDR*R - pow(deltamax_MDR,2.0)) + cA/M_PI);
A = 4.0*pY/Eeff*amax;
B = 2.0*amax;
const double deltae1Dmax = A/2.0; // maximum transformed elastic displacement
amax = sqrt((2.0*deltamax_MDR*R - pow(deltamax_MDR,2.0)) + cA/M_PI);
A = 4.0*pY/Eeff*amax;
B = 2.0*amax;
const double deltae1Dmax = A/2.0; // maximum transformed elastic displacement
const double Fmax = Eeff*(A*B/4.0)*(acos(1.0 - 2.0*deltae1Dmax/A) - (1.0 - 2.0*deltae1Dmax/A)*sqrt(4.0*deltae1Dmax/A - 4.0*pow(deltae1Dmax,2.0)/pow(A,2.0))); // force caused by full submersion of elliptical indenter to depth of A/2
const double zR = R - (deltamax_MDR - deltae1Dmax); // depth of particle center
deltaR = (Fmax*(2*pow(amax,2.0)*(-1 + nu) - (-1 + 2*nu)*zR*(-zR + sqrt(pow(amax,2.0) + pow(zR,2.0)))))/((M_PI*pow(amax,2.0))*2*G*sqrt(pow(amax,2.0) + pow(zR,2.0)));
deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR)/(1 + deltaR/deltae1Dmax); // transformed elastic displacement
deltaR = (Fmax*(2*pow(amax,2.0)*(-1 + nu) - (-1 + 2*nu)*zR*(-zR + sqrt(pow(amax,2.0) + pow(zR,2.0)))))/((M_PI*pow(amax,2.0))*2*G*sqrt(pow(amax,2.0) + pow(zR,2.0)));
deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR)/(1 + deltaR/deltae1Dmax); // transformed elastic displacement
// added for rigid flat placement
*deltap_offset = deltamax_MDR - (deltae1Dmax + deltaR);
//std::cout << *deltap_offset << std::endl;
@ -755,7 +755,7 @@ double GranSubModNormalMDR::calculate_forces()
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;
double aAdh = *aAdh_offset;
if (aAdh > a_fac*amax) aAdh = a_fac*amax;
//if (i_true == 4 && j_true == 52){
@ -764,9 +764,9 @@ double GranSubModNormalMDR::calculate_forces()
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)));
(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;
@ -775,7 +775,7 @@ double GranSubModNormalMDR::calculate_forces()
}
*aAdh_offset = a_fac*a_na;
} else {
const double lmax = sqrt(2.0*M_PI*aAdh*gamma/Eeff);
const double lmax = sqrt(2.0*M_PI*aAdh*gamma/Eeff);
g_aAdh = A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh,2.0));
const double acrit = (-((pow(B,2)*gamma*M_PI)/(pow(A,2)*Eeff)) + (pow(2,0.3333333333333333)*pow(B,4)*pow(gamma,2)*pow(M_PI,1.6666666666666667))/
(pow(A,2)*pow(Eeff,2)*pow((27*pow(A,4)*pow(B,4)*gamma)/Eeff - (2*pow(B,6)*pow(gamma,3)*pow(M_PI,2))/pow(Eeff,3) + (3*sqrt(3)*sqrt(27*pow(A,8)*pow(B,8)*pow(Eeff,2)*pow(gamma,2) -
@ -785,19 +785,19 @@ double GranSubModNormalMDR::calculate_forces()
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;
//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 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;
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;
} 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;
//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;
@ -807,31 +807,31 @@ double GranSubModNormalMDR::calculate_forces()
const double error = 1e-10;
const double error2 = 1e-16;
double aAdh_tmp = aAdh;
double fa;
double fa;
double fa2;
double dfda;
for (int lv1 = 0; lv1 < maxIterations; ++lv1) {
fa = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) );
if (abs(fa) < error) {
break;
}
}
dfda = -((aAdh_tmp*A)/(B*sqrt(-pow(aAdh_tmp,2.0) + pow(B,2.0)/4.0))) + (gamma*sqrt(M_PI/2.0))/(Eeff*sqrt((aAdh_tmp*gamma)/Eeff));
aAdh_tmp = aAdh_tmp - fa/dfda;
fa2 = deltae1D + sqrt(2.0*M_PI*aAdh_tmp*gamma/Eeff) - ( A/2 - A/B*sqrt(pow(B,2.0)/4 - pow(aAdh_tmp,2.0)) );
if (abs(fa-fa2) < error2) {
break;
}
}
if (lv1 == maxIterations-1){
aAdh_tmp = 0.0;
}
}
aAdh = aAdh_tmp;
g_aAdh = A/2.0 - A/B*sqrt(pow(B,2.0)/4.0 - pow(aAdh,2.0));
const double deltaeAdh = g_aAdh;
g_aAdh = A/2.0 - A/B*sqrt(pow(B,2.0)/4.0 - pow(aAdh,2.0));
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;
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;
}
*aAdh_offset = aAdh;
@ -839,13 +839,13 @@ double GranSubModNormalMDR::calculate_forces()
}
} else { // non-adhesive contact
*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)));
(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);
}
}
}
//std::cout << gm->i << ", " << gm->j << ", aAdh_offset: " << *aAdh_offset << ", aAdh: " << aAdh << ", a_na: " << a_na << std::endl;
@ -860,8 +860,8 @@ double GranSubModNormalMDR::calculate_forces()
//std::cout << gm->i << ", " << gm->j << ", " << gm->contact_type << ", " << *gm->xi[1] << ", " << *gm->xj[1] << std::endl;
// area related calculations
double Ac;
// area related calculations
double Ac;
(*Yflag_offset == 0.0) ? Ac = M_PI*delta*R : Ac = M_PI*((2.0*delta*R - pow(delta,2.0)) + cA/M_PI);
if (Ac < 0.0 ) Ac = 0.0;
Atot_sum[i] += wij*(Ac - 2.0*M_PI*R*(deltamax_MDR + delta_BULK));
@ -932,7 +932,7 @@ double GranSubModNormalMDR::calculate_forces()
//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 << "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;
@ -950,7 +950,7 @@ double GranSubModNormalMDR::calculate_forces()
// radius update scheme quantity calculation
Vcaps[i] += (M_PI/3.0)*pow(delta,2.0)*(3.0*R - delta);
const double Fntmp = wij*(F_MDR + F_BULK);
const double fx = Fntmp*gm->nx[0];
const double fy = Fntmp*gm->nx[1];
@ -960,7 +960,7 @@ double GranSubModNormalMDR::calculate_forces()
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;
//double **x = atom->x;
//const double xi = x[i_true][0];
//const double xj = x[j_true][0];
@ -977,7 +977,7 @@ double GranSubModNormalMDR::calculate_forces()
//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
if(delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0){
const double Vo = (4.0/3.0)*M_PI*pow(Ro,3.0);
@ -1020,7 +1020,7 @@ double GranSubModNormalMDR::calculate_forces()
//}
}
gm->i = i_true;
gm->j = j_true;
gm->radi = radi_true;
@ -1044,7 +1044,7 @@ double GranSubModNormalMDR::calculate_forces()
if (gm->contact_type != PAIR) {
F = wij*F0*wallForceMagnifer;
} else {
F = wij*(F0 + F1)/2.0;
F = wij*(F0 + F1)/2.0;
}
//std::cout << F << ", " << F0 << ", " << F1 << " | " << R0 << ", " << R1 << std::endl;
@ -1073,7 +1073,7 @@ double GranSubModNormalMDR::calculate_forces()
//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;
@ -1133,21 +1133,21 @@ void GranSubModNormalMDR::set_fncrit()
//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 << 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 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 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;
//double wall_contact_flag_offset;

18
src/GRANULAR/pair_granular.cpp
View File

@ -504,7 +504,7 @@ void PairGranular::init_style()
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") {
@ -521,22 +521,22 @@ void PairGranular::init_style()
index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THE INPUT SCRIPT
index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); // flag to check if history variables have been initialized
index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); // flag to check if history variables have been initialized
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_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
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;
@ -548,7 +548,7 @@ void PairGranular::init_style()
modify->add_fix("fix_mdr_mean_surf_disp all mdr/mean/surf/disp");
fix_flag = 1;
}
}
}
// check for FixFreeze and set freeze_group_bit

2
src/GRANULAR/pair_granular.h
View File

@ -74,7 +74,7 @@ class PairGranular : public Pair {
int index_Ro;
int index_Vcaps;
int index_Vgeo;
int index_Velas;
int index_Velas;
int index_eps_bar;
int index_dRnumerator;
int index_dRdenominator;

2
src/csv_writer.h
View File

@ -22,4 +22,4 @@ public:
private:
std::string filename_;
};
};

46
src/fix_mdr_mean_surf_disp.cpp
View File

@ -63,9 +63,9 @@ int FixMDRmeanSurfDisp::setmask()
void FixMDRmeanSurfDisp::setup_pre_force(int /*vflag*/)
{
int tmp1, tmp2;
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
Acon0 = atom->dvector[index_Acon0];
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
Acon0 = atom->dvector[index_Acon0];
ddelta_bar = atom->dvector[index_ddelta_bar];
pre_force(0);
@ -97,9 +97,9 @@ void FixMDRmeanSurfDisp::pre_force(int)
FixNeighHistory * fix_history = dynamic_cast<FixNeighHistory *>(modify->get_fix_by_id("NEIGH_HISTORY_GRANULAR"));
PairGranular * pair = dynamic_cast<PairGranular *>(force->pair_match("granular",1));
NeighList * list = pair->list;
const int size_history = pair->get_size_history();
//std::cout << " " << std::endl;
//std::cout << "New Step" << std::endl;
@ -135,14 +135,14 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double ytmp = x[i][1];
const double ztmp = x[i][2];
allhistory = firsthistory[i];
double radi = radius[i];
double radi = radius[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
double radj = radius[j];
double radj = radius[j];
const double delx_ij = x[j][0] - xtmp;
const double dely_ij = x[j][1] - ytmp;
const double delz_ij = x[j][2] - ztmp;
@ -153,7 +153,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double deltan_ij = radsum_ij - r_ij;
if (deltan_ij >= 0.0) {
for (int kk = jj; kk < jnum; kk++) {
k = jlist[kk];
k = jlist[kk];
k &= NEIGHMASK;
ktype = type[k];
if (kk != jj) {
@ -175,7 +175,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double radsum_jk = radj + radk;
const double deltan_jk = radsum_jk - r_jk;
if (deltan_ik >= 0.0 && deltan_jk >= 0.0) {
// pull ij history
history_ij = &allhistory[size_history * jj];
double * pij = &history_ij[22]; // penalty for contact i and j
@ -186,7 +186,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
// 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;
double * pjk = NULL;
int * const jklist = firstneigh[j];
const int jknum = numneigh[j];
for (int jk = 0; jk < jknum; jk++) {
@ -219,7 +219,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
//std::cout << "Print 198 pjk[0]: " << rank << ", " << pjk[0] << std::endl;
break;
}
}
}
}
//std::cout << "Print: " << __LINE__ << std::endl;
@ -234,8 +234,8 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double enx_kj = -delx_jk * rinv_jk;
const double eny_kj = -dely_jk * rinv_jk;
const double enz_kj = -delz_jk * rinv_jk;
const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj);
pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj);
pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
} else if (maxIndex == 1) { // the central particle is j
const double enx_ji = -delx_ij * rinv_ij;
const double eny_ji = -dely_ij * rinv_ij;
@ -243,7 +243,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double enx_jk = delx_jk * rinv_jk;
const double eny_jk = dely_jk * rinv_jk;
const double enz_jk = delz_jk * rinv_jk;
const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk);
const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk);
pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/M_PI - 1.0/2.0)) );
} else { // the central particle is i
if (j < atom->nlocal || k < atom->nlocal) {
@ -253,7 +253,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
const double enx_ik = delx_ik * rinv_ik;
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);
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/M_PI - 1.0/2.0)) );
//std::cout << "Print: " << __LINE__ << std::endl;
@ -264,7 +264,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
}
}
}
}
}
}
@ -286,7 +286,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
double *radius = atom->radius;
int nlocal = atom->nlocal;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -336,8 +336,8 @@ void FixMDRmeanSurfDisp::pre_force(int)
//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
const int delta_offset_1 = 1;
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;
@ -347,7 +347,7 @@ void FixMDRmeanSurfDisp::pre_force(int)
double deltamax = history[deltamax_offset_];
double deltap0 = history[deltap_offset_0];
double deltap1 = history[deltap_offset_1];
if (delta > deltamax) deltamax = delta;
double delta0old = history[delta_offset_0];
@ -511,13 +511,13 @@ void FixMDRmeanSurfDisp::pre_force(int)
if (Acon0[i] != 0.0) {
const double delta = model->radsum - model->r;
const double delta_offset0 = fix->history_many[i][fix->c2r[ic]][0];
const double ddelta = delta - delta_offset0;
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;
}
}
}
}
}
}
@ -536,4 +536,4 @@ void FixMDRmeanSurfDisp::pre_force(int)
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_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;
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;

6
src/fix_mdr_mean_surf_disp.h
View File

@ -25,10 +25,10 @@ FixStyle(mdr/mean/surf/disp,FixMDRmeanSurfDisp);
namespace LAMMPS_NS {
class FixMDRmeanSurfDisp : public Fix {
public:
double * Acon0;
public:
double * Acon0;
double * ddelta_bar;
FixMDRmeanSurfDisp(class LAMMPS *, int, char **);
int setmask() override;
void setup_pre_force(int) override;

96
src/fix_mdr_radius_update.cpp
View File

@ -35,7 +35,7 @@
#include "gran_sub_mod_normal.h"
#include "update.h"
#include "comm.h"
#include <iomanip>
#include <iomanip>
#include <sstream>
using namespace LAMMPS_NS;
@ -69,19 +69,19 @@ void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/)
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_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_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);
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];
@ -89,9 +89,9 @@ void FixMDRradiusUpdate::setup_pre_force(int /*vflag*/)
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
@ -116,19 +116,19 @@ void FixMDRradiusUpdate::setup(int /*vflag*/)
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_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_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);
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];
@ -136,9 +136,9 @@ void FixMDRradiusUpdate::setup(int /*vflag*/)
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
@ -170,20 +170,20 @@ void FixMDRradiusUpdate::pre_force(int)
//std::cout << "Preforce was called radius update" << std::endl;
// assign correct value to initially non-zero MDR particle history variables
// 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_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);
//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];
@ -196,13 +196,13 @@ void FixMDRradiusUpdate::pre_force(int)
//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++) {
for (int i = 0; i < nlocal; i++) {
if (history_setup_flag[i] < 1e-16) {
Ro[i] = radius[i];
Vgeo[i] = 4.0/3.0*M_PI*pow(Ro[i],3.0);
@ -232,7 +232,7 @@ int FixMDRradiusUpdate::pack_forward_comm(int n, int *list, double *buf, int /*p
buf[m++] = Vgeo[j]; // 2
buf[m++] = Velas[j]; // 3
buf[m++] = Vcaps[j]; // 4
buf[m++] = eps_bar[j]; // 5
buf[m++] = eps_bar[j]; // 5
buf[m++] = dRnumerator[j]; // 6
buf[m++] = dRdenominator[j]; // 7
buf[m++] = Acon0[j]; // 8
@ -261,7 +261,7 @@ void FixMDRradiusUpdate::unpack_forward_comm(int n, int first, double *buf)
Vgeo[i] = buf[m++]; // 2
Velas[i] = buf[m++]; // 3
Vcaps[i] = buf[m++]; // 4
eps_bar[i] = buf[m++]; // 5
eps_bar[i] = buf[m++]; // 5
dRnumerator[i] = buf[m++]; // 6
dRdenominator[i] = buf[m++]; // 7
Acon0[i] = buf[m++]; // 8
@ -282,7 +282,7 @@ 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;
@ -293,13 +293,13 @@ void FixMDRradiusUpdate::end_of_step()
int nlocal = atom->nlocal;
double sigmaxx_sum = 0.0;
double sigmayy_sum = 0.0;
double sigmazz_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++) {
const double R = radius[i];
Atot[i] = 4.0*M_PI*pow(R,2.0) + Atot_sum[i];
@ -316,7 +316,7 @@ void FixMDRradiusUpdate::end_of_step()
//}
const double dR = std::max(dRnumerator[i]/(dRdenominator[i] - 4.0*M_PI*pow(R,2.0)),0.0);
if (psi_b[i] < psi[i]) {
if (psi_b[i] < psi[i]) {
if ((radius[i] + dR) < (1.5*Ro[i])) radius[i] += dR;
//radius[i] += dR;
}
@ -327,10 +327,10 @@ void FixMDRradiusUpdate::end_of_step()
//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*M_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;
@ -339,7 +339,7 @@ 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;
//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
@ -381,4 +381,4 @@ void FixMDRradiusUpdate::end_of_step()
//std::cout << radius[i] << ", " << dR << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRdenominator[i] - 4.0*M_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;
//std::cout << Ro[0] << ", " << Vgeo[0] << ", " << Velas[0] << std::endl;

6
src/fix_mdr_radius_update.h
View File

@ -31,7 +31,7 @@ class FixMDRradiusUpdate : public Fix {
double * Velas;
double * Vcaps;
double * eps_bar;
double * dRnumerator;
double * dRnumerator;
double * dRdenominator;
double * Acon0;
double * Acon1;
@ -43,10 +43,10 @@ class FixMDRradiusUpdate : public Fix {
double * sigmaxx;
double * sigmayy;
double * sigmazz;
double * history_setup_flag;
double * history_setup_flag;
double * contacts;
double * adhesive_length;
FixMDRradiusUpdate(class LAMMPS *, int, char **);
int setmask() override;
void setup(int) override;

2
src/variable.cpp
View File

@ -3661,7 +3661,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree
strcmp(word,"logfreq") != 0 && strcmp(word,"logfreq2") != 0 &&
strcmp(word,"logfreq3") != 0 && strcmp(word,"stride") != 0 &&
strcmp(word,"stride2") != 0 && strcmp(word,"vdisplace") != 0 &&
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0)
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0)
return 0;
// parse contents for comma-separated args