#include <periDEMModel.h>

Inheritance diagram for PeriDEMModel:

Collaboration diagram for PeriDEMModel:

Public Member Functions
	PeriDEMModel (std::shared_ptr< inp::Input > &deck, std::string modelName="PeriDEMModel")

void	log (std::ostringstream &oss, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)

void	log (const std::string &str, int priority=0, bool check_condition=true, int override_priority=-1, bool screen_out=false)

void	run (std::shared_ptr< inp::Input > &deck)

void	restart (std::shared_ptr< inp::Input > &deck)

void	init ()

void	close ()

void	integrate ()

void	integrateStep ()

void	computeForces ()

void	computeExternalForces ()

void	applyDisplacementBC ()

void	applyInitialCondition ()

void	output ()

std::string	ppTwoParticleTest ()

std::string	ppCompressiveTest ()

void	checkStop ()

void	setContact (std::unique_ptr< contact::Contact > c)

void	setPostprocess (std::unique_ptr< postprocess::Postprocess > p)

Public Member Functions inherited from data::ModelData
	ModelData (std::shared_ptr< inp::Input > &deck)
	Constructor.

const particle::BaseParticle *	getParticleFromAllList (size_t i) const
	Get pointer to base particle.

particle::BaseParticle *&	getParticleFromAllList (size_t i)
	Get pointer to base particle.

const particle::BaseParticle *	getParticleFromParticleList (size_t i) const
	Get pointer to particle (excluding wall)

particle::BaseParticle *&	getParticleFromParticleList (size_t i)

const particle::BaseParticle *	getParticleFromWallList (size_t i) const
	Get pointer to wall.

particle::BaseParticle *&	getParticleFromWallList (size_t i)

double	getDensity (size_t i)
	Get density of particle.

double	getHorizon (size_t i)
	Get horizon of particle.

size_t &	getPtId (size_t i)
	Get particle id given the location in particle list.

const size_t &	getPtId (size_t i) const
	Get particle id given the location in particle list.

void	setPtId (size_t i, const size_t &id)
	Set particle id given the location in particle list.

double	getKeyData (std::string key, bool issue_err=false)
	Get data for a key.

void	appendKeyData (std::string key, double data, bool issue_err=false)
	Append value to data associated with key.

void	setKeyData (std::string key, double data, bool issue_err=false)
	Set value to data associated with key.

int	dimension () const

double	timeStep () const

double	currentDt () const

void	setCurrentDt (double dt)

void	advanceTime ()

size_t	numPdForceNodes () const

size_t	pdForceNode (size_t k) const

bool	isDofFree (size_t i, unsigned int dof) const

double	getVMag (size_t i) const

void	setVMag (size_t i, double mag)

size_t	currentStep () const

size_t	numTimeSteps () const

const std::string &	timeDiscretization () const

bool	performOutput () const

bool	shouldOutput () const

util::Point &	getXRef (size_t i)
	Get reference coordinate of the node.

const util::Point &	getXRef (size_t i) const
	Get reference coordinate of the node.

void	setXRef (size_t i, const util::Point &x)
	Set reference coordinate of the node.

void	addXRef (size_t i, const util::Point &x)
	Add reference coordinate of the node.

void	setXRef (size_t i, int dof, double x)
	Set specific reference coordinate of the node.

void	addXRef (size_t i, int dof, double x)
	Add specific reference coordinate of the node.

util::Point &	getX (size_t i)
	Get current coordinate of the node.

const util::Point &	getX (size_t i) const
	Get current coordinate of the node.

void	setX (size_t i, const util::Point &x)
	Set current coordinate of the node.

void	addX (size_t i, const util::Point &x)
	Add current coordinate of the node.

void	setX (size_t i, int dof, double x)
	Set specific current coordinate of the node.

void	addX (size_t i, int dof, double x)
	Add to specific current coordinate of the node.

util::Point &	getU (size_t i)
	Get displacement of the node.

const util::Point &	getU (size_t i) const
	Get displacement of the node.

void	setU (size_t i, const util::Point &u)
	Set displacement of the node.

void	addU (size_t i, const util::Point &u)
	Add to displacement of the node.

void	setU (size_t i, int dof, double u)
	Set displacement of the node.

void	addU (size_t i, int dof, double u)
	Add to displacement of the node.

util::Point &	getV (size_t i)
	Get velocity of the node.

const util::Point &	getV (size_t i) const
	Get velocity of the node.

void	setV (size_t i, const util::Point &v)
	Set velocity of the node.

void	addV (size_t i, const util::Point &v)
	Add to velocity of the node.

void	setV (size_t i, int dof, double v)
	Set velocity of the node.

void	addV (size_t i, int dof, double v)
	Add to velocity of the node.

util::Point &	getF (size_t i)
	Get force of the node.

const util::Point &	getF (size_t i) const
	Get force of the node.

void	setF (size_t i, const util::Point &f)
	Set force of the node.

void	addF (size_t i, const util::Point &f)
	Add to force of the node.

void	setF (size_t i, int dof, double f)
	Set force of the node.

void	addF (size_t i, int dof, double f)
	Add to force of the node.

double &	getVol (size_t i)
	Get volume of the node.

const double &	getVol (size_t i) const
	Get volume of the node.

void	setVol (size_t i, const double &vol)
	Set volume of the node.

void	addVol (size_t i, const double &vol)
	Add to volume of the node.

uint8_t &	getFix (size_t i)
	Get fixity of the node.

const uint8_t &	getFix (size_t i) const
	Get fixity of the node.

void	setFix (size_t i, const unsigned int &dof, const bool &flag)
	Set fixity of the node.

double &	getMx (size_t i)
	Get weighted-volume (mx) of the node.

const double &	getMx (size_t i) const
	Get weighted-volume (mx) of the node.

void	setMx (size_t i, const double &mx)
	Set weighted-volume (mx) of the node.

void	addMx (size_t i, const double &mx)
	Add to weighted-volume (mx) of the node.

double &	getThetax (size_t i)
	Get volumetric deformation (thetax) of the node.

const double &	getThetax (size_t i) const
	Get volumetric deformation (thetax) of the node.

void	setThetax (size_t i, const double &thetax)
	Set volumetric deformation (thetax) of the node.

void	addThetax (size_t i, const double &thetax)
	Add to volumetric deformation (thetax) of the node.

Data Fields
std::unique_ptr< contact::Contact >	d_contact_p

std::unique_ptr< postprocess::Postprocess >	d_postprocess_p

Data Fields inherited from data::ModelData
size_t	d_n
	Current time step.

double	d_time
	Current time.

double	d_currentDt
	Current timestep.

size_t	d_infoN
	Print log step interval.

std::string	d_name
	Model name for log prefixes.

std::vector< std::pair< double, std::string > >	d_pvdParticleEntries
	(simulation time, VTU filename) for output.pvd

std::vector< std::pair< double, std::string > >	d_pvdStrainEntries
	(simulation time, VTU filename) for output_strain.pvd

std::map< std::string, double >	d_dbgData
	Debug data.

std::ofstream	d_ppFile
	File stream to output information.

std::shared_ptr< inp::Input >	d_input_p
	Pointer to Input object.

std::shared_ptr< inp::ModelDeck >	d_modelDeck_p
	Pointer to deck holding problem related data.

std::shared_ptr< inp::OutputDeck >	d_outputDeck_p
	Pointer to deck holding output related data.

std::shared_ptr< inp::RestartDeck >	d_restartDeck_p
	Pointer to deck holding restart related data.

std::shared_ptr< inp::TestDeck >	d_testDeck_p
	Test deck.

std::shared_ptr< inp::BCDeck >	d_bcDeck_p
	Boundary condition deck.

std::shared_ptr< inp::ParticleDeck >	d_particleDeck_p
	Pointer to deck holding particle related data.

bool	d_stop
	flag to stop the simulation midway

double	d_hMax
	Minimum mesh over all particles and walls.

double	d_hMin
	Maximum mesh over all particles and walls.

double	d_maxContactR
	Maximum contact radius between over pairs of particles and walls.

size_t	d_contNeighUpdateInterval
	Neighborlist update interval.

size_t	d_contNeighTimestepCounter
	Contact neighborlist time step counter.

double	d_contNeighSearchRadius
	Neighborlist contact search radius (multiple of d_maxContactR). This variable will be updated during simulation based on maximum velocity.

std::vector< char >	d_mpiIncludeInContactCloud
	Per entry in d_particlesListTypeAll: include nodes in the MPI contact search cloud (owned grains, walls, and distance-limited ghosts). Empty / all-true when not using multi-rank particle parallel.

std::vector< std::vector< int > >	d_mpiGhostNeedFrom
	Cached ghost grain ids requested from each MPI rank (skin rebuild).

std::vector< std::vector< int > >	d_mpiGhostServeTo
	Cached grain ids this rank must send to each MPI rank.

bool	d_mpiGhostPlanValid
	True once a ghost communication plan has been built.

size_t	d_mpiGhostStepsSinceRebuild
	Force steps since last ghost-plan rebuild.

bool	d_pdDofMpi
	Nodal DOF-MPI active.

bool	d_pdGrainAligned
	DOF owners match whole-grain brick (same as particle-MPI). When true, every node of a grain lives on one rank — contact ghosts are the particle ghost plan, not a full nodal Allreduce. When false, Metis may split a grain across ranks — need syncAllOwned (or a PD halo large enough for contact) before contact reads remote u/v.

std::vector< size_t >	d_pdNodePartition
	Owner rank for each node when d_pdDofMpi is true.

std::vector< std::vector< int > >	d_pdGhostNeedFrom
	Ghost node ids requested from each rank (PD DOF halo).

std::vector< std::vector< int > >	d_pdGhostServeTo
	Local node ids this rank must send to each rank.

std::vector< std::shared_ptr< particle::RefParticle > >	d_referenceParticles
	Pointer to reference particle.

std::vector< particle::BaseParticle * >	d_particlesListTypeAll
	List of particles + walls.

std::vector< particle::BaseParticle * >	d_particlesListTypeParticle
	List of particles.

std::vector< particle::BaseParticle * >	d_particlesListTypeWall
	List of walls.

std::vector< inp::MatData >	d_particlesMatDataList
	List of particle material data. Only populated if needed for calculation of stress or other quantities.

std::vector< double >	d_maxVelocityParticlesListTypeAll
	Maximum velocity among all nodes in the particle for each particle.

double	d_maxVelocity
	Maximum velocity among all nodes.

std::vector< std::vector< size_t > >	d_zInfo
	Zone information of particles. For zone 0, d_zInfo[0] = [n1, n2] where n1 is the index at which the particle in this zone starts in d_particlesListTypeAll and n2 is the index + 1, where index is at which the particle in this zone ends.

std::unique_ptr< loading::ParticleULoading >	d_uLoading_p
	Pointer to displacement Loading object.

std::unique_ptr< loading::ParticleFLoading >	d_fLoading_p
	Pointer to force Loading object.

std::unique_ptr< geometry::Fracture >	d_fracture_p
	Fracture state of bonds.

std::unique_ptr< NSearch >	d_nsearch_p
	Pointer to nsearch.

std::vector< util::Point >	d_xRef
	reference positions of the nodes

std::vector< util::Point >	d_x
	Current positions of the nodes.

std::vector< util::Point >	d_u
	Displacement of the nodes.

std::vector< util::Point >	d_v
	Velocity of the nodes.

std::vector< double >	d_vMag
	Magnitude of velocity of the nodes.

std::vector< util::Point >	d_f
	Total force on the nodes.

std::vector< double >	d_vol
	Nodal volumes.

std::vector< size_t >	d_ptId
	Global node to particle id (walls are assigned id after last particle id)

std::vector< std::vector< size_t > >	d_neighC
	Neighbor data for contact forces.

std::vector< std::vector< size_t > >	d_neighPd
	Neighbor data for peridynamic forces.

std::vector< std::vector< float > >	d_neighPdSqdDist
	Square distance neighbor data for peridynamic forces.

std::vector< std::vector< std::vector< size_t > > >	d_neighWallNodes
	Neighbor data for contact between particle and walls.

std::vector< std::vector< std::vector< double > > >	d_neighWallNodesDistance
	Neighbor data (distance) for contact between particle and walls.

std::vector< std::vector< size_t > >	d_neighWallNodesCondensed
	Neighbor data for contact between particle and walls condensed into single vector for each particle.

std::vector< uint8_t >	d_fix
	Vector of fixity mask of each node.

std::vector< uint8_t >	d_forceFixity
	Vector of fixity mask of each node for force.

std::vector< double >	d_thetaX
	Dilation.

std::vector< double >	d_mX
	Weighted volume.

std::vector< size_t >	d_fPdCompNodes
	List of global nodes on which force (peridynamic/internal) is to be computed.

std::vector< size_t >	d_fContCompNodes
	List of global nodes on which force (contact) is to be computed.

std::vector< float >	d_Z
	Damage at nodes.

std::vector< float >	d_e
	Energy of the nodes.

std::vector< float >	d_w
	Work done on each of the nodes.

std::vector< float >	d_phi
	Damage function \( \phi \) at the nodes (volume-weighted, Silling 2000)

std::vector< float >	d_phiBond
	Damage as broken-bond count fraction (Bhattacharya & Lipton 2023)

std::vector< float >	d_eF
	Fracture energy of the nodes.

std::vector< float >	d_eFB
	Bond-based fracture energy of the nodes.

std::vector< util::Point >	d_xQuadCur
	Current position of quadrature points.

std::vector< util::SymMatrix3 >	d_strain
	Strain in elements (values at quadrature points)

std::vector< util::SymMatrix3 >	d_stress
	Stress in elements (values at quadrature points)

float	d_te
	Total internal energy.

float	d_tw
	Total work done.

float	d_tk
	Total kinetic energy.

float	d_teF
	Total fracture energy.

float	d_teFB
	Total bond-based fracture energy.

Detailed Description

Default driver. Hands the run to time_int::Integrator at integrate(). Physics lives here: forces, BCs, output, particles.

Definition at line 32 of file periDEMModel.h.

Constructor & Destructor Documentation

◆ PeriDEMModel()

PeriDEMModel::PeriDEMModel	(	std::shared_ptr< inp::Input > &	deck,
		std::string	modelName = `"PeriDEMModel"`
	)

explicit

Definition at line 58 of file periDEMModel.cpp.

  : data::ModelData(deck) {
 
  d_name = std::move(modelName);
 
  // Ensure output directory exists (Path from input; relative to cwd).
  {
    namespace fs = std::filesystem;
    fs::path out(d_outputDeck_p->d_path);
    if (!out.empty())
      fs::create_directories(out);
  }
 
  // initialize logger
  util::io::initLogger(d_outputDeck_p->d_debug,
                       d_outputDeck_p->d_path + "log.txt");
  if (!d_postprocess_p)
    d_postprocess_p = std::make_unique<postprocess::Postprocess>();
}

References data::ModelData::d_name, data::ModelData::d_outputDeck_p, d_postprocess_p, and util::io::initLogger().

Here is the call graph for this function:

Member Function Documentation

◆ applyDisplacementBC()

void PeriDEMModel::applyDisplacementBC ( )

Definition at line 512 of file periDEMModel.cpp.

                                       {
  log("    Computing external displacement bc \n", 3);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->apply(d_time, p); // applied in parallel
}

References data::ModelData::d_particlesListTypeAll, data::ModelData::d_time, data::ModelData::d_uLoading_p, and log().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ applyInitialCondition()

void PeriDEMModel::applyInitialCondition ( )

Definition at line 518 of file periDEMModel.cpp.

                                         {
  log("Applying initial condition \n", 3);
  for (auto &p : d_particlesListTypeAll)
    loading::applyIC(p, d_bcDeck_p->d_icDeck); // applied in parallel
}

References loading::applyIC(), data::ModelData::d_bcDeck_p, data::ModelData::d_particlesListTypeAll, and log().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ checkStop()

void PeriDEMModel::checkStop ( )

Definition at line 537 of file periDEMModel.cpp.

                             {
  if (d_testDeck_p->d_testName == "two_particle")
    log(ppTwoParticleTest(), 2, d_n % d_infoN == 0, 3);
  else if (d_testDeck_p->d_testName == "compressive_test")
    log(ppCompressiveTest(), 2, d_n % d_infoN == 0, 3);
  if (d_postprocess_p)
    d_postprocess_p->checkStop(*this);
}

References data::ModelData::d_infoN, data::ModelData::d_n, d_postprocess_p, data::ModelData::d_testDeck_p, log(), ppCompressiveTest(), and ppTwoParticleTest().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ close()

void PeriDEMModel::close ( )

Definition at line 124 of file periDEMModel.cpp.

                         {
  if (d_postprocess_p)
    d_postprocess_p->close(*this);
}

References d_postprocess_p.

Referenced by main(), run(), and anonymous_namespace{main.cpp}::runLaw().

Here is the caller graph for this function:

◆ computeExternalForces()

void PeriDEMModel::computeExternalForces ( )

Definition at line 476 of file periDEMModel.cpp.

                                         {
  log("    Computing external force \n", 3);
 
  auto gravity = d_bcDeck_p->d_gravity;
 
  if (gravity.length() > 1.0E-8) {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    // Only owned force nodes (walls on rank 0 for particle-MPI; Metis owners
    // for DOF-MPI). Applying gravity on every rank then Allreducing reaction
    // would multiply wall weight by mpiSize.
    const auto &nodes = d_fContCompNodes;
    taskflow.for_each_index((std::size_t) 0, nodes.size(), (std::size_t)1,
                            [this, gravity, &nodes](std::size_t II) {
                              const size_t i = nodes[II];
                              this->d_f[i] += this->getDensity(i) * gravity;
                            });
 
    executor.run(taskflow).get();
  }
 
  for (auto &p : d_particlesListTypeAll) {
    if (p->isWall()) {
      if (util::parallel::mpiRank() == 0)
        d_fLoading_p->apply(d_time, p);
    } else if (d_pdDofMpi) {
      // Nodal ownership: apply on every rank; loading loops all nodes of p,
      // but only owned force nodes are integrated.
      d_fLoading_p->apply(d_time, p);
    } else if (particle::isLocallyOwned(*p)) {
      d_fLoading_p->apply(d_time, p);
    }
  }
}

References data::ModelData::d_bcDeck_p, data::ModelData::d_fContCompNodes, data::ModelData::d_fLoading_p, data::ModelData::d_particlesListTypeAll, data::ModelData::d_pdDofMpi, data::ModelData::d_time, util::parallel::getNThreads(), particle::isLocallyOwned(), log(), and util::parallel::mpiRank().

Referenced by computeForces().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ computeForces()

void PeriDEMModel::computeForces ( )

Definition at line 322 of file periDEMModel.cpp.

                                 {
 
  bool dbg_condition = d_n % d_infoN == 0;
 
  log("  Compute forces \n", 2, dbg_condition, 3);
 
  // Refresh ghost grain kinematics from owners before contact / damping.
  {
    const auto t_ex0 = steady_clock::now();
    particle::exchangeGhostKinematics(*this);
    appendKeyData("mpi_exchange_wall_time",
                  util::methods::timeDiff(t_ex0, steady_clock::now()));
  }
 
  // Nodal DOF halo for PD neighbor reads (Single_Particle MPI).
  pd::exchangeGhostDisplacement(*this);
 
  // reset force
  auto t1 = steady_clock::now();
  float force_reset_time = 0;
  {
    tf::Executor executor(util::parallel::getNThreads());
    tf::Taskflow taskflow;
 
    taskflow.for_each_index(
      (std::size_t) 0, d_x.size(), (std::size_t) 1,
        [this](std::size_t i) { this->d_f[i] = util::Point(); }
    ); // for_each
 
    executor.run(taskflow).get();
    force_reset_time = util::methods::timeDiff(t1, steady_clock::now());
  }
 
  // compute peridynamic forces
  t1 = steady_clock::now();
  pd::computeForces(*this);
  auto pd_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("pd_compute_time", pd_time);
  appendKeyData("avg_peridynamics_force_time", pd_time/d_infoN);
 
  float current_contact_neigh_update_time = 0;
  float contact_time = 0;
  if (d_input_p->isMultiParticle() && d_contact_p) {
    // update contact neighborlist
    t1 = steady_clock::now();
    d_contact_p->updateNeighborlist(*this);
    current_contact_neigh_update_time = util::methods::timeDiff(t1,
                                                                     steady_clock::now());
    appendKeyData("contact_neigh_update_time",
                  current_contact_neigh_update_time);
    appendKeyData("avg_contact_neigh_update_time",
                  current_contact_neigh_update_time / d_infoN);
 
    // compute contact forces between particles
    t1 = steady_clock::now();
    d_contact_p->computeForces(*this);
    contact_time = util::methods::timeDiff(t1, steady_clock::now());
    appendKeyData("contact_compute_time", contact_time);
    appendKeyData("avg_contact_force_time", contact_time / d_infoN);
  }
 
  // Compute external forces
  t1 = steady_clock::now();
  computeExternalForces();
  auto extf_time = util::methods::timeDiff(t1, steady_clock::now());
  appendKeyData("extf_compute_time", extf_time);
  appendKeyData("avg_extf_compute_time", extf_time/d_infoN);
 
  // output avg time info
  if (dbg_condition) {
    if (d_input_p->isMultiParticle()) {
      log(std::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "tree update", size_t(getKeyData("avg_tree_update_time")),
                      "contact neigh update",
                      size_t(getKeyData("avg_contact_neigh_update_time")),
                      "contact force",
                      size_t(getKeyData("avg_contact_force_time")),
                      "total contact", size_t(getKeyData("avg_tree_update_time")
                                              + getKeyData(
                          "avg_contact_neigh_update_time")
                                              + getKeyData(
                          "avg_contact_force_time")),
                      "peridynamics force",
                      size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force",
                      size_t(getKeyData("avg_extf_compute_time") / d_infoN)),
          2, dbg_condition, 3);
 
      appendKeyData("avg_tree_update_time", 0.);
      appendKeyData("avg_contact_neigh_update_time", 0.);
      appendKeyData("avg_contact_force_time", 0.);
      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
    else {
      log(std::format("    Avg time (ms): \n"
                      "      {:48s} = {:8d}\n"
                      "      {:48s} = {:8d}\n",
                      "peridynamics force", size_t(getKeyData("avg_peridynamics_force_time")),
                      "external force", size_t(getKeyData("avg_extf_compute_time")/d_infoN)),
          2, dbg_condition, 3);
 
      appendKeyData("avg_peridynamics_force_time", 0.);
      appendKeyData("avg_extf_compute_time", 0.);
    }
  }
 
  log(std::format("    {:50s} = {:8d} \n",
                  "Force reset time (ms)",
                  size_t(force_reset_time)
      ),
      2, dbg_condition, 3);
 
  log(std::format("    {:50s} = {:8d} \n",
                  "External force time (ms)",
                  size_t(extf_time)
      ),
      2, dbg_condition, 3);
 
  log(std::format("    {:50s} = {:8d} \n",
                  "Peridynamics force time (ms)",
                  size_t(pd_time)
      ),
      2, dbg_condition, 3);
 
  if (d_input_p->isMultiParticle()) {
 
    log(std::format("    {:50s} = {:8d} \n",
                    "Point cloud update time (ms)",
                    size_t(getKeyData("pt_cloud_update_time"))
        ),
        2, dbg_condition, 3);
 
    log(std::format("    {:50s} = {:8d} \n",
                    "Contact neighborlist update time (ms)",
                    size_t(current_contact_neigh_update_time)
        ),
        2, dbg_condition, 3);
 
    log(std::format("    {:50s} = {:8d} \n",
                    "Contact force time (ms)",
                    size_t(contact_time)
        ),
        2, dbg_condition, 3);
  }
 
}

References data::ModelData::appendKeyData(), computeExternalForces(), pd::computeForces(), d_contact_p, data::ModelData::d_infoN, data::ModelData::d_input_p, data::ModelData::d_n, data::ModelData::d_x, pd::exchangeGhostDisplacement(), particle::exchangeGhostKinematics(), data::ModelData::getKeyData(), util::parallel::getNThreads(), log(), and util::methods::timeDiff().

Referenced by main(), and anonymous_namespace{main.cpp}::runLaw().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ init()

void PeriDEMModel::init ( )

Definition at line 129 of file periDEMModel.cpp.

                        {
 
  // init time step
  d_n = 0;
  d_time = 0.;
  if (d_outputDeck_p->d_dtOut < 1)
    throw std::runtime_error(
        "Output_Interval must be >= 1 (use 1 to write every step).");
  if (d_outputDeck_p->d_dtOutCriteria < 1)
    d_outputDeck_p->d_dtOutCriteria = d_outputDeck_p->d_dtOut;
  if (d_outputDeck_p->d_dtOutOld < 1)
    d_outputDeck_p->d_dtOutOld = d_outputDeck_p->d_dtOut;
  if (d_outputDeck_p->d_dtTestOut == 0)
    d_outputDeck_p->d_dtTestOut = std::max<size_t>(1, d_outputDeck_p->d_dtOut / 10);
  d_infoN = d_outputDeck_p->d_dtOut;
 
  // debug/information variables
  {
    appendKeyData("debug_once", -1);
    appendKeyData("update_contact_neigh_search_params_init_call_count", 0);
    appendKeyData("tree_compute_time", 0);
    appendKeyData("contact_compute_time", 0);
    appendKeyData("contact_neigh_update_time", 0);
    appendKeyData("peridynamics_neigh_update_time", 0);
    appendKeyData("pd_compute_time", 0);
    appendKeyData("extf_compute_time", 0);
    appendKeyData("integrate_compute_time", 0);
    appendKeyData("pt_cloud_update_time", 0);
    appendKeyData("avg_tree_update_time", 0);
    appendKeyData("avg_contact_neigh_update_time", 0);
    appendKeyData("avg_contact_force_time", 0);
    appendKeyData("avg_peridynamics_force_time", 0);
    appendKeyData("avg_extf_compute_time", 0);
    appendKeyData("pen_dist", 0);
    appendKeyData("max_y", 0);
    appendKeyData("contact_area_radius", 0);
  }
 
 
  auto t1 = steady_clock::now();
  auto t2 = steady_clock::now();
  log(d_name + ": Initializing objects.\n");
 
  // create particles
  log(d_name + ": Creating particles.\n");
  particle::createParticles(*this);
  particle::assignMpiOwners(*this);
 
  log(d_name + ": Creating maximum velocity data for particles.\n");
  d_maxVelocityParticlesListTypeAll
          = std::vector<double>(d_particlesListTypeAll.size(), 0.);
  d_maxVelocity = util::methods::max(d_maxVelocityParticlesListTypeAll);
 
  // setup contact
  if (d_input_p->isMultiParticle()) {
    if (!d_contact_p)
      d_contact_p = std::make_unique<contact::Contact>();
    log(d_name + ": Setting up contact.\n");
    d_contact_p->setup(*this);
  }
 
  // setup element-node connectivity data if needed
  log(d_name + ": Setting up element-node connectivity data for strain/stress.\n");
  data::setupQuadratureData(*this);
 
  // create search object
  log(d_name + ": Creating neighbor search tree.\n");
 
  // create tree object
  d_nsearch_p = std::make_unique<NSearch>(d_x, d_outputDeck_p->d_debug);
 
  // setup tree
  double set_tree_time = d_nsearch_p->setInputCloud();
  log(std::format("{}: Tree setup time (ms) = {}. \n", d_name, set_tree_time));
 
  // create neighborlists
  log(d_name + ": Creating neighborlist for peridynamics.\n");
  t1 = steady_clock::now();
  nsearch::updatePeridynamicNeighborlist(*this);
  t2 = steady_clock::now();
  appendKeyData("peridynamics_neigh_update_time", util::methods::timeDiff(t1, t2));
 
  if (d_input_p->isMultiParticle()) {
    log(d_name + ": Creating neighborlist for contact.\n");
    d_contNeighUpdateInterval = d_particleDeck_p->d_pNeighDeck.d_neighUpdateInterval;
    d_contNeighSearchRadius = d_particleDeck_p->d_pNeighDeck.d_sFactor * d_maxContactR;
    t1 = steady_clock::now();
    d_contact_p->updateNeighborlist(*this);
    t2 = steady_clock::now();
    appendKeyData("contact_neigh_update_time", util::methods::timeDiff(t1, t2));
  }
 
  // create peridynamic bonds
  log(d_name + ": Creating peridynamics bonds.\n");
  d_fracture_p = std::make_unique<geometry::Fracture>(&d_x, &d_neighPd);
 
  // compute quantities in state-based simulations
  log(d_name + ": Compute state-based peridynamic quantities.\n");
  material::computeStateMx(this, true);
 
  // initialize loading class
  log(d_name + ": Initializing displacement loading object.\n");
  d_uLoading_p =
      std::make_unique<loading::ParticleULoading>(d_bcDeck_p->d_dispDeck);
  for (auto &p : d_particlesListTypeAll)
    d_uLoading_p->setFixity(p);
 
  log(d_name + ": Initializing force loading object.\n");
  d_fLoading_p =
      std::make_unique<loading::ParticleFLoading>(d_bcDeck_p->d_forceDeck);
 
  // A particle with every degree of freedom fixed still has its force
  // computed, because the reaction force on it is wanted.
 
  // if this is a two-particle test, we set the force calculation off in
  // first particle
  if (d_testDeck_p->d_testName == "two_particle") {
    d_particlesListTypeAll[0]->d_computeForce = false;
  }
 
  log(std::format("{}: Total particles = {}. \n",
                  d_name, d_particlesListTypeAll.size()));
 
  for (const auto &p : d_particlesListTypeAll)
    if (!p->d_computeForce)
      log(std::format("{}: Force OFF in Particle i = {}. \n", d_name, p->getId()));
 
  // Single-particle DOF-MPI: Metis node owners + PD ghost plan (T10).
  pd::setupDofPartition(*this);
 
  log(d_name + ": Creating list of nodes on which force is to be computed.\n");
  // Selection is per particle or wall: a node enters d_fCompNodes if force is
  // computed on any node of its body. Selecting individual nodes is not
  // implemented.
  const int mpi_rank = util::parallel::mpiRank();
  for (size_t i = 0; i < d_x.size(); i++) {
    const auto &ptId = d_ptId[i];
    const auto &pi = getParticleFromAllList(ptId);
    bool node_owned = false;
    if (d_pdDofMpi)
      node_owned = (static_cast<int>(d_pdNodePartition[i]) == mpi_rank);
    else if (pi->isWall())
      // Walls are replicated; assemble wall contact/reaction on rank 0 only.
      node_owned = (mpi_rank == 0);
    else
      node_owned = particle::isLocallyOwned(*pi);
    if (pi->d_computeForce && node_owned) {
      d_fContCompNodes.push_back(i);
      // Walls keep contact (and reaction) but not peridynamic force. Treating a
      // container as a PD body on a thin/boolean mesh makes Damage_Z explode
      // and the neighbor search then allocates until the process is OOM-killed.
      if (!pi->isWall())
        d_fPdCompNodes.push_back(i);
    }
  }
 
  // initialize remaining fields (if any)
  d_Z = std::vector<float>(d_x.size(), 0.);
  // Damage: volume-weighted phi (Silling/Trask) and broken-bond count
  // fraction (Bhattacharya & Lipton).
  d_phi = std::vector<float>(d_x.size(), 0.);
  d_phiBond = std::vector<float>(d_x.size(), 0.);
 
  t2 = steady_clock::now();
  log(std::format("{}: Total setup time (ms) = {}. \n",
                  d_name, util::methods::timeDiff(t1, t2)));
 
  // compute complexity information
  size_t free_dofs = 0;
  for (const auto &f : d_fix) {
    for (size_t dof = 0; dof < 3; dof++)
      if (util::methods::isFree(f, dof))
        free_dofs++;
  }
  log(std::format("{}: Computational complexity information \n"
                  "  Total number of particles = {}, number of "
                  "particles = {}, number of walls = {}, \n"
                  "  number of dofs = {}, number of free dofs = {}. \n",
                  d_name, d_particlesListTypeAll.size(),
                  d_particlesListTypeParticle.size(),
                  d_particlesListTypeWall.size(),
                  3 * d_x.size(),
                  free_dofs));
}

Referenced by main(), run(), and anonymous_namespace{main.cpp}::runLaw().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ integrate()

void PeriDEMModel::integrate ( )

Definition at line 314 of file periDEMModel.cpp.

                             {
  time_int::Integrator().integrate(*this);
}

References time_int::Integrator::integrate().

Referenced by main(), and run().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ integrateStep()

void PeriDEMModel::integrateStep ( )

Definition at line 318 of file periDEMModel.cpp.

                                 {
  time_int::Integrator().step(*this);
}

References time_int::Integrator::step().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ log() [1/2]

void PeriDEMModel::log	(	const std::string &	str,
		int	priority = `0`,
		bool	check_condition = `true`,
		int	override_priority = `-1`,
		bool	screen_out = `false`
	)

Definition at line 83 of file periDEMModel.cpp.

                                           {
  util::io::log(priority, str, check_condition, override_priority, screen_out);
}

References util::io::log().

Here is the call graph for this function:

◆ log() [2/2]

void PeriDEMModel::log	(	std::ostringstream &	oss,
		int	priority = `0`,
		bool	check_condition = `true`,
		int	override_priority = `-1`,
		bool	screen_out = `false`
	)

Definition at line 78 of file periDEMModel.cpp.

                                           {
  util::io::log(priority, oss, check_condition, override_priority, screen_out);
}

References util::io::log().

Referenced by applyDisplacementBC(), applyInitialCondition(), checkStop(), computeExternalForces(), computeForces(), init(), output(), and restart().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ output()

void PeriDEMModel::output ( )

Definition at line 524 of file periDEMModel.cpp.

                          {
  if (currentStep() == 0)
    log(std::format("{}: Output step = {}, time = {:.6f} \n", d_name, d_n, d_time),
        2);
  rw::writeOutput(*this);
}

References data::ModelData::currentStep(), data::ModelData::d_n, data::ModelData::d_name, data::ModelData::d_time, log(), and rw::writeOutput().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ppCompressiveTest()

std::string PeriDEMModel::ppCompressiveTest ( )

Definition at line 546 of file periDEMModel.cpp.

                                          {
  if (!d_postprocess_p)
    return "";
  return d_postprocess_p->compressive(*this);
}

References d_postprocess_p.

Referenced by checkStop().

Here is the caller graph for this function:

◆ ppTwoParticleTest()

std::string PeriDEMModel::ppTwoParticleTest ( )

Definition at line 531 of file periDEMModel.cpp.

                                          {
  if (!d_postprocess_p)
    return "";
  return d_postprocess_p->twoParticle(*this);
}

References d_postprocess_p.

Referenced by checkStop().

Here is the caller graph for this function:

◆ restart()

void PeriDEMModel::restart ( std::shared_ptr< inp::Input > & deck )

Definition at line 104 of file periDEMModel.cpp.

                                                         {
 
  log(d_name + ": Restarting the simulation\n");
 
  // set time step to step specified in restart deck
  d_n = d_restartDeck_p->d_step;
  d_time = double(d_n) * d_modelDeck_p->d_dt;
  log(std::format("  Restart step = {}, time = {:.6f} \n", d_n, d_time));
 
  // get backup of reference configuration
  std::vector<util::Point> x_ref(d_x.size(), util::Point());
  for (auto &x : d_x)
    x_ref.push_back(x);
 
  // read displacement and velocity from restart file
  log("  Reading data from restart file = " + d_restartDeck_p->d_file + " \n");
  auto reader = rw::reader::VtkParticleReader(d_restartDeck_p->d_file);
  reader.readNodes(this);
}

References data::ModelData::d_modelDeck_p, data::ModelData::d_n, data::ModelData::d_name, data::ModelData::d_restartDeck_p, data::ModelData::d_time, data::ModelData::d_x, and log().

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ run()

void PeriDEMModel::run ( std::shared_ptr< inp::Input > & deck )

Definition at line 88 of file periDEMModel.cpp.

                                                     {
 
  // initialize data
  init();
 
  // check for restart
  if (d_modelDeck_p->d_isRestartActive)
    restart(deck);
 
  // integrate in time
  integrate();
 
  // close
  close();
}

References close(), data::ModelData::d_modelDeck_p, init(), integrate(), and restart().

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ setContact()

void PeriDEMModel::setContact ( std::unique_ptr< contact::Contact > c )

inline

Definition at line 69 of file periDEMModel.h.

69{ d_contact_p = std::move(c); }

References d_contact_p.

◆ setPostprocess()

void PeriDEMModel::setPostprocess ( std::unique_ptr< postprocess::Postprocess > p )

inline

Definition at line 71 of file periDEMModel.h.

                                                               {
    d_postprocess_p = std::move(p);
  }

References d_postprocess_p.

Referenced by main().

Here is the caller graph for this function:

Field Documentation

◆ d_contact_p

std::unique_ptr<contact::Contact> PeriDEMModel::d_contact_p

Definition at line 75 of file periDEMModel.h.

Referenced by anonymous_namespace{main.cpp}::checkWallSpringOnly(), computeForces(), init(), main(), and setContact().

◆ d_postprocess_p

std::unique_ptr<postprocess::Postprocess> PeriDEMModel::d_postprocess_p

Definition at line 76 of file periDEMModel.h.

Referenced by checkStop(), close(), PeriDEMModel(), ppCompressiveTest(), ppTwoParticleTest(), and setPostprocess().

The documentation for this class was generated from the following files:

/home/user/project/PeriDEM/periDEMModel.h
/home/user/project/PeriDEM/periDEMModel.cpp

Public Member Functions

Data Fields

Detailed Description

Constructor & Destructor Documentation

◆ PeriDEMModel()

Member Function Documentation

◆ applyDisplacementBC()

◆ applyInitialCondition()

◆ checkStop()

◆ close()

◆ computeExternalForces()

◆ computeForces()

◆ init()

◆ integrate()

◆ integrateStep()

◆ log() [1/2]

◆ log() [2/2]

◆ output()

◆ ppCompressiveTest()

◆ ppTwoParticleTest()

◆ restart()

◆ run()

◆ setContact()

◆ setPostprocess()

Field Documentation

◆ d_contact_p

◆ d_postprocess_p