main.cpp

Go to the documentation of this file.

/*
 * -------------------------------------------
 * Copyright (c) 2021 - 2026 Prashant K. Jha
 * -------------------------------------------
 * PeriDEM https://github.com/prashjha/PeriDEM
 *
 * Distributed under the Boost Software License, Version 1.0. (See accompanying
 * file LICENSE)
 */
 
#include "inp/deckIncludes.h"
#include "geom/geomObjects.h"
#include "mesh_gen/meshGenerator.h"
#include "util/io.h"
#include "util/function.h"
#include "material/materialUtil.h"
#include "periDEMModel.h"
#include <filesystem>
#include <format>
#include <fstream>
#include <cmath>
#include <memory>
#include <thread>
 
json getInputJson();
 
int main(int argc, char *argv[]) {
 
    // init parallel
  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));
  }
  // set number of threads
  util::parallel::initNThreads(nThreads);
  util::io::print(std::format("Number of threads = {}\n", util::parallel::getNThreads()));
  
  auto inputJson = getInputJson();
 
  // Create input deck
  auto deck = std::make_shared<inp::Input>(inputJson);
 
  // PeriDEM
  PeriDEMModel dem(deck);
  dem.run(deck);
 
  return EXIT_SUCCESS;
}
 
json getInputJson() {
 
  // Create output directory
  std::string output_dir = "./out/";
  std::string input_dir = "./inp/";
  std::filesystem::create_directories(output_dir);
  std::filesystem::create_directories(input_dir);
 
  // Simulation parameters
  const std::vector<double> center = {0.0, 0.0, 0.0};
  const double R1 = 0.001;  // Bottom particle radius
  const double R2 = 0.001;  // Top particle radius
  const double mesh_size = std::min(R1, R2) / 5.0;
  const double horizon = 3.0 * mesh_size;
  const double particle_dist = 0.001;
 
  // Material parameters
  const double poisson1 = 0.25;
  const double rho1 = 1200.0;
  const double K1 = 2.16e+7;
  const double E1 = material::toE(K1, poisson1);
  const double G1 = material::toGE(E1, poisson1);
  const double Gc1 = 50.0;
 
  const double poisson2 = 0.25;
  const double rho2 = 1200.0;
  const double K2 = 2.16e+7;
  const double E2 = material::toE(K2, poisson2);
  const double G2 = material::toGE(E2, poisson2);
  const double Gc2 = 50.0;
 
  // Contact parameters
  const double R_contact_factor = 0.95;
  const double Kn_11 = util::normalContactStiffness(K1, K1, horizon);
  const double Kn_22 = util::normalContactStiffness(K2, K2, horizon);
  const double Kn_12 = util::normalContactStiffness(K1, K2, horizon);
  const double beta_n_eps = 0.9;
  const double friction_coeff = 0.5;
  const double beta_n_factor = 100.0;
 
  // Generate meshes for both particles
  
  // Bottom particle (Zone 1)
  std::vector<double> p1_center = center;
  std::string mesh1_file_name = input_dir + "mesh_cir_1";
  auto mesh_geom_1 = std::make_shared<geom::Circle>(
      R1, util::Point(p1_center[0], p1_center[1], p1_center[2]));
  mesh_gen::generateBuiltinParticleMeshGmsh(mesh_geom_1, mesh_size, mesh1_file_name, false, true,
                                              nullptr, nullptr, nullptr);
 
  // Top particle (Zone 2)
  std::vector<double> p2_center = center;
  std::string mesh2_file_name = input_dir + "mesh_cir_2";
  auto mesh_geom_2 = std::make_shared<geom::Circle>(
      R2, util::Point(p2_center[0], p2_center[1], p2_center[2]));
  mesh_gen::generateBuiltinParticleMeshGmsh(mesh_geom_2, mesh_size, mesh2_file_name, false, true,
                                              nullptr, nullptr, nullptr);
 
  // Model deck
  const double final_time = 0.0001;
  const size_t num_steps = 1000;
  const size_t dt_out_n = num_steps / 4;
  auto modelDeckJson = inp::ModelDeck::getExampleJson(
      json{{"Dimension", 2},
           {"Final_Time", final_time},
           {"Time_Steps", num_steps},
           {"Discretization_Type", json{{"Spatial", "finite_difference"}, {"Time", "central_difference"}}},
           {"Populate_ElementNodeConnectivity", true},
           {"Quad_Approximation_Order", 2},
           {"Particle_Sim_Type", "Multi_Particle"},
           {"Seed", 0}});
 
  // Output deck
  auto outputDeckJson = inp::OutputDeck::getExampleJson(
      json{{"File_Format", "vtu"},
           {"Path", output_dir},
           {"Tags", std::vector<std::string>({"Displacement", "Velocity", "Force", "Damage_Z", "Damage", "Particle_ID"})},
           {"Output_Interval", dt_out_n},
           {"Debug", 2},
           {"Perform_FE_Out", true},
           {"Compress_Type", "zlib"},
           {"Perform_Out", true},
           {"Test_Output_Interval", 1},
           {"Tag_PP", ""}});
 
  // BC deck
  auto bcDeckJson = inp::BCDeck::getExampleJson(
      json{{"Displacement_BC_Sets", 1},
           {"IC_Sets", 1},
           {"Gravity", util::Point(0, -10, 0).toVec()}});
 
  // Displacement BC for fixing bottom particle
  bcDeckJson["Displacement_BC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "Displacement_BC"},
           {"Particle_List", {0}},
           {"Direction", {1, 2}},
           {"Zero_Displacement", true}});
 
  // Initial velocity for top particle
  const double free_fall_dist = particle_dist - horizon;
  const double free_fall_vel = -std::sqrt(2.0 * std::abs(-10.0) * free_fall_dist);
  bcDeckJson["IC"]["Set_1"] = inp::BCBaseDeck::getExampleJson(
      json{{"Type", "IC"},
           {"Particle_List", {1}},
           {"IC_Type", "Constant_Velocity"},
           {"IC_Vector", {0.0, free_fall_vel, 0.0}}});
 
  // Particle deck
  auto pDeckJson = json({});
 
  // Particle geometry
  std::vector<geom::GeomData> pGeomVec(2);
  
  // Bottom particle
  pGeomVec[0].d_geomName = "circle";
  pGeomVec[0].d_geomParams = {R1, p1_center[0], p1_center[1], p1_center[2]};
 
  // Top particle
  pGeomVec[1].d_geomName = "circle";
  pGeomVec[1].d_geomParams = {R2, p2_center[0], p2_center[1], p2_center[2]};
 
  pDeckJson["Particle"] = inp::ParticleDeck::getParticleGeomExampleJson(pGeomVec);
 
  // Mesh settings
  pDeckJson["Mesh"] = inp::ParticleDeck::getParticleMeshExampleJson({
      mesh1_file_name + ".msh",
      mesh2_file_name + ".msh"
  });
 
  // Material settings
  auto pMatJson = inp::ParticleDeck::getParticleMaterialExampleJson(2);
 
  // Material 1 (bottom particle)
  pMatJson["Set_1"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho1},
           {"K", K1},
           {"G", G1},
           {"Gc", Gc1},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
 
  // Material 2 (top particle)
  pMatJson["Set_2"] = inp::MaterialDeck::getExampleJson(
      json{{"Type", "PDState"},
           {"Is_Plane_Strain", false},
           {"Horizon", horizon},
           {"Density", rho2},
           {"K", K2},
           {"G", G2},
           {"Gc", Gc2},
           {"Compute_From_Classical", true},
           {"Influence_Function", json{{"Type", 0}}}});
 
  pDeckJson["Material"] = pMatJson;
 
  // Contact settings
  auto pContactJson = inp::ParticleDeck::getParticleContactExampleJson(2);
 
  // Base contact parameters
  json contact_base = inp::ContactPairDeck::getExampleJson(
          json{{"Contact_Radius_Factor", R_contact_factor},
               {"Kn", Kn_11},
               {"Damping_On", false},
               {"Epsilon", beta_n_eps},
               {"Friction_On", false},
               {"Friction_Coeff", friction_coeff},
               {"Kn_Factor", 1.0},
               {"Beta_n_Factor", beta_n_factor},
               {"K", 0.0}});
 
  // Contact pair 1-1 (bottom-bottom)
  pContactJson["Set_1_1"] = contact_base;
  pContactJson["Set_1_1"]["Kn"] = Kn_11;
 
  // Contact pair 1-2 (bottom-top)
  pContactJson["Set_1_2"] = contact_base;
  pContactJson["Set_1_2"]["Kn"] = Kn_12;
 
  // Contact pair 2-2 (top-top)
  pContactJson["Set_2_2"] = contact_base;
  pContactJson["Set_2_2"]["Kn"] = Kn_22;
 
  pDeckJson["Contact"] = pContactJson;
 
  // Neighbor settings
  pDeckJson["Neighbor"] = inp::PNeighborDeck::getExampleJson(
          json{{"Update_Criteria", "simple_all"},
               {"Search_Factor", 10.0},
               {"Search_Interval", 40},
               {"Near_Bd_Nodes_Tol", 0.5}});
 
  // Particle generation settings
  auto pGenJson = inp::PGenDeck::getExampleJson(json{{"Method", "From_File"}});
 
  // Add data that will be used to create particles
  pGenJson["Data"]["N"] = 2;  // two particles
 
  // Bottom particle
  pGenJson["Data"]["0"] = {
      {"x", R1}, {"y", R1}, {"z", 0.0},
      {"theta", 0.0}, {"s", 1.0},
      {"geom_id", 0}, {"mat_id", 0}, {"contact_id", 0}
  };
 
  // Top particle
  pGenJson["Data"]["1"] = {
      {"x", R1}, {"y", 2.0 * R1 + R2 + particle_dist}, {"z", 0.0},
      {"theta", M_PI * 0.5}, {"s", 1.0},
      {"geom_id", 1}, {"mat_id", 1}, {"contact_id", 1}
  };
 
  // Add to json
  pDeckJson["Particle_Generation"] = pGenJson;
 
  // Collect all decks into global input JSON
  auto inputJson = json({{"Model", modelDeckJson}, 
            {"Output", outputDeckJson}, 
            {"Force_BC", bcDeckJson["Force_BC"]},
            {"Displacement_BC", bcDeckJson["Displacement_BC"]},
            {"IC", bcDeckJson["IC"]}, 
            {"Particle", pDeckJson["Particle"]}, 
            {"Mesh", pDeckJson["Mesh"]}, 
            {"Material", pDeckJson["Material"]}, 
            {"Contact", pDeckJson["Contact"]},
            {"Neighbor", pDeckJson["Neighbor"]}, 
            {"Particle_Generation", pDeckJson["Particle_Generation"]}});
 
  // Write configuration to file
  std::cout << "\n\nPrinting global input deck json:\n";
  std::cout << inputJson.dump(2) << std::endl;
  
  // save to file
  std::ofstream f(input_dir + "/input.json");
  f << inputJson.dump(2);
  f.close();
 
  return inputJson;
}