diff --git a/examples/inertia_relief/inertia_relief_example.cpp b/examples/inertia_relief/inertia_relief_example.cpp index a1a4c3964..5aaba18fe 100644 --- a/examples/inertia_relief/inertia_relief_example.cpp +++ b/examples/inertia_relief/inertia_relief_example.cpp @@ -80,37 +80,6 @@ class ParaviewWriter { StateVecs dual_states; }; -auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector& states, - std::string output_name) -{ - if (output_name == "") { - output_name = "default"; - } - - ParaviewWriter::StateVecs output_states; - for (const auto& s : states) { - output_states.push_back(std::make_shared(s->space(), s->name())); - } - - auto non_const_mesh = const_cast(&mesh); - auto paraview_dc = std::make_unique(output_name, non_const_mesh); - int max_order_in_fields = 0; - - // Find the maximum polynomial order in the physics module's states - for (const auto& state : output_states) { - paraview_dc->RegisterField(state->name(), &state->gridFunction()); - max_order_in_fields = std::max(max_order_in_fields, state->space().GetOrder(0)); - } - - // Set the options for the paraview output files - paraview_dc->SetLevelsOfDetail(max_order_in_fields); - paraview_dc->SetHighOrderOutput(true); - paraview_dc->SetDataFormat(mfem::VTKFormat::BINARY); - paraview_dc->SetCompression(true); - - return ParaviewWriter(std::move(paraview_dc), output_states, {}); -} - /* Nonlinear problem of the form * F(X) = [ r(u) + (dc/du)^T l ] = [ 0 ] * [ -c(u) ] [ 0 ] @@ -302,7 +271,7 @@ int main(int argc, char* argv[]) return u; }); - auto writer = createParaviewOutput(mesh->mfemParMesh(), objective_states, "inertia_relief"); + auto writer = createParaviewWriter(mesh->mfemParMesh(), objective_states, "inertia_relief"); if (visualize) { writer.write(0, 0.0, objective_states); } diff --git a/src/smith/differentiable_numerics/CMakeLists.txt b/src/smith/differentiable_numerics/CMakeLists.txt index cfeb470ab..61d4e28b4 100644 --- a/src/smith/differentiable_numerics/CMakeLists.txt +++ b/src/smith/differentiable_numerics/CMakeLists.txt @@ -6,20 +6,33 @@ set(differentiable_numerics_sources field_state.cpp - differentiable_utils.cpp differentiable_physics.cpp lumped_mass_explicit_newmark_state_advancer.cpp + solid_mechanics_state_advancer.cpp + differentiable_solver.cpp + nonlinear_solve.cpp + evaluate_objective.cpp + dirichlet_boundary_conditions.cpp ) set(differentiable_numerics_headers field_state.hpp state_advancer.hpp + reaction.hpp + differentiable_solver.hpp + differentiable_physics.hpp timestep_estimator.hpp explicit_dynamic_solve.hpp lumped_mass_weak_form.hpp - differentiable_utils.hpp - differentiable_physics.hpp lumped_mass_explicit_newmark_state_advancer.hpp + solid_mechanics_state_advancer.hpp + nonlinear_solve.hpp + evaluate_objective.hpp + dirichlet_boundary_conditions.hpp + time_integration_rule.hpp + time_discretized_weak_form.hpp + differentiable_solid_mechanics.hpp + paraview_writer.hpp ) set(differentiable_numerics_depends diff --git a/src/smith/differentiable_numerics/differentiable_physics.cpp b/src/smith/differentiable_numerics/differentiable_physics.cpp index 2fca62825..6811b841d 100644 --- a/src/smith/differentiable_numerics/differentiable_physics.cpp +++ b/src/smith/differentiable_numerics/differentiable_physics.cpp @@ -1,13 +1,10 @@ -// 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) - +#include "gretl/data_store.hpp" #include "smith/differentiable_numerics/differentiable_physics.hpp" #include "smith/physics/weak_form.hpp" #include "smith/physics/mesh.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" #include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/differentiable_numerics/reaction.hpp" #include "gretl/data_store.hpp" namespace smith { @@ -36,10 +33,12 @@ gretl::State make_milestone(const std::vector& states) DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr mesh, std::shared_ptr graph, const FieldState& shape_disp, const std::vector& states, const std::vector& params, - std::shared_ptr advancer, std::string mech_name) + std::shared_ptr advancer, std::string mech_name, + const std::vector& reaction_names) : BasePhysics(mech_name, mesh, 0, 0.0, false), // the false is checkpoint_to_disk checkpointer_(graph), - advancer_(advancer) + advancer_(advancer), + reaction_names_(reaction_names) { SLIC_ERROR_IF(states.size() == 0, "Must have a least 1 state for a mechanics."); field_shape_displacement_ = std::make_unique(shape_disp); @@ -48,14 +47,18 @@ DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr mesh, std::sh field_states_.push_back(s); initial_field_states_.push_back(s); state_name_to_field_index_[s.get()->name()] = i; - state_names.push_back(s.get()->name()); + state_names_.push_back(s.get()->name()); } for (size_t i = 0; i < params.size(); ++i) { const auto& p = params[i]; field_params_.push_back(p); param_name_to_field_index_[p.get()->name()] = i; - param_names.push_back(p.get()->name()); + param_names_.push_back(p.get()->name()); + } + + for (size_t i = 0; i < reaction_names_.size(); ++i) { + reaction_name_to_reaction_index_[reaction_names_[i]] = i; } completeSetup(); @@ -84,9 +87,11 @@ void DifferentiablePhysics::resetAdjointStates() gretl_assert(checkpointer_->check_validity()); } -std::vector DifferentiablePhysics::stateNames() const { return state_names; } +std::vector DifferentiablePhysics::stateNames() const { return state_names_; } -std::vector DifferentiablePhysics::parameterNames() const { return param_names; } +std::vector DifferentiablePhysics::parameterNames() const { return param_names_; } + +std::vector DifferentiablePhysics::dualNames() const { return reaction_names_; } const FiniteElementState& DifferentiablePhysics::state([[maybe_unused]] const std::string& field_name) const { @@ -99,6 +104,21 @@ const FiniteElementState& DifferentiablePhysics::state([[maybe_unused]] const st return *field_states_[state_index].get(); } +const FiniteElementDual& DifferentiablePhysics::dual(const std::string& dual_name) const +{ + SLIC_ERROR_IF( + reaction_name_to_reaction_index_.find(dual_name) == reaction_name_to_reaction_index_.end(), + axom::fmt::format("Could not find dual named {0} in mesh with tag \"{1}\" to get", dual_name, mesh_->tag())); + size_t reaction_index = reaction_name_to_reaction_index_.at(dual_name); + SLIC_ERROR_IF( + reaction_index >= reaction_names_.size(), + "Dual reactions not correctly allocated yet, cannot get dual until after initializationStep is called."); + + TimeInfo time_info(time_old_, dt_old_, static_cast(cycle_old_)); + reaction_states_ = advancer_->computeReactions(time_info, *field_shape_displacement_, field_states_, field_params_); + return *reaction_states_[reaction_index].get(); +} + FiniteElementState DifferentiablePhysics::loadCheckpointedState(const std::string& state_name, int cycle) { SLIC_ERROR_IF(cycle != cycle_, @@ -154,14 +174,28 @@ void DifferentiablePhysics::setAdjointLoad( for (auto string_dual_pair : string_to_dual) { std::string field_name = string_dual_pair.first; const smith::FiniteElementDual& dual = string_dual_pair.second; - SLIC_ERROR_IF( - state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), - axom::fmt::format("Could not find dual named {0} in mesh with tag {1} to set", field_name, mesh_->tag())); + SLIC_ERROR_IF(state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), + axom::fmt::format("Could not find dual named {0} in mesh with tag {1}", field_name, mesh_->tag())); size_t state_index = state_name_to_field_index_.at(field_name); *field_states_[state_index].get_dual() += dual; } } +void DifferentiablePhysics::setDualAdjointBcs( + std::unordered_map string_to_bc) +{ + for (auto string_bc_pair : string_to_bc) { + std::string reaction_name = string_bc_pair.first; + const smith::FiniteElementState& reaction_dual = string_bc_pair.second; + SLIC_ERROR_IF( + reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(), + axom::fmt::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}", + reaction_name, mesh_->tag())); + size_t reaction_index = reaction_name_to_reaction_index_.at(reaction_name); + *reaction_states_[reaction_index].get_dual() += reaction_dual; + } +} + const FiniteElementState& DifferentiablePhysics::adjoint([[maybe_unused]] const std::string& adjoint_name) const { // MRT, not implemented @@ -176,8 +210,12 @@ void DifferentiablePhysics::advanceTimestep(double dt) milestones_.push_back(make_milestone(field_states_).step()); } + cycle_old_ = cycle_; + time_old_ = time_; + dt_old_ = dt; + TimeInfo time_info(time_, dt, static_cast(cycle_)); - field_states_ = advancer_->advanceState(*field_shape_displacement_, field_states_, field_params_, time_info); + field_states_ = advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_); cycle_++; time_ += dt; @@ -232,7 +270,7 @@ DifferentiablePhysics::computeInitialConditionSensitivity() const return map; } -std::vector DifferentiablePhysics::getAllFieldStates() const +std::vector DifferentiablePhysics::getFieldStatesAndParamStates() const { std::vector fields; fields.insert(fields.end(), field_states_.begin(), field_states_.end()); diff --git a/src/smith/differentiable_numerics/differentiable_physics.hpp b/src/smith/differentiable_numerics/differentiable_physics.hpp index 8806f4134..5ee99fa8c 100644 --- a/src/smith/differentiable_numerics/differentiable_physics.hpp +++ b/src/smith/differentiable_numerics/differentiable_physics.hpp @@ -26,6 +26,7 @@ class WeakForm; class DifferentiableSolver; class StateAdvancer; class TimestepEstimator; +class Reaction; /// @brief Implementation of BasePhysics which uses FieldStates and gretl to track the computational graph, dynamically /// checkpoint, and back-propagate sensitivities. @@ -35,7 +36,7 @@ class DifferentiablePhysics : public BasePhysics { DifferentiablePhysics(std::shared_ptr mesh, std::shared_ptr graph, const FieldState& shape_disp, const std::vector& states, const std::vector& params, std::shared_ptr advancer, - std::string mech_name); + std::string physics_name, const std::vector& reaction_names = {}); /// @brief destructor ~DifferentiablePhysics() {} @@ -54,11 +55,14 @@ class DifferentiablePhysics : public BasePhysics { /// @overload std::vector parameterNames() const override; + /// @overload + std::vector dualNames() const override; + /// @overload const FiniteElementState& state(const std::string& state_name) const override; /// @overload - FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) override; + const FiniteElementDual& dual(const std::string& dual_name) const override; /// @overload const FiniteElementState& shapeDisplacement() const override; @@ -69,6 +73,9 @@ class DifferentiablePhysics : public BasePhysics { /// @overload const FiniteElementState& parameter(const std::string& parameter_name) const override; + /// @overload + FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) override; + /// @overload void setState(const std::string& state_name, const FiniteElementState& s) override; @@ -81,6 +88,9 @@ class DifferentiablePhysics : public BasePhysics { /// @overload void setAdjointLoad(std::unordered_map string_to_dual) override; + /// @overload + void setDualAdjointBcs(std::unordered_map string_to_bc) override; + /// @overload const FiniteElementState& adjoint(const std::string& adjoint_name) const override; @@ -100,12 +110,27 @@ class DifferentiablePhysics : public BasePhysics { const std::unordered_map computeInitialConditionSensitivity() const override; + /// @brief Get all the initial FieldStates + std::vector getInitialFieldStates() const { return initial_field_states_; } + + /// @brief Get all the current FieldStates + std::vector getFieldStates() const { return field_states_; } + + /// @brief Get all the parameter FieldStates + std::vector getFieldParams() const { return field_params_; } + /// @brief Get all the FieldStates... states first, parameters next - std::vector getAllFieldStates() const; + std::vector getFieldStatesAndParamStates() const; /// @brief Get the shape displacement FieldState FieldState getShapeDispFieldState() const; + /// @brief Get the current reactionStates + std::vector getReactionStates() const { return reaction_states_; } + + /// @brief Get state advancer + std::shared_ptr getStateAdvancer() const { return advancer_; } + private: std::shared_ptr checkpointer_; ///< gretl data store manages dynamic checkpointing logic std::shared_ptr advancer_; ///< abstract interface for advancing state from one cycle to the next @@ -119,13 +144,21 @@ class DifferentiablePhysics : public BasePhysics { std::map state_name_to_field_index_; ///< map from state names to field index std::map param_name_to_field_index_; ///< map from param names to param index - std::vector state_names; ///< names of all the states in order - std::vector param_names; ///< names of all the states in order + std::vector state_names_; ///< names of all the states in order + std::vector param_names_; ///< names of all the states in order + + mutable std::vector reaction_states_; ///< all the reactions registered for the physics + std::map reaction_name_to_reaction_index_; ///< map from reaction names to reaction index + std::vector reaction_names_; ///< names for all the relevant reactions/reactions std::vector milestones_; ///< a record of the steps in the graph that represent the end conditions of ///< advanceTimestep(dt). this information is used to halt the gretl graph when ///< back-propagating to allow users of reverseAdjointTimestep to specify ///< adjoint loads and to retrieve timestep sensitivity information. + + double time_old_ = 0.0; + double dt_old_ = 0.0; + int cycle_old_ = 0; }; } // namespace smith diff --git a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp new file mode 100644 index 000000000..58a94a161 --- /dev/null +++ b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp @@ -0,0 +1,48 @@ +// Copyright (c) 2019-2024, 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 differentiable_solid_mechanics.hpp + * + */ + +#pragma once + +#include +#include "smith/differentiable_numerics/solid_mechanics_state_advancer.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" + +namespace smith { + +/// Helper function to generate the base-physics for solid mechanics +/// This will return a tuple of shared pointers to the: BasePhysics, WeakForm, and DirichetBoundaryConditions +/// Only the BasePhysics needs to stay in scope. The others are returned to the user so they can define the WeakForm +/// integrals, and to specify space and time varying boundary conditions +template +auto buildSolidMechanics(std::shared_ptr mesh, + std::shared_ptr d_solid_nonlinear_solver, + smith::SecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, + std::string physics_name, const std::vector& param_names = {}) +{ + auto graph = std::make_shared(num_checkpoints); + auto [shape_disp, states, params, time, solid_mechanics_weak_form] = + SolidMechanicsStateAdvancer::buildWeakFormAndStates( + mesh, graph, time_rule, physics_name, param_names); + + auto vector_bcs = std::make_shared( + mesh->mfemParMesh(), space(states[SolidMechanicsStateAdvancer::DISPLACEMENT])); + + auto state_advancer = std::make_shared(d_solid_nonlinear_solver, vector_bcs, + solid_mechanics_weak_form, time_rule); + + auto physics = std::make_shared(mesh, graph, shape_disp, states, params, state_advancer, + physics_name, std::vector{"reactions"}); + + return std::make_tuple(physics, solid_mechanics_weak_form, vector_bcs); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp new file mode 100644 index 000000000..7fabac39c --- /dev/null +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -0,0 +1,286 @@ +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/physics/state/finite_element_state.hpp" +#include "smith/physics/state/finite_element_dual.hpp" +#include "smith/numerics/stdfunction_operator.hpp" +#include "smith/numerics/equation_solver.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" +#include "smith/physics/mesh.hpp" +#include "mfem.hpp" + +namespace smith { + +using smith::FiniteElementDual; +using smith::FiniteElementState; + +/// @brief Utility to compute the matrix norm +inline double matrixNorm(std::unique_ptr& K) +{ + mfem::HypreParMatrix* H = K.get(); + hypre_ParCSRMatrix* Hhypre = static_cast(*H); + double Hfronorm; + hypre_ParCSRMatrixNormFro(Hhypre, &Hfronorm); + return Hfronorm; +} + +/// @brief Utility to compute 0.5*norm(K-K.T) +inline double skewMatrixNorm(std::unique_ptr& K) +{ + auto K_T = std::unique_ptr(K->Transpose()); + K_T->Add(-1.0, *K); + (*K_T) *= 0.5; + mfem::HypreParMatrix* H = K_T.get(); + hypre_ParCSRMatrix* Hhypre = static_cast(*H); + double Hfronorm; + hypre_ParCSRMatrixNormFro(Hhypre, &Hfronorm); + return Hfronorm; +} + +/// @brief Initialize mfem solver if near-nullspace is needed +void initializeSolver(mfem::Solver* mfem_solver, const smith::FiniteElementState& u) +{ + // If the user wants the AMG preconditioner with a linear solver, set the pfes + // to be the displacement + auto* amg_prec = dynamic_cast(mfem_solver); + if (amg_prec) { + // ZRA - Iterative refinement tends to be more expensive than it is worth + // We should add a flag allowing users to enable it + + // bool iterative_refinement = false; + // amg_prec->SetElasticityOptions(&displacement_.space(), iterative_refinement); + + // SetElasticityOptions only works with byVDIM ordering, some evidence that it is not often optimal + amg_prec->SetSystemsOptions(u.space().GetVDim(), smith::ordering == mfem::Ordering::byNODES); + } + +#ifdef SMITH_USE_PETSC + auto* space_dep_pc = dynamic_cast(mfem_solver); + if (space_dep_pc) { + // This call sets the displacement ParFiniteElementSpace used to get the spatial coordinates and to + // generate the near null space for the PCGAMG preconditioner + mfem::ParFiniteElementSpace* space = const_cast(&u.space()); + space_dep_pc->SetFESpace(space); + } +#endif +} + +LinearDifferentiableSolver::LinearDifferentiableSolver(std::unique_ptr s, std::unique_ptr p) + : mfem_solver(std::move(s)), mfem_preconditioner(std::move(p)) +{ +} + +void LinearDifferentiableSolver::completeSetup(const smith::FiniteElementState& u) +{ + initializeSolver(mfem_preconditioner.get(), u); +} + +std::shared_ptr LinearDifferentiableSolver::solve( + const FiniteElementState& u, // initial guess + std::function equation, + std::function(const FiniteElementState&)> jacobian) const +{ + SMITH_MARK_FUNCTION; + auto r = equation(u); + auto du = std::make_shared(u.space(), "u"); + *du = 0.0; + auto Jptr = jacobian(u); + mfem_solver->SetOperator(*Jptr); + mfem_solver->Mult(r, *du); + *du -= u; + *du *= -1.0; + return du; // return u - K^{-1}r +} + +std::shared_ptr LinearDifferentiableSolver::solveAdjoint( + const FiniteElementDual& u_bar, std::unique_ptr jacobian_transposed) const +{ + SMITH_MARK_FUNCTION; + + auto ds = std::make_shared(u_bar.space(), "ds"); + mfem_solver->SetOperator(*jacobian_transposed); + mfem_solver->Mult(u_bar, *ds); + return ds; +} + +NonlinearDifferentiableSolver::NonlinearDifferentiableSolver(std::unique_ptr s) + : nonlinear_solver_(std::move(s)) +{ +} + +void NonlinearDifferentiableSolver::completeSetup(const smith::FiniteElementState& u) +{ + initializeSolver(&nonlinear_solver_->preconditioner(), u); +} + +std::shared_ptr NonlinearDifferentiableSolver::solve( + const FiniteElementState& u_guess, // initial guess + std::function equation, + std::function(const FiniteElementState&)> jacobian) const +{ + SMITH_MARK_FUNCTION; + + auto u = std::make_shared(u_guess); + + auto residual_op_ = std::make_unique( + u->space().TrueVSize(), + + [&u, &equation](const mfem::Vector& u_, mfem::Vector& r_) { + FiniteElementState uu(u->space(), "uu"); + uu = u_; + r_ = equation(uu); + }, + + [&u, &jacobian, this](const mfem::Vector& u_) -> mfem::Operator& { + FiniteElementState uu(u->space(), "uu"); + uu = u_; + J_.reset(); + J_ = jacobian(uu); + return *J_; + }); + + nonlinear_solver_->setOperator(*residual_op_); + nonlinear_solver_->solve(*u); + + // std::cout << "solution norm = " << u->Norml2() << std::endl; + + return u; +} + +std::shared_ptr NonlinearDifferentiableSolver::solveAdjoint( + const FiniteElementDual& x_bar, std::unique_ptr jacobian_transposed) const +{ + SMITH_MARK_FUNCTION; + + auto ds = std::make_shared(x_bar.space(), "ds"); + auto& linear_solver = nonlinear_solver_->linearSolver(); + linear_solver.SetOperator(*jacobian_transposed); + linear_solver.Mult(x_bar, *ds); + + return ds; +} + +void NonlinearDifferentiableSolver::clearMemory() const { J_.reset(); } + +LinearDifferentiableBlockSolver::LinearDifferentiableBlockSolver(std::unique_ptr s, + std::unique_ptr p) + : mfem_solver(std::move(s)), mfem_preconditioner(std::move(p)) +{ +} + +void LinearDifferentiableBlockSolver::completeSetup(const std::vector& us) +{ + initializeSolver(mfem_preconditioner.get(), us[0]); +} + +std::vector LinearDifferentiableBlockSolver::solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residual_funcs, + std::function>(const std::vector&)> jacobian_funcs) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_guesses.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_guesses[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_du = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_du->GetBlock(row_i) = *u_guesses[static_cast(row_i)]; + } + + auto residuals = residual_funcs(u_guesses); + auto block_r = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_r->GetBlock(row_i) = residuals[static_cast(row_i)]; + } + + auto jacs = jacobian_funcs(u_guesses); + auto block_jac = std::make_unique(block_offsets); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + block_jac->SetBlock(i, j, jacs[static_cast(i)][static_cast(j)].get()); + } + } + + mfem_solver->SetOperator(*block_jac); + + mfem_solver->Mult(*block_r, *block_du); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] -= block_du->GetBlock(row_i); + } + *block_du = 0.0; + + return u_guesses; +} + +std::vector LinearDifferentiableBlockSolver::solveAdjoint( + const std::vector& u_bars, std::vector>& jacobian_transposed) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_bars.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + std::vector u_duals(static_cast(num_rows)); + for (int row_i = 0; row_i < num_rows; ++row_i) { + u_duals[static_cast(row_i)] = std::make_shared( + u_bars[static_cast(row_i)]->space(), "u_dual_" + std::to_string(row_i)); + } + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_bars[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_ds = std::make_unique(block_offsets); + *block_ds = 0.0; + + auto block_r = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_r->GetBlock(row_i) = *u_bars[static_cast(row_i)]; + } + + auto block_jac = std::make_unique(block_offsets); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + block_jac->SetBlock(i, j, jacobian_transposed[static_cast(i)][static_cast(j)].get()); + } + } + + mfem_solver->SetOperator(*block_jac); + + mfem_solver->Mult(*block_r, *block_ds); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_duals[static_cast(row_i)] = block_ds->GetBlock(row_i); + } + + return u_duals; +} + +std::shared_ptr buildDifferentiableLinearSolve(LinearSolverOptions linear_opts, + const smith::Mesh& mesh) +{ + auto [linear_solver, precond] = smith::buildLinearSolverAndPreconditioner(linear_opts, mesh.getComm()); + return std::make_shared(std::move(linear_solver), std::move(precond)); +} + +std::shared_ptr buildDifferentiableNonlinearSolve( + smith::NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) +{ + auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); + return std::make_shared(std::move(solid_solver)); +} + +} // namespace smith \ No newline at end of file diff --git a/src/smith/differentiable_numerics/differentiable_solver.hpp b/src/smith/differentiable_numerics/differentiable_solver.hpp new file mode 100644 index 000000000..e351d14b6 --- /dev/null +++ b/src/smith/differentiable_numerics/differentiable_solver.hpp @@ -0,0 +1,192 @@ +// 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 differentiable_solver.hpp + * + * @brief This file contains the declaration of the DifferentiableSolver interface + */ + +#pragma once + +#include +#include + +namespace mfem { +class Solver; +class Vector; +class HypreParMatrix; +} // namespace mfem + +namespace smith { + +class EquationSolver; +class BoundaryConditionManager; +class FiniteElementState; +class FiniteElementDual; +class Mesh; +struct NonlinearSolverOptions; +struct LinearSolverOptions; + +/// @brief Abstract interface to DifferentiableSolver inteface. Each dfferenriable solve should provide both its +/// forward solve and an adjoint solve +class DifferentiableSolver { + public: + /// @brief destructor + virtual ~DifferentiableSolver() {} + + /// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space + /// For these cases, it should be called before solve + virtual void completeSetup(const smith::FiniteElementState& u) = 0; + + /// @brief Solve a set of equations with a FiniteElementState as unknown + /// @param u_guess initial guess for solver + /// @param equation std::function for equation to be solved + /// @param jacobian std::function for evaluating the linearized Jacobian about the current solution + /// @return The solution FiniteElementState + virtual std::shared_ptr solve( + const smith::FiniteElementState& u_guess, std::function equation, + std::function(const smith::FiniteElementState&)> jacobian) const = 0; + + /// @brief Solve the (linear) adjoint set of equations with a FiniteElementState as unknown + /// @param u_bar rhs for the solve + /// @param jacobian_transposed the evaluated linearized adjoint space matrix + /// @return The adjoint solution field + virtual std::shared_ptr solveAdjoint( + const smith::FiniteElementDual& u_bar, std::unique_ptr jacobian_transposed) const = 0; + + /// @brief Interface option to clear memory between solves to avoid high-water mark memory usage. + virtual void clearMemory() const {} +}; + +/// @brief Implementation of the DifferentiableSolver interface for the special case of linear solves with linear +/// adjoint solves +class LinearDifferentiableSolver : public DifferentiableSolver { + public: + /// @brief Construct from a linear solver and linear precondition which may be used by the linear solver + LinearDifferentiableSolver(std::unique_ptr s, std::unique_ptr p); + + /// @overload + void completeSetup(const smith::FiniteElementState& u) override; + + /// @overload + std::shared_ptr solve( + const smith::FiniteElementState& u_guess, std::function equation, + std::function(const smith::FiniteElementState&)> jacobian) const override; + + /// @overload + std::shared_ptr solveAdjoint( + const smith::FiniteElementDual& u_bar, std::unique_ptr jacobian_transposed) const override; + + mutable std::unique_ptr mfem_solver; ///< linear solver + mutable std::unique_ptr mfem_preconditioner; ///< optionally used preconditioner +}; + +/// @brief Implementation of the DifferentiableSolver interface for the special case of nonlinear solves with linear +/// adjoint solves +class NonlinearDifferentiableSolver : public DifferentiableSolver { + public: + /// @brief Consruct from a smith nonlinear EquationSolver + NonlinearDifferentiableSolver(std::unique_ptr s); + + /// @overload + void completeSetup(const smith::FiniteElementState& u) override; + + /// @overload + std::shared_ptr solve( + const smith::FiniteElementState& u_guess, std::function equation, + std::function(const smith::FiniteElementState&)> jacobian) const override; + + /// @overload + std::shared_ptr solveAdjoint( + const smith::FiniteElementDual& u_bar, std::unique_ptr jacobian_transposed) const override; + + /// @overload + void clearMemory() const override; + + mutable std::unique_ptr J_; ///< stored linearized Jacobian matrix for memory reuse + mutable std::unique_ptr + nonlinear_solver_; ///< the nonlinear equation solver used for the forward pass +}; + +/// @brief Abstract interface to DifferentiableBlockSolver inteface. Each dfferenriable block solve should provide +/// both its forward solve and an adjoint solve +class DifferentiableBlockSolver { + public: + /// @brief destructor + virtual ~DifferentiableBlockSolver() {} + + using FieldT = FiniteElementState; ///< using + using FieldPtr = std::shared_ptr; ///< using + using FieldD = FiniteElementDual; ///< using + using DualPtr = std::shared_ptr; ///< using + using MatrixPtr = std::unique_ptr; ///< using + + /// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space + /// For these cases, it should be called before solve + virtual void completeSetup(const std::vector& us) = 0; + + /// @brief Solve a set of equations with a vector of FiniteElementState as unknown + /// @param u_guesses initial guess for solver + /// @param residuals std::vector for equations to be solved + /// @param jacobians std::vector> of std::function for evaluating the linearized Jacobians about the + /// current solution + /// @return std::vector of solution vectors (FiniteElementState) + virtual std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residuals, + std::function>(const std::vector&)> jacobians) const = 0; + + /// @brief Solve the (linear) adjoint set of equations with a vector of FiniteElementState as unknown + /// @param u_bars std::vector of right hand sides (rhs) for the solve + /// @param jacobian_transposed std::vector> of evaluated linearized adjoint space matrices + /// @return The adjoint vector of solution field + virtual std::vector solveAdjoint(const std::vector& u_bars, + std::vector>& jacobian_transposed) const = 0; + + /// @brief Interface option to clear memory between solves to avoid high-water mark memory usage. + virtual void clearMemory() const {} +}; + +/// @brief Implementation of the DifferentiableBlockSolver interface for the special case of linear solves with linear +/// adjoint solves +class LinearDifferentiableBlockSolver : public DifferentiableBlockSolver { + public: + /// @brief Construct from a linear solver and linear block precondition which may be used by the linear solver + LinearDifferentiableBlockSolver(std::unique_ptr s, std::unique_ptr p); + + /// @overload + void completeSetup(const std::vector& us) override; + + /// @overload + std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residuals, + std::function>(const std::vector&)> jacobians) const override; + + /// @overload + std::vector solveAdjoint(const std::vector& u_bars, + std::vector>& jacobian_transposed) const override; + + mutable std::unique_ptr mfem_solver; ///< stored mfem block solver + mutable std::unique_ptr mfem_preconditioner; ///< stored mfem block preconditioner +}; + +/// @brief Create a differentiable linear solver +/// @param linear_opts linear options struct +/// @param mesh mesh +std::shared_ptr buildDifferentiableLinearSolve(LinearSolverOptions linear_opts, + const smith::Mesh& mesh); + +/// @brief Create a differentiable nonlinear solver +/// @param nonlinear_opts nonlinear options struct +/// @param linear_opts linear options struct +/// @param mesh mesh +std::shared_ptr buildDifferentiableNonlinearSolve(NonlinearSolverOptions nonlinear_opts, + LinearSolverOptions linear_opts, + const smith::Mesh& mesh); + +} // namespace smith diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.cpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.cpp new file mode 100644 index 000000000..31005682f --- /dev/null +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.cpp @@ -0,0 +1,23 @@ +// 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) + +#include "smith/physics/mesh.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" + +namespace smith { + +DirichletBoundaryConditions::DirichletBoundaryConditions(const mfem::ParMesh& mfem_mesh, + mfem::ParFiniteElementSpace& space) + : bcs_(mfem_mesh), space_(space) +{ +} + +DirichletBoundaryConditions::DirichletBoundaryConditions(const Mesh& mesh, mfem::ParFiniteElementSpace& space) + : DirichletBoundaryConditions(mesh.mfemParMesh(), space) +{ +} + +} // namespace smith \ No newline at end of file diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp new file mode 100644 index 000000000..817a55898 --- /dev/null +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -0,0 +1,142 @@ +// 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 dirichlet_boundary_conditions.hpp + * + * @brief Contains DirichletBoundaryConditions class for interaction with the differentiable solve interfaces + */ + +#pragma once + +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" + +namespace smith { + +class Mesh; + +/// @brief A generic class for setting Dirichlet boundary conditions on arbitrary physics +class DirichletBoundaryConditions { + public: + /// @brief Construct from mfem::ParMesh + DirichletBoundaryConditions(const mfem::ParMesh& mfem_mesh, mfem::ParFiniteElementSpace& space); + + /// @brief Construct from smith::Mesh + DirichletBoundaryConditions(const Mesh& mesh, mfem::ParFiniteElementSpace& space); + + /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. + /// @param domain domain All dofs in this domain have boundary conditions applied to it. + /// @param components vectors of computents. The applied_displacement function returns the full vector, this + /// specifies which subset of those should have dirichlet boundary conditions applied. direction to apply boundary + /// condition to if the underlying field is a vector-field. + /// @param applied_displacement applied_displacement is a functor which takes time, and a + /// smith::tensor corresponding to the spatial coordinate. The functor must return a + /// smith::Tensor, where field_dim is the dimension of the vector space for the field. For example: + /// [](double t, smith::tensor X) { return smith::tensor{}; } + template + void setVectorBCs(const Domain& domain, std::vector components, AppliedDisplacementFunction applied_displacement) + { + int field_dim = space_.GetVDim(); + for (auto component : components) { + SLIC_ERROR_IF(component >= field_dim, + axom::fmt::format("Trying to set boundary conditions on a field with dim {}, using component {}", + field_dim, component)); + auto mfem_coefficient_function = [applied_displacement, component](const mfem::Vector& X_mfem, double t) { + auto X = make_tensor([&X_mfem](int k) { return X_mfem[k]; }); + return applied_displacement(t, X)[component]; + }; + + auto dof_list = domain.dof_list(&space_); + // scalar ldofs -> vector ldofs + space_.DofsToVDofs(static_cast(component), dof_list); + + auto component_disp_bdr_coef_ = std::make_shared(mfem_coefficient_function); + bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, static_cast(component)); + } + } + + /// @overload + template + void setVectorBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) + { + const int field_dim = space_.GetVDim(); + std::vector components(static_cast(field_dim)); + for (int component = 0; component < field_dim; ++component) { + components[static_cast(component)] = component; + } + setVectorBCs(domain, components, applied_displacement); + } + + /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. + /// @param domain domain All dofs in this domain have boundary conditions applied to it. + /// @param component component direction to apply boundary condition to if the underlying field is a vector-field. + /// @param applied_displacement applied_displacement is a functor which takes time, and a + /// smith::tensor corresponding to the spatial coordinate. The functor must return a double. For + /// example: [](double t, smith::tensor X) { return 1.0; } + template + void setScalarBCs(const Domain& domain, int component, AppliedDisplacementFunction applied_displacement) + { + const int field_dim = space_.GetVDim(); + SLIC_ERROR_IF(component >= field_dim, + axom::fmt::format("Trying to set boundary conditions on a field with dim {}, using component {}", + field_dim, component)); + auto mfem_coefficient_function = [applied_displacement](const mfem::Vector& X_mfem, double t) { + auto X = make_tensor([&X_mfem](int k) { return X_mfem[k]; }); + return applied_displacement(t, X); + }; + + auto dof_list = domain.dof_list(&space_); + // scalar ldofs -> vector ldofs + space_.DofsToVDofs(component, dof_list); + + auto component_disp_bdr_coef_ = std::make_shared(mfem_coefficient_function); + bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, component); + } + + /// @overload + template + void setScalarBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) + { + setScalarBCs(domain, 0, applied_displacement); + } + + /// @brief Constrain the dofs of a scalar field over a domain + template + void setFixedScalarBCs(const Domain& domain, int component = 0) + { + setScalarBCs(domain, component, [](auto, auto) { return 0.0; }); + } + + /// @brief Constrain the vector dofs over a domain corresponding to a subset of the vector components + template + void setFixedVectorBCs(const Domain& domain, std::vector components) + { + setVectorBCs(domain, components, [](auto, auto) { return smith::tensor{}; }); + } + + /// @brief Constrain all the vector dofs over a domain + template + void setFixedVectorBCs(const Domain& domain) + { + const int field_dim = space_.GetVDim(); + SLIC_ERROR_IF(field_dim != spatial_dim, + "Vector boundary conditions current only work if they match the spatial dimension"); + std::vector components(static_cast(field_dim)); + for (int component = 0; component < field_dim; ++component) { + components[static_cast(component)] = component; + } + setFixedVectorBCs(domain, components); + } + + /// @brief Return the smith BoundaryConditionManager + const smith::BoundaryConditionManager& getBoundaryConditionManager() const { return bcs_; } + + private: + smith::BoundaryConditionManager bcs_; ///< boundary condition manager that does the heavy lifting + mfem::ParFiniteElementSpace& space_; ///< save the space for the field which will be constrained +}; + +} // namespace smith \ No newline at end of file diff --git a/src/smith/differentiable_numerics/differentiable_utils.cpp b/src/smith/differentiable_numerics/evaluate_objective.cpp similarity index 97% rename from src/smith/differentiable_numerics/differentiable_utils.cpp rename to src/smith/differentiable_numerics/evaluate_objective.cpp index 6d0cd288d..cb97fcfad 100644 --- a/src/smith/differentiable_numerics/differentiable_utils.cpp +++ b/src/smith/differentiable_numerics/evaluate_objective.cpp @@ -5,7 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include "gretl/data_store.hpp" -#include "smith/differentiable_numerics/differentiable_utils.hpp" +#include "smith/differentiable_numerics/evaluate_objective.hpp" namespace smith { diff --git a/src/smith/differentiable_numerics/evaluate_objective.hpp b/src/smith/differentiable_numerics/evaluate_objective.hpp new file mode 100644 index 000000000..4f50f5050 --- /dev/null +++ b/src/smith/differentiable_numerics/evaluate_objective.hpp @@ -0,0 +1,30 @@ +// 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 evaluate_objective.hpp + * + * @brief Methods for evaluating objective functions and tracking these operations on the gretl graph with a custom vjp + */ + +#pragma once + +#include "gretl/double_state.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/physics/scalar_objective.hpp" + +namespace smith { + +/// @brief Evaluates a DoubleState using a provided ScalarObjective instance, and the input arguments to that objective. +/// This operation is tracked on the gretl graph. +DoubleState evaluateObjective(std::shared_ptr objective, const TimeInfo& time_info, + const FieldState& shape_disp, const std::vector& inputs); + +/// @brief operation is tracked on the gretl graph. +DoubleState evaluateObjective(std::shared_ptr objective, const FieldState& shape_disp, + const std::vector& inputs); + +} // namespace smith diff --git a/src/smith/differentiable_numerics/field_state.cpp b/src/smith/differentiable_numerics/field_state.cpp index 173d59f9d..9a18dceec 100644 --- a/src/smith/differentiable_numerics/field_state.cpp +++ b/src/smith/differentiable_numerics/field_state.cpp @@ -26,7 +26,7 @@ FieldState square(const FieldState& state) state); } -gretl::State inner_product(const FieldState& a, const FieldState& b) +gretl::State innerProduct(const FieldState& a, const FieldState& b) { return gretl::create_state( gretl::defaultInitializeZeroDual(), @@ -38,6 +38,18 @@ gretl::State inner_product(const FieldState& a, const FieldState& b) a, b); } +gretl::State innerProduct(const ReactionState& a, const ReactionState& b) +{ + return gretl::create_state( + gretl::defaultInitializeZeroDual(), + [](FEDualPtr A, FEDualPtr B) { return smith::innerProduct(*A, *B); }, + [](FEDualPtr A, FEDualPtr B, double, FEFieldPtr& A_, FEFieldPtr& B_, double product_) { + A_->Add(product_, *B); + B_->Add(product_, *A); + }, + a, b); +} + FieldState axpby(double a, const FieldState& x, double b, const FieldState& y) { auto z = x.clone({x, y}); diff --git a/src/smith/differentiable_numerics/field_state.hpp b/src/smith/differentiable_numerics/field_state.hpp index 3167080e8..cf0a70e6f 100644 --- a/src/smith/differentiable_numerics/field_state.hpp +++ b/src/smith/differentiable_numerics/field_state.hpp @@ -20,8 +20,8 @@ namespace smith { using FEFieldPtr = std::shared_ptr; ///< typedef using FEDualPtr = std::shared_ptr; ///< typedef using FieldState = gretl::State; ///< typedef +using ReactionState = gretl::State; ///< typedef using FieldVecState = gretl::State, std::vector>; ///< typedef -using ResultantState = gretl::State; ///< typedef using DoubleState = gretl::State; ///< typedef /// @brief functor which takes a std::shared_ptr, and returns a zero-valued @@ -59,26 +59,29 @@ FieldState create_field_state(gretl::DataStore& dataStore, function_space space, return create_field_state(dataStore, f); } -/// @brief initialize on the gretl::DataStore a ResultantState with values from s -inline ResultantState create_field_resultant(gretl::DataStore& dataStore, const smith::FEDualPtr& s) +/// @brief initialize on the gretl::DataStore a ReactionState with values from s +inline ReactionState create_reaction_state(gretl::DataStore& dataStore, const smith::FEDualPtr& s) { return dataStore.create_state(s, zero_state_from_dual()); } -/// @brief initialize on the gretl::DataStore a ResultantState from a FiniteElementDual of given space, name and mesh. +/// @brief initialize on the gretl::DataStore a ReactionState from a FiniteElementDual of given space, name and mesh. template -ResultantState create_field_resultant(gretl::DataStore& dataStore, function_space space, const std::string& name, - const std::string& mesh_tag) +ReactionState create_reaction_state(gretl::DataStore& dataStore, function_space space, const std::string& name, + const std::string& mesh_tag) { auto f = std::make_shared(StateManager::newDual(space, name, mesh_tag)); - return create_field_resultant(dataStore, f); + return create_reaction_state(dataStore, f); } /// @brief gretl-function to square (x^2) every component of the Field FieldState square(const FieldState& state); /// @brief gretl-function to compute the inner product (vector l2-norm) of a and b -gretl::State inner_product(const FieldState& a, const FieldState& b); +gretl::State innerProduct(const FieldState& a, const FieldState& b); + +/// @brief gretl-function to compute the inner product (vector l2-norm) of a and b +gretl::State innerProduct(const ReactionState& a, const ReactionState& b); /// @brief gretl-function to compute a*x + b*y FieldState axpby(double a, const FieldState& x, double b, const FieldState& y); diff --git a/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.cpp b/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.cpp index 2aa93f0f7..7441469a4 100644 --- a/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.cpp +++ b/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.cpp @@ -33,10 +33,9 @@ FieldState applyZeroBoundaryConditions(const FieldState& s, const BoundaryCondit return s_bc.finalize(); } -std::vector LumpedMassExplicitNewmarkStateAdvancer::advanceState(const FieldState& shape_disp, - const std::vector& states_in, - const std::vector& params, - const TimeInfo& time_info) const +std::vector LumpedMassExplicitNewmarkStateAdvancer::advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states_in, + const std::vector& params) const { SMITH_MARK_FUNCTION; SLIC_ERROR_IF(states_in.size() != 3, "ExplicitNewmark is a 2nd order time integrator requiring 3 states."); diff --git a/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.hpp b/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.hpp index 2fd7e512f..1dd292617 100644 --- a/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.hpp +++ b/src/smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.hpp @@ -36,8 +36,9 @@ class LumpedMassExplicitNewmarkStateAdvancer : public StateAdvancer { } /// @overload - std::vector advanceState(const FieldState& shape_disp, const std::vector& states, - const std::vector& params, const TimeInfo& time_info) const override; + std::vector advanceState(const TimeInfo& time_info, const FieldState& shape_disp, + const std::vector& states, + const std::vector& params) const override; private: const std::shared_ptr residual_eval_; ///< weak form to evaluate mechanical forces diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp new file mode 100644 index 000000000..7d2ed606a --- /dev/null +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -0,0 +1,574 @@ +#include "smith/physics/weak_form.hpp" +#include "smith/physics/field_types.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" + +namespace smith { + +/// @brief apply boundary conditions +void applyBoundaryConditions(double time, const smith::BoundaryConditionManager* bc_manager, + smith::FEFieldPtr& primal_field, const smith::FEFieldPtr& bc_field_ptr) +{ + auto constrained_dofs = bc_manager->allEssentialTrueDofs(); + if (bc_field_ptr) { + for (int i = 0; i < constrained_dofs.Size(); i++) { + int j = constrained_dofs[i]; + (*primal_field)[j] = (*bc_field_ptr)(j); + } + } else { + for (auto& bc : bc_manager->essentials()) { + bc.setDofs(*primal_field, time); + } + } +} + +/// @brief Solve a nonlinear system of equations as defined by the weak form +/// @param residual_eval The weak form which defines the equations to be solved +/// @param shape_disp The mesh-morphed shape displacement +/// @param states The time varying states as inputs to the weak form +/// @param params The fixed field parameters as inputs to the weak form +/// @param state_update_weights Specifies how to blend the arguments of the weak form into a unique unknown. The +/// primary unknown, p, is states[primal_solve_state_index]. However, other arguments to the weak form may also depend +/// on this unknown. The state_update_weights specify how this is done. For a nonzero weight, the field argument i is +/// assumed to change as states[i] + (p - p_initial) * state_update_weights[i]. +/// @param primal_solve_state_index Index specifying which of the states is the primary unknown. +/// @param dirichlet_state_index Index specifying which field has the Dirichlet boundary conditions applied to it. +/// Typically this will be the same as the primal_solve_state_index, but it can be different. An example is explicit +/// dynamics, where the unknown is the acceleration, and the boundary conditions are often applied directly to the +/// displacement field. +/// @param time_info Timestep information (time, dt, cycle) +/// @param solver The differentiable, potentially nonlinear, equation solver used to solve the system of equations +/// @param bc_manager Holds information about which degrees of freedom (DOFS) +/// @param bc_field A field which holds to desired values for the fixed degrees of freedom as specified by the +/// bc_manager +/// @return The field solution to the weak form +FieldState nonlinearSolve(const WeakForm* residual_eval, const FieldState& shape_disp, + const std::vector& states, const std::vector& params, + const std::vector& state_update_weights, size_t primal_solve_state_index, + size_t dirichlet_state_index, const TimeInfo& time_info, const DifferentiableSolver* solver, + const BoundaryConditionManager* bc_manager, const FieldState* bc_field = nullptr) +{ + SMITH_MARK_FUNCTION; + // there should be one less input state, as the higher time derivative term (e.g., acceleration), does not have a + // predictor + SLIC_ERROR_IF(states.size() != state_update_weights.size(), "State and state weight fields are inconsistent"); + SLIC_ERROR_IF(state_update_weights[primal_solve_state_index] != 1.0, "Primal state must have a weight of 1.0"); + + std::vector allFields; + for (auto& s : states) allFields.push_back(s); + for (auto& p : params) allFields.push_back(p); + allFields.push_back(shape_disp); + + bool have_bc_field = bc_field; + if (have_bc_field) { + allFields.push_back(*bc_field); + } + + FieldState sol = states[primal_solve_state_index].clone(allFields); + + sol.set_eval([=](const gretl::UpstreamStates& inputs, gretl::DownstreamState& output) { + SMITH_MARK_BEGIN("solve forward"); + + const size_t num_states = state_update_weights.size(); + + std::vector non_primal_to_state_index; + for (size_t i = 0; i < num_states; ++i) { + if (i != primal_solve_state_index) { + non_primal_to_state_index.push_back(i); + } + } + + const size_t num_extra_args = have_bc_field ? 2 : 1; + const size_t num_fields = inputs.size() - num_extra_args; + + std::vector corrected_fields(num_fields); + for (size_t field_index = 0; field_index < num_fields; ++field_index) { + if (field_index < state_update_weights.size() && state_update_weights[field_index] != 0.0) { + corrected_fields[field_index] = std::make_shared(*inputs[field_index].get()); + } else { + corrected_fields[field_index] = inputs[field_index].get(); + } + } + + const FEFieldPtr shape_disp_ptr = inputs[num_fields].get(); + + FEFieldPtr bc_field_ptr; + if (have_bc_field) { + bc_field_ptr = inputs[num_fields + num_extra_args - 1].get(); + } + + FEFieldPtr s0 = corrected_fields[primal_solve_state_index]; + FEFieldPtr s = std::make_shared(s0->space(), "s"); + + if (bc_manager && (dirichlet_state_index == primal_solve_state_index)) { + applyBoundaryConditions(time_info.time(), bc_manager, s0, bc_field_ptr); + } + + s = solver->solve( + *s0, // initial guess when solving for the primal index field + [=](const FiniteElementState& s_) { + FEFieldPtr primal_field = corrected_fields[primal_solve_state_index]; + *primal_field = s_; + + if (bc_manager && (dirichlet_state_index == primal_solve_state_index)) { + applyBoundaryConditions(time_info.time(), bc_manager, primal_field, bc_field_ptr); + } + + for (size_t corrected_field_index : non_primal_to_state_index) { + if (state_update_weights[corrected_field_index] != 0.0) { + *corrected_fields[corrected_field_index] = *inputs[corrected_field_index].get(); + corrected_fields[corrected_field_index]->Add(state_update_weights[corrected_field_index], *primal_field); + corrected_fields[corrected_field_index]->Add(-state_update_weights[corrected_field_index], *s0); + } + } + + auto r = residual_eval->residual(time_info, shape_disp_ptr.get(), getConstFieldPointers(corrected_fields)); + + if (bc_manager) { + if (dirichlet_state_index == primal_solve_state_index) { + auto constrained_dofs = bc_manager->allEssentialTrueDofs(); + r.SetSubVector(constrained_dofs, 0.0); + } + } + + return r; + }, + [=](const FiniteElementState& s_) { + FEFieldPtr primal_field = corrected_fields[primal_solve_state_index]; + *primal_field = s_; + + if (bc_manager && (dirichlet_state_index == primal_solve_state_index)) { + applyBoundaryConditions(time_info.time(), bc_manager, primal_field, bc_field_ptr); + } + + for (size_t corrected_field_index : non_primal_to_state_index) { + if (state_update_weights[corrected_field_index] != 0.0) { + *corrected_fields[corrected_field_index] = *inputs[corrected_field_index].get(); + corrected_fields[corrected_field_index]->Add(state_update_weights[corrected_field_index], *primal_field); + corrected_fields[corrected_field_index]->Add(-state_update_weights[corrected_field_index], *s0); + } + } + + auto J = residual_eval->jacobian(time_info, shape_disp_ptr.get(), getConstFieldPointers(corrected_fields), + state_update_weights); + + if (bc_manager) { + if (dirichlet_state_index == primal_solve_state_index) { + J->EliminateBC(bc_manager->allEssentialTrueDofs(), mfem::Operator::DiagonalPolicy::DIAG_ONE); + } + } + return J; + }); + + output.set(s); + + SMITH_MARK_END("solve forward"); + }); + + sol.set_vjp([=](gretl::UpstreamStates& inputs, const gretl::DownstreamState& output) { + SMITH_MARK_BEGIN("solve reverse"); + const FEFieldPtr s = output.get(); // get the final solution + const FEDualPtr s_dual = output.get_dual(); // get the dual load + + const size_t num_states = state_update_weights.size(); + + std::vector non_primal_to_state_index; + for (size_t i = 0; i < num_states; ++i) { + if (i != primal_solve_state_index) { + non_primal_to_state_index.push_back(i); + } + } + + const size_t num_extra_args = have_bc_field ? 2 : 1; + const size_t num_fields = inputs.size() - num_extra_args; + + std::vector corrected_fields(num_fields); + for (size_t field_index = 0; field_index < num_fields; ++field_index) { + if (field_index < state_update_weights.size() && state_update_weights[field_index] != 0.0) { + corrected_fields[field_index] = std::make_shared(*inputs[field_index].get()); + } else { + corrected_fields[field_index] = inputs[field_index].get(); + } + } + + const FEFieldPtr shape_disp_ptr = inputs[num_fields].get(); + + const FEFieldPtr s0 = inputs[primal_solve_state_index].get(); + + *corrected_fields[primal_solve_state_index] = *s; + for (size_t corrected_field_index : non_primal_to_state_index) { + if (state_update_weights[corrected_field_index] != 0.0) { + corrected_fields[corrected_field_index]->Add(state_update_weights[corrected_field_index], *s); + corrected_fields[corrected_field_index]->Add(-state_update_weights[corrected_field_index], *s0); + } + } + + solver->clearMemory(); + auto J = residual_eval->jacobian(time_info, shape_disp_ptr.get(), getConstFieldPointers(corrected_fields), + state_update_weights, {}); + + if (bc_manager) { + if (dirichlet_state_index == primal_solve_state_index) { + J->EliminateBC(bc_manager->allEssentialTrueDofs(), mfem::Operator::DiagonalPolicy::DIAG_ONE); + } + } + + auto J_T = std::unique_ptr(J->Transpose()); + J.reset(); + + auto s_adjoint_ptr = solver->solveAdjoint(*s_dual, std::move(J_T)); + + if (bc_manager) { + if (dirichlet_state_index == primal_solve_state_index) { + s_adjoint_ptr->SetSubVector(bc_manager->allEssentialTrueDofs(), 0.0); + } + } + + *s_adjoint_ptr *= -1.0; + + std::vector field_sensitivities(num_fields, nullptr); + FEDualPtr shape_disp_sensitivity = inputs[num_fields].get_dual(); + for (size_t state_index = 0; state_index < num_states; ++state_index) { + field_sensitivities[state_index] = inputs[state_index].get_dual().get(); + } + for (size_t param_index = num_states; param_index < num_fields; ++param_index) { + field_sensitivities[param_index] = inputs[param_index].get_dual().get(); + } + + auto primal_sensitivity = std::make_shared(*field_sensitivities[primal_solve_state_index]); + field_sensitivities[primal_solve_state_index] = primal_sensitivity.get(); + *field_sensitivities[primal_solve_state_index] = *s_dual; + + residual_eval->vjp(time_info, shape_disp_ptr.get(), getConstFieldPointers(corrected_fields), {}, + s_adjoint_ptr.get(), shape_disp_sensitivity.get(), field_sensitivities, {}); + + if (bc_manager && have_bc_field && dirichlet_state_index == primal_solve_state_index) { + auto bc_dual_ptr = inputs[num_fields + num_extra_args - 1].get_dual(); + field_sensitivities[primal_solve_state_index]->SetSubVectorComplement(bc_manager->allEssentialTrueDofs(), 0.0); + *bc_dual_ptr += *field_sensitivities[primal_solve_state_index]; + } + + SMITH_MARK_END("solve reverse"); + }); + + sol.finalize(); + + return sol; +} + +FieldState solve(const WeakForm& weak_form, const FieldState& shape_disp, const std::vector& states, + const std::vector& params, const TimeInfo& time_info, const DifferentiableSolver& solver, + const DirichletBoundaryConditions& bcs, size_t unknown_state_index) +{ + std::vector state_update_weights(states.size(), 0.0); + state_update_weights[unknown_state_index] = 1.0; + return nonlinearSolve(&weak_form, shape_disp, states, params, state_update_weights, unknown_state_index, + unknown_state_index, time_info, &solver, &bcs.getBoundaryConditionManager()); +} + +std::vector block_solve(const std::vector& residual_evals, + const std::vector> block_indices, const FieldState& shape_disp, + const std::vector>& states, + const std::vector>& params, const TimeInfo& time_info, + const DifferentiableBlockSolver* solver, + const std::vector bc_managers) +{ + SMITH_MARK_FUNCTION; + size_t num_rows_ = residual_evals.size(); + + SLIC_ERROR_IF(num_rows_ != block_indices.size(), "Block indices size not consistent with number of residual rows"); + SLIC_ERROR_IF(num_rows_ != states.size(), + "Number of state input vectors not consistent with number of residual rows"); + SLIC_ERROR_IF(num_rows_ != params.size(), + "Number of parameter input vectors not consistent with number of residual rows"); + SLIC_ERROR_IF(num_rows_ != bc_managers.size(), + "Number of boundary condition manager not consistent with number of residual rows"); + + for (size_t r = 0; r < num_rows_; ++r) { + SLIC_ERROR_IF(num_rows_ != block_indices[r].size(), "All block index rows must have the same number of columns"); + } + + std::vector num_state_inputs; + std::vector allFields; + for (auto& ss : states) { + num_state_inputs.push_back(ss.size()); + for (auto& s : ss) { + allFields.push_back(s); + } + } + std::vector num_param_inputs; + for (auto& ps : params) { + num_param_inputs.push_back(ps.size()); + for (auto& p : ps) { + allFields.push_back(p); + } + } + allFields.push_back(shape_disp); + + struct ZeroDualVectors { + std::vector operator()(const std::vector& fs) + { + std::vector ds(fs.size()); + for (size_t i = 0; i < fs.size(); ++i) { + ds[i] = std::make_shared(fs[i]->space(), fs[i]->name() + "_dual"); + } + return ds; + } + }; + + FieldVecState sol = + shape_disp.create_state, std::vector>(allFields, ZeroDualVectors()); + + sol.set_eval([=](const gretl::UpstreamStates& upstreams, gretl::DownstreamState& downstream) { + SMITH_MARK_BEGIN("solve forward"); + const size_t num_rows = num_state_inputs.size(); + std::vector> input_fields(num_rows); + SLIC_ERROR_IF(num_rows != num_param_inputs.size(), "row count for params and columns are inconsistent"); + + // The order of inputs in upstreams seems to be: + // states of residual 0 -> states of residual 1 -> params of residual 0 -> params of residual 1 + size_t field_count = 0; + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t state_i = 0; state_i < num_state_inputs[row_i]; ++state_i) { + input_fields[row_i].push_back(upstreams[field_count++].get()); + } + } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t param_i = 0; param_i < num_param_inputs[row_i]; ++param_i) { + input_fields[row_i].push_back(upstreams[field_count++].get()); + } + } + + std::vector diagonal_fields(num_rows); + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][block_indices[row_i][row_i]]); + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + } + } + + const FEFieldPtr shape_disp_ptr = upstreams[field_count].get(); + + auto eval_residuals = [=](const std::vector& unknowns) { + SLIC_ERROR_IF(unknowns.size() != num_rows, + "block solver unknown size must match the number or residuals in block_solve"); + std::vector residuals(num_rows); + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + FEFieldPtr primal_field_row_i = diagonal_fields[row_i]; + *primal_field_row_i = *unknowns[row_i]; + applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + residuals[row_i] = residual_evals[row_i]->residual(time_info, shape_disp_ptr.get(), + getConstFieldPointers(input_fields[row_i])); + residuals[row_i].SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); + } + + return residuals; + }; + + auto eval_jacobians = [=](const std::vector& unknowns) { + SLIC_ERROR_IF(unknowns.size() != num_rows, + "block solver unknown size must match the number or residuals in block_solve"); + std::vector>> jacobians(num_rows); + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + FEFieldPtr primal_field_row_i = diagonal_fields[row_i]; + *primal_field_row_i = *unknowns[row_i]; + applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + std::vector row_field_inputs = input_fields[row_i]; + std::vector tangent_weights(row_field_inputs.size(), 0.0); + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + size_t field_index_to_diff = block_indices[row_i][col_j]; + tangent_weights[field_index_to_diff] = 1.0; + auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), + getConstFieldPointers(row_field_inputs), tangent_weights); + jacobians[row_i].emplace_back(std::move(jac_ij)); + tangent_weights[field_index_to_diff] = 0.0; + } + } + + // Apply BCs to the block system + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + mfem::HypreParMatrix* Jii = + jacobians[row_i][row_i]->EliminateRowsCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jii; + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + if (col_j != row_i) { + jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); + mfem::HypreParMatrix* Jji = + jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jji; + } + } + } + + return jacobians; + }; + + diagonal_fields = solver->solve(diagonal_fields, eval_residuals, eval_jacobians); + + downstream.set, std::vector>(diagonal_fields); + + SMITH_MARK_END("solve forward"); + }); + + sol.set_vjp([=](gretl::UpstreamStates& upstreams, const gretl::DownstreamState& downstream) { + SMITH_MARK_BEGIN("solve reverse"); + const std::vector s = downstream.get>(); // get the final solution + const std::vector s_dual = + downstream.get_dual, std::vector>(); // get the dual load + + const size_t num_rows = num_state_inputs.size(); + SLIC_ERROR_IF(s_dual.size() != num_rows, + "block solver vjp downstream size must match the number or residuals in block_solve"); + + std::vector> input_fields(num_rows); + size_t field_count = 0; + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t state_i = 0; state_i < num_state_inputs[row_i]; ++state_i) { + input_fields[row_i].push_back(upstreams[field_count++].get()); + } + } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t param_i = 0; param_i < num_param_inputs[row_i]; ++param_i) { + input_fields[row_i].push_back(upstreams[field_count++].get()); + } + } + + // if the field is a primal variable we solved before, + // make a copy so we don't accidentally override the original copy + std::vector diagonal_fields(num_rows); + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][block_indices[row_i][row_i]]); + *diagonal_fields[row_i] = *s[row_i]; + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + } + } + + const FEFieldPtr shape_disp_ptr = upstreams[field_count].get(); + + // I'm not sure this will be the right timestamp to apply boundary condition during backward propagation + // Need to double check for time-dependent boundary conditions + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + FEFieldPtr primal_field_row_i = diagonal_fields[row_i]; + applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); + } + + solver->clearMemory(); + + std::vector>> jacobians(num_rows); + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + std::vector row_field_inputs = input_fields[row_i]; + std::vector tangent_weights(row_field_inputs.size(), 0.0); + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + size_t field_index_to_diff = block_indices[row_i][col_j]; + tangent_weights[field_index_to_diff] = 1.0; + auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), + getConstFieldPointers(row_field_inputs), tangent_weights); + jacobians[row_i].emplace_back(std::move(jac_ij)); + tangent_weights[field_index_to_diff] = 0.0; + } + } + + // Apply BCs to the block system + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + s_dual[row_i]->SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); + + mfem::HypreParMatrix* Jii = + jacobians[row_i][row_i]->EliminateRowsCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jii; + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + if (col_j != row_i) { + jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); + mfem::HypreParMatrix* Jji = + jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jji; + } + } + } + + // Take the transpose of the block system + std::vector>> jacobians_T(num_rows); + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + jacobians_T[col_j].emplace_back(std::unique_ptr(jacobians[row_i][col_j]->Transpose())); + } + } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + jacobians[row_i][col_j].reset(); + } + } + + std::vector adjoint_fields(num_rows); + adjoint_fields = solver->solveAdjoint(s_dual, jacobians_T); + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + *adjoint_fields[row_i] *= -1.0; + } + + // Update sensitivities + std::vector> field_sensitivities(num_rows); + FEDualPtr shape_disp_sensitivity = upstreams[field_count].get_dual(); + size_t dual_index = 0; + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t state_i = 0; state_i < num_state_inputs[row_i]; ++state_i) { + field_sensitivities[row_i].push_back(upstreams[dual_index++].get_dual()); + } + } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t param_i = 0; param_i < num_param_inputs[row_i]; ++param_i) { + field_sensitivities[row_i].push_back(upstreams[dual_index++].get_dual()); + } + } + SLIC_ERROR_IF(field_count != dual_index, "Number of sensitivities must equal to number of upstreams"); + + // No sensitivity needed for primal fields + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + for (size_t col_j = 0; col_j < num_rows; ++col_j) { + field_sensitivities[row_i][block_indices[row_i][col_j]] = nullptr; + } + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + residual_evals[row_i]->vjp(time_info, shape_disp_ptr.get(), getConstFieldPointers(input_fields[row_i]), {}, + adjoint_fields[row_i].get(), shape_disp_sensitivity.get(), + getFieldPointers(field_sensitivities[row_i]), {}); + } + + SMITH_MARK_END("solve reverse"); + }); + + sol.finalize(); + + std::vector results; + for (size_t i = 0; i < num_rows_; ++i) { + FieldState s = gretl::create_state( + zero_dual_from_state(), + [i](const std::vector& sols) { return std::make_shared(*sols[i]); }, + [i](const std::vector&, const FEFieldPtr&, std::vector& sols_, + const FEDualPtr& output_) { *sols_[i] += *output_; }, + sol); + + results.emplace_back(s); + } + + return results; +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp new file mode 100644 index 000000000..07027b274 --- /dev/null +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -0,0 +1,66 @@ +// Copyright (c) 2019-2024, 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 nonlinear_solve.hpp + * + * @brief Methods for solving systems of equations as given by WeakForms. Tracks these operations on the gretl graph + * with a custom vjp. + */ + +#pragma once + +#include +#include "smith/differentiable_numerics/field_state.hpp" + +namespace smith { + +class WeakForm; +class DifferentiableSolver; +class DifferentiableBlockSolver; +class BoundaryConditionManager; +class DirichletBoundaryConditions; + +/// @brief Solve a nonlinear system of equations as defined by the weak form, assuming that the field indexed by +/// unknown_index is the unknown field +/// @param residual_eval The weak form which defines the equations to be solved +/// @param shape_disp The mesh-morphed shape displacement +/// @param states The time varying states as inputs to the weak form +/// @param params All fixed fields pass to the weak form +/// @param time_info Timestep information (time, dt, cycle) +/// @param solver The differentiable, potentially nonlinear, equation solver used to solve the system of equations +/// @param bcs Holds information about which degrees of freedom (DOFS), and has the information about the time and space +/// varying values for the boundary conditions +/// @param unknown_state_index +/// @return The field solution to the weak form +FieldState solve(const WeakForm& residual_eval, const FieldState& shape_disp, const std::vector& states, + const std::vector& params, const TimeInfo& time_info, const DifferentiableSolver& solver, + const DirichletBoundaryConditions& bcs, size_t unknown_state_index = 0); + +/// @brief Solve a block nonlinear system of equations as defined by the vector of weak form +/// @param residual_evals Vector of weak forms which define the equations to be solved +/// @param block_indices Matrix of index arguments specifying where in each WeakForm the unknown fields are passed in. +/// Example: for a 2 weak-form system, with weak-forms, r1, r2 +/// r1(a,b,c) +/// r2(b,d,e,a) +// with unknowns (with respect to the solver) being a, and b. +// r1 has unknowns a,b in the ‘slots’ 0, 1 +// r2 has unknowns a,b, in the ‘slots’ 3,0 +/// @param shape_disp The mesh-morphed shape displacement +/// @param states The time varying states as inputs to the weak form +/// @param params The fixed field parameters as inputs to the weak form +/// @param time_info Timestep information (time, dt, cycle) +/// @param solver The differentiable, potentially nonlinear, equation solver used to solve the system of equations +/// @param bc_managers Holds information about which degrees of freedom (DOFS) +/// @return Vector of field solutions satisfying the weak forms +std::vector block_solve(const std::vector& residual_evals, + const std::vector> block_indices, const FieldState& shape_disp, + const std::vector>& states, + const std::vector>& params, const TimeInfo& time_info, + const DifferentiableBlockSolver* solver, + const std::vector bc_managers); + +} // namespace smith diff --git a/src/smith/differentiable_numerics/tests/paraview_helper.hpp b/src/smith/differentiable_numerics/paraview_writer.hpp similarity index 66% rename from src/smith/differentiable_numerics/tests/paraview_helper.hpp rename to src/smith/differentiable_numerics/paraview_writer.hpp index e2124789a..5f8081fd4 100644 --- a/src/smith/differentiable_numerics/tests/paraview_helper.hpp +++ b/src/smith/differentiable_numerics/paraview_writer.hpp @@ -5,7 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) /** - * @file paraview_helper.hpp + * @file paraview_writer.hpp * * @brief Simple paraview field output methods */ @@ -20,21 +20,21 @@ namespace smith { +/// @brief Class which interactions with ParaViewDataCollection to write arbitrary field results to disk. This allows +/// output independent of a particular BasePhysics. class ParaviewWriter { public: - using StateVecs = std::vector >; - using DualVecs = std::vector >; - - ParaviewWriter(std::unique_ptr pv_, const StateVecs& states_) - : pv(std::move(pv_)), states(states_) - { - } + using StateVecs = std::vector >; ///< using + /// Construct from ParaViewDataCollection, a vector of shared_ptr to FiniteElementState, and vector of shared_ptr to + /// FiniteElementState which dual fields will be copied into. ParaviewWriter(std::unique_ptr pv_, const StateVecs& states_, const StateVecs& duals_) : pv(std::move(pv_)), states(states_), dual_states(duals_) { } + /// @brief write paraview output from vector of finite element states. states must be passed in with a consistent + /// order as how the ParaviewWriter was constructed (consistent order of spaces) void write(size_t step, double time, const std::vector& current_states) { SMITH_MARK_FUNCTION; @@ -51,6 +51,8 @@ class ParaviewWriter { pv->Save(); } + /// @brief write paraview output from vector of finite element duals. duals must be passed in with a consistent order + /// as how the ParaviewWriter was constructed (consistent order of spaces) void write(size_t step, double time, const std::vector& current_duals) { SMITH_MARK_FUNCTION; @@ -67,6 +69,9 @@ class ParaviewWriter { pv->Save(); } + /// @brief write paraview output from vector of FieldState. These must be passed in with a consistent order as how the + /// ParaviewWriter was constructed (consistent order of spaces). Both the field, and its dual/reaction will be + /// written. void write(int step, double time, const std::vector& current_fields) { SMITH_MARK_FUNCTION; @@ -87,13 +92,21 @@ class ParaviewWriter { pv->Save(); } + /// overload + void write(size_t step, double time, const std::vector& current_fields) + { + write(static_cast(step), time, current_fields); + } + private: std::unique_ptr pv; StateVecs states; StateVecs dual_states; }; -inline auto createParaviewOutput(const smith::Mesh& mesh, const std::vector& states, +/// @brief Creates a ParaviewWriter from a mesh, vector of FieldState, and the name of the output paraview file. File +/// will be in directory filename/filename.pvd. +inline auto createParaviewWriter(const smith::Mesh& mesh, const std::vector& states, std::string output_name) { if (output_name == "") { @@ -105,7 +118,8 @@ inline auto createParaviewOutput(const smith::Mesh& mesh, const std::vector(&mesh.mfemParMesh()); auto paraview_dc = std::make_unique(output_name, non_const_mesh); - int max_order_in_fields = 0; + // visualization order has to be at least 1 for paraview (because there is no zero order mesh) + int max_order_in_fields = 1; // Find the maximum polynomial order in the physics module's states for (const auto& fstate : states) { @@ -129,7 +143,9 @@ inline auto createParaviewOutput(const smith::Mesh& mesh, const std::vector& states, +/// @brief Creates a ParaviewWriter from an mfem::ParMesh, vector of FiniteElementState pointers, and the name of the +/// output paraview file. File will be in directory filename/filename.pvd. +inline auto createParaviewWriter(const mfem::ParMesh& mesh, const std::vector& states, std::string output_name) { if (output_name == "") { @@ -143,7 +159,8 @@ inline auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector(&mesh); auto paraview_dc = std::make_unique(output_name, non_const_mesh); - int max_order_in_fields = 0; + // visualization order has to be at least 1 for paraview (because there is no zero order mesh) + int max_order_in_fields = 1; // Find the maximum polynomial order in the physics module's states for (const auto& state : output_states) { @@ -160,60 +177,4 @@ inline auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector& duals, - std::string output_name) -{ - if (output_name == "") { - output_name = "default"; - } - - auto non_const_mesh = const_cast(&mesh); - auto paraview_dc = std::make_unique(output_name, non_const_mesh); - int max_order_in_fields = 0; - - ParaviewWriter::StateVecs output_duals; - // Find the maximum polynomial order in the physics module's states - for (const auto& dual : duals) { - output_duals.push_back(std::make_shared(dual->space(), dual->name())); - paraview_dc->RegisterField(dual->name(), &output_duals.back()->gridFunction()); - max_order_in_fields = std::max(max_order_in_fields, dual->space().GetOrder(0)); - } - - // Set the options for the paraview output files - paraview_dc->SetLevelsOfDetail(max_order_in_fields); - paraview_dc->SetHighOrderOutput(true); - paraview_dc->SetDataFormat(mfem::VTKFormat::BINARY); - paraview_dc->SetCompression(true); - - return ParaviewWriter(std::move(paraview_dc), {}, output_duals); -} - -inline std::vector get_field_ptrs_for_output( - smith::FieldState shape_disp, const std::vector& states_in, - std::vector other_fields = {}) -{ - std::vector states_out{shape_disp.get().get()}; - for (auto& s : states_in) { - states_out.push_back(s.get().get()); - } - for (auto& s : other_fields) { - states_out.push_back(s.get().get()); - } - return states_out; -}; - -inline std::vector get_dual_ptrs_for_output( - smith::FieldState shape_disp, const std::vector& states_in, - std::vector other_fields = {}) -{ - std::vector duals_out{shape_disp.get_dual().get()}; - for (auto& s : states_in) { - duals_out.push_back(s.get_dual().get()); - } - for (auto& s : other_fields) { - duals_out.push_back(s.get_dual().get()); - } - return duals_out; -}; - } // namespace smith diff --git a/src/smith/differentiable_numerics/reaction.hpp b/src/smith/differentiable_numerics/reaction.hpp new file mode 100644 index 000000000..65d9ef69d --- /dev/null +++ b/src/smith/differentiable_numerics/reaction.hpp @@ -0,0 +1,89 @@ +// 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 reaction.hpp + * + * @brief Reaction class which is a names combination of a weak form and a set of dirichlet constrained nodes. + */ + +#pragma once + +#include +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" + +namespace smith { + +/// @brief gretl-function implementation which evaluates the residual force (which is minus the mechanical force) +/// given +/// shape displacement, states and params. The field_for_residual_space Field is only used to set the approriate size +/// (mfem::ParFiniteElementSpace) for the residual field so it can be returned as a ReactionState +inline auto evaluateWeakForm(const std::shared_ptr& weak_form, const TimeInfo& time_info, + FieldState shape_disp, const std::vector& field_states, + FieldState field_for_residual_space) +{ + std::vector all_state_bases{shape_disp}; + for (auto& f : field_states) all_state_bases.push_back(f); + all_state_bases.push_back(field_for_residual_space); + + auto z = shape_disp.create_state(all_state_bases, zero_state_from_dual()); + + z.set_eval([=](const gretl::UpstreamStates& inputs, gretl::DownstreamState& output) { + SMITH_MARK_FUNCTION; + + size_t num_fields = inputs.size() - 2; + ConstFieldPtr shape_disp_ = inputs[0].get().get(); + std::vector fields(num_fields); + for (size_t field_index = 0; field_index < num_fields; ++field_index) { + fields[field_index] = inputs[field_index + 1].get().get(); + } + ConstFieldPtr field_for_residual_space_ = inputs[num_fields + 1].get().get(); + + FEDualPtr R = std::make_shared(field_for_residual_space_->space(), + "residual"); // set up output pointer + // evaluate the residual with zero acceleration contribution + // std::cout << "time info = " << time_info.time() << std::endl; + // std::cout << "num fields = " << fields.size() << std::endl; + // std::cout << "shape disp name = " << inputs[0].get().get()->name() << std::endl; + // for (auto& f : fields) { + // std::cout << "name = " << f->name() << std::endl; + // } + *R = weak_form->residual(time_info, shape_disp_, fields); + output.set(R); + }); + + z.set_vjp([=](gretl::UpstreamStates& inputs, const gretl::DownstreamState& output) { + SMITH_MARK_FUNCTION; + + const FEDualPtr Z = output.get(); + const FEFieldPtr Z_dual = output.get_dual(); + FiniteElementState Z_dual_state(Z_dual->space(), Z_dual->name()); + Z_dual_state = *Z_dual; + + size_t num_fields = inputs.size() - 2; + std::vector fields(num_fields); + for (size_t field_index = 0; field_index < num_fields; ++field_index) { + fields[field_index] = inputs[field_index + 1].get().get(); + } + + std::vector field_sensitivities(num_fields); + for (size_t field_index = 0; field_index < num_fields; ++field_index) { + field_sensitivities[field_index] = inputs[field_index + 1].get_dual().get(); + } + + ConstFieldPtr shape_disp_ = inputs[0].get().get(); + DualFieldPtr shape_disp_sensitivity = inputs[0].get_dual().get(); + + // set the dual fields for each input, using the call to residual that pulls the derivative + weak_form->vjp(time_info, shape_disp_, fields, {}, &Z_dual_state, shape_disp_sensitivity, field_sensitivities, {}); + }); + + return z.finalize(); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp new file mode 100644 index 000000000..cdd1afa5e --- /dev/null +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -0,0 +1,65 @@ +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/solid_mechanics_state_advancer.hpp" +#include "smith/differentiable_numerics/reaction.hpp" + +namespace smith { + +SolidMechanicsStateAdvancer::SolidMechanicsStateAdvancer( + std::shared_ptr solver, std::shared_ptr vector_bcs, + std::shared_ptr solid_dynamic_weak_forms, + smith::SecondOrderTimeIntegrationRule time_rule) + : solver_(solver), + vector_bcs_(vector_bcs), + solid_dynamic_weak_forms_(solid_dynamic_weak_forms), + time_rule_(time_rule) +{ +} + +std::vector SolidMechanicsStateAdvancer::advanceState(const TimeInfo& time_info, + const FieldState& shape_disp, + const std::vector& states_old_, + const std::vector& params) const +{ + std::vector states_old = states_old_; + if (time_info.cycle() == 0) { + states_old[ACCELERATION] = solve(*solid_dynamic_weak_forms_->quasi_static_weak_form, shape_disp, states_old_, + params, time_info, *solver_, *vector_bcs_, ACCELERATION); + } + + double dt = time_info.dt(); + size_t cycle = time_info.cycle(); + double final_time = time_info.time() + dt; + TimeInfo final_time_info(final_time, dt, cycle); + + std::vector solid_inputs{states_old[DISPLACEMENT], states_old[DISPLACEMENT], states_old[VELOCITY], + states_old[ACCELERATION]}; + + auto displacement = solve(*solid_dynamic_weak_forms_->time_discretized_weak_form, shape_disp, solid_inputs, params, + final_time_info, *solver_, *vector_bcs_); + + std::vector states = states_old; + + states[DISPLACEMENT] = displacement; + states[VELOCITY] = time_rule_.derivative(final_time_info, displacement, states_old[DISPLACEMENT], + states_old[VELOCITY], states_old[ACCELERATION]); + states[ACCELERATION] = time_rule_.second_derivative(final_time_info, displacement, states_old[DISPLACEMENT], + states_old[VELOCITY], states_old[ACCELERATION]); + + return states; +} + +std::vector SolidMechanicsStateAdvancer::computeReactions(const TimeInfo& time_info, + const FieldState& shape_disp, + const std::vector& states, + const std::vector& params) const +{ + std::vector solid_inputs{states[DISPLACEMENT], states[VELOCITY], states[ACCELERATION]}; + solid_inputs.insert(solid_inputs.end(), params.begin(), params.end()); + return {evaluateWeakForm(solid_dynamic_weak_forms_->quasi_static_weak_form, time_info, shape_disp, solid_inputs, + states[DISPLACEMENT])}; +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp new file mode 100644 index 000000000..22dc42936 --- /dev/null +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp @@ -0,0 +1,118 @@ +// Copyright (c) 2019-2024, 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 solid_mechanics_state_advancer.hpp + * .hpp + * + * @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of + * equations + */ + +#pragma once + +#include "gretl/data_store.hpp" +#include "smith/smith_config.hpp" +#include "smith/physics/mesh.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" + +namespace smith { + +class DifferentiableSolver; +class DirichletBoundaryConditions; +class SecondOrderTimeDiscretizedWeakForms; + +/// @brief Implementation of the StateAdvancer interface for advancing the solution of solid mechanics problems +class SolidMechanicsStateAdvancer : public StateAdvancer { + public: + /// @brief Constructor + /// @param solid_solver differentiable solve + /// @param vector_bcs Dirichlet boundary conditions that can be applies to vector unknowns + /// @param solid_dynamic_weak_forms The weak-forms for time discretized solid mechanics equations + /// @param time_rule The specific time-integration rule, typically Implicit Newmark or Quasi-static + SolidMechanicsStateAdvancer(std::shared_ptr solid_solver, + std::shared_ptr vector_bcs, + std::shared_ptr solid_dynamic_weak_forms, + SecondOrderTimeIntegrationRule time_rule); + + /// State enum for indexing convenience + enum STATE + { + DISPLACEMENT, + VELOCITY, + ACCELERATION + }; + + /// @brief Recursive function for constructing parameter FieldStates of the appropriate space and name, register it on + /// the gretl graph. + template + static std::vector createParams(gretl::DataStore& graph, const std::string& name, + const std::vector& param_names, const std::string& tag, + size_t index = 0) + { + FieldState newParam = create_field_state(graph, FirstParamSpace{}, name + "_" + param_names[index], tag); + std::vector end_spaces{}; + if constexpr (sizeof...(ParamSpaces) > 0) { + end_spaces = createParams(graph, name, param_names, tag, ++index); + } + end_spaces.insert(end_spaces.begin(), newParam); + return end_spaces; + } + + /// @brief Utility function to consistently construct all the weak forms and FieldStates for a solid mechanics + /// application You will get back: shape_disp, states, params, time, and solid_mechanics_weak_form + template + static auto buildWeakFormAndStates(const std::shared_ptr& mesh, const std::shared_ptr& graph, + SecondOrderTimeIntegrationRule time_rule, std::string physics_name, + const std::vector& param_names, double initial_time = 0.0) + { + auto shape_disp = create_field_state(*graph, ShapeDispSpace{}, physics_name + "_shape_displacement", mesh->tag()); + auto disp = create_field_state(*graph, VectorSpace{}, physics_name + "_displacement", mesh->tag()); + auto velo = create_field_state(*graph, VectorSpace{}, physics_name + "_velocity", mesh->tag()); + auto acceleration = create_field_state(*graph, VectorSpace{}, physics_name + "_acceleration", mesh->tag()); + auto time = graph->create_state(initial_time); + std::vector params = + createParams(*graph, physics_name + "_param", param_names, mesh->tag()); + std::vector states{disp, velo, acceleration}; + + // weak form unknowns are disp, disp_old, velo_old, accel_old + using SolidWeakFormT = SecondOrderTimeDiscretizedWeakForm< + spatial_dim, VectorSpace, Parameters>; + auto input_spaces = spaces({states[DISPLACEMENT], states[DISPLACEMENT], states[VELOCITY], states[ACCELERATION]}); + auto param_spaces = spaces(params); + input_spaces.insert(input_spaces.end(), param_spaces.begin(), param_spaces.end()); + + auto solid_mechanics_weak_form = + std::make_shared(physics_name, mesh, time_rule, space(states[DISPLACEMENT]), input_spaces); + + return std::make_tuple(shape_disp, states, params, time, solid_mechanics_weak_form); + } + + /// @overload + std::vector advanceState(const TimeInfo& time_info, const FieldState& shape_disp, + const std::vector& states_old, + const std::vector& params) const override; + + /// @overload + std::vector computeReactions(const TimeInfo& time_info, const FieldState& shape_disp, + const std::vector& states, + const std::vector& params) const override; + + private: + std::shared_ptr solver_; ///< Differentiable solver + std::shared_ptr vector_bcs_; ///< Dirichlet boundary conditions on a vector-field + std::shared_ptr + solid_dynamic_weak_forms_; ///< Solid mechanics time discretized weak forms, user must setup the appropriate + ///< integrals. Has both the time discretized and the undiscretized weak forms. + SecondOrderTimeIntegrationRule time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, + ///< given the current predicted u and the previous u, u_dot, u_dot_dot +}; + +} // namespace smith diff --git a/src/smith/differentiable_numerics/state_advancer.hpp b/src/smith/differentiable_numerics/state_advancer.hpp index 6e4d53c9b..4c17147e0 100644 --- a/src/smith/differentiable_numerics/state_advancer.hpp +++ b/src/smith/differentiable_numerics/state_advancer.hpp @@ -28,9 +28,18 @@ class StateAdvancer { /// @brief interface method to advance the states from a given cycle and time, to the next cycle (cycle+1) and time /// (time+dt). shape_disp and params are assumed to be fixed in this advance. Time and time increment (dt) are /// gretl::State in order to record the duals on the reverse pass - virtual std::vector advanceState(const FieldState& shape_disp, const std::vector& states, - const std::vector& params, - const TimeInfo& time_info) const = 0; + virtual std::vector advanceState(const TimeInfo& time_info, const FieldState& shape_disp, + const std::vector& states, + const std::vector& params) const = 0; + + /// @brief interface method to compute reactions given previous, current states and + /// parameters. + virtual std::vector computeReactions(const TimeInfo& /*time_info*/, const FieldState& /*shape_disp*/, + const std::vector& /*states*/, + const std::vector& /*params*/) const + { + return std::vector{}; + } }; /// @brief creates a time info struct from gretl::State diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index 2ef0f624d..76e85d8b3 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -4,14 +4,23 @@ # # SPDX-License-Identifier: (BSD-3-Clause) -set(differentiable_numerics_test_depends smith_physics smith_differentiable_numerics gtest) +set(differentiable_numerics_test_headers + differentiable_test_utils.hpp + ) -set(differentiable_numerics_tests - test_explicit_dynamics.cpp +smith_add_library( + NAME differentiable_numerics_test_utils + HEADERS ${differentiable_numerics_test_headers} + ) + +set(differentiable_numerics_test_depends smith_physics smith_differentiable_numerics differentiable_numerics_test_utils gtest) + +set(differentiable_numerics_test_source test_field_state.cpp - paraview_helper.hpp + test_explicit_dynamics.cpp + test_solid_mechanics_state_advancer.cpp ) -smith_add_tests( SOURCES ${differentiable_numerics_tests} +smith_add_tests( SOURCES ${differentiable_numerics_test_source} DEPENDS_ON ${differentiable_numerics_test_depends} NUM_MPI_TASKS 1) diff --git a/src/smith/differentiable_numerics/differentiable_utils.hpp b/src/smith/differentiable_numerics/tests/differentiable_test_utils.hpp similarity index 89% rename from src/smith/differentiable_numerics/differentiable_utils.hpp rename to src/smith/differentiable_numerics/tests/differentiable_test_utils.hpp index 90a3be53f..761172d3c 100644 --- a/src/smith/differentiable_numerics/differentiable_utils.hpp +++ b/src/smith/differentiable_numerics/tests/differentiable_test_utils.hpp @@ -5,7 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) /** - * @file differentiable_utils.hpp + * @file differentiable_test_utils.hpp * * @brief Utility functions for testing. */ @@ -18,15 +18,6 @@ namespace smith { -/// @brief Evaluates a DoubleState using a provided ScalarObjective instance, and the input arguments to that objective. -/// This operation is tracked on the gretl graph. -DoubleState evaluateObjective(std::shared_ptr objective, const TimeInfo& time_info, - const FieldState& shape_disp, const std::vector& inputs); - -/// @brief operation is tracked on the gretl graph. -DoubleState evaluateObjective(std::shared_ptr objective, const FieldState& shape_disp, - const std::vector& inputs); - /// @brief Utility function to construct a smith::functional which evaluates the total kinetic energy template auto createKineticEnergyIntegrator(smith::Domain& domain, const mfem::ParFiniteElementSpace& velocity_space, @@ -125,9 +116,11 @@ inline auto checkGradients(const gretl::State& objectiveState, gretl::St /// @brief Testing utility function which runs a gretl graph num_fd_steps (with increasingly smaller finite difference /// steps) to check if the computed graph gradients are converging to the finite differenced gradients at the expected /// rate -inline double checkGradWrt(const gretl::State& objective, smith::FieldState& input, gretl::DataStore& graph, - double eps, size_t num_fd_steps = 4, bool printmore = false) +inline double checkGradWrt(const gretl::State& objective, smith::FieldState& input, double eps, + size_t num_fd_steps = 4, bool printmore = false) { + auto& graph = objective.data_store(); + // reset each time, just to be sure graph.reset(); @@ -165,9 +158,11 @@ inline double checkGradWrt(const gretl::State& objective, smith::FieldSt /// @brief Testing utility function which runs a gretl graph num_fd_steps (with increasingly smaller finite difference /// steps) to check if the computed graph gradients are converging to the finite differenced gradients at the expected /// rate -inline double checkGradWrt(const gretl::State& objective, gretl::State& input, - gretl::DataStore& graph, double eps, size_t num_fd_steps = 4, bool printmore = false) +inline double checkGradWrt(const gretl::State& objective, gretl::State& input, double eps, + size_t num_fd_steps = 4, bool printmore = false) { + auto& graph = objective.data_store(); + // reset each time, just to be sure graph.reset(); diff --git a/src/smith/differentiable_numerics/tests/test_explicit_dynamics.cpp b/src/smith/differentiable_numerics/tests/test_explicit_dynamics.cpp index e8de7df4d..06462ef4f 100644 --- a/src/smith/differentiable_numerics/tests/test_explicit_dynamics.cpp +++ b/src/smith/differentiable_numerics/tests/test_explicit_dynamics.cpp @@ -16,10 +16,11 @@ #include "smith/differentiable_numerics/lumped_mass_explicit_newmark_state_advancer.hpp" #include "smith/differentiable_numerics/lumped_mass_weak_form.hpp" -#include "smith/differentiable_numerics/tests/paraview_helper.hpp" -#include "smith/differentiable_numerics/differentiable_utils.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" #include "smith/differentiable_numerics/timestep_estimator.hpp" #include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/differentiable_numerics/evaluate_objective.hpp" +#include "smith/differentiable_numerics/tests/differentiable_test_utils.hpp" // This tests the interface between the new smith::WeakForm with gretl and its conformity to the existing base_physics // interface @@ -112,7 +113,7 @@ struct MeshFixture : public testing::Test { auto mfem_shape = mfem::Element::QUADRILATERAL; // mfem::Element::TRIANGLE; double length = 0.5; double width = 2.0; - mesh_ = std::make_shared(mfem::Mesh::MakeCartesian2D(5, 5, mfem_shape, true, length, width), MESHTAG, + mesh_ = std::make_shared(mfem::Mesh::MakeCartesian2D(5, 4, mfem_shape, true, length, width), MESHTAG, 0, 0); // checkpointing graph checkpointer_ = std::make_shared(200); @@ -161,7 +162,8 @@ struct MeshFixture : public testing::Test { // specify dirichlet bcs bc_manager_ = std::make_shared(mesh_->mfemParMesh()); - auto dt_estimator = std::make_shared(dt / double(dt_reduction)); + auto dt_estimator = + std::make_shared(dt / 10.0); // reduce the timestep a bit, so it subcycles std::shared_ptr time_integrator = std::make_shared(solid_mechanics_residual, solid_mass_residual, dt_estimator, bc_manager_); @@ -204,14 +206,14 @@ struct MeshFixture : public testing::Test { double integrateForward() { resetAndApplyInitialConditions(); - double lido_qoi = 0.0; + double base_physics_qoi = 0.0; for (size_t m = 0; m < num_steps; ++m) { physics_->advanceTimestep(dt); - lido_qoi += (*kinetic_energy_integrator_)(physics_->time(), physics_->shapeDisplacement(), - physics_->state(velo_name), physics_->parameter(DENSITY)); + base_physics_qoi += (*kinetic_energy_integrator_)(physics_->time(), physics_->shapeDisplacement(), + physics_->state(velo_name), physics_->parameter(DENSITY)); } - return lido_qoi; + return base_physics_qoi; } void adjointBackward(smith::FiniteElementDual& shape_sensitivity, @@ -263,12 +265,12 @@ struct MeshFixture : public testing::Test { std::shared_ptr bc_manager_; - const double dt = 1e-2; - const size_t dt_reduction = 10; - const size_t num_steps = 4; + static constexpr double total_simulation_time = 0.005; + static constexpr size_t num_steps = 10; + static constexpr double dt = total_simulation_time / num_steps; }; -TEST_F(MeshFixture, TRANSIENT_DYNAMICS_LIDO) +TEST_F(MeshFixture, TransientDynamicsBasePhysics) { SMITH_MARK_FUNCTION; @@ -300,7 +302,7 @@ TEST_F(MeshFixture, TRANSIENT_DYNAMICS_LIDO) } } -TEST_F(MeshFixture, TRANSIENT_DYNAMICS_GRETL) +TEST_F(MeshFixture, TransientDynamicsGretl) { SMITH_MARK_FUNCTION; @@ -309,19 +311,17 @@ TEST_F(MeshFixture, TRANSIENT_DYNAMICS_GRETL) resetAndApplyInitialConditions(); - auto all_fields = mechanics_->getAllFieldStates(); + auto all_fields = mechanics_->getFieldStatesAndParamStates(); gretl::State gretl_qoi = 0.0 * smith::evaluateObjective(objective_, *shape_disp_, {all_fields[F_VELO], all_fields[F_DENSITY]}); std::string pv_dir = std::string("paraview_") + mechanics_->name(); - auto pv_writer = smith::createParaviewOutput(*mesh_, all_fields, pv_dir); + auto pv_writer = smith::createParaviewWriter(*mesh_, all_fields, pv_dir); pv_writer.write(mechanics_->cycle(), mechanics_->time(), all_fields); for (size_t m = 0; m < num_steps; ++m) { - for (size_t n = 0; n < dt_reduction; ++n) { - mechanics_->advanceTimestep(dt / double(dt_reduction)); - } - all_fields = mechanics_->getAllFieldStates(); + mechanics_->advanceTimestep(dt); + all_fields = mechanics_->getFieldStatesAndParamStates(); gretl_qoi = gretl_qoi + smith::evaluateObjective(objective_, *shape_disp_, {all_fields[F_VELO], all_fields[F_DENSITY]}); pv_writer.write(mechanics_->cycle(), mechanics_->time(), all_fields); @@ -344,29 +344,28 @@ TEST_F(MeshFixture, TRANSIENT_DYNAMICS_GRETL) << std::endl; } - EXPECT_GT(smith::checkGradWrt(gretl_qoi, *shape_disp_, *checkpointer_, 0.01, 4, true), 0.8); - EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[DISP], *checkpointer_, 0.01, 4, true), 0.8); - EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[VELO], *checkpointer_, 0.01, 4, true), 0.8); - // EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states[ACCEL], *checkpointer_, 1.0, 4, true), 0.8); - EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[DENSITY], *checkpointer_, 0.01, 4, true), 0.8); + EXPECT_GT(smith::checkGradWrt(gretl_qoi, *shape_disp_, 0.01, 4, true), 0.8); + EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[DISP], 0.01, 4, true), 0.8); + EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[VELO], 0.01, 4, true), 0.8); + EXPECT_GT(smith::checkGradWrt(gretl_qoi, initial_states_[DENSITY], 0.01, 4, true), 0.8); } -TEST_F(MeshFixture, TRANSIENT_CONSTANT_GRAVITY) +TEST_F(MeshFixture, TransientConstantGravity) { SMITH_MARK_FUNCTION; mechanics_->resetStates(); - auto all_fields = mechanics_->getAllFieldStates(); + auto all_fields = mechanics_->getFieldStatesAndParamStates(); std::string pv_dir = std::string("paraview_") + mechanics_->name(); std::cout << "Writing output to " << pv_dir << std::endl; - auto pv_writer = smith::createParaviewOutput(*mesh_, all_fields, pv_dir); + auto pv_writer = smith::createParaviewWriter(*mesh_, all_fields, pv_dir); pv_writer.write(mechanics_->cycle(), mechanics_->time(), all_fields); double time = 0.0; - for (size_t m = 0; m < dt_reduction * num_steps; ++m) { - double timestep = dt / double(dt_reduction); + for (size_t m = 0; m < num_steps; ++m) { + double timestep = dt / double(num_steps); mechanics_->advanceTimestep(timestep); - all_fields = mechanics_->getAllFieldStates(); + all_fields = mechanics_->getFieldStatesAndParamStates(); pv_writer.write(mechanics_->cycle(), mechanics_->time(), all_fields); time += timestep; } diff --git a/src/smith/differentiable_numerics/tests/test_field_state.cpp b/src/smith/differentiable_numerics/tests/test_field_state.cpp index 49bfa4fd6..115a34a96 100644 --- a/src/smith/differentiable_numerics/tests/test_field_state.cpp +++ b/src/smith/differentiable_numerics/tests/test_field_state.cpp @@ -11,7 +11,7 @@ #include "gretl/data_store.hpp" #include "smith/differentiable_numerics/field_state.hpp" -#include "smith/differentiable_numerics/differentiable_utils.hpp" +#include "smith/differentiable_numerics/tests/differentiable_test_utils.hpp" // This tests the interface between the new smith::WeakForm with gretl and its conformity to the existing base_physics // interface @@ -86,17 +86,15 @@ TEST_F(MeshFixture, FieldStateWithDifferentiable_axpby) auto u = axpby(dt, disp, dt, velo); auto u_exact = axpby(dt_f, disp, dt_f, velo); - auto uu_exact = smith::inner_product(u_exact, u_exact); - auto uu = smith::inner_product(u, u); + auto uu_exact = smith::innerProduct(u_exact, u_exact); + auto uu = smith::innerProduct(u, u); gretl::set_as_objective(uu); EXPECT_EQ(uu.get(), uu_exact.get()); - checkpointer_->back_prop(); - - EXPECT_GT(smith::checkGradWrt(uu, disp, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, velo, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, dt, *checkpointer_, 1e-7, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, disp, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, velo, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, dt, 1e-7, 4, true), 0.95); } TEST_F(MeshFixture, FieldStateDifferentiablyWeightedSum_WithOperators) @@ -133,8 +131,8 @@ TEST_F(MeshFixture, FieldStateDifferentiablyWeightedSum_WithOperators) u_exact = axpby(2.0, velo, -2.0, u_exact); u_exact = axpby(initial_dt * initial_h, u_exact, 0.0, u_exact); - auto uu_exact = smith::inner_product(u_exact, u_exact); - auto uu = smith::inner_product(u, u); + auto uu_exact = smith::innerProduct(u_exact, u_exact); + auto uu = smith::innerProduct(u, u); gretl::set_as_objective(uu); @@ -142,11 +140,11 @@ TEST_F(MeshFixture, FieldStateDifferentiablyWeightedSum_WithOperators) checkpointer_->back_prop(); - EXPECT_GT(smith::checkGradWrt(uu, disp, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, velo, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, accel, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, dt, *checkpointer_, 1e-5, 4, true), 0.95); - EXPECT_GT(smith::checkGradWrt(uu, h, *checkpointer_, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, disp, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, velo, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, accel, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, dt, 1e-5, 4, true), 0.95); + EXPECT_GT(smith::checkGradWrt(uu, h, 1e-5, 4, true), 0.95); } int main(int argc, char* argv[]) diff --git a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp new file mode 100644 index 000000000..5ec9bb3a4 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp @@ -0,0 +1,399 @@ +// 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) + +#include + +#include "gretl/data_store.hpp" + +#include "smith/smith_config.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/equation_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/mesh_utils/mesh_utils.hpp" + +#include "smith/physics/state/state_manager.hpp" +#include "smith/physics/functional_objective.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" +#include "smith/physics/materials/parameterized_solid_material.hpp" + +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/differentiable_solid_mechanics.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/differentiable_numerics/tests/differentiable_test_utils.hpp" + +namespace smith { + +smith::LinearSolverOptions solid_linear_options{.linear_solver = smith::LinearSolver::CG, + .preconditioner = smith::Preconditioner::HypreJacobi, + .relative_tol = 1e-11, + .absolute_tol = 1e-11, + .max_iterations = 10000, + .print_level = 0}; + +smith::NonlinearSolverOptions solid_nonlinear_opts{.nonlin_solver = NonlinearSolver::TrustRegion, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-10, + .max_iterations = 500, + .print_level = 0}; + +static constexpr int dim = 3; +static constexpr int order = 1; + +using ShapeDispSpace = H1<1, dim>; +using VectorSpace = H1; +using ScalarParameterSpace = L2<0>; + +struct SolidMechanicsMeshFixture : public testing::Test { + double length = 1.0; + double width = 0.04; + int num_elements_x = 12; + int num_elements_y = 2; + int num_elements_z = 2; + double elem_size = length / num_elements_x; + + void SetUp() + { + smith::StateManager::initialize(datastore, "solid"); + auto mfem_shape = mfem::Element::QUADRILATERAL; + mesh = std::make_shared( + mfem::Mesh::MakeCartesian3D(num_elements_x, num_elements_y, num_elements_z, mfem_shape, length, width, width), + "mesh", 0, 0); + mesh->addDomainOfBoundaryElements("left", smith::by_attr(3)); + mesh->addDomainOfBoundaryElements("right", smith::by_attr(5)); + } + + static constexpr double total_simulation_time_ = 1.1; + static constexpr size_t num_steps_ = 4; + static constexpr double dt_ = total_simulation_time_ / num_steps_; + + axom::sidre::DataStore datastore; + std::shared_ptr mesh; +}; + +TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) +{ + SMITH_MARK_FUNCTION; + + enum STATE + { + DISP, + VELO, + ACCEL + }; + + std::string physics_name = "solid"; + + std::shared_ptr d_solid_nonlinear_solver = + buildDifferentiableNonlinearSolve(solid_nonlinear_opts, solid_linear_options, *mesh); + + smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + + auto [physics, solid_weak_form, bcs] = + buildSolidMechanics( + mesh, d_solid_nonlinear_solver, time_rule, 100, physics_name, {"bulk", "shear"}); + + static constexpr double gravity = -9.0; + + double E = 100.0; + double nu = 0.25; + auto K = E / (3.0 * (1.0 - 2.0 * nu)); + auto G = E / (2.0 * (1.0 + nu)); + using MaterialType = solid_mechanics::ParameterizedNeoHookeanSolid; + MaterialType material{.density = 1.0, .K0 = K, .G0 = G}; + + solid_weak_form->addBodyIntegral( + smith::DependsOn<0, 1>{}, mesh->entireBodyName(), + [material](const auto& /*time_info*/, auto /*X*/, auto u, auto /*v*/, auto a, auto bulk, auto shear) { + MaterialType::State state; + auto pk_stress = material(state, get(u), bulk, shear); + smith::tensor b{}; + b[1] = gravity; + return smith::tuple{get(a) * material.density - b, pk_stress}; + }); + + auto shape_disp = physics->getShapeDispFieldState(); + auto params = physics->getFieldParams(); + auto states = physics->getInitialFieldStates(); + + params[0].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.K0; + }); + + params[1].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.G0; + }); + + physics->resetStates(); + auto all_fields = physics->getFieldStatesAndParamStates(); + + std::string pv_dir = std::string("paraview_") + physics->name(); + auto pv_writer = createParaviewWriter(*mesh, all_fields, pv_dir); + pv_writer.write(physics->cycle(), physics->time(), all_fields); + + for (size_t m = 0; m < num_steps_; ++m) { + physics->advanceTimestep(dt_); + all_fields = physics->getFieldStatesAndParamStates(); + pv_writer.write(physics->cycle(), physics->time(), all_fields); + } + + double a_exact = gravity; + double v_exact = gravity * total_simulation_time_; + double u_exact = 0.5 * gravity * total_simulation_time_ * total_simulation_time_; + + TimeInfo end_time_info(physics->time(), dt_, static_cast(physics->cycle())); + + FunctionalObjective> accel_error("accel_error", mesh, spaces({all_fields[ACCEL]})); + accel_error.addBodyIntegral(DependsOn<0>{}, mesh->entireBodyName(), [a_exact](auto /*t*/, auto /*X*/, auto A) { + auto a = get(A); + auto da0 = a[0]; + auto da1 = a[1] - a_exact; + return da0 * da0 + da1 * da1; + }); + double a_err = + accel_error.evaluate(end_time_info, shape_disp.get().get(), getConstFieldPointers({all_fields[ACCEL]})); + EXPECT_NEAR(0.0, a_err, 1e-14); + + FunctionalObjective> velo_error("velo_error", mesh, spaces({all_fields[VELO]})); + velo_error.addBodyIntegral(DependsOn<0>{}, mesh->entireBodyName(), [v_exact](auto /*t*/, auto /*X*/, auto V) { + auto v = get(V); + auto dv0 = v[0]; + auto dv1 = v[1] - v_exact; + return dv0 * dv0 + dv1 * dv1; + }); + double v_err = + velo_error.evaluate(TimeInfo(0.0, 1.0, 0), shape_disp.get().get(), getConstFieldPointers({all_fields[VELO]})); + EXPECT_NEAR(0.0, v_err, 1e-14); + + FunctionalObjective> disp_error("disp_error", mesh, spaces({all_fields[DISP]})); + disp_error.addBodyIntegral(DependsOn<0>{}, mesh->entireBodyName(), [u_exact](auto /*t*/, auto /*X*/, auto U) { + auto u = get(U); + auto du0 = u[0]; + auto du1 = u[1] - u_exact; + return du0 * du0 + du1 * du1; + }); + double u_err = + disp_error.evaluate(TimeInfo(0.0, 1.0, 0), shape_disp.get().get(), getConstFieldPointers({all_fields[DISP]})); + EXPECT_NEAR(0.0, u_err, 1e-14); +} + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) +{ + SMITH_MARK_FUNCTION; + + std::string physics_name = "solid"; + + std::shared_ptr d_solid_nonlinear_solver = + buildDifferentiableNonlinearSolve(solid_nonlinear_opts, solid_linear_options, *mesh); + + smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + + auto [physics, solid_weak_form, bcs] = + buildSolidMechanics( + mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); + + bcs->setFixedVectorBCs(mesh->domain("right")); + bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { + auto bc = 0.0 * X; + bc[0] = 0.01 * t; + bc[1] = -0.05 * t; + return bc; + }); + + double E = 100.0; + double nu = 0.25; + auto K = E / (3.0 * (1.0 - 2 * nu)); + auto G = E / (2.0 * (1.0 + nu)); + using MaterialType = solid_mechanics::ParameterizedNeoHookeanSolid; + MaterialType material{.density = 10.0, .K0 = K, .G0 = G}; + + solid_weak_form->addBodyIntegral( + smith::DependsOn<0, 1>{}, mesh->entireBodyName(), + [material](const auto& /*time_info*/, auto /*X*/, auto u, auto /*v*/, auto a, auto bulk, auto shear) { + MaterialType::State state; + auto pk_stress = material(state, get(u), bulk, shear); + return smith::tuple{get(a) * material.density, pk_stress}; + }); + + auto shape_disp = physics->getShapeDispFieldState(); + auto params = physics->getFieldParams(); + auto initial_states = physics->getInitialFieldStates(); + + params[0].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.K0; + }); + + params[1].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.G0; + }); + + physics->resetStates(); + + auto pv_writer = smith::createParaviewWriter(*mesh, physics->getFieldStatesAndParamStates(), physics_name); + pv_writer.write(0, physics->time(), physics->getFieldStatesAndParamStates()); + for (size_t m = 0; m < num_steps_; ++m) { + physics->advanceTimestep(dt_); + pv_writer.write(m + 1, physics->time(), physics->getFieldStatesAndParamStates()); + } + + TimeInfo time_info(physics->time(), dt_); + + auto state_advancer = physics->getStateAdvancer(); + auto reactions = state_advancer->computeReactions(time_info, shape_disp, physics->getFieldStates(), params); + + auto reaction_squared = 0.5 * innerProduct(reactions[0], reactions[0]); + + gretl::set_as_objective(reaction_squared); + std::cout << "final residual norm2 = " << reaction_squared.get() << std::endl; + + EXPECT_GT(checkGradWrt(reaction_squared, shape_disp, 1.1e-2, 4, true), 0.7); + EXPECT_GT(checkGradWrt(reaction_squared, params[0], 6.2e-1, 4, true), 0.7); + EXPECT_GT(checkGradWrt(reaction_squared, params[1], 6.2e-1, 4, true), 0.7); + + // re-evaluate the final objective value, and backpropagate again + reaction_squared.get(); + gretl::set_as_objective(reaction_squared); + reaction_squared.data_store().back_prop(); + + for (auto s : initial_states) { + std::cout << s.get()->name() << " " << s.get()->Norml2() << " " << s.get_dual()->Norml2() << std::endl; + } + + std::cout << shape_disp.get()->name() << " " << shape_disp.get()->Norml2() << " " << shape_disp.get_dual()->Norml2() + << std::endl; + + for (size_t p = 0; p < params.size(); ++p) { + std::cout << params[p].get()->name() << " " << params[p].get()->Norml2() << " " << params[p].get_dual()->Norml2() + << std::endl; + } +} + +// these functions mimic the BasePhysics style of running smith + +void resetAndApplyInitialConditions(std::shared_ptr physics) { physics->resetStates(); } + +double integrateForward(std::shared_ptr physics, size_t num_steps, double dt) +{ + resetAndApplyInitialConditions(physics); + for (size_t m = 0; m < num_steps; ++m) { + physics->advanceTimestep(dt); + } + FiniteElementDual reaction = physics->dual("reactions"); + + return 0.5 * innerProduct(reaction, reaction); +} + +void adjointBackward(std::shared_ptr physics, smith::FiniteElementDual& shape_sensitivity, + std::vector& parameter_sensitivities) +{ + smith::FiniteElementDual reaction = physics->dual("reactions"); + smith::FiniteElementState reaction_dual(reaction.space(), "reactions_dual"); + reaction_dual = reaction; + + physics->resetAdjointStates(); + + physics->setDualAdjointBcs({{"reactions", reaction_dual}}); + + while (physics->cycle() > 0) { + physics->reverseAdjointTimestep(); + shape_sensitivity += physics->computeTimestepShapeSensitivity(); + for (size_t param_index = 0; param_index < parameter_sensitivities.size(); ++param_index) { + parameter_sensitivities[param_index] += physics->computeTimestepSensitivity(param_index); + } + } +} + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) +{ + SMITH_MARK_FUNCTION; + + std::string physics_name = "solid"; + + std::shared_ptr d_solid_nonlinear_solver = + buildDifferentiableNonlinearSolve(solid_nonlinear_opts, solid_linear_options, *mesh); + + smith::SecondOrderTimeIntegrationRule time_rule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + + auto [physics, solid_weak_form, bcs] = + buildSolidMechanics( + mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); + + bcs->setFixedVectorBCs(mesh->domain("right")); + bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { + auto bc = 0.0 * X; + bc[0] = 0.01 * t; + bc[1] = -0.05 * t; + return bc; + }); + + double E = 100.0; + double nu = 0.25; + auto K = E / (3.0 * (1.0 - 2 * nu)); + auto G = E / (2.0 * (1.0 + nu)); + using MaterialType = solid_mechanics::ParameterizedNeoHookeanSolid; + MaterialType material{.density = 10.0, .K0 = K, .G0 = G}; + + solid_weak_form->addBodyIntegral( + smith::DependsOn<0, 1>{}, mesh->entireBodyName(), + [material](const auto& /*time_info*/, auto /*X*/, auto u, auto /*v*/, auto a, auto bulk, auto shear) { + MaterialType::State state; + auto pk_stress = material(state, get(u), bulk, shear); + return smith::tuple{get(a) * material.density, pk_stress}; + }); + + auto shape_disp = physics->getShapeDispFieldState(); + auto params = physics->getFieldParams(); + auto states = physics->getInitialFieldStates(); + + params[0].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.K0; + }); + + params[1].get()->setFromFieldFunction([=](smith::tensor) { + double scaling = 1.0; + return scaling * material.G0; + }); + + physics->resetStates(); + + double qoi = integrateForward(physics, num_steps_, dt_); + std::cout << "qoi = " << qoi << std::endl; + + size_t num_params = physics->parameterNames().size(); + + smith::FiniteElementDual shape_sensitivity(*shape_disp.get_dual()); + std::vector parameter_sensitivities; + for (size_t p = 0; p < num_params; ++p) { + parameter_sensitivities.emplace_back(*params[p].get_dual()); + } + + adjointBackward(physics, shape_sensitivity, parameter_sensitivities); + + auto state_sensitivities = physics->computeInitialConditionSensitivity(); + for (auto name_and_state_sensitivity : state_sensitivities) { + std::cout << name_and_state_sensitivity.first << " " << name_and_state_sensitivity.second.Norml2() << std::endl; + } + + std::cout << shape_sensitivity.name() << " " << shape_sensitivity.Norml2() << std::endl; + + for (size_t p = 0; p < num_params; ++p) { + std::cout << parameter_sensitivities[p].name() << " " << parameter_sensitivities[p].Norml2() << std::endl; + } +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp new file mode 100644 index 000000000..fbafedab3 --- /dev/null +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -0,0 +1,154 @@ +// Copyright (c) 2019-2024, 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 time_discretized_weak_form.hpp + * + * @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of + * equations + */ + +#pragma once + +#include "smith/physics/functional_weak_form.hpp" +#include "smith/physics/mesh.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" + +namespace smith { + +template > +class TimeDiscretizedWeakForm; + +/// @brief A time discretized weakform gets a TimeInfo object passed as arguments to q-function (lambdas which are +/// integrated over quadrature points) so users can have access to time increments, and timestep cycle. These +/// quantities are often valuable for time integrated PDEs. +/// @tparam OutputSpace The output residual for the weak form (test-space) +/// @tparam ...InputSpaces All the input FiniteElementState fields (trial-spaces) +/// @tparam spatial_dim The spatial dimension for the problem +template +class TimeDiscretizedWeakForm> + : public FunctionalWeakForm> { + public: + using WeakFormT = FunctionalWeakForm>; ///< using + + /// Constructor + TimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr mesh, + const mfem::ParFiniteElementSpace& output_mfem_space, + const typename WeakFormT::SpacesT& input_mfem_spaces) + : WeakFormT(physics_name, mesh, output_mfem_space, input_mfem_spaces) + { + } + + /// @overload + template + void addBodyIntegral(DependsOn depends_on, std::string body_name, BodyIntegralType integrand) + { + const double* dt = &this->dt_; + const size_t* cycle = &this->cycle_; + WeakFormT::addBodyIntegral(depends_on, body_name, [dt, cycle, integrand](double t, auto X, auto... inputs) { + TimeInfo time_info(t, *dt, *cycle); + return integrand(time_info, X, inputs...); + }); + } + + /// @overload + template + void addBodyIntegral(std::string body_name, BodyForceType body_integral) + { + addBodyIntegral(DependsOn<>{}, body_name, body_integral); + } +}; + +/// @brief A container holding the two types of weak forms useful for solving time discretized second order (in time) +/// systems of equations +class SecondOrderTimeDiscretizedWeakForms { + public: + std::shared_ptr time_discretized_weak_form; ///< this publically available abstract weak form is a + ///< functions of the current u, u_old, v_old, and a_old, + std::shared_ptr quasi_static_weak_form; ///< this publically available abstract weak form is structly a + ///< function of the current u, v, and a (no time discretization) +}; + +template > +class SecondOrderTimeDiscretizedWeakForm; + +/// @brief Useful for time-discretized PDEs of second order (involves for first and second derivatives of time). Users +/// write q-functions in terns of u, u_dot, u_dot_dot, and the weak form is transformed by the +/// SecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, u_dot_dot_old, with u as +/// the distinct unknown for the time discretized system. +/// @tparam spatial_dim Spatial dimension, 2 or 3. +/// @tparam OutputSpace The space corresponding to the output residual for the weak form (test-space). +/// @tparam TrialInputSpace The space corresponding to the predicted solution u, i.e., the trial solution, the unique +/// unknown of the time discretized equation. +/// @tparam ...InputSpaces Spaces for all the remaining input FiniteElementState fields. +/// @tparam spatial_dim The spatial dimension for the problem. +template +class SecondOrderTimeDiscretizedWeakForm> + : public SecondOrderTimeDiscretizedWeakForms { + public: + static constexpr int NUM_STATE_VARS = 4; ///< u, u_old, v_old, a_old + + using TimeDiscretizedWeakFormT = + TimeDiscretizedWeakForm>; ///< using + using QuasiStaticWeakFormT = + TimeDiscretizedWeakForm>; ///< using + + /// @brief Constructor + SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr mesh, + SecondOrderTimeIntegrationRule time_rule, + const mfem::ParFiniteElementSpace& output_mfem_space, + const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces) + : time_rule_(time_rule) + { + time_discretized_weak_form_ = + std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces); + time_discretized_weak_form = time_discretized_weak_form_; + + typename TimeDiscretizedWeakFormT::SpacesT input_mfem_spaces_trial_removed(std::next(input_mfem_spaces.begin()), + input_mfem_spaces.end()); + quasi_static_weak_form_ = + std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); + quasi_static_weak_form = quasi_static_weak_form_; + } + + /// @overload + template + void addBodyIntegral(DependsOn /*depends_on*/, std::string body_name, + BodyIntegralType integrand) + { + auto time_rule = time_rule_; + time_discretized_weak_form_->addBodyIntegral( + DependsOn<0, 1, 2, 3, NUM_STATE_VARS + active_parameters...>{}, body_name, + [integrand, time_rule](const TimeInfo& t, auto X, auto U, auto U_old, auto U_dot_old, auto U_dot_dot_old, + auto... inputs) { + return integrand(t, X, time_rule.value(t, U, U_old, U_dot_old, U_dot_dot_old), + time_rule.derivative(t, U, U_old, U_dot_old, U_dot_dot_old), + time_rule.second_derivative(t, U, U_old, U_dot_old, U_dot_dot_old), inputs...); + }); + quasi_static_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, body_name, + integrand); + } + + /// @overload + template + void addBodyIntegral(std::string body_name, BodyForceType body_integral) + { + addBodyIntegral(DependsOn<>{}, body_name, body_integral); + } + + private: + std::shared_ptr + time_discretized_weak_form_; ///< fully templated time discretized weak form (with time integration rule injected + ///< into the q-function) + std::shared_ptr + quasi_static_weak_form_; ///< fully template underlying weak form (no time integration included, a function of + ///< current u, v, and a) + + SecondOrderTimeIntegrationRule time_rule_; ///< encodes the time integration rule +}; + +} // namespace smith diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp new file mode 100644 index 000000000..edef366e0 --- /dev/null +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -0,0 +1,100 @@ +// 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 time_integration_rule.hpp + * + * @brief Provides templated implementations for discretizing values, velocities and accelerations from current and + * previous states + */ + +#pragma once + +#include "smith/physics/common.hpp" + +namespace smith { + +/// @brief encodes rules for time discretizing first order odes (involving first time derivatives). +/// When solving f(u, u_dot, t) = 0 +/// this class provides the current discrete approximation for u and u_dot as a function of +/// (u^{n+1}, u^n). +struct BackwardEulerFirstOrderTimeIntegrationRule { + /// @brief Constructor + BackwardEulerFirstOrderTimeIntegrationRule() {} + + /// @brief evaluate value of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto value(const TimeInfo& /*t*/, const T1& field_new, const T2& /*field_old*/) const + { + return field_new; + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto derivative(const TimeInfo& t, const T1& field_new, const T2& field_old) const + { + return (1.0 / t.dt()) * (field_new - field_old); + } +}; + +/// @brief Options for second order time integration methods +enum class SecondOrderTimeIntegrationMethod +{ + IMPLICIT_NEWMARK, /// implicit newmark discretization + QUASI_STATIC /// quasi-static, specifies current field, velocity is central difference (for quasi-static artificial + /// viscosity), and acceleration is lagged for cases where it is set to some fixed value in time. +}; + +/// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). +/// When solving f(u, u_dot, u_dot_dot, t) = 0 +/// this class provides the current discrete approximation for u, u_dot, and u_dot_dot as a function of +/// (u^{n+1},u^n,u_dot^n,u_dot_dot^n). +struct SecondOrderTimeIntegrationRule { + /// @brief Constructor + SecondOrderTimeIntegrationRule(SecondOrderTimeIntegrationMethod method) : method_(method) {} + + /// @brief evaluate value of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto value([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return field_new; + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return (2.0 / t.dt()) * (field_new - field_old) - velo_old; + // if (method_ == SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK) { + // return (2.0 / t.dt()) * (field_new - field_old) - velo_old; + // } else { + // return (1.0 / t.dt()) * (field_new - field_old) - 0.0 * velo_old; + // } + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto second_derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + auto dt = t.dt(); + return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; + // auto accel = (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; + // if (method_ == SecondOrderTimeIntegrationMethod::QUASI_STATIC) { + // accel = accel_old + 0.0 * accel; + // } + // return accel; + } + + SecondOrderTimeIntegrationMethod method_; ///< method specifying time integration rule to inject into the q-function. +}; + +} // namespace smith diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index fde6b24e8..6c24646dc 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -110,7 +110,7 @@ class NewtonSolver : public mfem::NewtonSolver { int it = 0; for (; true; it++) { MFEM_ASSERT(mfem::IsFinite(norm), "norm = " << norm); - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Newton iteration " << std::setw(3) << it << " : ||r|| = " << std::setw(13) << norm; if (it > 0) { mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm); @@ -190,10 +190,10 @@ class NewtonSolver : public mfem::NewtonSolver { } if (ls_iter_sum) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Number of line search steps taken = " << ls_iter_sum << std::endl; } - if (print_level == 2 && (ls_iter_sum == 2 * max_ls_iters + 1)) { + if (print_level >= 2 && (ls_iter_sum == 2 * max_ls_iters + 1)) { mfem::out << "The maximum number of line search cut back have occurred, the resulting residual may not have " "decreased. " << std::endl; @@ -400,7 +400,7 @@ class TrustRegion : public mfem::NewtonSolver { try { std::tie(directions, H_directions) = removeDependentDirections(directions, H_directions); } catch (const std::exception& e) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "remove dependent directions failed with " << e.what() << std::endl; } return; @@ -432,7 +432,7 @@ class TrustRegion : public mfem::NewtonSolver { double base_energy = computeEnergy(g, hess_vec_func, z); double subspace_energy = computeEnergy(g, hess_vec_func, sol); - if (print_level == 2) { + if (print_level >= 2) { double leftval = leftvals.size() ? leftvals[0] : 1.0; mfem::out << "Energy using subspace solver from: " << base_energy << ", to: " << subspace_energy << " / " << energy_change << ". Min eig: " << leftval << std::endl; @@ -468,7 +468,7 @@ class TrustRegion : public mfem::NewtonSolver { if (cc >= tt) { add(s, std::sqrt(tt / cc), cp, s); } else if (cc > nn) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "cp outside newton, preconditioner likely inaccurate\n"; } add(s, 1.0, cp, s); @@ -513,7 +513,7 @@ class TrustRegion : public mfem::NewtonSolver { const double cg_tol_squared = settings.cg_tol * settings.cg_tol; if (Dot(r0, r0) <= cg_tol_squared && settings.min_cg_iterations == 0) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Trust region solution state within tolerance on first iteration." << "\n"; } @@ -566,7 +566,7 @@ class TrustRegion : public mfem::NewtonSolver { z = zPred; if (results.interior_status == TrustRegionResults::Status::NonDescentDirection) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Found a non descent direction\n"; } return; @@ -635,6 +635,9 @@ class TrustRegion : public mfem::NewtonSolver { MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator)."); MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver)."); + print_level = print_options.iterations ? 1 : print_level; + print_level = print_options.summary ? 2 : print_level; + using real_t = mfem::real_t; num_hess_vecs = 0; @@ -678,7 +681,7 @@ class TrustRegion : public mfem::NewtonSolver { int it = 0; for (; true; it++) { MFEM_ASSERT(mfem::IsFinite(norm), "norm = " << norm); - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "TrustRegion iteration " << std::setw(3) << it << " : ||r|| = " << std::setw(13) << norm; if (it > 0) { mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm); @@ -726,14 +729,14 @@ class TrustRegion : public mfem::NewtonSolver { } else { const double alphaTr = -tr_size / std::sqrt(Dot(r, r)); add(trResults.cauchy_point, alphaTr, r, trResults.cauchy_point); - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Negative curvature un-preconditioned cauchy point direction found." << "\n"; } } if (cauchyPointNormSquared >= tr_size * tr_size) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Un-preconditioned gradient cauchy point outside trust region, step size = " << std::sqrt(cauchyPointNormSquared) << "\n"; } @@ -806,7 +809,7 @@ class TrustRegion : public mfem::NewtonSolver { X = x_pred; r = r_pred; norm = normPred; - if (print_options.iterations) { + if (print_level >= 2) { printTrustRegionInfo(realObjective, modelObjective, trResults.cg_iterations_count, tr_size, true); trResults.cg_iterations_count = 0; // zero this output so it doesn't look like the linesearch is doing cg iterations @@ -819,7 +822,7 @@ class TrustRegion : public mfem::NewtonSolver { double rho = realImprove / modelImprove; if (modelObjective > 0) { - if (print_level == 2) { + if (print_level >= 2) { mfem::out << "Found a positive model objective increase. Debug if you see this.\n"; } rho = realImprove / -modelImprove; @@ -846,7 +849,7 @@ class TrustRegion : public mfem::NewtonSolver { // realResNorm = np.linalg.norm(gy) bool willAccept = rho >= settings.eta1 && rho <= settings.eta4; // or (rho >= -0 and realResNorm <= gNorm) - if (print_options.iterations) { + if (print_level >= 2) { printTrustRegionInfo(realObjective, modelObjective, trResults.cg_iterations_count, tr_size, willAccept); trResults.cg_iterations_count = 0; // zero this output so it doesn't look like the linesearch is doing cg iterations @@ -873,7 +876,7 @@ class TrustRegion : public mfem::NewtonSolver { mfem::out << "TrustRegion: No convergence!\n"; } - if (false && print_level == 2) { + if (false && print_level >= 2) { mfem::out << "num hess vecs = " << num_hess_vecs << "\n"; mfem::out << "num preconds = " << num_preconds << "\n"; mfem::out << "num residuals = " << num_residuals << "\n"; diff --git a/src/smith/physics/tests/solid_dynamics_patch.cpp b/src/smith/physics/tests/solid_dynamics_patch.cpp index 5a2a2dba5..5e10ee984 100644 --- a/src/smith/physics/tests/solid_dynamics_patch.cpp +++ b/src/smith/physics/tests/solid_dynamics_patch.cpp @@ -318,23 +318,6 @@ double solution_error(solution_type exact_solution, PatchBoundaryCondition bc) // Integrate in time for (int i = 0; i < 3; i++) { solid.advanceTimestep(1.0); - - // Output solution for debugging - // solid.outputStateToDisk("paraview_output"); - // std::cout << "cycle " << i << std::endl; - // std::cout << "time = " << solid.time() << std::endl; - // std::cout << "displacement =\n"; - // solid.displacement().Print(std::cout); - // std::cout << "forces =\n"; - // solid.reactions().Print(); - // tensor resultant = make_tensor([&](int j) { - // double y = 0; - // for (int n = 0; n < solid.reactions().Size()/dim; n++) { - // y += solid.reactions()(dim*n + j); - // } - // return y; - // }); - // std::cout << "resultant = " << resultant << std::endl; } // Compute norm of error