Tupek/system solver extended

examples/CMakeLists.txt

@@ -23,6 +23,12 @@ install(
 # Add the conduction examples
 add_subdirectory(conduction)
 
+# Add the solid mechanics examples
+add_subdirectory(solid_mechanics)
+
+# Add the thermo-mechanics examples
+add_subdirectory(thermo_mechanics)
+
 # Add the contact examples
 add_subdirectory(contact)
 

examples/contact/homotopy/two_blocks.cpp

@@ -30,9 +30,9 @@ using DensitySpace = smith::L2<disp_order - 1>;
 
 using SolidMaterial = smith::solid_mechanics::NeoHookeanWithFieldDensity;
 using SolidWeakFormT =
-    smith::TimeDiscretizedWeakForm<dim, smith::H1<disp_order, dim>,
-                                   smith::Parameters<smith::H1<disp_order, dim>, smith::H1<disp_order, dim>,
-                                                     smith::H1<disp_order, dim>, DensitySpace>>;
+    smith::FunctionalWeakForm<dim, smith::H1<disp_order, dim>,
+                              smith::Parameters<smith::H1<disp_order, dim>, smith::H1<disp_order, dim>,
+                                                smith::H1<disp_order, dim>, DensitySpace>>;
 
 enum FIELD
 {

examples/inertia_relief/inertia_relief_example.cpp

@@ -38,9 +38,9 @@ using DensitySpace = smith::L2<disp_order - 1>;
 
 using SolidMaterial = smith::solid_mechanics::NeoHookeanWithFieldDensity;
 using SolidWeakFormT =
-    smith::TimeDiscretizedWeakForm<dim, smith::H1<disp_order, dim>,
-                                   smith::Parameters<smith::H1<disp_order, dim>, smith::H1<disp_order, dim>,
-                                                     smith::H1<disp_order, dim>, DensitySpace>>;
+    smith::FunctionalWeakForm<dim, smith::H1<disp_order, dim>,
+                              smith::Parameters<smith::H1<disp_order, dim>, smith::H1<disp_order, dim>,
+                                                smith::H1<disp_order, dim>, DensitySpace>>;
 
 enum FIELD
 {
@@ -212,20 +212,17 @@ int main(int argc, char* argv[])
 
   ObjectiveT mass_objective("mass constraining", mesh, param_space_ptrs);
 
-  mass_objective.addBodyIntegral(smith::DependsOn<1>{}, mesh->entireBodyName(),
-                                 [](double /*t*/, auto /*X*/, auto RHO) { return get<smith::VALUE>(RHO); });
+  mass_objective.addBodyIntegral(
+      mesh->entireBodyName(), [](auto /*t_info*/, auto /*X*/, auto /*U*/, auto RHO) { return get<smith::VALUE>(RHO); });
   double mass = mass_objective.evaluate(time_info, shape_disp.get(), objective_states);
 
   smith::tensor<double, dim> initial_cg;  // center of gravity
 
   for (int i = 0; i < dim; ++i) {
     auto cg_objective = std::make_shared<ObjectiveT>("translation " + std::to_string(i), mesh, param_space_ptrs);
-    cg_objective->addBodyIntegral(smith::DependsOn<0, 1>{}, mesh->entireBodyName(),
-                                  [i](double
-                                      /*time*/,
-                                      auto X, auto U, auto RHO) {
-                                    return (get<smith::VALUE>(X)[i] + get<smith::VALUE>(U)[i]) * get<smith::VALUE>(RHO);
-                                  });
+    cg_objective->addBodyIntegral(mesh->entireBodyName(), [i](auto /*t_info*/, auto X, auto U, auto RHO) {
+      return (get<smith::VALUE>(X)[i] + get<smith::VALUE>(U)[i]) * get<smith::VALUE>(RHO);
+    });
     initial_cg[i] = cg_objective->evaluate(time_info, shape_disp.get(), objective_states) / mass;
 
     constraints.push_back(cg_objective);
@@ -234,8 +231,8 @@ int main(int argc, char* argv[])
   for (int i = 0; i < dim; ++i) {
     auto center_rotation_objective =
         std::make_shared<ObjectiveT>("rotation" + std::to_string(i), mesh, param_space_ptrs);
-    center_rotation_objective->addBodyIntegral(smith::DependsOn<0, 1>{}, mesh->entireBodyName(),
-                                               [i, initial_cg](double /*time*/, auto X, auto U, auto RHO) {
+    center_rotation_objective->addBodyIntegral(mesh->entireBodyName(),
+                                               [i, initial_cg](auto /*t_info*/, auto X, auto U, auto RHO) {
                                                  auto u = get<smith::VALUE>(U);
                                                  auto x = get<smith::VALUE>(X) + u;
                                                  auto dx = x - initial_cg;

examples/solid_mechanics/CMakeLists.txt

@@ -0,0 +1,24 @@
+# Copyright (c) Lawrence Livermore National Security, LLC and
+# other Smith Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+smith_add_executable( NAME        composable_solid_mechanics
+                      SOURCES     composable_solid_mechanics.cpp
+                      OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
+                      DEPENDS_ON  smith
+                      )
+
+install(
+    FILES
+        composable_solid_mechanics.cpp
+    DESTINATION
+        examples/smith/solid_mechanics
+)
+
+if(SMITH_ENABLE_TESTS)
+    blt_add_test(NAME          composable_solid_mechanics
+                 COMMAND       composable_solid_mechanics
+                 NUM_MPI_TASKS 1 )
+endif()

examples/solid_mechanics/composable_solid_mechanics.cpp

@@ -0,0 +1,205 @@
+// Copyright (c) Lawrence Livermore National Security, LLC and
+// other Smith Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/**
+ * @file composable_solid_mechanics.cpp
+ * @brief Dynamic solid-mechanics example using composable differentiable numerics systems.
+ */
+
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+// _includes_start
+#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"
+// _includes_end
+
+namespace {
+
+struct YoungsModulusNeoHookeanWithTimeInfo {
+  using State = smith::solid_mechanics::NeoHookean::State;
+
+  double density;
+  double nu;
+
+  template <typename T, int dim, typename GradVType, typename YoungsType>
+  auto operator()(const smith::TimeInfo&, [[maybe_unused]] State& state, const smith::tensor<T, dim, dim>& grad_u,
+                  const GradVType&, const YoungsType& youngs_modulus) const
+  {
+    using std::log1p;
+    constexpr auto I = smith::Identity<dim>();
+    auto E = smith::get<0>(youngs_modulus);
+    auto G = E / (2.0 * (1.0 + nu));
+    auto K = E / (3.0 * (1.0 - 2.0 * nu));
+    auto lambda = K - (2.0 / dim) * G;
+    auto B_minus_I = grad_u * transpose(grad_u) + transpose(grad_u) + grad_u;
+    auto logJ = log1p(detApIm1(grad_u));
+    auto TK = lambda * logJ * I + G * B_minus_I;
+    auto F = grad_u + I;
+    return dot(TK, inv(transpose(F)));
+  }
+};
+
+}  // namespace
+
+int main(int argc, char* argv[])
+{
+  // _init_start
+  smith::ApplicationManager application_manager(argc, argv);
+  axom::sidre::DataStore datastore;
+  smith::StateManager::initialize(datastore, "composable_solid_mechanics");
+  // _init_end
+
+  // _mesh_start
+  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));
+  // _mesh_end
+
+  // _solver_start
+  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};
+  // _solver_end
+
+  // _build_start
+  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());
+  // _build_end
+
+  // _bc_start
+  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;
+  });
+  // _bc_end
+
+  // _ic_start
+  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>{};
+  });
+  // _ic_end
+
+  // _run_start
+  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';
+  // _run_end
+
+  // _output_start
+  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";
+  // _output_end
+
+  return 0;
+}

examples/thermo_mechanics/CMakeLists.txt

@@ -0,0 +1,34 @@
+# Copyright (c) Lawrence Livermore National Security, LLC and
+# other Smith Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+smith_add_executable( NAME        composable_thermo_mechanics
+                      SOURCES     composable_thermo_mechanics.cpp
+                      OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
+                      DEPENDS_ON  smith
+                      )
+
+smith_add_executable( NAME        composable_thermo_mechanics_advanced
+                      SOURCES     composable_thermo_mechanics_advanced.cpp
+                      OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
+                      DEPENDS_ON  smith
+                      )
+
+install(
+    FILES
+        composable_thermo_mechanics.cpp
+        composable_thermo_mechanics_advanced.cpp
+    DESTINATION
+        examples/smith/thermo_mechanics
+)
+
+if(SMITH_ENABLE_TESTS)
+    blt_add_test(NAME          composable_thermo_mechanics
+                 COMMAND       composable_thermo_mechanics
+                 NUM_MPI_TASKS 1 )
+    blt_add_test(NAME          composable_thermo_mechanics_advanced
+                 COMMAND       composable_thermo_mechanics_advanced
+                 NUM_MPI_TASKS 1 )
+endif()