From 9d01ac2caf1c0946399c438bfd717f2d653b647c Mon Sep 17 00:00:00 2001 From: William Zunker Date: Wed, 26 Mar 2025 10:43:51 -0400 Subject: [PATCH] Stablized radius update, added WFM, added outputs --- src/GRANULAR/fix_granular_mdr.cpp | 42 +++++- src/GRANULAR/fix_granular_mdr.h | 3 +- src/GRANULAR/gran_sub_mod_damping.cpp | 10 +- src/GRANULAR/gran_sub_mod_normal.cpp | 189 ++++++++++++++++++++++---- src/GRANULAR/gran_sub_mod_normal.h | 3 +- src/GRANULAR/granular_model.cpp | 4 +- src/csv_writer.h | 25 ++++ 7 files changed, 236 insertions(+), 40 deletions(-) create mode 100644 src/csv_writer.h diff --git a/src/GRANULAR/fix_granular_mdr.cpp b/src/GRANULAR/fix_granular_mdr.cpp index 9efabdf465..712d4ad1ff 100644 --- a/src/GRANULAR/fix_granular_mdr.cpp +++ b/src/GRANULAR/fix_granular_mdr.cpp @@ -38,6 +38,10 @@ #include "update.h" #include "variable.h" +// wzunker +#include "csv_writer.h" +#include + using namespace LAMMPS_NS; using namespace Granular_NS; using namespace Granular_MDR_NS; @@ -45,7 +49,7 @@ using namespace FixConst; using MathConst::MY_PI; static constexpr double EPSILON = 1e-16; -static constexpr double OVERLAP_LIMIT = 0.75; +static constexpr double OVERLAP_LIMIT = 0.95; enum { COMM_1, COMM_2 }; @@ -85,7 +89,7 @@ void FixGranularMDR::post_constructor() modify->add_fix( fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator " "d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi " - "d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz ghost yes", + "d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_dRavg ghost yes", id_fix)); index_Ro = atom->find_custom("Ro", tmp1, tmp2); @@ -105,6 +109,7 @@ void FixGranularMDR::post_constructor() index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); + index_dRavg = atom->find_custom("dRavg", tmp1, tmp2); } /* ---------------------------------------------------------------------- */ @@ -208,6 +213,7 @@ void FixGranularMDR::pre_force(int) double *sigmayy = atom->dvector[index_sigmayy]; double *sigmazz = atom->dvector[index_sigmazz]; double *history_setup_flag = atom->dvector[index_history_setup_flag]; + double *dRavg = atom->dvector[index_dRavg]; int new_atom; int nlocal = atom->nlocal; @@ -250,10 +256,35 @@ void FixGranularMDR::pre_force(int) psi[i] = (Atot[i] - Acon1[i]) / Atot[i]; if (psi_b_coeff < psi[i]) { - const double dR = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)), 0.0); - if ((radius[i] + dR) < (1.5 * Ro[i])) radius[i] += dR; + double w_confinement; + ( psi[i] > 0.1 ) ? w_confinement = 1.0/(1.0 + exp(-75.0*(psi[i]-0.2))) : w_confinement = 0.0; + const double dR = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0))*w_confinement, 0.0); + + const double N_window = 10.0; + if (dR > 0.0) dRavg[i] += (dR - dRavg[i]) / N_window; + + if (((radius[i] + dR) < (1.5 * Ro[i])) && (dR > 0.0)) radius[i] += dRavg[i]; + + //(dR + dRavg[i])/2.0; + //dRavg[i] = (dR + dRavg[i])/2.0; + + //// wzunker + //const double dRdenominatorTrue = (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)); + //if (i == 0) { + // CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters.csv"); + // std::stringstream rowDataStream; + // rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision + // rowDataStream << dRnumerator[i] << ", " << dRdenominator[i] << ", " << dRavg[i] << ", " << R << ", " << dRdenominatorTrue << ", " << lmp->update->ntimestep; + // std::string rowData = rowDataStream.str(); + // csvWriter.writeRow(rowData); + //} + } Acon0[i] = Acon1[i]; + + // wzunker + //const double dR_tmp = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)), 0.0); + //printf("i = %d, psi_b = %f, psi = %f, Atot = %f, Acon = %f, dRnumerator = %f, dRdenominator = %f, R = %f, dR = %f \n", i, psi_b_coeff, psi[i], Atot[i], Acon1[i], dRnumerator[i], dRdenominator[i],radius[i],dR_tmp); } comm_stage = COMM_1; @@ -510,6 +541,9 @@ void FixGranularMDR::calculate_contact_penalty() "(e.g. increase the skin distance)."); pjk[0] += 1.0 / (1.0 + std::exp(-50.0 * (alpha / MY_PI - 0.5))); + + // wzunker + //printf("i = %d, j = %d, k = %d, pij = %f, pik = %f, pjk = %f \n", i, j, k, pij[0], pik[0], pjk[0]); } } } diff --git a/src/GRANULAR/fix_granular_mdr.h b/src/GRANULAR/fix_granular_mdr.h index f0ba76d155..56c576d79b 100644 --- a/src/GRANULAR/fix_granular_mdr.h +++ b/src/GRANULAR/fix_granular_mdr.h @@ -97,8 +97,7 @@ class FixGranularMDR : public Fix { int index_sigmayy; // yy-component of the stress tensor, not necessary forforce calculation int index_sigmazz; // zz-component of the stress tensor, not necessary forforce calculation int index_history_setup_flag; // flag to check if history variables have beeninitialized - int index_contacts; // total contacts on particle - int index_adhesive_length; // total length of adhesive contact on a particle + int index_dRavg; // total contacts on particle }; } // namespace LAMMPS_NS diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp index 2fdd7e1f82..4072944d5e 100644 --- a/src/GRANULAR/gran_sub_mod_damping.cpp +++ b/src/GRANULAR/gran_sub_mod_damping.cpp @@ -18,6 +18,8 @@ #include "math_special.h" #include "math_const.h" +#include "style_gran_sub_mod.h" // IWYU pragma: keep + #include using namespace LAMMPS_NS; @@ -74,8 +76,12 @@ GranSubModDampingVelocity::GranSubModDampingVelocity(GranularModel *gm, LAMMPS * double GranSubModDampingVelocity::calculate_forces() { - damp_prefactor = damp; - return -damp_prefactor * gm->vnnr; + if (gm->normal_model->name == "mdr") { + return 0.0; + } else { + damp_prefactor = damp; + return -damp_prefactor * gm->vnnr; + } } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index 226d26d86b..6511f75c30 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -26,6 +26,9 @@ #include #include +// wzunker +#include "csv_writer.h" + using namespace LAMMPS_NS; using namespace Granular_NS; using namespace MathConst; @@ -47,7 +50,7 @@ static constexpr int MDR_MAX_IT = 100; // Newton-Raphs static constexpr double MDR_EPSILON1 = 1e-10; // Newton-Raphson for MDR static constexpr double MDR_EPSILON2 = 1e-16; // Newton-Raphson for MDR static constexpr double MDR_EPSILON3 = 1e-20; // For precision checks -static constexpr double MDR_OVERLAP_LIMIT = 0.75; // Maximum contact overlap for MDR +static constexpr double MDR_OVERLAP_LIMIT = 0.95; // Maximum contact overlap for MDR static const char cite_mdr[] = "MDR contact model command: (i) https://doi.org/10.1016/j.jmps.2023.105492 || (ii) https://doi.org/10.1016/j.jmps.2023.105493 || (iii) https://doi.org/10.31224/4289\n\n" @@ -472,14 +475,16 @@ void GranSubModNormalMDR::coeffs_to_local() Y = coeffs[2]; // yield stress gamma = coeffs[3]; // effective surface energy psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total} - CoR = coeffs[5]; // coefficent of restitution + //CoR = coeffs[5]; // coefficent of restitution + damp = coeffs[5]; // coefficent of restitution if (E <= 0.0) error->all(FLERR, "Illegal MDR normal model, Young's modulus must be greater than 0"); if (nu < 0.0 || nu > 0.5) error->all(FLERR, "Illegal MDR normal model, Poisson's ratio must be between 0 and 0.5"); if (Y < 0.0) error->all(FLERR, "Illegal MDR normal model, yield stress must be greater than or equal to 0"); if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0"); if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0"); - if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0"); + //if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0"); + if (damp < 0.0) error->all(FLERR, "Illegal MDR normal model, damping coefficent must be greater than 1"); G = E / (2.0 * (1.0 + nu)); // shear modulus kappa = E / (3.0 * (1.0 - 2.0 * nu)); // bulk modulus @@ -490,6 +495,7 @@ void GranSubModNormalMDR::coeffs_to_local() Eeffinv = 1.0 / Eeff; Eeffsq = Eeff * Eeff; Eeffsqinv = Eeffinv * Eeffinv; + Eeff2particle = 0.5 * Eeff; gammasq = gamma * gamma; gamma3 = gammasq * gamma; @@ -528,6 +534,7 @@ void GranSubModNormalMDR::init() index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); // xx-component of the stress tensor, not necessary for force calculation index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); // yy-component of the stress tensor, not necessary for force calculation index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); // zz-component of the stress tensor, not necessary for force calculation + index_dRavg = atom->find_custom("dRavg", tmp1, tmp2); // radius update increment } /* ---------------------------------------------------------------------- */ @@ -566,6 +573,7 @@ double GranSubModNormalMDR::calculate_forces() double *sigmaxx = atom->dvector[index_sigmaxx]; double *sigmayy = atom->dvector[index_sigmayy]; double *sigmazz = atom->dvector[index_sigmazz]; + double *dRavg = atom->dvector[index_dRavg]; const int itag_true = atom->tag[gm->i]; // true i particle tag const int jtag_true = atom->tag[gm->j]; // true j particle tag @@ -579,6 +587,7 @@ double GranSubModNormalMDR::calculate_forces() double F1 = 0.0; // force on contact side 1 double delta = gm->delta; // apparent overlap double Ac_avg = 0.0; // average contact area across both sides + double a_damp = 0.0; // damping contact radius double *history = & gm->history[history_index]; // load in all history variables int history_update = gm->history_update; @@ -596,6 +605,11 @@ double GranSubModNormalMDR::calculate_forces() if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta; double deltamax = *deltamax_offset; + //wzunker + double F_MDR0; + double F_BULK0; + double F_MDR1; + double F_BULK1; for (int contactSide = 0; contactSide < 2; contactSide++) { @@ -607,7 +621,7 @@ double GranSubModNormalMDR::calculate_forces() // displacement partitioning only necessary for particle-particle contact // itag and jtag persist after neighbor list builds, use tags to compare to match - // contact history variables consistently across steps for a particle pair. + // contact history variables consistently across steps for a particle pair. if ((contactSide == 0 && itag_true > jtag_true) || (contactSide != 0 && itag_true < jtag_true)) { gm->i = i_true; gm->j = j_true; @@ -804,8 +818,10 @@ double GranSubModNormalMDR::calculate_forces() F_MDR = calculate_nonadhesive_mdr_force(deltae1D, Ainv, Eeff, A, B); } - if (std::isnan(F_MDR)) + if (std::isnan(F_MDR)) { + printf("itag = %d, jtag = %d \n", itag_true, jtag_true); error->one(FLERR, "F_MDR is NaN, case 1: no tensile springs"); + } if (history_update) *aAdh_offset = a_fac * a_na; } else { @@ -861,15 +877,20 @@ double GranSubModNormalMDR::calculate_forces() } aAdh = aAdh_tmp; - g_aAdh = A * 0.5 - A * Binv * sqrt(Bsq * 0.25 - pow(aAdh, 2)); - g_aAdh = round_up_negative_epsilon(g_aAdh); + if (aAdh < acrit) { + aAdh = 0.0; + F_MDR = 0.0; + } else { + g_aAdh = A * 0.5 - A * Binv * sqrt(Bsq * 0.25 - pow(aAdh, 2)); + g_aAdh = round_up_negative_epsilon(g_aAdh); - const double deltaeAdh = g_aAdh; - const double F_na = calculate_nonadhesive_mdr_force(deltaeAdh, Ainv, Eeff, A, B); - const double F_Adhes = 2.0 * Eeff * (deltae1D - deltaeAdh) * aAdh; - F_MDR = F_na + F_Adhes; - if (std::isnan(F_MDR)) - error->one(FLERR, "F_MDR is NaN, case 3: tensile springs exceed critical length"); + const double deltaeAdh = g_aAdh; + const double F_na = calculate_nonadhesive_mdr_force(deltaeAdh, Ainv, Eeff, A, B); + const double F_Adhes = 2.0 * Eeff * (deltae1D - deltaeAdh) * aAdh; + F_MDR = F_na + F_Adhes; + if (std::isnan(F_MDR)) + error->one(FLERR, "F_MDR is NaN, case 3: tensile springs exceed critical length"); + } } if (history_update) *aAdh_offset = aAdh; } @@ -891,13 +912,22 @@ double GranSubModNormalMDR::calculate_forces() } Ac_avg += wij * Ac; + // wzunker damping contact area related things + (gamma > 0.0) ? a_damp += aAdh : a_damp += a_na; + // bulk force calculation double F_BULK; (delta_BULK <= 0.0) ? F_BULK = 0.0 : F_BULK = (1.0 / Vgeo[i]) * Acon0[i] * delta_BULK * kappa * Ac; + // force-magnifier + // total force calculation (contactSide == 0) ? F0 = F_MDR + F_BULK : F1 = F_MDR + F_BULK; + // wzunker + (contactSide == 0) ? F_MDR0 = F_MDR : F_MDR1 = F_MDR; + (contactSide == 0) ? F_BULK0 = F_BULK : F_BULK1 = F_BULK; + if (history_update) { // mean surface displacement calculation *Ac_offset = wij * Ac; @@ -921,6 +951,32 @@ double GranSubModNormalMDR::calculate_forces() dRnumerator[i] -= Vo * (eps_bar_contact - *eps_bar_offset); dRnumerator[i] -= wij * MY_PI * ddeltao * (2 * deltao * Ro - pow(deltao, 2) + pow(R, 2) - pow(Ro, 2)); dRdenominator[i] += wij * 2.0 * MY_PI * R * (deltao + R - Ro); + + + // wzunker + //if (gm->contact_type == PAIR && i == 0) { + // CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters_0.csv"); + // std::stringstream rowDataStream; + // rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision + // rowDataStream << wij << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << eps_bar_contact << ", " << *eps_bar_offset << ", " << ddeltao << ", " << deltao << ", " << R << ", " << Fntmp << ", " << lmp->update->ntimestep; + // std::string rowData = rowDataStream.str(); + // csvWriter.writeRow(rowData); + //} +// + //if (gm->contact_type == PAIR && i == 1) { + // CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/compression_sleeve/dR_parameters_1.csv"); + // std::stringstream rowDataStream; + // rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision + // rowDataStream << wij << ", " << dRnumerator[i] << ", " << dRdenominator[i] << ", " << eps_bar_contact << ", " << *eps_bar_offset << ", " << ddeltao << ", " << deltao << ", " << R << ", " << Fntmp << ", " << lmp->update->ntimestep; + // std::string rowData = rowDataStream.str(); + // csvWriter.writeRow(rowData); + //} + + // wzunker + //if (gm->contact_type == PAIR) { + // printf("i = %d, j = %d, contact_type = %d, dRnumerator = %f, dRdenominator = %f, eps_bar_contact = %f, eps_bar_offset = %f, wij = %f, ddelto = %f, deltao = %f, R = %f, Fntmp = %f \n", i, gm->j, gm->contact_type, dRnumerator[i], dRdenominator[i], eps_bar_contact, *eps_bar_offset, wij, ddeltao, deltao, R, Fntmp); + //} + } if (history_update) { @@ -944,31 +1000,106 @@ double GranSubModNormalMDR::calculate_forces() const double wij = MAX(1.0 - pij, 0.0); // assign final force + double damp_prefactor; if (gm->contact_type != PAIR) { - F = wij * F0; + a_damp = a_damp/2.0; + damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff2particle * a_damp); + double *deltao_offset = &history[DELTAO_0]; + const double wfm = std::exp(10.7*(*deltao_offset)/Rinitial[gm->i] - 10.0) + 1.0; // wall force magnifier + //const double wfm = 1.0; + F = wij * F0 * wfm; } else { + damp_prefactor = 0.5*sqrt(gm->meff * 2.0 * Eeff * a_damp); F = wij * (F0 + F1) * 0.5; } // calculate damping force - if (F > 0.0) { - double Eeff2; - double Reff2; - if (gm->contact_type == PAIR) { - Eeff2 = E / (2.0 * (1.0 - pow(nu, 2))); - Reff2 = 1.0 / ((1.0 / gm->radi + 1.0 / gm->radj)); - } else { - Eeff2 = E / (1.0 - pow(nu, 2)); - Reff2 = gm->radi; - } - const double kn = Eeff2 * Reff2; - const double beta = -log(CoR) / sqrt(pow(log(CoR), 2) + PISQ); - const double damp_prefactor = beta * sqrt(gm->meff * kn); - const double F_DAMP = -damp_prefactor * gm->vnnr; + //if (F > 0.0) { + // double Eeff2; + // double Reff2; + // if (gm->contact_type == PAIR) { + // Eeff2 = E / (2.0 * (1.0 - pow(nu, 2))); + // Reff2 = 1.0 / ((1.0 / gm->radi + 1.0 / gm->radj)); + // } else { + // Eeff2 = E / (1.0 - pow(nu, 2)); + // Reff2 = gm->radi; + // } + // const double kn = Eeff2 * Reff2; + // const double beta = -log(CoR) / sqrt(pow(log(CoR), 2) + PISQ); + // const double damp_prefactor = beta * sqrt(gm->meff * kn); + // const double F_DAMP = -damp_prefactor * gm->vnnr; + // + // F += wij * F_DAMP; + //} - F += wij * F_DAMP; + // calculate damping force + if (F > 0.0) F += -wij * damp * gm->vnnr; + + //F += -wij * damp_prefactor * gm->vnnr; + + //printf("i = %d, j = %d, wij = %f, damp = %f, vnnr = %e, F_DAMP = %e \n", gm->i, gm->j, wij, damp, gm->vnnr, -wij * damp * gm->vnnr); + + // wzunker + double *delta_offset_0 = &history[DELTA_0]; + double *delta_offset_1 = &history[DELTA_0 + 1]; + const double delta0 = *delta_offset_0; + const double delta1 = *delta_offset_1; + double *delta_BULK_offset_0 = &history[DELTA_BULK_0]; + double *delta_BULK_offset_1 = &history[DELTA_BULK_0 + 1]; + const double deltaBULK0 = *delta_BULK_offset_0; + const double deltaBULK1 = *delta_BULK_offset_1; + double *aAdh_offset_0 = &history[AADH_0]; + double *aAdh_offset_1 = &history[AADH_0 + 1]; + const double aAdh0 = *aAdh_offset_0; + const double aAdh1 = *aAdh_offset_1; + double *Ac_offset_0 = &history[AC_0]; + double *Ac_offset_1 = &history[AC_0 + 1]; + const double Ac0 = *Ac_offset_0; + const double Ac1 = *Ac_offset_1; + double F_DAMP = -wij * damp * gm->vnnr; + + const double i_0 = 0; + const double i_1 = 1; + double psi_0; + double psi_1; + double Acon_0; + double Acon_1; + double Vgeo_0; + double Vgeo_1; + if (itag_true == i_0) { + psi_0 = psi[gm->i]; + psi_1 = psi[gm->j]; + Acon_0 = Acon0[gm->i]; + Acon_1 = Acon0[gm->j]; + Vgeo_0 = Vgeo[gm->i]; + Vgeo_1 = Vgeo[gm->j]; + } else { + psi_0 = psi[gm->j]; + psi_1 = psi[gm->i]; + Acon_0 = Acon0[gm->j]; + Acon_1 = Acon0[gm->i]; + Vgeo_0 = Vgeo[gm->j]; + Vgeo_1 = Vgeo[gm->i]; } + //if ( (itag_true == i_0 && jtag_true == i_1) || (itag_true == i_1 && jtag_true == i_0) ) { + // CSVWriter csvWriter("/Users/willzunker/simulations/lammps/avicel_tableting/avicel_tableting_pair_info.csv"); + // std::stringstream rowDataStream; + // rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision + // rowDataStream << itag_true << ", " << jtag_true << ", " << delta0 << ", " << delta1 << ", " << F0 << ", " << F1 << ", " << F << ", " << gm->radi << ", " << gm->radj << ", " << wij << ", " << psi_0 << ", " << psi_1 << ", " << aAdh0 << ", " << aAdh1 << ", " << Ac0 << ", " << Ac1 << ", " << F_MDR0 << ", " << F_MDR1 << ", " << F_BULK0 << ", " << F_BULK1 << ", " << deltaBULK0 << ", " << deltaBULK1 << ", " << Acon_0 << ", " << Acon_1 << ", " << Vgeo_0 << ", " << Vgeo_1 << ", " << kappa << ", " << lmp->update->ntimestep; + // std::string rowData = rowDataStream.str(); + // csvWriter.writeRow(rowData); + //} + + //if ( (gm->contact_type == PAIR) ) { + // CSVWriter csvWriter("/Users/willzunker/simulations/lammps/bulk_response/sticky/pair_info.csv"); + // std::stringstream rowDataStream; + // rowDataStream << std::scientific << std::setprecision(8); // Set the format and precision + // rowDataStream << itag_true << ", " << jtag_true << ", " << delta0 << ", " << delta1 << ", " << F0 << ", " << F1 << ", " << F << ", " << gm->radi << ", " << gm->radj << ", " << wij << ", " << psi_0 << ", " << psi_1 << ", " << aAdh0 << ", " << aAdh1 << ", " << Ac0 << ", " << Ac1 << ", " << F_MDR0 << ", " << F_MDR1 << ", " << F_BULK0 << ", " << F_BULK1 << ", " << deltaBULK0 << ", " << deltaBULK1 << ", " << Acon_0 << ", " << Acon_1 << ", " << Vgeo_0 << ", " << Vgeo_1 << ", " << kappa << ", " << F_DAMP << ", " << gm->vnnr << ", " << delta << ", " << lmp->update->ntimestep; + // std::string rowData = rowDataStream.str(); + // csvWriter.writeRow(rowData); + //} + return F; } diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index db96227f13..b096e03ca8 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -147,7 +147,7 @@ namespace Granular_NS { protected: double G, kappa, Eeff; // derived coeffs - double Eeffsq, Eeffinv, Eeffsqinv; + double Eeffsq, Eeffinv, Eeffsqinv, Eeff2particle; double gammasq, gamma3, gamma4; int warn_flag; @@ -155,6 +155,7 @@ namespace Granular_NS { int index_Ro, index_Vgeo, index_Velas, index_Vcaps, index_eps_bar, index_dRnumerator; int index_dRdenominator, index_Acon0, index_Acon1, index_Atot, index_Atot_sum, index_ddelta_bar; int index_psi, index_sigmaxx, index_sigmayy, index_sigmazz, index_contacts, index_adhesive_length; + int index_dRavg; int fix_mdr_flag; char *id_fix; diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index e96debd59e..448281cfa6 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -251,8 +251,8 @@ void GranularModel::init() // Must have valid normal, damping, and tangential models if (normal_model->name == "none") error->all(FLERR, "Must specify normal granular model"); if (normal_model->name == "mdr") { - if (damping_model->name != "none") - error->all(FLERR, "MDR require 'none' damping model. To damp, specify a coefficient of restitution < 1."); + //if (damping_model->name != "none") + //error->all(FLERR, "MDR require 'none' damping model. To damp, specify a coefficient of restitution < 1."); } else { if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); } diff --git a/src/csv_writer.h b/src/csv_writer.h new file mode 100644 index 0000000000..0a01604393 --- /dev/null +++ b/src/csv_writer.h @@ -0,0 +1,25 @@ +#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_; +}; \ No newline at end of file