main.cpp File Reference

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

Include dependency graph for main.cpp:

Go to the source code of this file.

Data Structures
struct	anonymous_namespace{main.cpp}::PackedParticle
class	anonymous_namespace{main.cpp}::GrainContactProbe

Namespaces
namespace	anonymous_namespace{main.cpp}

Functions
std::string	anonymous_namespace{main.cpp}::directoryPathWithTrailingSep (const std::filesystem::path &dir)
std::vector< PackedParticle >	anonymous_namespace{main.cpp}::generateCircularGrid (int ncols, int nrows, double x_lo, double y_lo, double R, double gap)
json	anonymous_namespace{main.cpp}::buildInputJson (const std::string &output_path_for_deck, const std::filesystem::path &mesh_cir, const std::filesystem::path &mesh_fixed, const std::filesystem::path &mesh_moving, double final_time, size_t num_steps, bool file_mesh, bool write_meshes, size_t search_interval, int ncols, int nrows, const std::string &mpi_strategy="auto")
void	anonymous_namespace{main.cpp}::writeLocations (const std::filesystem::path &csv, const std::vector< PackedParticle > &packed)
double	anonymous_namespace{main.cpp}::maxAbsForceFromCsv (const std::filesystem::path &csv)
int	main (int argc, char *argv[])

Function Documentation

main()

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

Definition at line 460 of file main.cpp.

                                 {
 
  util::parallel::initMpi(argc, argv);
  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::parallel::initNThreads(nThreads);
  util::io::print(std::format("Number of threads = {}\n", util::parallel::getNThreads()));
 
  // Close ~1.08 Rc at vy=-0.06, then squeeze. 0.004 s × 0.06 m/s = 0.24 mm.
  double final_time = 0.004;
  size_t num_steps = 20000;
  if (input.cmdOptionExists("-finalTime"))
    final_time = std::stod(input.getCmdOption("-finalTime"));
  if (input.cmdOptionExists("-numSteps"))
    num_steps = std::stoul(input.getCmdOption("-numSteps"));
  size_t search_interval = 40;
  if (input.cmdOptionExists("-searchInterval"))
    search_interval = std::stoul(input.getCmdOption("-searchInterval"));
  int ncols = 4;
  int nrows = 3;
  if (input.cmdOptionExists("-nCols"))
    ncols = std::stoi(input.getCmdOption("-nCols"));
  if (input.cmdOptionExists("-nRows"))
    nrows = std::stoi(input.getCmdOption("-nRows"));
  std::string mpi_strategy = "auto";
  if (input.cmdOptionExists("-mpiStrategy"))
    mpi_strategy = input.getCmdOption("-mpiStrategy");
  util::io::print(std::format("MPI_Strategy = {}, pack = {}x{} ({} grains)\n",
                              mpi_strategy, ncols, nrows, ncols * nrows));
 
  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"))
    inp_dir = out_dir.parent_path() / "inp";
 
  fs::create_directories(out_dir);
  fs::create_directories(inp_dir);
 
  const fs::path mesh_cir = inp_dir / "mesh_cir.msh";
  const fs::path mesh_fixed = inp_dir / "mesh_fixed_container.msh";
  const fs::path mesh_moving = inp_dir / "mesh_moving_container.msh";
  const bool file_mesh = input.cmdOptionExists("-fileMesh");
  const bool inbuilt = !file_mesh;
  const bool write_meshes = input.cmdOptionExists("-writeMeshes") || inbuilt;
  if (file_mesh) {
#ifndef JHA2021_COMP_MESH_DIR
    throw std::runtime_error("JHA2021_COMP_MESH_DIR is not set.");
#else
    const fs::path src(JHA2021_COMP_MESH_DIR);
    for (const char *name :
         {"mesh_cir.msh", "mesh_fixed_container.msh", "mesh_moving_container.msh"}) {
      const fs::path from = src / name;
      if (!fs::exists(from))
        throw std::runtime_error("Missing frozen mesh " + from.string());
      fs::copy_file(from, inp_dir / name, fs::copy_options::overwrite_existing);
    }
    util::io::print(std::format("Using frozen meshes from {}\n", src.string()));
#endif
  }
 
  const std::string output_path_for_deck = directoryPathWithTrailingSep(out_dir);
  auto inputJson = buildInputJson(output_path_for_deck, mesh_cir, mesh_fixed, mesh_moving,
                                  final_time, num_steps, file_mesh, write_meshes,
                                  search_interval, ncols, nrows, mpi_strategy);
 
  {
    std::ofstream os(inp_dir / "input.json");
    os << inputJson.dump(2);
  }
 
  {
    const double R = 0.001;
    const double mesh_size = R / 5.0;
    const double pad = 1.15 * 0.95 * 0.7 * mesh_size;
    writeLocations(inp_dir / "particle_locations.csv",
                   generateCircularGrid(ncols, nrows, 0.0, 0.0, R, pad));
  }
 
  auto deck = std::make_shared<inp::Input>(inputJson);
  PeriDEMModel dem(deck);
  const size_t ts_every = std::max<size_t>(1, num_steps / 100);
  auto probe = std::make_unique<GrainContactProbe>(out_dir, ts_every);
  GrainContactProbe *probe_p = probe.get();
  dem.setPostprocess(std::move(probe));
  dem.run(deck);
 
  // Under particle-MPI, ranks only refresh ghost centers; far grains stay stale.
  // Reduce contact probe metrics so requireContact is rank-consistent.
  int max_pairs = probe_p->maxPairs();
  double min_gap = probe_p->minGap();
  double gap0 = probe_p->gap0();
  if (util::parallel::mpiSize() > 1) {
    int max_pairs_g = 0;
    double min_gap_g = 0., gap0_g = 0.;
    MPI_Allreduce(&max_pairs, &max_pairs_g, 1, MPI_INT, MPI_MAX,
                  util::parallel::mpiComm());
    MPI_Allreduce(&min_gap, &min_gap_g, 1, MPI_DOUBLE, MPI_MIN,
                  util::parallel::mpiComm());
    MPI_Allreduce(&gap0, &gap0_g, 1, MPI_DOUBLE, MPI_MAX,
                  util::parallel::mpiComm());
    max_pairs = max_pairs_g;
    min_gap = min_gap_g;
    gap0 = gap0_g;
  }
 
  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(
      "grain contact: gap0={}, min_gap={} (t={}), Rc={}, pairs_in_Rc={}, max Damage_Z={}\n",
      gap0, min_gap, probe_p->tMin(), probe_p->Rc(), max_pairs, zmax));
 
  const bool require_contact = !input.cmdOptionExists("-noRequireContact");
  if (require_contact) {
    const bool grains_touched = max_pairs > 0 && min_gap < gap0 - 1.0e-8;
    if (!grains_touched) {
      util::io::print("requireContact: no grain–grain pair entered the contact radius.\n");
      util::parallel::finalizeMpi();
      return EXIT_FAILURE;
    }
    if (!(gap0 > probe_p->Rc())) {
      util::io::print("requireContact: grains already in Rc at t=0; packing is too tight.\n");
      util::parallel::finalizeMpi();
      return EXIT_FAILURE;
    }
  }
 
  if (input.cmdOptionExists("-assertForce")) {
    // Rank 0 owns the reaction CSV; broadcast pass/fail.
    int force_ok = 1;
    if (util::parallel::mpiRank() == 0) {
      const fs::path csv = out_dir / "pp_compressive_test_0.csv";
      const double fmax = maxAbsForceFromCsv(csv);
      util::io::print(std::format("assertForce: max |plate reaction| = {} from {}\n", fmax,
                                  csv.string()));
      if (fmax <= 0.) {
        util::io::print("assertForce: plate reaction is zero.\n");
        force_ok = 0;
      }
    }
    if (util::parallel::mpiSize() > 1)
      MPI_Bcast(&force_ok, 1, MPI_INT, 0, util::parallel::mpiComm());
    if (!force_ok) {
      util::parallel::finalizeMpi();
      return EXIT_FAILURE;
    }
  }
 
  // Metric for serial vs MPI checks (rank 0). Allreduce max|u| for DOF-MPI.
  {
    double max_u = 0.;
    for (const auto &u : dem.d_u)
      max_u = std::max(max_u, u.length());
    double com0x = 0., com0y = 0.;
    if (!dem.d_particlesListTypeParticle.empty()) {
      const auto c = dem.d_particlesListTypeParticle[0]->getXCenter();
      com0x = c.d_x;
      com0y = c.d_y;
    }
    if (util::parallel::mpiSize() > 1)
      MPI_Allreduce(MPI_IN_PLACE, &max_u, 1, MPI_DOUBLE, MPI_MAX,
                    util::parallel::mpiComm());
    if (util::parallel::mpiRank() == 0) {
      std::ofstream os(out_dir / "mpi_metric.txt");
      os << std::format("{:.12e} {:.12e} {:.12e}\n", max_u, com0x, com0y);
      util::io::print(std::format("mpi_metric: max|u|={:.12e} grain0_com=({:.12e},{:.12e})\n",
                                  max_u, com0x, com0y));
    }
  }
 
  util::parallel::finalizeMpi();
  return EXIT_SUCCESS;
}

References anonymous_namespace{main.cpp}::buildInputJson(), util::io::InputParser::cmdOptionExists(), data::ModelData::d_particlesListTypeParticle, data::ModelData::d_u, data::ModelData::d_Z, util::parallel::finalizeMpi(), anonymous_namespace{main.cpp}::GrainContactProbe::gap0(), anonymous_namespace{main.cpp}::generateCircularGrid(), util::io::InputParser::getCmdOption(), util::parallel::getNThreads(), util::parallel::initMpi(), util::parallel::initNThreads(), anonymous_namespace{main.cpp}::maxAbsForceFromCsv(), util::parallel::mpiComm(), util::parallel::mpiRank(), util::parallel::mpiSize(), util::io::print(), PeriDEMModel::run(), PeriDEMModel::setPostprocess(), and anonymous_namespace{main.cpp}::writeLocations().

Here is the call graph for this function: