Composable Solid Mechanics Demo¶
This demo builds one composable solid mechanics system. It registers fields, adds a Young's modulus parameter, applies a Neo-Hookean material, advances a dynamic solve, checks sensitivities, and writes output.
The full source code lives in examples/solid_mechanics/composable_solid_mechanics.cpp.
Includes and Initialization¶
#include "smith/infrastructure/application_manager.hpp"
#include "smith/physics/state/state_manager.hpp"
#include "smith/physics/mesh.hpp"
#include "smith/numerics/solver_config.hpp"
#include "smith/physics/functional_objective.hpp"
#include "smith/differentiable_numerics/nonlinear_block_solver.hpp"
#include "smith/differentiable_numerics/system_solver.hpp"
#include "smith/differentiable_numerics/solid_mechanics_system.hpp"
#include "smith/differentiable_numerics/differentiable_physics.hpp"
#include "smith/differentiable_numerics/evaluate_objective.hpp"
#include "smith/differentiable_numerics/differentiable_test_utils.hpp"
#include "smith/differentiable_numerics/paraview_writer.hpp"
#include "smith/physics/materials/solid_material.hpp"
smith::ApplicationManager application_manager(argc, argv);
axom::sidre::DataStore datastore;
smith::StateManager::initialize(datastore, "composable_solid_mechanics");
Mesh Construction¶
constexpr int dim = 3;
constexpr int order = 1;
auto mesh = std::make_shared<smith::Mesh>(
mfem::Mesh::MakeCartesian3D(8, 2, 2, mfem::Element::HEXAHEDRON, 1.0, 0.1, 0.1), "mesh", 0, 0);
mesh->addDomainOfBoundaryElements("left", smith::by_attr<dim>(3));
mesh->addDomainOfBoundaryElements("right", smith::by_attr<dim>(5));
Solver and Field Registration¶
Registration declares the displacement fields and the Young's modulus parameter
on the shared FieldStore. Register parameter fields with
registerParameterFields(field_store, ...) before building the system.
smith::LinearSolverOptions linear_options{.linear_solver = smith::LinearSolver::CG,
.preconditioner = smith::Preconditioner::HypreAMG,
.relative_tol = 1e-6,
.absolute_tol = 1e-10,
.max_iterations = 80,
.print_level = 0};
smith::NonlinearSolverOptions nonlinear_options{.nonlin_solver = smith::NonlinearSolver::TrustRegion,
.relative_tol = 1e-7,
.absolute_tol = 1e-8,
.max_iterations = 15,
.print_level = 0};
smith::SolidMechanicsOptions output_options{.enable_stress_output = true, .output_cauchy_stress = true};
System Build and Material Setup¶
The build step consumes the registered field pack and parameter bundle. The
material reads the parameter field and uses the TimeInfo material interface.
auto solid_system = smith::buildSolidMechanicsSystem<dim, order>(
nonlinear_options, linear_options, output_options, mesh, smith::FieldType<smith::L2<0>>("youngs_modulus"));
constexpr double E = 100.0;
constexpr double nu = 0.25;
solid_system->setMaterial(YoungsModulusNeoHookeanWithTimeInfo{.density = 1.0, .nu = nu}, mesh->entireBodyName());
Boundary Conditions and Loads¶
The left boundary fixes selected displacement components. The body force and traction provide the applied loads.
solid_system->setDisplacementBC(mesh->domain("left"), std::vector<int>{0, 2});
solid_system->addBodyForce(mesh->entireBodyName(), [](double, auto X, auto, auto, auto, auto... /*args*/) {
auto body_force = 0.0 * X;
body_force[1] = -0.02;
return body_force;
});
solid_system->addTraction("right", [](double, auto X, auto, auto, auto, auto, auto... /*args*/) {
auto traction = 0.0 * X;
traction[0] = -0.01;
return traction;
});
Physics Creation and Initial Conditions¶
The differentiable physics object owns the timestep loop. The example sets the parameter field and seeds non-zero initial displacement and velocity fields.
auto physics = smith::makeDifferentiablePhysics(solid_system, "composable_solid_mechanics");
if (solid_system->cycle_zero_systems.empty()) {
throw std::runtime_error("Expected cycle-zero solve for implicit dynamics.");
}
physics->getFieldParam("param_youngs_modulus").get()->setFromFieldFunction([=](smith::tensor<double, dim>) {
return E;
});
auto initial_displacement = [](smith::tensor<double, dim> X) {
auto displacement = 0.0 * X;
displacement[0] = 1.0e-3 * X[0];
return displacement;
};
auto initial_velocity = [](smith::tensor<double, dim> X) {
auto velocity = 0.0 * X;
velocity[1] = 2.0e-2 * X[0];
return velocity;
};
physics->getInitialFieldState("displacement_solve_state").get()->setFromFieldFunction(initial_displacement);
physics->getInitialFieldState("displacement").get()->setFromFieldFunction(initial_displacement);
physics->getInitialFieldState("velocity").get()->setFromFieldFunction(initial_velocity);
physics->getInitialFieldState("acceleration").get()->setFromFieldFunction([](smith::tensor<double, dim>) {
return smith::tensor<double, dim>{};
});
Advance, Sensitivities, and Reactions¶
using DispSpace = smith::H1<order, dim>;
const auto qoi_fields = std::vector<smith::FieldState>{physics->getInitialFieldState("displacement"),
physics->getInitialFieldState("velocity")};
smith::FunctionalObjective<dim, smith::Parameters<DispSpace, DispSpace>> qoi("solid_dynamic_energy_proxy", mesh,
smith::spaces(qoi_fields));
qoi.addBodyIntegral(mesh->entireBodyName(), [](const smith::TimeInfo&, auto, auto U, auto V) {
auto u = smith::get<smith::VALUE>(U);
auto v = smith::get<smith::VALUE>(V);
return 0.5 * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]) + 0.05 * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
});
constexpr double dt = 0.25;
constexpr int num_steps = 3;
auto current_qoi_fields = [&]() {
return std::vector<smith::FieldState>{physics->getFieldState("displacement"), physics->getFieldState("velocity")};
};
auto qoi_state =
0.0 * smith::evaluateObjective(qoi, physics->getShapeDispFieldState(), qoi_fields,
smith::TimeInfo(physics->time(), dt, static_cast<size_t>(physics->cycle())));
for (int step = 0; step < num_steps; ++step) {
physics->advanceTimestep(dt);
qoi_state = qoi_state +
smith::evaluateObjective(qoi, physics->getShapeDispFieldState(), current_qoi_fields(),
smith::TimeInfo(physics->time(), dt, static_cast<size_t>(physics->cycle())));
}
std::cout << "reaction norm: " << physics->getReactionStates().front().get()->Norml2() << '\n';
gretl::set_as_objective(qoi_state);
std::cout << "QoI value: " << qoi_state.get() << '\n';
qoi_state.data_store().back_prop();
auto shape_displacement = physics->getShapeDispFieldState();
auto initial_displacement_state = physics->getInitialFieldState("displacement");
auto initial_velocity_state = physics->getInitialFieldState("velocity");
auto youngs_modulus_state = physics->getFieldParam("param_youngs_modulus");
std::cout << "dQoI/d(shape) norm: " << shape_displacement.get_dual()->Norml2() << '\n';
std::cout << "dQoI/d(youngs_modulus) norm: " << youngs_modulus_state.get_dual()->Norml2() << '\n';
std::cout << "dQoI/d(initial displacement) norm: " << initial_displacement_state.get_dual()->Norml2() << '\n';
std::cout << "dQoI/d(initial velocity) norm: " << initial_velocity_state.get_dual()->Norml2() << '\n';
std::cout << "shape FD rate: \n" << smith::checkGradWrt(qoi_state, shape_displacement, 1.0e-2, 4, false) << '\n';
std::cout << "youngs_modulus FD rate: \n"
<< smith::checkGradWrt(qoi_state, youngs_modulus_state, 5.0e-2, 4, false) << '\n';
std::cout << "initial displacement FD rate: \n"
<< smith::checkGradWrt(qoi_state, initial_displacement_state, 5.0e-3, 4, false) << '\n';
std::cout << "initial velocity FD rate: \n"
<< smith::checkGradWrt(qoi_state, initial_velocity_state, 5.0e-3, 4, false) << '\n';
Write ParaView Output¶
auto writer =
smith::createParaviewWriter(*mesh, physics->getFieldStatesAndParamStates(), "paraview_composable_solid_mechanics",
smith::ParaviewWriter::Options{.write_duals = false});
writer.write(physics->cycle(), physics->time(), physics->getFieldStatesAndParamStates());
std::cout << "ParaView output: paraview_composable_solid_mechanics\n";