Composable Thermo-Mechanics Advanced Example¶
This demo extends the minimal thermo-mechanics setup. It adds a parameter field, a staged solver, a differentiable quantity of interest, a finite-difference check, and ParaView output.
The full source code lives in examples/thermo_mechanics/composable_thermo_mechanics_advanced.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/thermal_system.hpp"
#include "smith/differentiable_numerics/thermo_mechanics_system.hpp"
#include "smith/differentiable_numerics/combined_system.hpp"
#include "smith/differentiable_numerics/differentiable_physics.hpp"
#include "smith/differentiable_numerics/paraview_writer.hpp"
#include "smith/differentiable_numerics/evaluate_objective.hpp"
#include "smith/differentiable_numerics/differentiable_test_utils.hpp"
#include "smith/differentiable_numerics/make_time_info_material.hpp"
#include "smith/physics/materials/green_saint_venant_thermoelastic.hpp"
smith::ApplicationManager application_manager(argc, argv);
axom::sidre::DataStore datastore;
smith::StateManager::initialize(datastore, "composable_thermo_mechanics_advanced");
Mesh and Field Setup¶
constexpr int dim = 3;
constexpr int order = 1;
using DispSpace = smith::H1<order, dim>;
using TempSpace = smith::H1<order>;
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 Config and Field Registration¶
Registration declares solid, thermal, and parameter fields on one shared
FieldStore. The thermal-expansion parameter is registered with
registerParameterFields(field_store, ...) before either system is built.
smith::LinearSolverOptions coupled_linear{.linear_solver = smith::LinearSolver::SuperLU,
.relative_tol = 1e-8,
.absolute_tol = 1e-10,
.max_iterations = 200,
.print_level = 0};
smith::NonlinearSolverOptions coupled_nonlin{.nonlin_solver = smith::NonlinearSolver::NewtonLineSearch,
.relative_tol = 1e-8,
.absolute_tol = 1e-8,
.max_iterations = 12,
.max_line_search_iterations = 6,
.print_level = 0};
auto field_store = std::make_shared<smith::FieldStore>(mesh, 200);
smith::SolidMechanicsOptions solid_options{.enable_stress_output = true, .output_cauchy_stress = true};
auto solid_fields = smith::registerSolidMechanicsFields<dim, order, smith::QuasiStaticSecondOrderTimeIntegrationRule>(
field_store, solid_options);
auto thermal_fields =
smith::registerThermalFields<dim, order, smith::BackwardEulerFirstOrderTimeIntegrationRule>(field_store);
auto param_fields =
smith::registerParameterFields(field_store, smith::FieldType<smith::L2<0>>("thermal_expansion_scaling"));
System Build and Coupling¶
The build step creates solid and thermal systems from the registered field
packs. combineSystems(...) attaches them to one staged solver.
auto solid_system = smith::buildSolidMechanicsSystem(nullptr, solid_options, solid_fields,
smith::couplingFields(thermal_fields), param_fields);
auto thermal_system = smith::buildThermalSystem(nullptr, smith::ThermalOptions{}, thermal_fields,
smith::couplingFields(solid_fields), param_fields);
auto material =
smith::makeTimeInfoMaterial(smith::thermomechanics::ParameterizedGreenSaintVenantThermoelasticMaterial{
1.0, 100.0, 0.25, 1.0, 0.0025, 0.0, 0.05});
smith::setCoupledThermoMechanicsMaterial(solid_system, thermal_system, material, mesh->entireBodyName());
auto coupled_solver = std::make_shared<smith::SystemSolver>(10);
coupled_solver->addSubsystemSolver({0, 1}, smith::buildNonlinearBlockSolver(coupled_nonlin, coupled_linear, *mesh),
1.0);
auto coupled_system = smith::combineSystems(coupled_solver, solid_system, thermal_system);
std::string physics_name = "composable_thermo_mechanics_advanced";
auto physics = smith::makeDifferentiablePhysics(coupled_system, physics_name);
auto output_states = field_store->getOutputFieldStates();
auto output_writer =
smith::createParaviewWriter(*mesh, output_states, "paraview_composable_thermo_mechanics_advanced",
smith::ParaviewWriter::Options{.write_duals = false});
Boundary Conditions and Loads¶
Boundary conditions are applied on the left and right boundaries. Loads are added through the solid and thermal systems before timestepping.
field_store->getParameterFields()[0].get()->setFromFieldFunction([](smith::tensor<double, dim>) { return 1.0; });
solid_system->setDisplacementBC(mesh->domain("left"));
thermal_system->setTemperatureBC(mesh->domain("left"), [](auto, auto) { return 1.0; });
thermal_system->setTemperatureBC(mesh->domain("right"), [](auto, auto) { return 0.0; });
solid_system->addTraction("right",
[](double, auto X, auto /*n*/, auto /*u*/, auto /*v*/, auto /*a*/, auto... /*args*/) {
auto traction = 0.0 * X;
traction[0] = -0.005;
return traction;
});
thermal_system->addHeatSource(mesh->entireBodyName(), [](auto, auto... /*unused*/) { return 0.1; });
QoI Definition and Timestep Advance¶
auto fetch_qoi_fields = [&]() {
return std::vector<smith::FieldState>{field_store->getField("displacement"), field_store->getField("temperature")};
};
smith::FunctionalObjective<dim, smith::Parameters<DispSpace, TempSpace>> qoi("thermo_mechanical_energy_proxy", mesh,
smith::spaces(fetch_qoi_fields()));
qoi.addBodyIntegral(mesh->entireBodyName(), [](auto, auto /*X*/, auto U, auto Theta) {
auto u = smith::get<smith::VALUE>(U);
auto theta = smith::get<smith::VALUE>(Theta);
return 0.5 * u[0] * u[0] + 0.05 * theta * theta;
});
auto qoi_state =
0.0 * smith::evaluateObjective(qoi, physics->getShapeDispFieldState(), fetch_qoi_fields(),
smith::TimeInfo(physics->time(), 1.0, static_cast<size_t>(physics->cycle())));
constexpr double dt = 0.5;
constexpr int qoi_steps = 1;
for (int step = 0; step < qoi_steps; ++step) {
physics->advanceTimestep(dt);
qoi_state = qoi_state +
smith::evaluateObjective(qoi, physics->getShapeDispFieldState(), fetch_qoi_fields(),
smith::TimeInfo(physics->time(), dt, static_cast<size_t>(physics->cycle())));
}
ParaView Output¶
output_writer.write(physics->cycle(), physics->time(), field_store->getOutputFieldStates());
std::cout << "ParaView output: paraview_composable_thermo_mechanics_advanced\n";
Sensitivity and Finite-Difference Check¶
gretl::set_as_objective(qoi_state);
auto qoi_value = qoi_state.get();
std::cout << "QoI value: " << qoi_value << '\n';
qoi_state.data_store().back_prop();
auto parameter_state = field_store->getParameterFields()[0];
auto parameter_sensitivity = parameter_state.get_dual()->Norml2();
std::cout << "dQoI/d(thermal_expansion_scaling) norm: " << parameter_sensitivity << '\n';
SLIC_ERROR_ROOT_IF(parameter_sensitivity <= 0.0, "Expected non-zero QoI sensitivity.");
auto fd_order = smith::checkGradWrt(qoi_state, parameter_state, 5.0e-2, 4, true);
std::cout << "finite-difference convergence rate: " << fd_order << '\n';
SLIC_ERROR_ROOT_IF(fd_order < 0.7, "Finite-difference check did not converge.");