#include "inp/deckIncludes.h"
#include "util/io.h"
#include "util/function.h"
#include "util/parallelUtil.h"
#include "geom/geomObjectsUtil.h"
#include "material/materialUtil.h"
#include "periDEMModel.h"
#include "particle/baseParticle.h"
#include "particle/particleMpi.h"
#include "postprocess/postprocess.h"
#include <mpi.h>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <filesystem>
#include <format>
#include <fstream>
#include <memory>
#include <stdexcept>
#include <thread>
#include <utility>
#include <vector>

Include dependency graph for main.cpp:

Data Structures
struct	Jha2021Table2

class	RestitutionProbe

class	LateralProbe

class	MpiMetricTsProbe

Namespaces
namespace	anonymous_namespace{main.cpp}

Functions
Jha2021Table2	jha2021Table2 (int test_id)

std::string	anonymous_namespace{main.cpp}::directoryPathWithTrailingSep (const std::filesystem::path &dir)

json	anonymous_namespace{main.cpp}::buildInputJson (const std::string &output_path_for_deck, const std::filesystem::path &mesh_file_1, const std::filesystem::path &mesh_file_2, double final_time, size_t num_steps, bool zero_ic, double mesh_size_in, double horizon_in, bool damping_on, double eps_n, bool two_particle_test, bool file_mesh, double beta_n_factor, const std::string &damping_law="com_and_node", const std::string &friction_law="coulomb_simple", bool friction_on=false, double friction_mu=-1., double ic_vx=0., bool bottom_patch_bc=false, const std::string &mpi_strategy="auto")

int	main (int argc, char *argv[])

Function Documentation

◆ jha2021Table2()

Jha2021Table2 jha2021Table2 ( int test_id )

Definition at line 75 of file main.cpp.

                                         {
  static const Jha2021Table2 k[] = {
      {1.0, 0.996}, {0.95, 0.925}, {0.9, 0.861}, {0.85, 0.801}, {0.8, 0.744}};
  if (test_id < 1 || test_id > 5)
    throw std::runtime_error("jha2021Table2Test must be 1..5");
  return k[test_id - 1];
}

Referenced by main().

Here is the caller graph for this function:

◆ main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 549 of file main.cpp.

                                 {
 
  util::parallel::initMpi(argc, argv);
  int mpiSize = util::parallel::mpiSize(), mpiRank = util::parallel::mpiRank();
  util::io::print(std::format("Initialized MPI. MPI size = {}, MPI rank = {}\n", mpiSize, mpiRank));
  util::io::print(util::parallel::getMpiStatus()->printStr());
 
  util::io::InputParser input(argc, argv);
 
  unsigned int nThreads;
  if (input.cmdOptionExists("-nThreads"))
    nThreads = std::stoi(input.getCmdOption("-nThreads"));
  else {
    nThreads = std::thread::hardware_concurrency();
    util::io::print(std::format("Running test with default number of threads = {}\n", nThreads));
  }
  util::parallel::initNThreads(nThreads);
  util::io::print(std::format("Number of threads = {}\n", util::parallel::getNThreads()));
 
#ifndef TWOP_CONTACT_EXAMPLE
  double final_time = 0.002;
  size_t num_steps = 6000;
#else
  double final_time = 0.012;
  size_t num_steps = 36000;
#endif
  bool zero_ic = false;
  double mesh_size_in = -1.;
  double horizon_in = -1.;
  bool assert_cr = false;
  double cr_ref = 1.;
  double cr_tol = 0.1;
  bool damping_on = false;
  double eps_n = 0.9;
  double beta_n_factor = 100.0;
  std::string damping_law = "com_and_node";
  std::string friction_law = "coulomb_simple";
  bool friction_on = false;
  double friction_mu = -1.;
  double ic_vx = 0.;
  bool assert_lat = false;
  double lat_min = 0.;
  bool bottom_patch_bc = false;
  std::string mpi_strategy = "auto";
  bool write_mpi_metric = false;
  int jha_test = 0;
  if (input.cmdOptionExists("-jha2021Table2Test"))
    jha_test = std::stoi(input.getCmdOption("-jha2021Table2Test"));
  else if (input.cmdOptionExists("-jha2021Table2"))
    jha_test = 1;
  if (jha_test > 0) {
    const auto c = jha2021Table2(jha_test);
    // Paper §4.1 / v0.1.0 circ_damp: T=0.04 s, Δt=0.2 μs, horizon 0.6 mm.
    // Mesh is the frozen v0.1.0 .msh (hmin = 0.1423 mm), not in-process Gmsh.
    final_time = 0.04;
    num_steps = 200000;
    horizon_in = 0.0006;
    // v0.1.0 circ_damp: not drop-from-rest. Gap is shortened to the horizon
    // and the leftover 0.4 mm of the 1 mm drop is an impact velocity.
    zero_ic = false;
    assert_cr = true;
    cr_ref = c.cr;
    cr_tol = 0.05;
    eps_n = c.eps;
    damping_on = c.eps < 1.0 - 1.0e-12;
  }
  if (input.cmdOptionExists("-finalTime"))
    final_time = std::stod(input.getCmdOption("-finalTime"));
  if (input.cmdOptionExists("-numSteps"))
    num_steps = std::stoul(input.getCmdOption("-numSteps"));
  if (input.cmdOptionExists("-zeroIC"))
    zero_ic = true;
  if (input.cmdOptionExists("-meshSize"))
    mesh_size_in = std::stod(input.getCmdOption("-meshSize"));
  if (input.cmdOptionExists("-horizon"))
    horizon_in = std::stod(input.getCmdOption("-horizon"));
  if (input.cmdOptionExists("-epsN")) {
    eps_n = std::stod(input.getCmdOption("-epsN"));
    damping_on = eps_n < 1.0 - 1.0e-12;
  }
  if (input.cmdOptionExists("-betaNFactor"))
    beta_n_factor = std::stod(input.getCmdOption("-betaNFactor"));
  if (input.cmdOptionExists("-dampingLaw"))
    damping_law = input.getCmdOption("-dampingLaw");
  if (input.cmdOptionExists("-frictionLaw"))
    friction_law = input.getCmdOption("-frictionLaw");
  if (input.cmdOptionExists("-enableFriction"))
    friction_on = true;
  if (input.cmdOptionExists("-mu"))
    friction_mu = std::stod(input.getCmdOption("-mu"));
  if (input.cmdOptionExists("-icVx"))
    ic_vx = std::stod(input.getCmdOption("-icVx"));
  if (input.cmdOptionExists("-assertLateralVx")) {
    assert_lat = true;
    lat_min = std::stod(input.getCmdOption("-assertLateralVx"));
    assert_cr = false; // skip CR gate from table2 defaults for this lateral check
  }
  if (input.cmdOptionExists("-assertCR")) {
    assert_cr = true;
    cr_ref = std::stod(input.getCmdOption("-assertCR"));
  }
  if (input.cmdOptionExists("-crTol"))
    cr_tol = std::stod(input.getCmdOption("-crTol"));
  if (input.cmdOptionExists("-bottomPatchBC"))
    bottom_patch_bc = true;
  if (input.cmdOptionExists("-mpiStrategy"))
    mpi_strategy = input.getCmdOption("-mpiStrategy");
  if (input.cmdOptionExists("-writeMpiMetric"))
    write_mpi_metric = true;
  util::io::print(std::format(
      "final_time = {}, num_steps = {}, zero_ic = {}, eps_n = {}, C_bar = {}, "
      "damping = {}, Damping_Law = {}, Friction_Law = {}, bottomPatchBC = {}, "
      "MPI_Strategy = {}\n",
      final_time, num_steps, zero_ic, eps_n, beta_n_factor, damping_on,
      damping_law, friction_law, bottom_patch_bc, mpi_strategy));
 
  namespace fs = std::filesystem;
  const fs::path cwd = fs::current_path();
 
  fs::path out_dir = cwd / "out";
  fs::path inp_dir = cwd / "inp";
 
  if (input.cmdOptionExists("-outputDir")) {
    fs::path p = input.getCmdOption("-outputDir");
    out_dir = p.is_absolute() ? std::move(p) : cwd / p;
  }
  if (input.cmdOptionExists("-inputDir")) {
    fs::path p = input.getCmdOption("-inputDir");
    inp_dir = p.is_absolute() ? std::move(p) : cwd / p;
  } else if (input.cmdOptionExists("-outputDir")) {
    // Place inp next to out: .../run/out -> .../run/inp
    inp_dir = out_dir.parent_path() / "inp";
  }
 
  fs::create_directories(out_dir);
  fs::create_directories(inp_dir);
 
  const std::string output_path_for_deck = directoryPathWithTrailingSep(out_dir);
  const fs::path mesh1 = inp_dir / "mesh_cir_1.msh";
  const fs::path mesh2 = inp_dir / "mesh_cir_2.msh";
  const bool inbuilt_mesh = input.cmdOptionExists("-inbuiltMesh");
  const bool file_mesh = jha_test > 0 && !inbuilt_mesh;
  if (file_mesh) {
#ifndef JHA2021_MESH_DIR
    throw std::runtime_error("JHA2021_MESH_DIR is not set (v0.1.0 circ_damp meshes).");
#else
    const fs::path src(JHA2021_MESH_DIR);
    fs::copy_file(src / "mesh_cir_1.msh", mesh1, fs::copy_options::overwrite_existing);
    fs::copy_file(src / "mesh_cir_2.msh", mesh2, fs::copy_options::overwrite_existing);
    util::io::print(std::format("Using v0.1.0 circ_damp meshes from {}\n", src.string()));
#endif
  } else if (jha_test > 0)
    util::io::print("Using in-process Gmsh (gmsh_builtin_mesh) instead of frozen v0.1.0 .msh\n");
 
  util::io::print(std::format("Output directory (VTU, log.txt): {}\n", fs::absolute(out_dir).string()));
  util::io::print(std::format("Input directory (input.json, meshes): {}\n", fs::absolute(inp_dir).string()));
 
  auto inputJson = buildInputJson(output_path_for_deck, mesh1, mesh2, final_time,
                                  num_steps, zero_ic, mesh_size_in, horizon_in,
                                  damping_on, eps_n, jha_test > 0, file_mesh,
                                  beta_n_factor, damping_law,
                                  friction_law, friction_on, friction_mu,
                                  ic_vx, bottom_patch_bc, mpi_strategy);
 
  const fs::path input_json_path = inp_dir / "input.json";
  {
    std::ofstream os(input_json_path);
    if (!os)
      throw std::runtime_error("Failed to open " + input_json_path.string() + " for writing.");
    os << inputJson.dump(2);
  }
  util::io::print(std::format("Wrote deck to {}\n", fs::absolute(input_json_path).string()));
 
  // -deckOnly stops after the deck is written. The Python comparison reads
  // the deck and does not need this driver to run the simulation.
  if (input.cmdOptionExists("-deckOnly")) {
    util::io::print("deck written; -deckOnly, not running\n");
    util::parallel::finalizeMpi();
    return 0;
  }
  if (jha_test > 0) {
    const auto &p1 = inputJson["Particle_Generation"]["Data"]["1"];
    const auto &v = inputJson["IC"]["Set_1"]["Constant_Velocity"]["Velocity_Vector"];
    const double y1 = p1.at("y").get<double>();
    const double gap0 = y1 - 0.001 - 0.001 - 0.001; // y_top - R_bottom_center - R_bot - R_top
    util::io::print(std::format(
        "Table 2 kinematics: top y = {}, surface gap = {}, v_y = {}\n",
        y1, gap0, v.at(1).get<double>()));
  }
 
  auto deck = std::make_shared<inp::Input>(inputJson);
 
  PeriDEMModel dem(deck);
  RestitutionProbe *probe = nullptr;
  LateralProbe *lat_probe = nullptr;
  if (assert_cr) {
    auto p = std::make_unique<RestitutionProbe>(0.001);
    probe = p.get();
    dem.setPostprocess(std::move(p));
  } else if (assert_lat) {
    auto p = std::make_unique<LateralProbe>(ic_vx);
    lat_probe = p.get();
    dem.setPostprocess(std::move(p));
  } else if (write_mpi_metric) {
    const size_t ts_every = std::max<size_t>(1, num_steps / 100);
    dem.setPostprocess(
        std::make_unique<MpiMetricTsProbe>(out_dir, ts_every));
  }
  dem.run(deck);
 
  if (input.cmdOptionExists("-requireContact")) {
    const float zmax =
        dem.d_Z.empty() ? 0.f : *std::max_element(dem.d_Z.begin(), dem.d_Z.end());
    util::io::print(std::format("requireContact: max Damage_Z = {}\n", zmax));
    if (zmax <= 0.f) {
      util::io::print("requireContact: no damage; particles did not contact.\n");
      return EXIT_FAILURE;
    }
  }
 
  if (assert_cr) {
    {
      const fs::path gap_csv = out_dir / "gap.csv";
      std::ofstream gs(gap_csv);
      gs << "t,gap\n";
      for (const auto &s : probe->samples())
        gs << s.first << "," << s.second << "\n";
      util::io::print(std::format("Wrote {}\n", fs::absolute(gap_csv).string()));
    }
    const double cr = probe->CR();
    util::io::print(std::format(
        "assertCR: contacted={}, min_gap={} (t={}), H1={} (t={}), CR={}, ref={}, tol={}\n",
        probe->contacted(), probe->minGap(), probe->tMin(), probe->H1(),
        probe->tH1(), cr, cr_ref, cr_tol));
    if (cr > 1.02) {
      util::io::print(
          "assertCR: CR > 1.02 on the 1 mm energy-equivalent drop "
          "(surface gap = horizon, leftover height already in v_y). "
          "Unphysical energy gain.\n");
      return EXIT_FAILURE;
    }
    if (!probe->contacted() || std::abs(cr - cr_ref) > cr_tol) {
      util::io::print("assertCR: coefficient of restitution out of range.\n");
      return EXIT_FAILURE;
    }
  }
 
  if (assert_lat) {
    const double vmin = lat_probe->minAbsVxAfterContact();
    const double vmax = lat_probe->maxAbsVx();
    const double ic = lat_probe->icAbsVx();
    util::io::print(std::format(
        "assertLateralVx: contacted={}, max|vx|={}, min|vx|_after_contact={}, "
        "ic|vx|={}, slow_factor_max={}\n",
        lat_probe->contacted(), vmax, vmin, ic, lat_min));
    if (!lat_probe->contacted()) {
      util::io::print("assertLateralVx: particles never contacted.\n");
      return EXIT_FAILURE;
    }
    // lat_min is the maximum allowed ratio min_vx/ic_vx after contact (e.g. 0.95).
    if (!(ic > 0.) || !(vmin < lat_min * ic)) {
      util::io::print(
          "assertLateralVx: friction did not slow lateral speed enough "
          "(or kinematics missed).\n");
      return EXIT_FAILURE;
    }
  }
 
  if (write_mpi_metric) {
    double max_u = 0.;
    for (const auto &u : dem.d_u)
      max_u = std::max(max_u, u.length());
    double com1x = 0., com1y = 0.;
    if (dem.d_particlesListTypeParticle.size() > 1) {
      const auto c = dem.d_particlesListTypeParticle[1]->getXCenter();
      com1x = c.d_x;
      com1y = c.d_y;
    }
    if (util::parallel::isMpiEnabled()) {
      MPI_Allreduce(MPI_IN_PLACE, &max_u, 1, MPI_DOUBLE, MPI_MAX,
                    util::parallel::mpiComm());
      // Center node is unique; take from owning rank via MAX of abs then
      // broadcast of values — simpler: Allreduce SUM after zeroing non-owners
      // is wrong for COM. After DOF sync every rank has the same center x.
    }
    if (util::parallel::mpiRank() == 0) {
      std::ofstream os(out_dir / "mpi_metric.txt");
      os << std::format("{:.12e} {:.12e} {:.12e}\n", max_u, com1x, com1y);
      util::io::print(std::format(
          "mpi_metric: max|u|={:.12e} grain1_com=({:.12e},{:.12e})\n", max_u,
          com1x, com1y));
    }
  }
 
  util::parallel::finalizeMpi();
  return EXIT_SUCCESS;
}

Here is the call graph for this function:

Data Structures

Namespaces

Functions

Function Documentation

◆ jha2021Table2()

◆ main()