Apply ground motion input as nodal acceleration constraint

[thiagordonho]

Describe the PR
This implementation is based on the issue presented in the RFC #692 . It adds a feature to read and apply input ground motion in terms of acceleration-time history.

Related Issues/PRs
Continuation of PR #701.

Additional context
This design involves creating a class for acceleration constraints (AccelerationConstraint) based on the already existing class for velocity constraints (VelocityConstraint). However, a math function can be associated with the acceleration constraint defined from an acceleration-time history input. The acceleration constraints are applied at the nodes with the implemented Node::apply_acceleration_constraint function within the Node::compute_acceleration_velocity function as follows:

//! Compute acceleration and velocity
template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
bool mpm::Node<Tdim, Tdof, Tnphases>::compute_acceleration_velocity(
    unsigned phase, double dt) noexcept {
  bool status = false;
  const double tolerance = 1.0E-15;
  if (mass_(phase) > tolerance) {
    // acceleration = (unbalaced force / mass)
    this->acceleration_.col(phase) =
        (this->external_force_.col(phase) + this->internal_force_.col(phase)) /
        this->mass_(phase);

    // Apply friction constraints
    this->apply_friction_constraints(dt);

    // Apply acceleration constraints
    this->apply_acceleration_constraints();

    // Velocity += acceleration * dt
    this->velocity_.col(phase) += this->acceleration_.col(phase) * dt;
    // Apply velocity constraints, which also sets acceleration to 0,
    // when velocity is set.
    this->apply_velocity_constraints();

    // Set a threshold
    for (unsigned i = 0; i < Tdim; ++i)
      if (std::abs(velocity_.col(phase)(i)) < tolerance)
        velocity_.col(phase)(i) = 0.;
    for (unsigned i = 0; i < Tdim; ++i)
      if (std::abs(acceleration_.col(phase)(i)) < tolerance)
        acceleration_.col(phase)(i) = 0.;
    status = true;
  }
  return status;
}

The acceleration-time input, as mentioned, is done with math functions. An additional way to initialize a ”Linear” math function is to read a .csv file with the time and respective acceleration entries. This procedure is herein extended from the one presented in PR #701, where the CSVReader header file is used for reading the .csv file. In this implementation, reading the math function .csv file was moved to the io/io_mesh_ascii.tcc within the newly created function called IOMeshAscii::read_math_functions:

// Return array with math function entries
template <unsigned Tdim>
std::array<std::vector<double>, 2> mpm::IOMeshAscii<Tdim>::read_math_functions(
    const std::string& math_file) {
  // Initialise vector with 2 empty vectors
  std::array<std::vector<double>, 2> xfx_values;

  // Read from csv file
  try {
    io::CSVReader<2> in(math_file);
    double x_value, fx_value;
    while (in.read_row(x_value, fx_value)) {
      xfx_values[0].push_back(x_value);
      xfx_values[1].push_back(fx_value);
    }
  } catch (std::exception& exception) {
    console_->error("Read math functions: {}", exception.what());
  }

  return xfx_values;
}

Benchmark tests

To test the implemented feature, a time-dependent acceleration was imposed to a Linear elastic, 2D block. The acceleration was imposed on the nodes at the bottom boundary of a 3 x 2 (width x height) mesh. The mesh has a spacing of 0.25 m. The block is 1 meter high and 1 meter wide. A detail of the model is shown in Figure 1. No gravity was considered in this model. Four particles per cell and direction were considered in this model. The following are the elastic characteristics of the block:

• density: 3000 kg/m3
• Young's modulus; 7 MPa
• Poisson's ratio: 0.38

Figure 1; Schematic of the model in its starting position.

Two acceleration-time series were considered for two separate simulations. The charts below show the imposed nodal acceleration at the boundary (Figure 2) and the results of the center of mass of the block (Figure 3 and 4).

Figure 2: Imposed acceleration-time history at the nodes of the bottom boundary.

Figure 3: Horizontal velocity of the center of mass of the block.

Figure 4: Horizontal position of the center of mass of the block.

The velocity of the center of mass shows a slight deviation of the expected results, but the horizontal position of the center of mass is very close to the expected values.

[kks32]

Suggested change

	// ID
	// Particle and cell ids

[kks32]

inline this function

[kks32]

Is there a general constraint class from which this could be derived?

[thiagordonho]

Currently, no. But I do agree that we should implement a general constraint from which this one and others (VelocityConstraint, FrictionConstraint, Traction) should be derived. They all have a setid_, dir_ and their main value (acceleration_, velocity_, traction_, or friction_). They only differ on the way these main values are named and if they can or cannot be time-dependent -- i.e. if they can have a pointer to a math function. Currently, only traction and acceleration constraints can be time-dependent.

[kks32]

@thiagordonho is there an elastic wave propagation benchmark that demonstrates this implementation?

[thiagordonho]

@thiagordonho is there an elastic wave propagation benchmark that demonstrates this implementation?

@kks32 I ran two simple simulations of a linear elastic block, no gravity, positioned on top of the bottom boundary. I applied the acceleration-time constraint to the nodes at the bottom boundary. The first simulation considers a sinusoidal math function an the second considers a piecewise linear math function. I will attach the obtained results to the body of the PR description.

[kks32]

Do you interpolate between timesteps? Is that linear?

[kks32]

n is usually reserved for number not nodal, use nodal_acceleration

[codecov]

Codecov Report

Merging #715 (34c848d) into develop (86ba10e) will decrease coverage by 0.05%.
The diff coverage is 95.02%.

@@             Coverage Diff             @@
##           develop     #715      +/-   ##
===========================================
- Coverage    96.78%   96.73%   -0.05%     
===========================================
  Files          130      132       +2     
  Lines        25899    26456     +557     
===========================================
+ Hits         25065    25591     +526     
- Misses         834      865      +31     

Impacted Files	Coverage Δ
include/io/io_mesh.h	100.00% <ø> (ø)
include/io/io_mesh_ascii.h	100.00% <ø> (ø)
include/loads_bcs/constraints.h	100.00% <ø> (ø)
include/mesh.h	100.00% <ø> (ø)
include/node.h	100.00% <ø> (ø)
include/node_base.h	100.00% <ø> (ø)
include/solvers/mpm_base.h	50.00% <ø> (ø)
include/solvers/mpm_base.tcc	74.39% <37.21%> (-2.65%)	⬇️
include/io/io_mesh_ascii.tcc	76.33% <90.00%> (+1.52%)	⬆️
include/loads_bcs/constraints.tcc	96.39% <96.30%> (-0.04%)	⬇️
... and 13 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 86ba10e...34c848d. Read the comment docs.