diff --git a/examples/granular/in.mpfem.triaxial12particles b/examples/granular/in.mpfem.triaxial12particles index 8c417086bc..98b5621c1f 100644 --- a/examples/granular/in.mpfem.triaxial12particles +++ b/examples/granular/in.mpfem.triaxial12particles @@ -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 ############################### diff --git a/examples/granular/spheres.data b/examples/granular/spheres12.data similarity index 100% rename from examples/granular/spheres.data rename to examples/granular/spheres12.data diff --git a/examples/granular/spheresSTLgeneration.data b/examples/granular/spheresSTLgeneration.data deleted file mode 100644 index 960c7d357b..0000000000 --- a/examples/granular/spheresSTLgeneration.data +++ /dev/null @@ -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 diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 41db281568..3336a8a4d7 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -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(); diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 16b1b228ab..e62ab4e38e 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -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; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index ba8cd467ec..cc1ac3e04b 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -29,7 +29,6 @@ #include "region.h" #include "update.h" #include "variable.h" -#include 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]; diff --git a/src/GRANULAR/fix_wall_gran_region.h b/src/GRANULAR/fix_wall_gran_region.h index 6f9866d8ae..50cd3236b5 100644 --- a/src/GRANULAR/fix_wall_gran_region.h +++ b/src/GRANULAR/fix_wall_gran_region.h @@ -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 diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index 169756b41f..36bbd18982 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -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 -#include #include #include @@ -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 @@ -457,19 +470,15 @@ double GranSubModNormalMDR::calculate_forces() 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 - double F0 = 0.0; // force on contact side 0 - double F1 = 0.0; // force on contact side 1 + 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; double R1 = 0.0; 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; @@ -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; diff --git a/src/GRANULAR/gran_sub_mod_rolling.cpp b/src/GRANULAR/gran_sub_mod_rolling.cpp index b236b9a142..109bbee65d 100644 --- a/src/GRANULAR/gran_sub_mod_rolling.cpp +++ b/src/GRANULAR/gran_sub_mod_rolling.cpp @@ -17,7 +17,6 @@ #include "gran_sub_mod_normal.h" #include "granular_model.h" #include "math_extra.h" -#include #include @@ -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); diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index d4fb6f66a8..546c7a3d4e 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -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 diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 230cba446b..aee341aef1 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -39,7 +39,6 @@ #include "gran_sub_mod_normal.h" #include -#include 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"); diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 10b581968a..199431afe3 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -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; diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp index 054ef83349..9b153a1c28 100644 --- a/src/SRD/fix_srd.cpp +++ b/src/SRD/fix_srd.cpp @@ -40,7 +40,6 @@ #include #include -#include 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); } diff --git a/src/csv_writer.h b/src/csv_writer.h deleted file mode 100644 index 29bb44e419..0000000000 --- a/src/csv_writer.h +++ /dev/null @@ -1,25 +0,0 @@ -#include -#include -#include - -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_; -}; diff --git a/src/fix_mdr_mean_surf_disp.cpp b/src/fix_mdr_mean_surf_disp.cpp index 91804e51fb..324f9fbd1b 100644 --- a/src/fix_mdr_mean_surf_disp.cpp +++ b/src/fix_mdr_mean_surf_disp.cpp @@ -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 #include 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 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; diff --git a/src/fix_mdr_radius_update.cpp b/src/fix_mdr_radius_update.cpp index 36ada70f55..7ff883a0cf 100644 --- a/src/fix_mdr_radius_update.cpp +++ b/src/fix_mdr_radius_update.cpp @@ -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 #include #include @@ -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(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(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;