main.cpp File Reference

#include "contact/contact.h"
#include "geom/geomIncludes.h"
#include "inp/deckIncludes.h"
#include "material/materialUtil.h"
#include "periDEMModel.h"
#include "util/io.h"
#include "util/parallelUtil.h"
#include "util/point.h"
#include <cmath>
#include <cstdlib>
#include <filesystem>
#include <format>
#include <fstream>
#include <memory>
#include <stdexcept>
#include <string>

Include dependency graph for main.cpp:

Go to the source code of this file.

Data Structures
struct	anonymous_namespace{main.cpp}::Probe

Namespaces
namespace	anonymous_namespace{main.cpp}

Functions
std::string	anonymous_namespace{main.cpp}::directoryPathWithTrailingSep (const std::filesystem::path &dir)
json	anonymous_namespace{main.cpp}::buildInputJson (const std::string &output_path, const std::filesystem::path &mesh_cir, const std::filesystem::path &mesh_wall, double R, double mesh_size, double horizon, double Rc_factor, double Kn, double final_time, size_t num_steps, bool policy_combo=false)
Probe	anonymous_namespace{main.cpp}::checkWallSpringOnly (PeriDEMModel &dem, geom::GeomObject &wall_geom)
double	anonymous_namespace{main.cpp}::minGrainWallGap (PeriDEMModel &dem, geom::GeomObject &wall_geom)
int	main (int argc, char *argv[])

Function Documentation

main()

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

Definition at line 261 of file main.cpp.

                                 {
  util::parallel::initMpi(argc, argv);
  util::io::InputParser input(argc, argv);
 
  unsigned n_threads = 2;
  if (input.cmdOptionExists("-nThreads"))
    n_threads = static_cast<unsigned>(std::stoi(input.getCmdOption("-nThreads")));
  util::parallel::initNThreads(n_threads);
 
  namespace fs = std::filesystem;
  fs::path base = fs::current_path() / "wall_analytical_run";
  if (input.cmdOptionExists("-outputDir"))
    base = input.getCmdOption("-outputDir");
  const fs::path out_dir = base / "out";
  const fs::path inp_dir = base / "inp";
  fs::create_directories(out_dir);
  fs::create_directories(inp_dir);
 
  const double R = 0.001;
  const double mesh_size = R / 5.0;
  const double horizon = 3.0 * mesh_size;
  const double Rc_factor = 0.95;
  const double Kn = 1.0e11;
  const double final_time = 5.0e-4;
  const size_t num_steps = 5000;
 
  const bool policy_combo = input.cmdOptionExists("-policyCombo");
  auto input_json =
      buildInputJson(directoryPathWithTrailingSep(out_dir), inp_dir / "mesh_cir.msh",
                     inp_dir / "mesh_wall.msh", R, mesh_size, horizon, Rc_factor, Kn,
                     final_time, num_steps, policy_combo);
  {
    std::ofstream os(inp_dir / "input.json");
    os << input_json.dump(2);
  }
 
  // -deckOnly stops after the deck is written. The deck 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;
  }
 
  auto deck = std::make_shared<inp::Input>(input_json);
  PeriDEMModel dem(deck);
  dem.init();
 
  if (dem.d_modelDeck_p->d_wallContact != "analytical_plane")
    throw std::runtime_error("analytical wall: Wall_Contact not analytical_plane");
  if (policy_combo) {
    if (dem.d_modelDeck_p->d_selfContact != "reference_gap")
      throw std::runtime_error("alternate deck: Self_Contact is not reference_gap");
    if (dem.d_modelDeck_p->d_bondBreak != "absolute_stretch")
      throw std::runtime_error("alternate deck: Bond_Break is not absolute_stretch");
    if (dem.d_particleDeck_p->d_contactDeck.d_frictionLaw != "stick_slip")
      throw std::runtime_error("alternate deck: Friction_Law is not stick_slip");
    if (dem.d_particleDeck_p->d_contactDeck.d_dampingLaw != "node")
      throw std::runtime_error("alternate deck: Damping_Law is not node");
    if (!dem.d_contact_p || !dem.d_contact_p->d_pairForce ||
        dem.d_contact_p->d_useNodeDamping != true)
      throw std::runtime_error("alternate deck: node damping not active on contact");
    util::io::print("alternate deck fields applied\n");
  }
  if (!dem.d_contact_p || !dem.d_contact_p->d_wallContact ||
      !dem.d_contact_p->d_wallContact->skipsMeshedGrainWall())
    throw std::runtime_error("analytical wall: must skip meshed grain-wall");
 
  const double Rc = dem.d_particleDeck_p->d_contactDeck.getContact(0, 1).d_contactR;
  // Floor top face is at y=0 (matches rectangle placement in the deck).
  geom::Plane floor_plane(util::Point(0., 1., 0.), util::Point(0., 0., 0.));
 
  // Wall spring-only force density vs closed form.
  {
    const double push = 0.5 * Rc;
    for (size_t i = 0; i < dem.d_x.size(); ++i) {
      if (dem.d_ptId[i] != 0)
        continue;
      dem.d_u[i].d_y -= push;
      dem.d_x[i] = dem.d_xRef[i] + dem.d_u[i];
      dem.d_v[i] = util::Point();
    }
    const auto frozen = checkWallSpringOnly(dem, floor_plane);
    util::io::print(std::format(
        "analytical wall formula: n_contact={}, min_gap={}, Rc={}, max|f_ref|={}, max|f-f_ref|={}\n",
        frozen.n_contact_nodes, frozen.min_gap, Rc, frozen.max_force_ref,
        frozen.max_force_err));
    if (frozen.n_contact_nodes == 0)
      throw std::runtime_error("analytical wall: no grain nodes in contact");
    if (!(frozen.max_force_ref > 0.))
      throw std::runtime_error("analytical wall: expected nonzero force");
    if (!(frozen.max_force_err <= 1.e-6 * (1. + frozen.max_force_ref)))
      throw std::runtime_error(std::format(
          "analytical wall: force density mismatch vs Kn*(gap-Rc)*voli*(-n): err={}",
          frozen.max_force_err));
    if (!(frozen.min_gap > -Rc))
      throw std::runtime_error(std::format(
          "analytical wall: deep penetration past Rc in contact config: min_gap={}, Rc={}",
          frozen.min_gap, Rc));
  }
 
  // Full contact assembly (PD + analytical wall) must keep the same gap bound.
  dem.computeForces();
  const double ming = minGrainWallGap(dem, floor_plane);
  util::io::print(std::format("analytical wall after computeForces: min_gap={}, Rc={}\n", ming, Rc));
  if (!(ming < Rc))
    throw std::runtime_error(
        "analytical wall: expected grain-wall contact (min_gap < Rc)");
  if (!(ming > -Rc))
    throw std::runtime_error(std::format(
        "analytical wall: deep penetration past Rc: min_gap={}, Rc={}", ming, Rc));
 
  dem.close();
  util::io::print("analytical wall OK\n");
  return 0;
}

References anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::checkWallSpringOnly(), PeriDEMModel::close(), util::io::InputParser::cmdOptionExists(), PeriDEMModel::computeForces(), PeriDEMModel::d_contact_p, data::ModelData::d_modelDeck_p, data::ModelData::d_particleDeck_p, data::ModelData::d_ptId, data::ModelData::d_u, data::ModelData::d_v, data::ModelData::d_x, data::ModelData::d_xRef, util::parallel::finalizeMpi(), util::io::InputParser::getCmdOption(), PeriDEMModel::init(), util::parallel::initMpi(), util::parallel::initNThreads(), main(), anonymous_namespace{main.cpp}::minGrainWallGap(), and util::io::print().

Here is the call graph for this function: