diff --git a/src/smith/differentiable_numerics/CMakeLists.txt b/src/smith/differentiable_numerics/CMakeLists.txt index 961f480717..e1d526bf8d 100644 --- a/src/smith/differentiable_numerics/CMakeLists.txt +++ b/src/smith/differentiable_numerics/CMakeLists.txt @@ -6,17 +6,20 @@ set(differentiable_numerics_sources field_state.cpp + field_store.cpp differentiable_physics.cpp lumped_mass_explicit_newmark_state_advancer.cpp - solid_mechanics_state_advancer.cpp + solid_mechanics_time_integrator.cpp differentiable_solver.cpp nonlinear_solve.cpp evaluate_objective.cpp dirichlet_boundary_conditions.cpp + multiphysics_time_integrator.cpp ) set(differentiable_numerics_headers field_state.hpp + field_store.hpp state_advancer.hpp reaction.hpp differentiable_solver.hpp @@ -25,15 +28,17 @@ set(differentiable_numerics_headers explicit_dynamic_solve.hpp lumped_mass_weak_form.hpp lumped_mass_explicit_newmark_state_advancer.hpp - solid_mechanics_state_advancer.hpp + solid_mechanics_time_integrator.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 - differentiable_test_utils.hpp + multiphysics_time_integrator.hpp + solid_dynamics_system.hpp + thermo_mechanics_system.hpp + system_base.hpp ) set(differentiable_numerics_depends diff --git a/src/smith/differentiable_numerics/differentiable_physics.cpp b/src/smith/differentiable_numerics/differentiable_physics.cpp index ae091e0aed..7b8aa09c3a 100644 --- a/src/smith/differentiable_numerics/differentiable_physics.cpp +++ b/src/smith/differentiable_numerics/differentiable_physics.cpp @@ -120,12 +120,14 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& dual_nam 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_states_.empty() && !reaction_names_.empty(), + "Reactions were not computed during advanceState, but were requested."); + SLIC_ERROR_IF( - reaction_index >= reaction_names_.size(), + reaction_index >= reaction_states_.size(), "Dual reactions not correctly allocated yet, cannot get dual until after initializationStep is called."); - TimeInfo time_info(time_prev_, dt_prev_, static_cast(cycle_prev_)); - reaction_states_ = advancer_->computeReactions(time_info, *field_shape_displacement_, field_states_, field_params_); return *reaction_states_[reaction_index].get(); } @@ -225,7 +227,10 @@ void DifferentiablePhysics::advanceTimestep(double dt) dt_prev_ = dt; TimeInfo time_info(time_, dt, static_cast(cycle_)); - field_states_ = advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_); + auto [states, reactions] = + advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_); + field_states_ = states; + reaction_states_ = reactions; cycle_++; time_ += dt; diff --git a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp deleted file mode 100644 index 952705d97c..0000000000 --- a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp +++ /dev/null @@ -1,60 +0,0 @@ -// 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_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" -#include "gretl/strumm_walther_checkpoint_strategy.hpp" - -namespace smith { - -/// @brief Helper function to generate the base-physics for solid mechanics -/// @tparam ShapeDispSpace Space for shape displacement, must be H1<1, dim> in most cases -/// @tparam VectorSpace Space for displacement, velocity, acceleration field, typically H1 -/// @tparam ...ParamSpaces Additional parameter spaces, either H1 or L1 -/// @tparam dim Spatial dimension -/// @param mesh smith::Mesh -/// @param d_solid_nonlinear_solver Abstract differentiable solver -/// @param time_rule Time integration rule for second order systems. Likely either quasi-static or implicit Newmark -/// @param num_checkpoints Number of checkpointed states for gretl to store for reverse mode derivatives -/// @param physics_name Name of the physics/WeakForm -/// @param param_names Names for the parameter fields with a one-to-one correspondence with the templated ParamSpaces -/// @return tuple of shared pointers to the: BasePhysics, WeakForm, and DirichetBoundaryConditions -/// @note 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::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, - std::string physics_name, const std::vector& param_names = {}) -{ - auto graph = - std::make_shared(std::make_unique(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 index ca13e1722a..1914894f71 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -232,6 +232,13 @@ std::vector LinearDifferentiableBlockSolver u_bars[static_cast(row_i)]->space(), "u_dual_" + std::to_string(row_i)); } + // If it's a 1x1 system, pass the operator directly to avoid potential BlockOperator issues with smoothers + if (num_rows == 1) { + mfem_solver->SetOperator(*jacobian_transposed[0][0]); + mfem_solver->Mult(*u_bars[0], *u_duals[0]); + return u_duals; + } + mfem::Array block_offsets; block_offsets.SetSize(num_rows + 1); block_offsets[0] = 0; @@ -303,7 +310,7 @@ std::vector NonlinearDifferentiableBlockSol auto residual_op_ = std::make_unique( block_u->Size(), - [&residual_funcs, num_rows, &u_guesses, &block_r](const mfem::Vector& u_, mfem::Vector& r_) { + [&residual_funcs, num_rows, &u_guesses, &block_r, &block_offsets](const mfem::Vector& u_, mfem::Vector& r_) { const mfem::BlockVector* u = dynamic_cast(&u_); SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a block vector"); for (int row_i = 0; row_i < num_rows; ++row_i) { @@ -319,7 +326,11 @@ std::vector NonlinearDifferentiableBlockSol }, [this, &block_offsets, &u_guesses, jacobian_funcs, num_rows](const mfem::Vector& u_) -> mfem::Operator& { const mfem::BlockVector* u = dynamic_cast(&u_); - SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a block vector"); + mfem::BlockVector u_block_wrapper; + if (!u) { + u_block_wrapper.Update(const_cast(u_.GetData()), block_offsets); + u = &u_block_wrapper; + } for (int row_i = 0; row_i < num_rows; ++row_i) { *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); } @@ -359,6 +370,15 @@ std::vector NonlinearDifferentiableBlockSol u_bars[static_cast(row_i)]->space(), "u_dual_" + std::to_string(row_i)); } + auto& linear_solver = nonlinear_solver_->linearSolver(); + + // If it's a 1x1 system, pass the operator directly to avoid potential BlockOperator issues with smoothers + if (num_rows == 1) { + linear_solver.SetOperator(*jacobian_transposed[0][0]); + linear_solver.Mult(*u_bars[0], *u_duals[0]); + return u_duals; + } + mfem::Array block_offsets; block_offsets.SetSize(num_rows + 1); block_offsets[0] = 0; @@ -382,7 +402,6 @@ std::vector NonlinearDifferentiableBlockSol } } - auto& linear_solver = nonlinear_solver_->linearSolver(); linear_solver.SetOperator(*block_jac); linear_solver.Mult(*block_r, *block_ds); diff --git a/src/smith/differentiable_numerics/differentiable_test_utils.hpp b/src/smith/differentiable_numerics/differentiable_test_utils.hpp index 761172d3c2..2f21aacf86 100644 --- a/src/smith/differentiable_numerics/differentiable_test_utils.hpp +++ b/src/smith/differentiable_numerics/differentiable_test_utils.hpp @@ -15,9 +15,28 @@ #include "gretl/double_state.hpp" #include "smith/differentiable_numerics/field_state.hpp" #include "smith/physics/scalar_objective.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" namespace smith { +/** + * @brief Verify that reaction forces are zero at non-Dirichlet DOFs. + * @param reaction The reaction field to check. + * @param bc_manager Boundary condition manager to identify Dirichlet DOFs. + * @param tolerance Absolute tolerance for zero check. + */ +inline void checkUnconstrainedReactions(const FiniteElementDual& reaction, const BoundaryConditionManager& bc_manager, + double tolerance = 1e-8) +{ + FiniteElementState unconstrained_reactions(reaction.space(), "unconstrained_reactions"); + unconstrained_reactions = reaction; + unconstrained_reactions.SetSubVector(bc_manager.allEssentialTrueDofs(), 0.0); + + double max_unconstrained = unconstrained_reactions.Normlinf(); + EXPECT_LT(max_unconstrained, tolerance) + << "Reaction forces should be zero at non-Dirichlet DOFs. Max violation: " << max_unconstrained; +} + /// @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, diff --git a/src/smith/differentiable_numerics/field_store.cpp b/src/smith/differentiable_numerics/field_store.cpp new file mode 100644 index 0000000000..38d6b9d1b9 --- /dev/null +++ b/src/smith/differentiable_numerics/field_store.cpp @@ -0,0 +1,222 @@ +// 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/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "gretl/wang_checkpoint_strategy.hpp" + +namespace smith { + +FieldStore::FieldStore(std::shared_ptr mesh, size_t storage_size) + : mesh_(mesh), + graph_(std::make_shared(std::make_unique(storage_size))) +{ +} + +std::shared_ptr FieldStore::addBoundaryConditions(FEFieldPtr field) +{ + boundary_conditions_.push_back(std::make_shared(mesh_->mfemParMesh(), field->space())); + return boundary_conditions_.back(); +} + +void FieldStore::addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index) +{ + FieldLabel argument_name_and_index{.field_name = argument_name, .field_index_in_residual = argument_index}; + if (weak_form_name_to_unknown_name_index_.count(weak_form_name)) { + weak_form_name_to_unknown_name_index_.at(weak_form_name).push_back(argument_name_and_index); + } else { + weak_form_name_to_unknown_name_index_[weak_form_name] = std::vector{argument_name_and_index}; + } +} + +void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) +{ + // Store the field name instead of index to avoid confusion between states_ and params_ indices + if (weak_form_name_to_field_names_.count(weak_form_name)) { + weak_form_name_to_field_names_.at(weak_form_name).push_back(argument_name); + } else { + weak_form_name_to_field_names_[weak_form_name] = std::vector{argument_name}; + } + SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_names_.at(weak_form_name).size(), + "Invalid order for adding weak form arguments."); +} + +void FieldStore::printMap() +{ + for (auto& keyval : weak_form_name_to_unknown_name_index_) { + std::cout << "for residual: " << keyval.first << " "; + for (auto& name_index : keyval.second) { + std::cout << "arg " << name_index.field_name << " at " << name_index.field_index_in_residual << ", "; + } + std::cout << std::endl; + } +} + +std::vector> FieldStore::indexMap(const std::vector& residual_names) const +{ + std::vector> block_indices(residual_names.size()); + + for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { + std::vector& res_indices = block_indices[res_i]; + res_indices = std::vector(num_unknowns_, invalid_block_index); + const std::string& res_name = residual_names[res_i]; + const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); + + for (const auto& field_name_and_arg_index : arg_info) { + const std::string field_name = field_name_and_arg_index.field_name; + size_t unknown_index = to_unknown_index_.at(field_name); + SLIC_ASSERT(unknown_index < num_unknowns_); + res_indices[unknown_index] = field_name_and_arg_index.field_index_in_residual; + } + } + + return block_indices; +} + +std::vector FieldStore::getBoundaryConditionManagers() const +{ + std::vector bcs; + for (auto& bc : boundary_conditions_) { + bcs.push_back(&bc->getBoundaryConditionManager()); + } + return bcs; +} + +std::shared_ptr FieldStore::getBoundaryConditions(size_t unknown_index) const +{ + SLIC_ERROR_IF(unknown_index >= boundary_conditions_.size(), "Unknown index out of bounds for boundary conditions"); + return boundary_conditions_[unknown_index]; +} + +size_t FieldStore::getFieldIndex(const std::string& field_name) const +{ + if (to_states_index_.count(field_name)) { + return to_states_index_.at(field_name); + } + if (to_params_index_.count(field_name)) { + return to_params_index_.at(field_name); + } + SLIC_ERROR("Field or parameter '" << field_name << "' not found in getFieldIndex"); + return 0; // unreachable +} + +FieldState FieldStore::getField(const std::string& field_name) const +{ + // Check if it's a state field + if (to_states_index_.count(field_name)) { + size_t field_index = to_states_index_.at(field_name); + return states_[field_index]; + } + // Otherwise check if it's a parameter + if (to_params_index_.count(field_name)) { + size_t param_index = to_params_index_.at(field_name); + return params_[param_index]; + } + SLIC_ERROR("Field or parameter '" << field_name << "' not found"); + return states_[0]; // unreachable, but needed for compilation +} + +FieldState FieldStore::getParameter(const std::string& param_name) const +{ + size_t param_index = to_params_index_.at(param_name); + return params_[param_index]; +} + +void FieldStore::setField(const std::string& field_name, FieldState updated_field) +{ + size_t field_index = getFieldIndex(field_name); + states_[field_index] = updated_field; +} + +FieldState FieldStore::getShapeDisp() const { return shape_disp_[0]; } + +const std::vector& FieldStore::getAllFields() const { return states_; } + +std::vector FieldStore::getStates(const std::string& weak_form_name) const +{ + // Validate that weak form is registered + SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0, + axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addResidual()?", + weak_form_name)); + + auto field_names = weak_form_name_to_field_names_.at(weak_form_name); + std::vector fields_for_residual; + for (auto& name : field_names) { + // Validate that field exists + SLIC_ERROR_ROOT_IF( + to_states_index_.count(name) == 0 && to_params_index_.count(name) == 0, + axom::fmt::format("Field '{}' (required by weak form '{}') not found in FieldStore", name, weak_form_name)); + + // Only include state fields, not parameters + // Parameters are passed separately to avoid duplication in block_solve + if (to_states_index_.count(name)) { + fields_for_residual.push_back(getField(name)); + } + } + return fields_for_residual; +} + +std::vector FieldStore::getStatesFromVectors(const std::string& weak_form_name, + const std::vector& state_fields, + const std::vector& param_fields) const +{ + // Validate that weak form is registered + SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0, + axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addResidual()?", + weak_form_name)); + + auto field_names = weak_form_name_to_field_names_.at(weak_form_name); + std::vector fields_for_residual; + for (auto& name : field_names) { + // Check if it's a state field + if (to_states_index_.count(name)) { + size_t idx = to_states_index_.at(name); + SLIC_ERROR_ROOT_IF(idx >= state_fields.size(), + axom::fmt::format("State field index {} out of bounds (size={}) for field '{}'", idx, + state_fields.size(), name)); + fields_for_residual.push_back(state_fields[idx]); + } + // Otherwise check if it's a parameter + else if (to_params_index_.count(name)) { + size_t idx = to_params_index_.at(name); + SLIC_ERROR_ROOT_IF(idx >= param_fields.size(), + axom::fmt::format("Parameter field index {} out of bounds (size={}) for field '{}'", idx, + param_fields.size(), name)); + fields_for_residual.push_back(param_fields[idx]); + } else { + SLIC_ERROR_ROOT(axom::fmt::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore", + name, weak_form_name)); + } + } + return fields_for_residual; +} + +const std::shared_ptr& FieldStore::getMesh() const { return mesh_; } + +const std::shared_ptr& FieldStore::graph() const { return graph_; } + +const std::vector, FieldStore::TimeIntegrationMapping>>& +FieldStore::getTimeIntegrationRules() const +{ + return time_integration_rules_; +} + +size_t FieldStore::getUnknownIndex(const std::string& field_name) const { return to_unknown_index_.at(field_name); } + +void FieldStore::setField(size_t index, FieldState updated_field) { states_[index] = updated_field; } + +void FieldStore::addWeakFormReaction(std::string weak_form_name, std::string field_name) +{ + weak_form_to_test_field_[weak_form_name] = field_name; +} + +std::string FieldStore::getWeakFormReaction(const std::string& weak_form_name) const +{ + return weak_form_to_test_field_.at(weak_form_name); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp new file mode 100644 index 0000000000..ec103a4fd4 --- /dev/null +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -0,0 +1,434 @@ +// 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) + +#pragma once + +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/physics/mesh.hpp" + +#include +#include +#include +#include + +namespace smith { + +class DirichletBoundaryConditions; +class BoundaryConditionManager; + +/** + * @brief Representation of a field type with a name and an optional unknown index. + * @tparam Space The finite element space type. + * @tparam Time The time integration type (unused by default). + */ +template +struct FieldType { + /** + * @brief Construct a new FieldType object. + * @param n Name of the field. + * @param unknown_index_ Index of the unknown in the solver (default: -1). + */ + FieldType(std::string n, int unknown_index_ = -1) : name(n), unknown_index(unknown_index_) {} + std::string name; ///< Name of the field. + int unknown_index; ///< Index of the unknown in the solver. +}; + +/** + * @brief Manages storage and metadata for fields, parameters, and weak forms. + */ +struct FieldStore { + /** + * @brief Construct a new FieldStore object. + * @param mesh The mesh associated with the fields. + * @param storage_size Initial storage size for fields (default: 50). + */ + FieldStore(std::shared_ptr mesh, size_t storage_size = 50); + + /** + * @brief Enum for different types of time derivatives. + */ + enum class TimeDerivative + { + VAL, //< The value of the field. + DOT, ///< The first time derivative. + DDOT, ///< The second time derivative. + DDDOT ///< The third time derivative. + }; + + /** + * @brief Add a shape displacement field to the store. + * @tparam Space The finite element space type. + * @param type The field type specification. + */ + template + void addShapeDisp(FieldType type) + { + shape_disp_.push_back(smith::createFieldState(*graph_, Space{}, type.name, mesh_->tag())); + } + + /** + * @brief Add a parameter field to the store. + * @tparam Space The finite element space type. + * @param type The field type specification. + */ + template + void addParameter(FieldType type) + { + to_params_index_[type.name] = params_.size(); + params_.push_back(smith::createFieldState(*graph_, Space{}, type.name, mesh_->tag())); + } + + /** + * @brief Add an independent field (a solver unknown) to the store. + * + * Registers the field as an unknown and assigns it an index in the solver's block structure. + * The @p type argument is mutated in-place: its @c unknown_index is set to the assigned index, + * so the same @c FieldType object can later be passed to @c createSpaces to tell the + * weak form that this argument is an active unknown (i.e. the Jacobian should be computed + * with respect to it). + * + * @tparam Space The finite element space type. + * @param type The field type specification; @c type.unknown_index is set on return. + * @param time_rule The time integration rule governing how this unknown and its dependents + * are related across time steps. + * @return std::shared_ptr The boundary conditions for this field. + */ + template + std::shared_ptr addIndependent(FieldType& type, + std::shared_ptr time_rule) + { + type.unknown_index = static_cast(num_unknowns_); + to_states_index_[type.name] = states_.size(); + to_unknown_index_[type.name] = num_unknowns_; + FieldState new_field = smith::createFieldState(*graph_, Space{}, type.name, mesh_->tag()); + states_.push_back(new_field); + auto latest_bc = addBoundaryConditions(new_field.get()); + ++num_unknowns_; + SLIC_ERROR_IF(num_unknowns_ != boundary_conditions_.size(), + "Inconcistency between num unknowns and boundary condition size"); + + SLIC_ERROR_IF(!time_rule, "Invalid time_rule"); + + TimeIntegrationMapping mapping; + mapping.primary_name = type.name; + independent_name_to_rule_index_[type.name] = time_integration_rules_.size(); + time_integration_rules_.push_back({time_rule, mapping}); + + return latest_bc; + } + + /** + * @brief Add a dependent field (history value, velocity, or acceleration) to the store. + * + * Creates and registers a new field that carries the previous time-step value of a particular + * time derivative of an independent field. The relationship is recorded in the + * @c TimeIntegrationMapping for the parent independent field so that, at evaluation time, the + * time integration rule can reconstruct the current rate from the pair + * (predicted_value, stored_old_value). + * + * **Return value:** a @c FieldType whose @c name is the name of the newly registered + * field. This object is intentionally returned (rather than discarded) so that callers can + * pass it directly to @c createSpaces when assembling the weak-form argument list. The + * position at which a @c FieldType appears in that argument list determines which slot of the + * lambda the field occupies at quadrature-point evaluation time, which is how the time + * integration rules reconstruct quantities such as + * @f$ \dot{\alpha} = (\alpha_\text{predicted} - \alpha_\text{old}) / \Delta t @f$. + * + * The field name is derived automatically from the independent field's name plus a suffix that + * reflects the derivative level (@c _old for VALUE, @c _dot_old for DOT, + * @c _ddot_old for DDOT), unless @p name_override is supplied. + * + * @tparam Space The finite element space type (must match the independent field). + * @param independent_field The @c FieldType of the independent (predicted) field. + * @param derivative Which time-derivative level this history field stores. + * @param name_override If non-empty, use this as the field name instead of the auto-generated one. + * @return FieldType Type descriptor for the newly created dependent field; pass this to + * @c createSpaces to register it as a weak-form argument. + */ + template + auto addDependent(FieldType independent_field, TimeDerivative derivative, std::string name_override = "") + { + std::string suffix; + if (derivative == TimeDerivative::VAL) { + suffix = "_old"; + } else if (derivative == TimeDerivative::DOT) { + suffix = "_dot_old"; + } else if (derivative == TimeDerivative::DDOT) { + suffix = "_ddot_old"; + } else { + SLIC_ERROR("Unsupported TimeDerivative"); + } + + std::string name = name_override.empty() ? independent_field.name + suffix : name_override; + + if (independent_name_to_rule_index_.count(independent_field.name)) { + size_t rule_idx = independent_name_to_rule_index_.at(independent_field.name); + auto& mapping = time_integration_rules_[rule_idx].second; + if (derivative == TimeDerivative::VAL) { + mapping.history_name = name; + } else if (derivative == TimeDerivative::DOT) { + mapping.dot_name = name; + } else if (derivative == TimeDerivative::DDOT) { + mapping.ddot_name = name; + } + } else { + SLIC_WARNING("Adding dependent time integration field for independent field '" + << independent_field.name << "' which has no registered TimeIntegrationRule."); + } + + to_states_index_[name] = states_.size(); + states_.push_back(smith::createFieldState(*graph_, Space{}, name, mesh_->tag())); + return FieldType(name); + } + + /** + * @brief Register an argument to a weak form as an unknown. + * @param weak_form_name Name of the weak form. + * @param argument_name Name of the argument field. + * @param argument_index Index of the argument in the weak form's argument list. + */ + void addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index); + + /** + * @brief Register an argument to a weak form. + * @param weak_form_name Name of the weak form. + * @param argument_name Name of the argument field. + * @param argument_index Index of the argument in the weak form's argument list. + */ + void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index); + + /** + * @brief Register the reaction (test) field for a weak form. + * + * The reaction field is the field whose test function space the weak form integrates against. + * It determines which field's degrees of freedom the assembled residual is "returned to" + * (i.e. the field whose force/flux vector is populated). + * + * @param weak_form_name Name of the weak form. + * @param field_name Name of the reaction field. + */ + void addWeakFormReaction(std::string weak_form_name, std::string field_name); + + /** + * @brief Get the name of the reaction (test) field for a weak form. + * @param weak_form_name Name of the weak form. + * @return std::string Name of the reaction field. + */ + std::string getWeakFormReaction(const std::string& weak_form_name) const; + + /** + * @brief Register all input fields for a weak form and return their FE spaces. + * + * This is the primary setup method for constructing a weak form. It: + * 1. Registers @p reaction_field_name as the reaction/test field via @c addWeakFormReaction. + * 2. Iterates over every @c FieldType in @p types (in order), registering each as an input + * argument to the weak form and recording whether it is an active unknown. + * 3. Returns the ordered vector of finite element spaces, which can be passed directly to + * the @c TimeDiscretizedWeakForm constructor without creating a named temporary. + * + * @param weak_form_name Name of the weak form being constructed. + * @param reaction_field_name Name of the test/reaction field (may differ from the first input). + * @param types Ordered list of @c FieldType descriptors for every input argument. + * @return std::vector Ordered input FE spaces. + */ + template + std::vector createSpaces(const std::string& weak_form_name, + const std::string& reaction_field_name, + FieldTypes... types) + { + addWeakFormReaction(weak_form_name, reaction_field_name); + std::vector spaces; + size_t arg_num = 0; + auto register_field = [&](auto type) { + spaces.push_back(&getField(type.name).get()->space()); + addWeakFormArg(weak_form_name, type.name, arg_num); + if (type.unknown_index >= 0) { + addWeakFormUnknownArg(weak_form_name, type.name, arg_num); + } + ++arg_num; + }; + (register_field(types), ...); + return spaces; + } + + /** + * @brief Mapping between primary and history/derivative fields for time integration. + */ + struct TimeIntegrationMapping { + std::string primary_name; ///< Primary unknown field name. + std::string history_name; ///< Previous time step value field name. + std::string dot_name; ///< First time derivative field name. + std::string ddot_name; ///< Second time derivative field name. + }; + + /** + * @brief Get all registered time integration rules and their mappings. + * @return const std::vector, TimeIntegrationMapping>>& List of rules + * and mappings. + */ + const std::vector, TimeIntegrationMapping>>& getTimeIntegrationRules() + const; + + /** + * @brief Print the internal field maps for debugging. + */ + void printMap(); + + /** + * @brief Generate an index map for the residuals. + * @param residual_names Names of the residuals. + * @return std::vector> The index map. + */ + std::vector> indexMap(const std::vector& residual_names) const; + + /** + * @brief Get the boundary condition managers for all independent fields. + * @return std::vector List of boundary condition managers. + */ + std::vector getBoundaryConditionManagers() const; + + /** + * @brief Get the Dirichlet boundary conditions for an independent field by its unknown index. + * @param unknown_index The unknown index of the independent field. + * @return std::shared_ptr The boundary conditions. + */ + std::shared_ptr getBoundaryConditions(size_t unknown_index) const; + + /** + * @brief Get the internal index of a field by name. + * @param field_name Name of the field. + * @return size_t Index of the field. + */ + size_t getFieldIndex(const std::string& field_name) const; + + /** + * @brief Get the unknown index of a field by name. + * @param field_name Name of the field. + * @return size_t Unknown index of the field. + */ + size_t getUnknownIndex(const std::string& field_name) const; + + /** + * @brief Get a FieldState by name. + * @param field_name Name of the field. + * @return FieldState The field state. + */ + FieldState getField(const std::string& field_name) const; + + /** + * @brief Get a parameter field by name. + * @param param_name Name of the parameter. + * @return FieldState The parameter field state. + */ + FieldState getParameter(const std::string& param_name) const; + + /** + * @brief Update a field in the store by name. + * @param field_name Name of the field. + * @param updated_field The new field state. + */ + void setField(const std::string& field_name, FieldState updated_field); + + /** + * @brief Update a field in the store by index. + * @param index Index of the field. + * @param updated_field The new field state. + */ + void setField(size_t index, FieldState updated_field); + + /** + * @brief Get the shape displacement field. + * @return FieldState The shape displacement field. + */ + FieldState getShapeDisp() const; + + /** + * @brief Get all fields stored in the FieldStore. + * @return const std::vector& List of all fields. + */ + const std::vector& getAllFields() const; + + /** + * @brief Get the state fields associated with a weak form. + * @param weak_form_name Name of the weak form. + * @return std::vector List of state fields. + */ + std::vector getStates(const std::string& weak_form_name) const; + + /** + * @brief Extract state fields for a weak form from provided state and parameter vectors. + * @param weak_form_name Name of the weak form. + * @param state_fields Vector of all state fields. + * @param param_fields Vector of all parameter fields. + * @return std::vector Subset of fields relevant to the weak form. + */ + std::vector getStatesFromVectors(const std::string& weak_form_name, + const std::vector& state_fields, + const std::vector& param_fields) const; + + /** + * @brief Get the associated mesh. + * @return const std::shared_ptr& The mesh. + */ + const std::shared_ptr& getMesh() const; + + /** + * @brief Get the associated data store graph. + * @return const std::shared_ptr& The graph. + */ + const std::shared_ptr& graph() const; + + private: + std::shared_ptr mesh_; + std::shared_ptr graph_; + + std::vector shape_disp_; + std::vector params_; + std::vector states_; + + std::map to_states_index_; + std::map to_params_index_; + + size_t num_unknowns_ = 0; + std::map to_unknown_index_; + std::vector> boundary_conditions_; + + struct FieldLabel { + std::string field_name; + size_t field_index_in_residual; + }; + + std::shared_ptr addBoundaryConditions(FEFieldPtr field); + + std::map> weak_form_name_to_unknown_name_index_; + + std::map> weak_form_name_to_field_indices_; + std::map> weak_form_name_to_field_names_; + + std::map weak_form_to_test_field_; + + std::vector, TimeIntegrationMapping>> time_integration_rules_; + std::map independent_name_to_rule_index_; +}; + +/** + * @brief Create a TimeDiscretizedWeakForm and register its fields in the FieldStore. + * + * Thin convenience wrapper: registers @p test_type as the reaction field, registers all + * @p field_types as input arguments, and constructs the weak form in one call. + */ +template +auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, + FieldType... field_types) +{ + return std::make_shared>>( + name, field_store.getMesh(), field_store.getField(test_type.name).get()->space(), + field_store.createSpaces(name, test_type.name, field_types...)); +} + +} // namespace smith 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 61b9a6f80e..35ab7f5c87 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,7 +33,7 @@ FieldState applyZeroBoundaryConditions(const FieldState& s, const BoundaryCondit return s_bc.finalize(); } -std::vector LumpedMassExplicitNewmarkStateAdvancer::advanceState( +std::pair, std::vector> LumpedMassExplicitNewmarkStateAdvancer::advanceState( const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states_in, const std::vector& params) const { @@ -110,7 +110,7 @@ std::vector LumpedMassExplicitNewmarkStateAdvancer::advanceState( } // place all solved updated states into the output - return states; + return {states, std::vector{}}; } } // namespace smith 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 1dd292617b..d2b9585505 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,9 +36,9 @@ class LumpedMassExplicitNewmarkStateAdvancer : public StateAdvancer { } /// @overload - std::vector advanceState(const TimeInfo& time_info, const FieldState& shape_disp, - const std::vector& states, - const std::vector& params) const override; + std::pair, 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/multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp new file mode 100644 index 0000000000..560a0f3972 --- /dev/null +++ b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp @@ -0,0 +1,126 @@ +// 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/differentiable_numerics/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/reaction.hpp" + +namespace smith { + +MultiphysicsTimeIntegrator::MultiphysicsTimeIntegrator(std::shared_ptr field_store, + const std::vector>& weak_forms, + std::shared_ptr solver) + : field_store_(field_store), weak_forms_(weak_forms), solver_(solver) +{ +} + +std::pair, std::vector> MultiphysicsTimeIntegrator::advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states, + const std::vector& params) const +{ + // Sync FieldStore with input states + for (size_t i = 0; i < states.size(); ++i) { + field_store_->setField(i, states[i]); + } + + std::vector primary_unknowns = solve(weak_forms_, *field_store_, solver_.get(), time_info, params); + + // Create states for reaction computation: newly solved primary unknowns + original input states + std::vector states_for_reactions = states; + for (const auto& [rule, mapping] : field_store_->getTimeIntegrationRules()) { + size_t u_idx = field_store_->getFieldIndex(mapping.primary_name); + size_t unknown_idx = field_store_->getUnknownIndex(mapping.primary_name); + FieldState u_new = primary_unknowns[unknown_idx]; + states_for_reactions[u_idx] = u_new; + } + + // Compute reactions using newly solved unknowns but BEFORE time integration state updates + std::vector reactions; + for (const auto& wf : weak_forms_) { + std::vector wf_fields = field_store_->getStatesFromVectors(wf->name(), states_for_reactions, params); + std::string test_field_name = field_store_->getWeakFormReaction(wf->name()); + size_t test_field_idx = field_store_->getFieldIndex(test_field_name); + FieldState test_field = states_for_reactions[test_field_idx]; + reactions.push_back(smith::evaluateWeakForm(wf, time_info, shape_disp, wf_fields, test_field)); + } + + // Now do time integration to compute corrected velocities/accelerations and update all states + const auto& all_current_states = field_store_->getAllFields(); + std::vector new_states = states; + for (size_t i = 0; i < states.size(); ++i) { + new_states[i] = all_current_states[i]; + } + + for (const auto& [rule, mapping] : field_store_->getTimeIntegrationRules()) { + size_t u_idx = field_store_->getFieldIndex(mapping.primary_name); + size_t unknown_idx = field_store_->getUnknownIndex(mapping.primary_name); + FieldState u_new = primary_unknowns[unknown_idx]; + new_states[u_idx] = u_new; + + std::vector rule_inputs; + rule_inputs.push_back(u_new); // u_{n+1} + if (rule->num_args() >= 2) { + rule_inputs.push_back(states[u_idx]); // u_n + } + + if (rule->num_args() >= 3 && !mapping.dot_name.empty()) { + size_t v_idx = field_store_->getFieldIndex(mapping.dot_name); + rule_inputs.push_back(states[v_idx]); + } + + if (rule->num_args() >= 4 && !mapping.ddot_name.empty()) { + size_t a_idx = field_store_->getFieldIndex(mapping.ddot_name); + rule_inputs.push_back(states[a_idx]); + } + + if (!mapping.dot_name.empty()) { + size_t v_idx = field_store_->getFieldIndex(mapping.dot_name); + new_states[v_idx] = rule->corrected_dot(time_info, rule_inputs); + } + + if (!mapping.ddot_name.empty()) { + size_t a_idx = field_store_->getFieldIndex(mapping.ddot_name); + new_states[a_idx] = rule->corrected_ddot(time_info, rule_inputs); + } + + if (!mapping.history_name.empty()) { + size_t hist_idx = field_store_->getFieldIndex(mapping.history_name); + new_states[hist_idx] = u_new; + } + } + + return {new_states, reactions}; +} + +std::vector solve(const std::vector>& weak_forms, const FieldStore& field_store, + const DifferentiableBlockSolver* solver, const TimeInfo& time_info, + const std::vector& params) +{ + std::vector weak_form_names; + for (const auto& wf : weak_forms) { + weak_form_names.push_back(wf->name()); + } + std::vector> index_map = field_store.indexMap(weak_form_names); + + std::vector> inputs; + for (size_t i = 0; i < weak_forms.size(); ++i) { + std::string wf_name = weak_forms[i]->name(); + std::vector fields_for_wk = field_store.getStates(wf_name); + inputs.push_back(fields_for_wk); + } + std::vector> wk_params(weak_forms.size(), params); + + std::vector weak_form_ptrs; + for (auto& p : weak_forms) { + weak_form_ptrs.push_back(p.get()); + } + return block_solve(weak_form_ptrs, index_map, field_store.getShapeDisp(), inputs, wk_params, time_info, solver, + field_store.getBoundaryConditionManagers()); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/multiphysics_time_integrator.hpp b/src/smith/differentiable_numerics/multiphysics_time_integrator.hpp new file mode 100644 index 0000000000..6d85dcf449 --- /dev/null +++ b/src/smith/differentiable_numerics/multiphysics_time_integrator.hpp @@ -0,0 +1,70 @@ +// 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) + +#pragma once + +#include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/physics/mesh.hpp" +#include "smith/differentiable_numerics/field_store.hpp" + +namespace smith { + +class DifferentiableBlockSolver; +class DirichletBoundaryConditions; +class BoundaryConditionManager; + +/** + * @brief Solve a set of weak forms. + * @param weak_forms List of weak forms to solve. + * @param field_store Field store containing the fields. + * @param solver The solver to use. + * @param time_info Current time information. + * @param params Optional parameter fields. + * @return std::vector The updated state fields. + */ +std::vector solve(const std::vector>& weak_forms, const FieldStore& field_store, + const DifferentiableBlockSolver* solver, const TimeInfo& time_info, + const std::vector& params = {}); + +/** + * @brief Time integrator for multiphysics problems, coordinating multiple weak forms. + */ +class MultiphysicsTimeIntegrator : public StateAdvancer { + public: + /** + * @brief Construct a new MultiphysicsTimeIntegrator object. + * @param field_store Field store containing the fields. + * @param weak_forms List of weak forms to coordinate. + * @param solver The block solver to use. + */ + MultiphysicsTimeIntegrator(std::shared_ptr field_store, + const std::vector>& weak_forms, + std::shared_ptr solver); + + /** + * @brief Advance the multiphysics state by one time step. + * @param time_info Current time information. + * @param shape_disp Shape displacement field. + * @param states Current state fields. + * @param params Parameter fields. + * @return std::pair, std::vector> Updated states and reactions. + */ + std::pair, std::vector> advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states, + const std::vector& params) const override; + + private: + std::shared_ptr field_store_; + std::vector> weak_forms_; + std::shared_ptr solver_; +}; + +} // namespace smith diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 37777a1734..d634897658 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -289,6 +289,22 @@ std::vector block_solve(const std::vector& residual_evals SLIC_ERROR_IF(num_rows_ != block_indices[r].size(), "All block index rows must have the same number of columns"); } + // Validate block_indices bounds against states sizes + for (size_t row = 0; row < num_rows_; ++row) { + for (size_t col = 0; col < num_rows_; ++col) { + size_t idx = block_indices[row][col]; + if (idx != invalid_block_index) { + SLIC_ERROR_IF(idx >= states[row].size(), + axom::fmt::format("block_indices[{}][{}] = {} is out of bounds (states[{}].size() = {})", row, + col, idx, row, states[row].size())); + } + } + // Validate that diagonal entry is not invalid + SLIC_ERROR_IF( + block_indices[row][row] == invalid_block_index, + axom::fmt::format("block_indices[{}][{}] (diagonal entry) must not be invalid_block_index", row, row)); + } + std::vector num_state_inputs; std::vector allFields; for (auto& ss : states) { @@ -520,7 +536,11 @@ std::vector block_solve(const std::vector& residual_evals 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())); + if (jacobians[row_i][col_j]) { + jacobians_T[col_j].emplace_back(std::unique_ptr(jacobians[row_i][col_j]->Transpose())); + } else { + jacobians_T[col_j].emplace_back(nullptr); + } } } for (size_t row_i = 0; row_i < num_rows; ++row_i) { diff --git a/src/smith/differentiable_numerics/solid_dynamics_system.hpp b/src/smith/differentiable_numerics/solid_dynamics_system.hpp new file mode 100644 index 0000000000..bc6a2cbbe9 --- /dev/null +++ b/src/smith/differentiable_numerics/solid_dynamics_system.hpp @@ -0,0 +1,372 @@ +// 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 solid_dynamics_system.hpp + * @brief Defines the SolidDynamicsSystem struct and its factory function + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/differentiable_numerics/solid_mechanics_time_integrator.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { + +/** + * @brief System struct for solid dynamics with second-order time integration. + * @tparam dim Spatial dimension. + * @tparam order Polynomial order for displacement field. + * @tparam parameter_space Parameter spaces for material properties. + */ +template +struct SolidDynamicsSystem : public SystemBase { + /// using + using SolidWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, H1, H1, H1, parameter_space...>>; + + /// using + using CycleZeroWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, H1, H1, H1, parameter_space...>>; + + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + std::shared_ptr + cycle_zero_weak_form; ///< Cycle-zero weak form for initial acceleration solve. + std::shared_ptr disp_bc; ///< Displacement boundary conditions. + std::shared_ptr time_rule; ///< Time integration rule. + + /** + * @brief Get the list of all state fields (displacement, displacement_old, velocity_old, acceleration_old). + * @return std::vector List of state fields. + */ + std::vector getStateFields() const + { + return {field_store->getField(prefix("displacement_predicted")), field_store->getField(prefix("displacement")), + field_store->getField(prefix("velocity")), field_store->getField(prefix("acceleration"))}; + } + + /** + * @brief Create a DifferentiablePhysics object for this system. + * @param physics_name The name of the physics. + * @return std::shared_ptr The differentiable physics object. + */ + std::shared_ptr createDifferentiablePhysics(std::string physics_name) + { + return std::make_shared( + field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), + getParameterFields(), advancer, physics_name, std::vector{prefix("reactions")}); + } + + /** + * @brief Set the material model for a domain, defining integrals for the solid weak form. + * @tparam MaterialType The material model type. + * @param material The material model instance. + * @param domain_name The name of the domain to apply the material to. + */ + template + void setMaterial(const MaterialType& material, const std::string& domain_name) + { + // Add to solid weak form (inputs: u, u_old, v_old, a_old, params...) + // Manually apply time integration rule to compute current state + auto captured_rule = time_rule; + solid_weak_form->addBodyIntegral( + domain_name, [=](auto t_info, auto /*X*/, auto u, auto u_old, auto v_old, auto a_old, auto... params) { + // Apply time integration rule to compute current state from history + auto u_current = captured_rule->value(t_info, u, u_old, v_old, a_old); + auto a_current = captured_rule->ddot(t_info, u, u_old, v_old, a_old); + + typename MaterialType::State state; + auto pk_stress = material(state, get(u_current), params...); + + return smith::tuple{get(a_current) * material.density, pk_stress}; + }); + + // Add to cycle-zero weak form (inputs: u, u_old, v_old, a, params...) + // At cycle 0, we directly use u, v, a (no time integration needed) + cycle_zero_weak_form->addBodyIntegral( + domain_name, [=](auto /*t_info*/, auto /*X*/, auto u, auto /*u_old*/, auto /*v_old*/, auto a, auto... params) { + typename MaterialType::State state; + auto pk_stress = material(state, get(u), params...); + + return smith::tuple{get(a) * material.density, pk_stress}; + }); + } + + /** + * @brief Add a body force to the system (with DependsOn). + * @tparam active_parameters Indices of fields this force depends on. + * @tparam BodyForceType The body force function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the domain to apply the force to. + * @param force_function The force function (t, X, selected time-integrated inputs...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addBodyForce(DependsOn depends_on, const std::string& domain_name, + BodyForceType force_function) + { + auto captured_rule = time_rule; + solid_weak_form->addBodySource( + depends_on, domain_name, [=](auto t_info, auto X, auto u, auto u_old, auto v_old, auto a_old, auto... params) { + // Apply time integration rule to get current state + auto u_current = captured_rule->value(t_info, u, u_old, v_old, a_old); + auto v_current = captured_rule->dot(t_info, u, u_old, v_old, a_old); + auto a_current = captured_rule->ddot(t_info, u, u_old, v_old, a_old); + + return force_function(t_info.time(), X, u_current, v_current, a_current, params...); + }); + + cycle_zero_weak_form->addBodySource( + depends_on, domain_name, [=](auto t_info, auto X, auto u, auto /*u_old*/, auto v_old, auto a, auto... params) { + return force_function(t_info.time(), X, u, v_old, a, params...); + }); + } + + /** + * @brief Add a body force to the system. + * @tparam BodyForceType The body force function type. + * @param domain_name The name of the domain to apply the force to. + * @param force_function The force function (t, X, u, v, a, params...). + */ + template + void addBodyForce(const std::string& domain_name, BodyForceType force_function) + { + addBodyForceAllParams(domain_name, force_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a surface traction (flux) to the system (with DependsOn). + * @tparam active_parameters Indices of fields this traction depends on. + * @tparam TractionType The traction function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the boundary domain to apply the traction to. + * @param traction_function The traction function (t, X, n, selected time-integrated inputs...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addTraction(DependsOn depends_on, const std::string& domain_name, + TractionType traction_function) + { + auto captured_rule = time_rule; + solid_weak_form->addBoundaryFlux( + depends_on, domain_name, + [=](auto t_info, auto X, auto n, auto u, auto u_old, auto v_old, auto a_old, auto... params) { + // Apply time integration rule to get current state + auto u_current = captured_rule->value(t_info, u, u_old, v_old, a_old); + auto v_current = captured_rule->dot(t_info, u, u_old, v_old, a_old); + auto a_current = captured_rule->ddot(t_info, u, u_old, v_old, a_old); + + return traction_function(t_info.time(), X, n, u_current, v_current, a_current, params...); + }); + + cycle_zero_weak_form->addBoundaryFlux( + depends_on, domain_name, + [=](auto t_info, auto X, auto n, auto u, auto /*u_old*/, auto v_old, auto a, auto... params) { + return traction_function(t_info.time(), X, n, u, v_old, a, params...); + }); + } + + /** + * @brief Add a surface traction (flux) to the system. + * @tparam TractionType The traction function type. + * @param domain_name The name of the boundary domain to apply the traction to. + * @param traction_function The traction function (t, X, n, u, v, a, params...). + */ + template + void addTraction(const std::string& domain_name, TractionType traction_function) + { + addTractionAllParams(domain_name, traction_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a pressure boundary condition (follower force) (with DependsOn). + * @tparam active_parameters Indices of fields this pressure depends on. + * @tparam PressureType The pressure function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the boundary domain. + * @param pressure_function The pressure function (t, X, selected time-integrated inputs...). + * @note Pressure is applied in the current configuration: P * n_deformed. + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addPressure(DependsOn depends_on, const std::string& domain_name, + PressureType pressure_function) + { + auto captured_rule = time_rule; + solid_weak_form->addBoundaryIntegral( + depends_on, domain_name, [=](auto t_info, auto X, auto u, auto u_old, auto v_old, auto a_old, auto... params) { + // Apply time integration rule to get current state + auto u_current = captured_rule->value(t_info, u, u_old, v_old, a_old); + + // Compute deformed normal and apply correction for reference configuration integration + auto x_current = X + u_current; + auto n_deformed = cross(get(x_current)); + auto n_shape_norm = norm(cross(get(X))); + + auto pressure = pressure_function(t_info.time(), get(X), get(params)...); + + return pressure * n_deformed * (1.0 / n_shape_norm); + }); + + cycle_zero_weak_form->addBoundaryIntegral( + depends_on, domain_name, + [=](auto t_info, auto X, auto u, auto /*u_old*/, auto /*v_old*/, auto /*a*/, auto... params) { + // At cycle 0, u is the current displacement + auto u_current = u; + + // Compute deformed normal and apply correction for reference configuration integration + auto x_current = X + u_current; + auto n_deformed = cross(get(x_current)); + auto n_shape_norm = norm(cross(get(X))); + + auto pressure = pressure_function(t_info.time(), get(X), get(params)...); + + return pressure * n_deformed * (1.0 / n_shape_norm); + }); + } + + /** + * @brief Add a pressure boundary condition (follower force). + * @tparam PressureType The pressure function type. + * @param domain_name The name of the boundary domain. + * @param pressure_function The pressure function (t, X, params...). + * @note Pressure is applied in the current configuration: P * n_deformed. + */ + template + void addPressure(const std::string& domain_name, PressureType pressure_function) + { + addPressureAllParams(domain_name, pressure_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + private: + // Helper functions to forward non-DependsOn calls to DependsOn versions with all parameters + template + void addBodyForceAllParams(const std::string& domain_name, BodyForceType force_function, std::index_sequence) + { + addBodyForce(DependsOn(Is)...>{}, domain_name, force_function); + } + + template + void addTractionAllParams(const std::string& domain_name, TractionType traction_function, std::index_sequence) + { + addTraction(DependsOn(Is)...>{}, domain_name, traction_function); + } + + template + void addPressureAllParams(const std::string& domain_name, PressureType pressure_function, std::index_sequence) + { + addPressure(DependsOn(Is)...>{}, domain_name, pressure_function); + } +}; + +/** + * @brief Factory function to build a solid dynamics system with second-order time integration. + * @tparam dim Spatial dimension. + * @tparam order Polynomial order for displacement field. + * @tparam parameter_space Parameter spaces for material properties. + * @param mesh The mesh. + * @param solver The differentiable block solver. + * @param time_rule The time integration rule. + * @param prepend_name The name of the physics (used as field prefix). + * @param parameter_types Parameter field types. + * @return SolidDynamicsSystem with all components initialized. + */ +template +SolidDynamicsSystem buildSolidDynamicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, std::string prepend_name = "", + FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100); + + auto prefix = [&](const std::string& name) { + if (prepend_name.empty()) { + return name; + } + return prepend_name + "_" + name; + }; + + // Add shape displacement + FieldType> shape_disp_type(prefix("shape_displacement")); + field_store->addShapeDisp(shape_disp_type); + + // Add displacement as independent (unknown) with time integration rule + auto time_rule_ptr = std::make_shared(time_rule); + FieldType> disp_type(prefix("displacement_predicted")); + auto disp_bc = field_store->addIndependent(disp_type, time_rule_ptr); + + // Add dependent fields for time integration history + auto disp_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::VAL, prefix("displacement")); + auto velo_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::DOT, prefix("velocity")); + auto accel_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::DDOT, prefix("acceleration")); + + // Add parameters + std::vector parameter_fields; + (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); + (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + + // Create solid mechanics weak form (u, u_old, v_old, a_old) + std::string force_name = prefix("solid_force"); + auto solid_weak_form = + std::make_shared::SolidWeakFormType>( + force_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), + field_store->createSpaces(force_name, disp_type.name, disp_type, disp_old_type, velo_old_type, accel_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Create cycle-zero weak form (u, u_old, v_old, a) for initial acceleration solve + // Note: We solve R(u_0, v_0, a_0) = 0 for a_0. + // We include history terms to keep trial space indices consistent with solid_weak_form. + std::string cycle_zero_name = prefix("solid_reaction"); + auto cycle_zero_weak_form = + std::make_shared::CycleZeroWeakFormType>( + cycle_zero_name, field_store->getMesh(), field_store->getField(accel_old_type.name).get()->space(), + field_store->createSpaces(cycle_zero_name, accel_old_type.name, disp_type, disp_old_type, velo_old_type, + accel_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Build advancer using SolidMechanicsTimeIntegrator which wraps MultiphysicsTimeIntegrator + // and handles the initial acceleration solve at cycle=0 + auto advancer = + std::make_shared(field_store, solid_weak_form, cycle_zero_weak_form, solver); + + return SolidDynamicsSystem{ + {field_store, solver, advancer, parameter_fields, prepend_name}, + solid_weak_form, + cycle_zero_weak_form, + disp_bc, + time_rule_ptr}; +} + +/** + * @brief Factory function to build a solid dynamics system (without physics name). + * @tparam dim Spatial dimension. + * @tparam order Polynomial order for displacement field. + * @tparam parameter_space Parameter spaces for material properties. + * @param mesh The mesh. + * @param solver The differentiable block solver. + * @param time_rule The time integration rule. + * @param parameter_types Parameter field types. + * @return SolidDynamicsSystem with all components initialized. + */ +template +SolidDynamicsSystem buildSolidDynamicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, FieldType... parameter_types) +{ + return buildSolidDynamicsSystem(mesh, solver, time_rule, "", parameter_types...); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp deleted file mode 100644 index 67e2b66a97..0000000000 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ /dev/null @@ -1,66 +0,0 @@ -// 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/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::ImplicitNewmarkSecondOrderTimeIntegrationRule 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_->final_reaction_weak_form, shape_disp, states_old_, - params, time_info, *solver_, *vector_bcs_, ACCELERATION); - } - - 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, - time_info, *solver_, *vector_bcs_); - - std::vector states = states_old; - - states[DISPLACEMENT] = displacement; - states[VELOCITY] = - time_rule_.dot(time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]); - states[ACCELERATION] = time_rule_.ddot(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_->final_reaction_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 deleted file mode 100644 index 1218724912..0000000000 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp +++ /dev/null @@ -1,119 +0,0 @@ -// 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 solid_mechanics_state_advancer.hpp - * .hpp - * - * @brief Specifies parameterized 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, - ImplicitNewmarkSecondOrderTimeIntegrationRule 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 = createFieldState(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, - ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, std::string physics_name, - const std::vector& param_names, double initial_time = 0.0) - { - auto shape_disp = createFieldState(*graph, ShapeDispSpace{}, physics_name + "_shape_displacement", mesh->tag()); - auto disp = createFieldState(*graph, VectorSpace{}, physics_name + "_displacement", mesh->tag()); - auto velo = createFieldState(*graph, VectorSpace{}, physics_name + "_velocity", mesh->tag()); - auto acceleration = createFieldState(*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. - ImplicitNewmarkSecondOrderTimeIntegrationRule - 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/solid_mechanics_time_integrator.cpp b/src/smith/differentiable_numerics/solid_mechanics_time_integrator.cpp new file mode 100644 index 0000000000..1e9dc4b655 --- /dev/null +++ b/src/smith/differentiable_numerics/solid_mechanics_time_integrator.cpp @@ -0,0 +1,86 @@ +// 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/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_time_integrator.hpp" +#include "smith/differentiable_numerics/reaction.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" + +namespace smith { + +SolidMechanicsTimeIntegrator::SolidMechanicsTimeIntegrator(std::shared_ptr field_store, + std::shared_ptr solid_weak_form, + std::shared_ptr cycle_zero_weak_form, + std::shared_ptr solver) + : field_store_(field_store), cycle_zero_weak_form_(cycle_zero_weak_form), solver_(solver) +{ + std::vector> weak_forms = {solid_weak_form}; + integrator_ = std::make_shared(field_store, weak_forms, solver); +} + +std::pair, std::vector> SolidMechanicsTimeIntegrator::advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states, + const std::vector& params) const +{ + // Handle initial acceleration solve at cycle 0 + std::vector current_states = states; + if (time_info.cycle() == 0 && cycle_zero_weak_form_) { + // Sync FieldStore with input states + for (size_t i = 0; i < states.size(); ++i) { + field_store_->setField(i, states[i]); + } + + // The test field of cycle_zero_weak_form is the field we need to solve for (acceleration) + std::string test_field_name = field_store_->getWeakFormReaction(cycle_zero_weak_form_->name()); + + // Get the weak form's state fields + std::vector wf_fields = field_store_->getStates(cycle_zero_weak_form_->name()); + + // Find which argument index corresponds to the test field (the unknown we're solving for) + FieldState test_field = field_store_->getField(test_field_name); + size_t test_field_idx_in_wf = invalid_block_index; + for (size_t j = 0; j < wf_fields.size(); ++j) { + if (wf_fields[j].get() == test_field.get()) { + test_field_idx_in_wf = j; + break; + } + } + SLIC_ERROR_IF(test_field_idx_in_wf == invalid_block_index, "Test field '" << test_field_name + << "' not found in weak form '" + << cycle_zero_weak_form_->name() << "'"); + + // Set up block solve for this single unknown + std::vector wf_ptrs = {cycle_zero_weak_form_.get()}; + std::vector> block_indices = {{test_field_idx_in_wf}}; + + // Get boundary conditions (assume first BC manager for now) + std::vector bcs; + auto all_bcs = field_store_->getBoundaryConditionManagers(); + if (!all_bcs.empty()) { + bcs.push_back(all_bcs[0]); + } + + std::vector> states_vec = {wf_fields}; + std::vector> params_vec = {params}; + + auto result = + block_solve(wf_ptrs, block_indices, shape_disp, states_vec, params_vec, time_info, solver_.get(), bcs); + + // Update the acceleration field in our current states + size_t test_field_state_idx = field_store_->getFieldIndex(test_field_name); + current_states[test_field_state_idx] = result[0]; + } + + // Now perform the regular time step + auto [new_states, reactions] = integrator_->advanceState(time_info, shape_disp, current_states, params); + + return {new_states, reactions}; +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_time_integrator.hpp b/src/smith/differentiable_numerics/solid_mechanics_time_integrator.hpp new file mode 100644 index 0000000000..2982f62e47 --- /dev/null +++ b/src/smith/differentiable_numerics/solid_mechanics_time_integrator.hpp @@ -0,0 +1,111 @@ +// 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 solid_mechanics_time_integrator.hpp + * + * @brief Specifies parameterized residuals and various linearized evaluations for arbitrary nonlinear systems of + * equations + */ + +#pragma once + +#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/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" + +namespace smith { + +class DifferentiableBlockSolver; +class DirichletBoundaryConditions; + +/// @brief Implementation of the StateAdvancer interface for advancing the solution of solid mechanics problems +class SolidMechanicsTimeIntegrator : public StateAdvancer { + public: + /** + * @brief Construct a new SolidMechanicsTimeIntegrator object. + * @param field_store Field store containing the fields. + * @param solid_weak_form Primary solid mechanics weak form. + * @param cycle_zero_weak_form Weak form for initial acceleration solve at cycle=0. + * @param solver The block solver to use. + */ + SolidMechanicsTimeIntegrator(std::shared_ptr field_store, std::shared_ptr solid_weak_form, + std::shared_ptr cycle_zero_weak_form, + std::shared_ptr solver); + + /// State enum for indexing convenience (deprecated, use FieldType instead) + enum STATE + { + DISPLACEMENT, + VELOCITY, + ACCELERATION + }; + + /** + * @brief Utility function to consistently construct all the weak forms and FieldStates for a solid mechanics + * application. + */ + template + static auto buildWeakFormAndStates([[maybe_unused]] std::shared_ptr mesh, + std::shared_ptr field_store, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, std::string physics_name, + FieldType... parameter_types) + { + // Add shape displacement + FieldType shape_disp_type(physics_name + "_shape_displacement"); + field_store->addShapeDisp(shape_disp_type); + + // Add displacement as independent (unknown) with time integration rule + auto time_rule_ptr = std::make_shared(time_rule); + FieldType disp_type(physics_name + "_displacement_predicted"); + field_store->addIndependent(disp_type, time_rule_ptr); + + // Add dependent fields for time integration history + auto disp_old_type = + field_store->addDependent(disp_type, FieldStore::TimeDerivative::VAL, physics_name + "_displacement"); + auto velo_old_type = + field_store->addDependent(disp_type, FieldStore::TimeDerivative::DOT, physics_name + "_velocity"); + auto accel_old_type = + field_store->addDependent(disp_type, FieldStore::TimeDerivative::DDOT, physics_name + "_acceleration"); + + // Add parameters + (field_store->addParameter(parameter_types), ...); + + // Create solid mechanics weak form (u, u_old, v_old, a_old) + using SolidWeakFormType = + TimeDiscretizedWeakForm>; + + auto solid_weak_form = std::make_shared( + physics_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), + field_store->createSpaces(physics_name, disp_type.name, disp_type, disp_old_type, velo_old_type, accel_old_type, + parameter_types...)); + + // Create cycle-zero weak form (u, v, a) for initial acceleration solve at cycle=0 + auto cycle_zero_weak_form = + createWeakForm(physics_name + "_reaction", accel_old_type, *field_store, disp_type, velo_old_type, + accel_old_type, parameter_types...); + + return std::make_tuple(solid_weak_form, cycle_zero_weak_form); + } + + /// @overload + std::pair, std::vector> advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states, + const std::vector& params) const override; + + private: + std::shared_ptr field_store_; + std::shared_ptr cycle_zero_weak_form_; + std::shared_ptr solver_; + std::shared_ptr integrator_; +}; + +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_statics_with_internal_vars_system.hpp b/src/smith/differentiable_numerics/solid_statics_with_internal_vars_system.hpp new file mode 100644 index 0000000000..486e969cc7 --- /dev/null +++ b/src/smith/differentiable_numerics/solid_statics_with_internal_vars_system.hpp @@ -0,0 +1,220 @@ +// 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 solid_statics_with_internal_vars_system.hpp + * @brief Defines the SolidStaticsWithInternalVarsSystem struct and its factory function + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/state_advancer.hpp" +#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/numerics/functional/tuple.hpp" +#include "smith/numerics/functional/tensor.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { + +/** + * @brief System struct for solid statics with an additional L2 state variable. + * @tparam dim Spatial dimension. + * @tparam disp_order Polynomial order for displacement field. + * @tparam StateSpace Finite element space for the state variable (e.g., L2). + * @tparam parameter_space Parameter spaces for material properties. + */ +template +struct SolidStaticsWithInternalVarsSystem : public SystemBase { + // Primary weak form: residual for displacement (u). + // Inputs: u, u_old, alpha, alpha_old, params... + /// using + using SolidWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, H1, StateSpace, StateSpace, parameter_space...>>; + + // State weak form: residual for state variable (alpha). + // Inputs: alpha, alpha_old, u, u_old, params... + /// using + using StateWeakFormType = TimeDiscretizedWeakForm< + dim, StateSpace, + Parameters, H1, parameter_space...>>; + + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + std::shared_ptr state_weak_form; ///< State variable weak form. + std::shared_ptr disp_bc; ///< Displacement boundary conditions. + std::shared_ptr state_bc; ///< State variable boundary conditions. + std::shared_ptr disp_time_rule; ///< Time integration for displacement. + std::shared_ptr state_time_rule; ///< Time integration for state. + + /** + * @brief Get the list of all state fields. + * @return std::vector List of state fields. + */ + std::vector getStateFields() const + { + return {field_store->getField(prefix("displacement_predicted")), field_store->getField(prefix("displacement")), + field_store->getField(prefix("state_predicted")), field_store->getField(prefix("state"))}; + } + + /** + * @brief Create a DifferentiablePhysics object for this system. + * @param physics_name The name of the physics. + * @return std::shared_ptr The differentiable physics object. + */ + std::shared_ptr createDifferentiablePhysics(std::string physics_name) + { + return std::make_shared( + field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), + getParameterFields(), advancer, physics_name, + std::vector{prefix("solid_residual"), prefix("state_residual")}); + } + + /** + * @brief Set the material model for the solid mechanics part. + * @tparam MaterialType The material model type. + * @param material The material model instance. + * @param domain_name The name of the domain to apply the material to. + * @note The material model should accept (state, deform_grad, alpha, params...). + */ + template + void setMaterial(const MaterialType& material, const std::string& domain_name) + { + auto captured_disp_rule = disp_time_rule; + auto captured_state_rule = state_time_rule; + + solid_weak_form->addBodyIntegral( + domain_name, [=](auto t_info, auto /*X*/, auto u, auto u_old, auto alpha, auto alpha_old, auto... params) { + // Apply time integration + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto alpha_current = captured_state_rule->value(t_info, alpha, alpha_old); + + typename MaterialType::State state; + // Material model typically needs deformation gradient and internal state variable + auto pk_stress = material(state, get(u_current), get(alpha_current), params...); + + // Return {source, flux} + // Source is body force (zero here, internal forces only) + // Flux is stress + tensor source{}; + return smith::tuple{source, pk_stress}; + }); + } + + /** + * @brief Add the evolution law for the state variable. + * @tparam EvolutionType The evolution law function type. + * @param domain_name The name of the domain. + * @param evolution_law Function with signature (t_info, alpha, alpha_dot, grad_u, params...) returning + * the residual of the ODE: alpha_dot - f(alpha, grad_u, params...). + * Time integration is applied so alpha and alpha_dot are the current predicted values. + */ + template + void addStateEvolution(const std::string& domain_name, EvolutionType evolution_law) + { + auto captured_disp_rule = disp_time_rule; + auto captured_state_rule = state_time_rule; + + state_weak_form->addBodyIntegral( + domain_name, [=](auto t_info, auto /*X*/, auto alpha, auto alpha_old, auto u, auto u_old, auto... params) { + // Apply time integration to get current state and rate + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto alpha_current = captured_state_rule->value(t_info, alpha, alpha_old); + auto alpha_dot = captured_state_rule->dot(t_info, alpha, alpha_old); + + // The evolution law is in the form: alpha_dot = f(alpha, grad_u, params...) + // Residual: alpha_dot - f(alpha, grad_u, params...) = 0 + // Pass only the scalar VALUE of alpha/alpha_dot (not the gradient) since this is a + // local pointwise ODE. Dual numbers live in the VALUE part so Jacobians are preserved. + auto residual_val = evolution_law(t_info, get(alpha_current), get(alpha_dot), + get(u_current), params...); + + // Flux is zero for a local (pointwise) ODE + tensor flux{}; + return smith::tuple{residual_val, flux}; + }); + } + + // Add other methods (body forces, etc.) as needed... +}; + +/** + * @brief Factory function to build a solid statics system with L2 state variable. + */ +template +SolidStaticsWithInternalVarsSystem buildSolidStaticsWithL2StateSystem( + std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", + FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100); + + auto prefix = [&](const std::string& name) { + if (prepend_name.empty()) { + return name; + } + return prepend_name + "_" + name; + }; + + // Add shape displacement + FieldType> shape_disp_type(prefix("shape_displacement")); + field_store->addShapeDisp(shape_disp_type); + + // 1. Displacement fields + auto disp_time_rule = std::make_shared(); + FieldType> disp_type(prefix("displacement_predicted")); + auto disp_bc = field_store->addIndependent(disp_type, disp_time_rule); + auto disp_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::VAL, prefix("displacement")); + + // 2. State variable fields + auto state_time_rule = std::make_shared(); + FieldType state_type(prefix("state_predicted")); + auto state_bc = field_store->addIndependent(state_type, state_time_rule); + auto state_old_type = field_store->addDependent(state_type, FieldStore::TimeDerivative::VAL, prefix("state")); + + // 3. Parameters + std::vector parameter_fields; + (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); + (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + + // 4. Solid Weak Form (Residual for u) + // Inputs: u, u_old, alpha, alpha_old, params... + std::string solid_res_name = prefix("solid_residual"); + auto solid_weak_form = std::make_shared< + typename SolidStaticsWithInternalVarsSystem::SolidWeakFormType>( + solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), + field_store->createSpaces(solid_res_name, disp_type.name, disp_type, disp_old_type, state_type, state_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // 5. State Weak Form (Residual for alpha) + // Inputs: alpha, alpha_old, u, u_old, params... + std::string state_res_name = prefix("state_residual"); + auto state_weak_form = std::make_shared< + typename SolidStaticsWithInternalVarsSystem::StateWeakFormType>( + state_res_name, field_store->getMesh(), field_store->getField(state_type.name).get()->space(), + field_store->createSpaces(state_res_name, state_type.name, state_type, state_old_type, disp_type, disp_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // 6. Solver and Advancer + std::vector> weak_forms{solid_weak_form, state_weak_form}; + auto advancer = std::make_shared(field_store, weak_forms, solver); + + return SolidStaticsWithInternalVarsSystem{ + {field_store, solver, advancer, parameter_fields, prepend_name}, + solid_weak_form, + state_weak_form, + disp_bc, + state_bc, + disp_time_rule, + state_time_rule}; +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/state_advancer.hpp b/src/smith/differentiable_numerics/state_advancer.hpp index 4c17147e05..2bf25138bc 100644 --- a/src/smith/differentiable_numerics/state_advancer.hpp +++ b/src/smith/differentiable_numerics/state_advancer.hpp @@ -28,18 +28,9 @@ 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 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{}; - } + virtual std::pair, std::vector> advanceState( + const TimeInfo& time_info, const FieldState& shape_disp, const std::vector& states, + const std::vector& params) const = 0; }; /// @brief creates a time info struct from gretl::State diff --git a/src/smith/differentiable_numerics/system_base.hpp b/src/smith/differentiable_numerics/system_base.hpp new file mode 100644 index 0000000000..2ad3dcde82 --- /dev/null +++ b/src/smith/differentiable_numerics/system_base.hpp @@ -0,0 +1,54 @@ +// 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 system_base.hpp + * @brief Defines the SystemBase struct for common system functionality. + */ + +#pragma once + +#include +#include +#include +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/state_advancer.hpp" + +namespace smith { + +/** + * @brief Base struct for physics systems containing common members and helper functions. + */ +struct SystemBase { + std::shared_ptr field_store; ///< Field store managing the system's fields. + std::shared_ptr solver; ///< The solver for the system. + std::shared_ptr advancer; ///< The state advancer. + std::vector parameter_fields; ///< Optional parameter fields. + std::string prepend_name_; ///< Optional prepended name for all fields. + + /** + * @brief Get the list of all parameter fields. + * @return const std::vector& List of parameter fields. + */ + const std::vector& getParameterFields() const { return parameter_fields; } + + /** + * @brief Helper function to prepend the physics name to a string. + * @param name The name to prepend to. + * @return std::string The prepended name. + */ + std::string prefix(const std::string& name) const + { + if (prepend_name_.empty()) { + return name; + } + return prepend_name_ + "_" + name; + } +}; + +} // namespace smith diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index 06681ec856..679eef93cd 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -8,9 +8,10 @@ set(differentiable_numerics_test_depends smith_physics smith_differentiable_nume set(differentiable_numerics_test_source test_field_state.cpp - test_explicit_dynamics.cpp - test_solid_mechanics_state_advancer.cpp + test_solid_dynamics.cpp test_thermo_mechanics.cpp + test_thermal_static.cpp + test_solid_static_with_internal_vars.cpp ) smith_add_tests( SOURCES ${differentiable_numerics_test_source} diff --git a/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp new file mode 100644 index 0000000000..2c71902234 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp @@ -0,0 +1,436 @@ +// 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 "gtest/gtest.h" + +#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/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/differentiable_numerics/differentiable_test_utils.hpp" +#include "smith/differentiable_numerics/solid_dynamics_system.hpp" + +namespace smith { + +LinearSolverOptions solid_linear_options{.linear_solver = LinearSolver::CG, + .preconditioner = Preconditioner::HypreJacobi, + .relative_tol = 1e-11, + .absolute_tol = 1e-11, + .max_iterations = 10000, + .print_level = 0}; + +NonlinearSolverOptions solid_nonlinear_opts{.nonlin_solver = NonlinearSolver::TrustRegion, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-10, + .max_iterations = 500, + .print_level = 1, + .force_monolithic = true}; + +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() + { + StateManager::initialize(datastore, "solid"); + auto mfem_shape = mfem::Element::HEXAHEDRON; + 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", by_attr(3)); + mesh->addDomainOfBoundaryElements("right", 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; + + auto d_solid_nonlinear_solver = + buildDifferentiableNonlinearBlockSolver(solid_nonlinear_opts, solid_linear_options, *mesh); + + auto system = buildSolidDynamicsSystem( + mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule{}, + FieldType("bulk"), FieldType("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}; + + // Set parameters + auto params = system.getParameterFields(); + params[0].get()->setFromFieldFunction([=](tensor) { return material.K0; }); + params[1].get()->setFromFieldFunction([=](tensor) { return material.G0; }); + + system.setMaterial(material, mesh->entireBodyName()); + + // Add gravity body force + system.addBodyForce(mesh->entireBodyName(), + [](double /*time*/, auto /*X*/, auto /*u*/, auto /*v*/, auto /*a*/, auto... /*params*/) { + tensor b{}; + b[1] = gravity; + return b; + }); + + // Add dummy traction to test compilation + system.addTraction("right", [](double /*time*/, auto /*X*/, auto /*n*/, auto /*u*/, auto /*v*/, auto /*a*/, + auto... /*params*/) { return tensor{}; }); + + auto shape_disp = system.field_store->getShapeDisp(); + auto states = system.getStateFields(); + + std::string pv_dir = "paraview_solid"; + auto pv_writer = createParaviewWriter(*mesh, {states[0], params[0], params[1]}, pv_dir); + pv_writer.write(0, 0.0, {states[0], params[0], params[1]}); + + double time = 0.0; + size_t cycle = 0; + std::vector reactions; + + for (size_t m = 0; m < num_steps_; ++m) { + TimeInfo t_info(time, dt_, cycle); + std::tie(states, reactions) = system.advancer->advanceState(t_info, shape_disp, states, params); + time += dt_; + cycle++; + pv_writer.write(m + 1, time, {states[0], params[0], params[1]}); + } + + double a_exact = gravity; + double v_exact = gravity * total_simulation_time_; + double u_exact = 0.5 * gravity * total_simulation_time_ * total_simulation_time_; + + TimeInfo endTimeInfo(time, dt_, cycle); + + // Test acceleration (states[3] is acceleration) + FunctionalObjective> accel_error("accel_error", mesh, spaces({states[3]})); + 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(endTimeInfo, shape_disp.get().get(), getConstFieldPointers({states[3]})); + EXPECT_NEAR(0.0, a_err, 1e-14); + + // Test velocity (states[2] is velocity) + FunctionalObjective> velo_error("velo_error", mesh, spaces({states[2]})); + 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({states[2]})); + EXPECT_NEAR(0.0, v_err, 1e-14); + + // Test displacement (states[1] is displacement) + FunctionalObjective> disp_error("disp_error", mesh, spaces({states[1]})); + 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({states[0]})); + EXPECT_NEAR(0.0, u_err, 1e-14); +} + +auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) +{ + std::shared_ptr d_solid_nonlinear_solver = + buildDifferentiableNonlinearBlockSolver(solid_nonlinear_opts, solid_linear_options, *mesh); + + auto time_rule = ImplicitNewmarkSecondOrderTimeIntegrationRule(); + + auto system = buildSolidDynamicsSystem(mesh, d_solid_nonlinear_solver, time_rule, physics_name, + FieldType("bulk"), + FieldType("shear")); + + auto physics = system.createDifferentiablePhysics(physics_name); + auto bcs = system.disp_bc; + + bcs->setFixedVectorBCs(mesh->domain("right")); + bcs->setVectorBCs(mesh->domain("left"), [](double t, 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}; + + system.setMaterial(material, mesh->entireBodyName()); + + auto shape_disp = physics->getShapeDispFieldState(); + auto params = physics->getFieldParams(); + auto states = physics->getInitialFieldStates(); + + params[0].get()->setFromFieldFunction([=](tensor) { + double scaling = 1.0; + return scaling * material.K0; + }); + + params[1].get()->setFromFieldFunction([=](tensor) { + double scaling = 1.0; + return scaling * material.G0; + }); + + physics->resetStates(); + + return std::make_tuple(physics, shape_disp, states, params, bcs); +} + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) +{ + SMITH_MARK_FUNCTION; + std::string physics_name = "solid"; + auto [physics, shape_disp, initial_states, params, bcs] = createSolidMechanicsBasePhysics(physics_name, mesh); + + 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()); + } + + auto reactions = physics->getReactionStates(); + + // Check that reaction forces are zero away from Dirichlet DOFs + checkUnconstrainedReactions(*reactions[0].get(), bcs->getBoundaryConditionManager()); + + auto reaction_squared = 0.5 * innerProduct(reactions[0], reactions[0]); + + gretl::set_as_objective(reaction_squared); + + 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) { + SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", s.get()->name(), s.get()->Norml2(), s.get_dual()->Norml2())); + } + + SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", shape_disp.get()->name(), shape_disp.get()->Norml2(), + shape_disp.get_dual()->Norml2())); + + for (size_t p = 0; p < params.size(); ++p) { + SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", params[p].get()->name(), params[p].get()->Norml2(), + params[p].get_dual()->Norml2())); + } +} + +// 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, std::string reaction_name) +{ + resetAndApplyInitialConditions(physics); + for (size_t m = 0; m < num_steps; ++m) { + physics->advanceTimestep(dt); + } + FiniteElementDual reaction = physics->dual(reaction_name); + + return 0.5 * innerProduct(reaction, reaction); +} + +void adjointBackward(std::shared_ptr physics, smith::FiniteElementDual& shape_sensitivity, + std::vector& parameter_sensitivities, std::string reaction_name) +{ + smith::FiniteElementDual reaction = physics->dual(reaction_name); + smith::FiniteElementState reaction_dual(reaction.space(), reaction_name + "_dual"); + reaction_dual = reaction; + + physics->resetAdjointStates(); + + physics->setDualAdjointBcs({{reaction_name, 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"; + auto [physics, shape_disp, initial_states, params, bcs] = createSolidMechanicsBasePhysics(physics_name, mesh); + + double qoi = integrateForward(physics, num_steps_, dt_, physics_name + "_reactions"); + SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); + + // Check that reaction forces are zero away from Dirichlet DOFs + auto reactions = physics->getReactionStates(); + checkUnconstrainedReactions(*reactions[0].get(), bcs->getBoundaryConditionManager()); + + 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, physics_name + "_reactions"); + + auto state_sensitivities = physics->computeInitialConditionSensitivity(); + for (auto name_and_state_sensitivity : state_sensitivities) { + SLIC_INFO_ROOT( + axom::fmt::format("{} {}", name_and_state_sensitivity.first, name_and_state_sensitivity.second.Norml2())); + } + + SLIC_INFO_ROOT(axom::fmt::format("{} {}", shape_sensitivity.name(), shape_sensitivity.Norml2())); + + for (size_t p = 0; p < num_params; ++p) { + SLIC_INFO_ROOT(axom::fmt::format("{} {}", parameter_sensitivities[p].name(), parameter_sensitivities[p].Norml2())); + } +} + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) +{ + SMITH_MARK_FUNCTION; + std::string physics_name = "solid"; + + // 1. Calculate sensitivities using Gretl + auto [physicsGretl, shape_dispG, initial_statesG, paramsG, bcsG] = + createSolidMechanicsBasePhysics(physics_name + "_gretl", mesh); + + // Forward pass + for (size_t m = 0; m < num_steps_; ++m) { + physicsGretl->advanceTimestep(dt_); + } + + auto reactionsG = physicsGretl->getReactionStates(); + auto reaction_squaredG = 0.5 * innerProduct(reactionsG[0], reactionsG[0]); + + // Backprop + gretl::set_as_objective(reaction_squaredG); + reaction_squaredG.data_store().back_prop(); + + // 2. Calculate sensitivities using BasePhysics manual adjoint + auto [physicsBase, shape_dispB, initial_statesB, paramsB, bcsB] = + createSolidMechanicsBasePhysics(physics_name + "_base", mesh); + + // Forward pass + double qoiB = integrateForward(physicsBase, num_steps_, dt_, physics_name + "_base_reactions"); + + // Adjoint pass + size_t num_params = physicsBase->parameterNames().size(); + smith::FiniteElementDual shape_sensitivityB(*shape_dispB.get_dual()); + shape_sensitivityB = 0.0; + std::vector parameter_sensitivitiesB; + for (size_t p = 0; p < num_params; ++p) { + parameter_sensitivitiesB.emplace_back(*paramsB[p].get_dual()); + parameter_sensitivitiesB.back() = 0.0; + } + + adjointBackward(physicsBase, shape_sensitivityB, parameter_sensitivitiesB, physics_name + "_base_reactions"); + auto initial_condition_sensitivitiesB = physicsBase->computeInitialConditionSensitivity(); + + // 3. Compare sensitivities + double tol = 1e-12; + + // Compare objective values + EXPECT_NEAR(reaction_squaredG.get(), qoiB, tol); + + auto diff_norm = [](const mfem::Vector& a, const mfem::Vector& b) { + mfem::Vector diff = a; + diff -= b; + return diff.Norml2(); + }; + + // Compare shape sensitivities + double shape_diff = diff_norm(*shape_dispG.get_dual(), shape_sensitivityB); + SLIC_INFO_ROOT(axom::fmt::format("Shape sensitivity difference: {:.6e}", shape_diff)); + EXPECT_LT(shape_diff, tol); + + // Compare parameter sensitivities + for (size_t p = 0; p < num_params; ++p) { + double param_diff = diff_norm(*paramsG[p].get_dual(), parameter_sensitivitiesB[p]); + SLIC_INFO_ROOT(axom::fmt::format("Parameter {} sensitivity difference: {:.6e}", p, param_diff)); + EXPECT_LT(param_diff, tol); + } + + // Compare initial condition sensitivities + std::vector state_suffixes = {"displacement_predicted", "displacement", "velocity", "acceleration"}; + for (const auto& suffix : state_suffixes) { + std::string nameG = physics_name + "_gretl_" + suffix; + std::string nameB = physics_name + "_base_" + suffix; + SLIC_INFO_ROOT(axom::fmt::format("Comparing sensitivity for {}: {} vs {}", suffix, nameG, nameB)); + + // Find Gretl dual + const FiniteElementDual* dualG = nullptr; + for (auto const& s : initial_statesG) { + if (s.get()->name() == nameG) { + dualG = s.get_dual().get(); + break; + } + } + ASSERT_NE(dualG, nullptr) << "Could not find Gretl dual for " << nameG; + + double state_diff = diff_norm(*dualG, initial_condition_sensitivitiesB.at(nameB)); + SLIC_INFO_ROOT(axom::fmt::format("Initial state {} sensitivity difference: {:.6e}", suffix, state_diff)); + EXPECT_LT(state_diff, tol); + } +} + +} // 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/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp deleted file mode 100644 index e22e123fbb..0000000000 --- a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp +++ /dev/null @@ -1,355 +0,0 @@ -// 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 "gtest/gtest.h" - -#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/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 = - buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - - smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule; - - 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 endTimeInfo(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(endTimeInfo, 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); -} - -auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) -{ - size_t num_checkpoints = 200; - std::shared_ptr d_solid_nonlinear_solver = - buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - - auto [physics, solid_weak_form, bcs] = - buildSolidMechanics( - mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule(), num_checkpoints, - 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(); - - return std::make_tuple(physics, shape_disp, states, params); -} - -TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) -{ - SMITH_MARK_FUNCTION; - std::string physics_name = "solid"; - auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); - - 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); - - 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) { - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", s.get()->name(), s.get()->Norml2(), s.get_dual()->Norml2())); - } - - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", shape_disp.get()->name(), shape_disp.get()->Norml2(), - shape_disp.get_dual()->Norml2())); - - for (size_t p = 0; p < params.size(); ++p) { - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", params[p].get()->name(), params[p].get()->Norml2(), - params[p].get_dual()->Norml2())); - } -} - -// 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"; - auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); - - double qoi = integrateForward(physics, num_steps_, dt_); - SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); - - 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) { - SLIC_INFO_ROOT( - axom::fmt::format("{} {}", name_and_state_sensitivity.first, name_and_state_sensitivity.second.Norml2())); - } - - SLIC_INFO_ROOT(axom::fmt::format("{} {}", shape_sensitivity.name(), shape_sensitivity.Norml2())); - - for (size_t p = 0; p < num_params; ++p) { - SLIC_INFO_ROOT(axom::fmt::format("{} {}", parameter_sensitivities[p].name(), parameter_sensitivities[p].Norml2())); - } -} - -} // 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/tests/test_solid_static_with_internal_vars.cpp b/src/smith/differentiable_numerics/tests/test_solid_static_with_internal_vars.cpp new file mode 100644 index 0000000000..107e20ffac --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_solid_static_with_internal_vars.cpp @@ -0,0 +1,146 @@ +// 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 "gtest/gtest.h" +#include "smith/smith_config.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/differentiable_numerics/solid_statics_with_internal_vars_system.hpp" +#include "smith/differentiable_numerics/differentiable_test_utils.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" + +namespace smith { + +LinearSolverOptions solid_linear_options{.linear_solver = LinearSolver::Strumpack, + .preconditioner = Preconditioner::HypreJacobi, + .relative_tol = 1e-12, + .absolute_tol = 1e-12, + .max_iterations = 2000, + .print_level = 0}; + +NonlinearSolverOptions solid_nonlinear_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-10, + .max_iterations = 100, + .max_line_search_iterations = 50, + .print_level = 1, + .force_monolithic = true}; + +static constexpr int dim = 3; +static constexpr int disp_order = 1; +static constexpr int state_order = 0; + +using StateSpace = L2; + +struct SolidStaticWithInternalVarsFixture : public testing::Test { + void SetUp() override + { + StateManager::initialize(datastore, "solid_static_with_internal_vars"); + mesh = std::make_shared(mfem::Mesh::MakeCartesian3D(4, 4, 4, mfem::Element::HEXAHEDRON, 1.0, 1.0, 1.0), + "mesh", 0, 0); + mesh->addDomainOfBoundaryElements("bottom", by_attr(1)); // z=0 + mesh->addDomainOfBoundaryElements("top", by_attr(6)); // z=1 + } + + axom::sidre::DataStore datastore; + std::shared_ptr mesh; +}; + +// Simple damage-like material: stiffness decreases with internal state variable (isv) +struct DamageMaterial { + using State = smith::QOI; + + double E = 100.0; + double nu = 0.3; + + template + SMITH_HOST_DEVICE auto operator()(StateType /*state*/, DerivType deriv_u, ISVType isv, Params... /*params*/) const + { + auto epsilon = sym(deriv_u); + auto tr_eps = tr(epsilon); + auto I = Identity(); + + double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + double mu = E / (2.0 * (1.0 + nu)); + + // Stiffness degradation: (1 - isv) * Stress_elastic + auto damage = isv; + if (damage > 1.0) damage = 1.0; + if (damage < 0.0) damage = 0.0; + + auto factor = 1.0 - damage; + auto sigma = factor * (lambda * tr_eps * I + 2.0 * mu * epsilon); + + return sigma; + } +}; + +// Evolution law: isv_dot = (1 - isv) * eps_norm +// +// This ODE drives isv toward 1 as eps_norm increases. +// For constant eps_norm, the solution is: isv(t) = 1 - (1 - isv_0) * exp(-eps_norm * t) +// which asymptotes to 1 as eps_norm -> infinity (fast convergence) or t -> infinity. +// +// Residual form returned: isv_dot - (1 - isv) * eps_norm = 0 +struct StrainNormEvolution { + template + SMITH_HOST_DEVICE auto operator()(TimeInfo /*t_info*/, ISVType isv, ISVDotType isv_dot, DerivType deriv_u, + Params... /*params*/) const + { + using std::sqrt; + auto epsilon = sym(deriv_u); + auto eps_norm = sqrt(inner(epsilon, epsilon)); + + // ODE: isv_dot = (1 - isv) * eps_norm + return isv_dot - (1.0 - isv) * eps_norm; + } +}; + +TEST_F(SolidStaticWithInternalVarsFixture, CoupledSolve) +{ + auto solver = buildDifferentiableNonlinearBlockSolver(solid_nonlinear_opts, solid_linear_options, *mesh); + + auto system = + buildSolidStaticsWithL2StateSystem(mesh, solver, "solid_static_with_internal_vars"); + + // Material and Evolution + system.setMaterial(DamageMaterial{}, mesh->entireBodyName()); + system.addStateEvolution(mesh->entireBodyName(), StrainNormEvolution{}); + + // Boundary Conditions + + // Fix bottom face + system.disp_bc->setFixedVectorBCs(mesh->domain("bottom")); + + // Pull top face + double pull_rate = 0.05; + system.disp_bc->setVectorBCs(mesh->domain("top"), [pull_rate](double t, tensor /*X*/) { + tensor u{}; + u[2] = pull_rate * t; + return u; + }); + + auto physics = system.createDifferentiablePhysics("physics"); + + // Create ParaView writer + auto writer = createParaviewWriter(*mesh, physics->getFieldStates(), "solid_state_output"); + writer.write(0, 0.0, physics->getFieldStates()); + // Advance multiple steps + for (int step = 1; step <= 5; ++step) { + physics->advanceTimestep(1.0); + writer.write(step, step * 1.0, physics->getFieldStates()); + SLIC_INFO("Completed step " << step); + } +} + +} // 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/tests/test_thermal_static.cpp b/src/smith/differentiable_numerics/tests/test_thermal_static.cpp new file mode 100644 index 0000000000..fa59f8550b --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_thermal_static.cpp @@ -0,0 +1,137 @@ +// 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 +#include "smith/differentiable_numerics/thermal_system.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/physics/mesh.hpp" +#include "smith/physics/common.hpp" +#include "smith/physics/state/state_manager.hpp" +#include "smith/numerics/functional/functional.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "axom/slic.hpp" +#include "axom/fmt.hpp" +#include "axom/sidre.hpp" + +using namespace smith; + +struct ThermalStaticFixture : public testing::Test { + axom::sidre::DataStore datastore; + std::shared_ptr mesh; + + void setupMesh(int n_elements, std::string tag) + { + StateManager::reset(); + StateManager::initialize(datastore, "thermal_static_" + tag); + auto mfem_mesh = mfem::Mesh::MakeCartesian2D(n_elements, n_elements, mfem::Element::QUADRILATERAL, true, 1.0, 1.0); + mesh = std::make_shared(std::move(mfem_mesh), "thermal_mesh_" + tag); + } + + void SetUp() override + { + // Default mesh for single-mesh tests + setupMesh(32, "default"); + } + + void TearDown() override { StateManager::reset(); } + + template + double run_thermal_solve() + { + auto solver_options = NonlinearSolverOptions(); + solver_options.relative_tol = 1e-12; + solver_options.force_monolithic = true; + auto linear_options = LinearSolverOptions(); + auto block_solver = buildDifferentiableNonlinearBlockSolver(solver_options, linear_options, *mesh); + + auto thermal_system = buildThermalSystem<2, temp_order>(mesh, block_solver); + + double k = 1.0; + thermal_system.setThermalIntegrand("entire_body", [=](auto /*t_info*/, auto /*X*/, auto T) { + auto gradT = smith::get(T); + return smith::tuple{smith::zero{}, k * gradT}; + }); + + thermal_system.addBodyHeatSource("entire_body", [=](auto /*t*/, auto X, auto /*T*/) { + auto x = X[0]; + auto y = X[1]; + double pi = 3.14159265358979323846; + return 2.0 * k * pi * pi * sin(pi * x) * sin(pi * y); + }); + + thermal_system.temperature_bc->template setScalarBCs<2>(mesh->entireBoundary(), + [](double /*t*/, tensor /*X*/) { return 0.0; }); + + TimeInfo t_info(0.0, 1.0); + auto [new_states, reactions] = thermal_system.advancer->advanceState( + t_info, thermal_system.field_store->getShapeDisp(), thermal_system.field_store->getAllFields(), + thermal_system.getParameterFields()); + + for (size_t i = 0; i < new_states.size(); ++i) { + thermal_system.field_store->setField(i, new_states[i]); + } + + auto temperature = thermal_system.field_store->getField(thermal_system.prefix("temperature")); + + auto exact_sol_func = [](const mfem::Vector& X, mfem::Vector& T) { + double x = X(0); + double y = X(1); + double pi = 3.14159265358979323846; + T(0) = std::sin(pi * x) * std::sin(pi * y); + }; + mfem::VectorFunctionCoefficient exact_sol_coeff(1, exact_sol_func); + + return computeL2Error(*temperature.get(), exact_sol_coeff); + } +}; + +TEST_F(ThermalStaticFixture, ManufacturedSolutionOrder1) +{ + double error = run_thermal_solve<1>(); + SLIC_INFO("L2 Error (Order 1, h=1/32): " << error); + EXPECT_LT(error, 1e-3); +} + +TEST_F(ThermalStaticFixture, ManufacturedSolutionOrder2) +{ + double error = run_thermal_solve<2>(); + SLIC_INFO("L2 Error (Order 2, h=1/32): " << error); + EXPECT_LT(error, 1e-5); +} + +TEST_F(ThermalStaticFixture, ConvergenceOrder1) +{ + setupMesh(16, "conv1_h16"); + double e1 = run_thermal_solve<1>(); + setupMesh(32, "conv1_h32"); + double e2 = run_thermal_solve<1>(); + + double rate = std::log2(e1 / e2); + SLIC_INFO("Convergence Rate (Order 1): " << rate); + EXPECT_NEAR(rate, 2.0, 0.1); +} + +TEST_F(ThermalStaticFixture, ConvergenceOrder2) +{ + setupMesh(8, "conv2_h8"); + double e1 = run_thermal_solve<2>(); + setupMesh(16, "conv2_h16"); + double e2 = run_thermal_solve<2>(); + + double rate = std::log2(e1 / e2); + SLIC_INFO("Convergence Rate (Order 2): " << rate); + EXPECT_NEAR(rate, 3.0, 0.1); +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + ApplicationManager app(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index f5fa7d773a..8cc61ba2ed 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -9,18 +9,17 @@ #include "smith/smith_config.hpp" #include "smith/infrastructure/application_manager.hpp" #include "smith/numerics/solver_config.hpp" - #include "smith/physics/state/state_manager.hpp" +#include "smith/physics/mesh.hpp" -#include "smith/differentiable_numerics/field_state.hpp" #include "smith/differentiable_numerics/differentiable_solver.hpp" -#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" -#include "smith/differentiable_numerics/time_integration_rule.hpp" -#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/thermo_mechanics_system.hpp" #include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" #include "smith/differentiable_numerics/differentiable_test_utils.hpp" #include "smith/differentiable_numerics/nonlinear_solve.hpp" -#include "gretl/strumm_walther_checkpoint_strategy.hpp" +#include "smith/physics/functional_objective.hpp" +#include "gretl/wang_checkpoint_strategy.hpp" namespace smith { @@ -39,7 +38,7 @@ auto greenStrain(const tensor& grad_u) /// autodifferentiation struct GreenSaintVenantThermoelasticMaterial { double density; ///< density - double E; ///< Young's modulus + double E0; ///< base Young's modulus double nu; ///< Poisson's ratio double C_v; ///< volumetric heat capacity double alpha; ///< thermal expansion coefficient @@ -66,12 +65,13 @@ struct GreenSaintVenantThermoelasticMaterial { * step (units of energy), and the referential heat flux (units of * energy per unit time and per unit area). */ - template + template auto operator()(double, State&, const tensor& grad_u, const tensor& grad_v, T3 theta, - const tensor& grad_theta) const + const tensor& grad_theta, const T5& E_param) const { - const double K = E / (3.0 * (1.0 - 2.0 * nu)); - const double G = 0.5 * E / (1.0 + nu); + auto E = E0 + get<0>(E_param); + const auto K = E / (3.0 * (1.0 - 2.0 * nu)); + const auto G = 0.5 * E / (1.0 + nu); const auto Eg = greenStrain(grad_u); const auto trEg = tr(Eg); @@ -84,13 +84,15 @@ struct GreenSaintVenantThermoelasticMaterial { // internal heat power auto greenStrainRate = 0.5 * (grad_v + transpose(grad_v) + dot(transpose(grad_v), grad_u) + dot(transpose(grad_u), grad_v)); - const auto s0 = -dim * K * alpha * (theta + 273.1) * tr(greenStrainRate); + const auto s0 = -dim * K * alpha * (theta + 273.1) * tr(greenStrainRate) + 0.0 * E; // heat flux const auto q0 = -kappa * grad_theta; return smith::tuple{Piola, C_v, s0, q0}; } + + static constexpr int numParameters() { return 1; } }; smith::LinearSolverOptions linear_options{.linear_solver = smith::LinearSolver::Strumpack, @@ -100,16 +102,18 @@ smith::LinearSolverOptions linear_options{.linear_solver = smith::LinearSolver:: .print_level = 0}; smith::NonlinearSolverOptions nonlinear_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, - .relative_tol = 1.9e-6, + .relative_tol = 1.0e-10, .absolute_tol = 1.0e-10, .max_iterations = 500, .max_line_search_iterations = 50, .print_level = 2}; static constexpr int dim = 3; -static constexpr int order = 1; +static constexpr int vdim = dim; +static constexpr int displacement_order = 1; +static constexpr int temperature_order = 1; -struct SolidMechanicsMeshFixture : public testing::Test { +struct ThermoMechanicsMeshFixture : public testing::Test { double length = 1.0; double width = 0.04; int num_elements_x = 12; @@ -128,304 +132,311 @@ struct SolidMechanicsMeshFixture : public testing::Test { 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_; }; -template -struct FieldType { - FieldType(std::string n, int unknown_index_ = -1) : name(n), unknown_index(unknown_index_) {} - std::string name; - int unknown_index; -}; - -struct FieldStore { - FieldStore(std::shared_ptr mesh, size_t storage_size = 50) - : mesh_(mesh), - data_store_( - std::make_shared(std::make_unique(storage_size))) - { - } - - template - void addShapeDisp(FieldType type) - { - shape_disp_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); - } - - template - std::shared_ptr& addUnknown(FieldType& type) - { - type.unknown_index = static_cast(num_unknowns_); - to_fields_index_[type.name] = fields_.size(); - to_unknown_index_[type.name] = num_unknowns_; - FieldState new_field = smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag()); - fields_.push_back(new_field); - ++num_unknowns_; - boundary_conditions_.push_back( - std::make_shared(mesh_->mfemParMesh(), new_field.get()->space())); - SLIC_ERROR_IF(num_unknowns_ != boundary_conditions_.size(), - "Inconsistency between num unknowns and boundary condition size"); - return boundary_conditions_.back(); - } - - template - auto addDerived(FieldType, std::string name) - { - to_fields_index_[name] = fields_.size(); - fields_.push_back(smith::createFieldState(*data_store_, Space{}, name, mesh_->tag())); - return FieldType(name); - } - - void addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index) - { - FieldLabel argument_name_and_index{.field_name = argument_name, .field_index_in_residual = argument_index}; - if (weak_form_name_to_unknown_name_index_.count(weak_form_name)) { - weak_form_name_to_unknown_name_index_.at(weak_form_name).push_back(argument_name_and_index); - } else { - weak_form_name_to_unknown_name_index_[weak_form_name] = std::vector{argument_name_and_index}; - } - } +TEST_F(ThermoMechanicsMeshFixture, RunThermoMechanicalCoupled) +{ + SMITH_MARK_FUNCTION; + double rho = 1.0; + double E0 = 100.0; + double nu = 0.25; + double specific_heat = 1.0; + double kappa = 0.1; + GreenSaintVenantThermoelasticMaterial material{rho, E0, nu, specific_heat, 0.0, 1.0, kappa}; - void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) - { - size_t field_index = to_fields_index_.at(argument_name); - if (weak_form_name_to_field_indices_.count(weak_form_name)) { - weak_form_name_to_field_indices_.at(weak_form_name).push_back(field_index); - } else { - weak_form_name_to_field_indices_[weak_form_name] = std::vector{field_index}; - } - SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_indices_.at(weak_form_name).size(), - "Invalid order for adding weak form arguments."); - } + auto solver = buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); - void printMap() - { - for (auto& keyval : weak_form_name_to_unknown_name_index_) { - std::cout << "for residual: " << keyval.first << " "; - for (auto& name_index : keyval.second) { - std::cout << "arg " << name_index.field_name << " at " << name_index.field_index_in_residual << ", "; - } - std::cout << std::endl; - } - } + FieldType> youngs_modulus("youngs_modulus"); + auto system = buildThermoMechanicsSystem(mesh_, solver, youngs_modulus); + system.setMaterial(material, mesh_->entireBodyName()); - std::vector> indexMap(const std::vector& residual_names) const - { - std::vector> block_indices(residual_names.size()); - - for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { - std::vector& res_indices = block_indices[res_i]; - res_indices = std::vector(num_unknowns_, invalid_block_index); - const std::string& res_name = residual_names[res_i]; - const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); - - for (const auto& field_name_and_arg_index : arg_info) { - const std::string field_name = field_name_and_arg_index.field_name; - size_t unknown_index = to_unknown_index_.at(field_name); - SLIC_ASSERT(unknown_index < num_unknowns_); - res_indices[unknown_index] = field_name_and_arg_index.field_index_in_residual; - } - } - - return block_indices; - } + system.parameter_fields[0].get()->setFromFieldFunction([=](smith::tensor) { return E0; }); - std::vector getBoundaryConditionManagers() const - { - std::vector bcs; - for (auto& bc : boundary_conditions_) { - bcs.push_back(&bc->getBoundaryConditionManager()); - } - return bcs; + system.disp_bc->setVectorBCs(mesh_->domain("left"), [](double t, smith::tensor X) { + auto bc = 0.0 * X; + bc[0] = 0.01 * t; + return bc; + }); + system.disp_bc->setFixedVectorBCs(mesh_->domain("right")); + system.temperature_bc->setFixedScalarBCs(mesh_->domain("left")); + system.temperature_bc->setFixedScalarBCs(mesh_->domain("right")); + + system.addThermalHeatSource(mesh_->entireBodyName(), [](auto /*t*/, auto /*x*/, auto /*u*/, auto /*v*/, auto /*T*/, + auto /*T_dot*/, auto /*E_param*/) { return 100.0; }); + + std::string pv_dir = "paraview_thermo_mechanics"; + auto pv_writer = smith::createParaviewWriter(*mesh_, system.getStateFields(), pv_dir); + + double dt = 0.001; + double time = 0.0; + int cycle = 0; + + auto shape_disp = system.field_store->getShapeDisp(); + auto states = system.getStateFields(); + auto params = system.getParameterFields(); + std::vector reactions; + + pv_writer.write(cycle, time, states); + for (size_t step = 0; step < 2; ++step) { + TimeInfo t_info(time, dt, step); + std::tie(states, reactions) = system.advancer->advanceState(t_info, shape_disp, states, params); + time += dt; + cycle++; + pv_writer.write(cycle, time, states); } - size_t getFieldIndex(const std::string& field_name) const { return to_fields_index_.at(field_name); } - - const FieldState& getField(const std::string& field_name) const - { - size_t field_index = getFieldIndex(field_name); - return fields_[field_index]; - } + // Check that reactions are zero for unconstrained DOFs (should be within solver tolerance) + checkUnconstrainedReactions(*reactions[0].get(), system.disp_bc->getBoundaryConditionManager()); + checkUnconstrainedReactions(*reactions[1].get(), system.temperature_bc->getBoundaryConditionManager()); - void setField(const std::string& field_name, FieldState updated_field) - { - size_t field_index = getFieldIndex(field_name); - fields_[field_index] = updated_field; - } + auto reaction_squared = innerProduct(reactions[0], reactions[0]); + gretl::set_as_objective(reaction_squared); - const FieldState& getShapeDisp() const { return shape_disp_[0]; } + 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); +} - const std::vector& getAllFields() const { return fields_; } +TEST_F(ThermoMechanicsMeshFixture, TransientHeatEquationAnalytic) +{ + SMITH_MARK_FUNCTION; - std::vector getFields(const std::string& weak_form_name) const - { - auto unknown_field_indices = weak_form_name_to_field_indices_.at(weak_form_name); - std::vector fields_for_residual; - for (auto& i : unknown_field_indices) { - fields_for_residual.push_back(fields_[i]); - } - return fields_for_residual; + double rho = 1.0; + double E0 = 100.0; + double nu = 0.25; + double specific_heat = 1.0; + double kappa = 0.1; + GreenSaintVenantThermoelasticMaterial material{rho, E0, nu, specific_heat, 0.0, 0.0, kappa}; + + auto solver = buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + FieldType> youngs_modulus("youngs_modulus"); + auto system = buildThermoMechanicsSystem(mesh_, solver, youngs_modulus); + system.setMaterial(material, mesh_->entireBodyName()); + + system.parameter_fields[0].get()->setFromFieldFunction([=](smith::tensor) { return E0; }); + + system.disp_bc->setFixedVectorBCs(mesh_->domain("left")); + system.disp_bc->setFixedVectorBCs(mesh_->domain("right")); + system.temperature_bc->setScalarBCs(mesh_->domain("left"), [](double, auto) { return 100.0; }); + system.temperature_bc->setScalarBCs(mesh_->domain("right"), [](double, auto) { return 100.0; }); + + auto& temp_field = const_cast(*system.field_store->getField("temperature_predicted").get()); + temp_field.setFromFieldFunction([](tensor x) { + using std::sin; + return 100.0 + sin(M_PI * x[0]); + }); + const_cast(*system.field_store->getField("temperature").get()) = temp_field; + + double dt = 0.01; + double time = 0.0; + auto shape_disp = system.field_store->getShapeDisp(); + auto states = system.getStateFields(); + auto params = system.getParameterFields(); + std::vector reactions; + + double diffusivity = kappa / (rho * specific_heat); + + size_t num_steps = 2; + for (size_t step = 0; step < num_steps; ++step) { + TimeInfo t_info(time, dt, step); + std::tie(states, reactions) = system.advancer->advanceState(t_info, shape_disp, states, params); + time += dt; } - const std::shared_ptr& getMesh() const { return mesh_; } - - private: - std::shared_ptr mesh_; - std::shared_ptr data_store_; + // Check that reactions are zero for unconstrained DOFs + checkUnconstrainedReactions(*reactions[0].get(), system.disp_bc->getBoundaryConditionManager()); + checkUnconstrainedReactions(*reactions[1].get(), system.temperature_bc->getBoundaryConditionManager()); - std::vector shape_disp_; - std::vector fields_; - std::map to_fields_index_; + // Check nodal error against exact solution: T(x,t) = 100 + sin(pi*x) * exp(-diffusivity * pi^2 * t) + auto temp_idx = system.field_store->getFieldIndex("temperature_predicted"); + FieldState final_temp = states[temp_idx]; - size_t num_unknowns_ = 0; - std::map to_unknown_index_; - std::vector> boundary_conditions_; - - struct FieldLabel { - std::string field_name; - size_t field_index_in_residual; - }; - - std::map> weak_form_name_to_unknown_name_index_; - - std::map> weak_form_name_to_field_indices_; -}; + FiniteElementState exact_temp(final_temp.get()->space(), "exact_temp"); + exact_temp.setFromFieldFunction([&](tensor x) { + return 100.0 + std::sin(M_PI * x[0]) * std::exp(-diffusivity * M_PI * M_PI * time); + }); + mfem::Vector diff(final_temp.get()->Size()); + subtract(*final_temp.get(), exact_temp, diff); + double max_nodal_error = diff.Normlinf(); -template -void createSpaces(const std::string& weak_form_name, FieldStore& field_store, - std::vector& spaces, size_t arg_num, FirstType type, - Types... types) -{ - SLIC_ERROR_IF(spaces.size() != arg_num, "Error creating spaces recursively"); - spaces.push_back(&field_store.getField(type.name).get()->space()); - field_store.addWeakFormArg(weak_form_name, type.name, arg_num); - if (type.unknown_index >= 0) { - field_store.addWeakFormUnknownArg(weak_form_name, type.name, arg_num); - } - if constexpr (sizeof...(types) > 0) { - createSpaces(weak_form_name, field_store, spaces, arg_num + 1, types...); - } + std::cout << "Transient Heat max nodal error: " << max_nodal_error << std::endl; + EXPECT_LT(max_nodal_error, 1e-4); } -template -auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, - FieldType... field_types) +TEST_F(ThermoMechanicsMeshFixture, StaticElasticityAnalytic) { - const mfem::ParFiniteElementSpace& test_space = field_store.getField(test_type.name).get()->space(); - std::vector input_spaces; - createSpaces(name, field_store, input_spaces, 0, field_types...); - return std::make_shared>>( - name, field_store.getMesh(), test_space, input_spaces); -} + SMITH_MARK_FUNCTION; -std::vector solve(const std::vector& weak_forms, const FieldStore& field_store, - const DifferentiableBlockSolver* solver, const TimeInfo& time_info) -{ - std::vector weak_form_names; - for (const auto& wf : weak_forms) { - weak_form_names.push_back(wf->name()); - } - std::vector> index_map = field_store.indexMap(weak_form_names); + double rho = 1.0; + double E0 = 100.0; + double nu = 0.25; + double specific_heat = 1.0; + double kappa = 0.1; + GreenSaintVenantThermoelasticMaterial material{rho, E0, nu, specific_heat, 0.0, 0.0, kappa}; + + auto solver = buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + FieldType> youngs_modulus("youngs_modulus"); + auto system = buildThermoMechanicsSystem(mesh_, solver, youngs_modulus); + system.setMaterial(material, mesh_->entireBodyName()); + + system.parameter_fields[0].get()->setFromFieldFunction([=](smith::tensor) { return E0; }); + + // Arbitrary affine displacement: u(X) = G * X, where G is a constant displacement gradient + // Choose a small uniform deformation with both normal and shear components + tensor G; + G[0][0] = 0.02; + G[0][1] = 0.01; + G[0][2] = 0.005; + G[1][0] = 0.01; + G[1][1] = 0.03; + G[1][2] = 0.002; + G[2][0] = 0.005; + G[2][1] = 0.002; + G[2][2] = 0.015; + + auto u_exact_func = [G](auto X) { return dot(G, X); }; + + system.disp_bc->setVectorBCs(mesh_->entireBoundary(), + [=](double, tensor X) { return u_exact_func(X); }); + system.temperature_bc->setFixedScalarBCs(mesh_->entireBoundary()); + + double dt = 1.0; + double time = 0.0; + auto shape_disp = system.field_store->getShapeDisp(); + auto states = system.getStateFields(); + auto params = system.getParameterFields(); + std::vector reactions; + + // Run 1 step + TimeInfo t_info(time, dt, 0); + auto states_and_reactions = system.advancer->advanceState(t_info, shape_disp, states, params); + states = states_and_reactions.first; + reactions = states_and_reactions.second; + + // Check that reactions are zero for unconstrained DOFs + checkUnconstrainedReactions(*reactions[0].get(), system.disp_bc->getBoundaryConditionManager()); + checkUnconstrainedReactions(*reactions[1].get(), system.temperature_bc->getBoundaryConditionManager()); + + // Check error - for affine displacement, FEM solution should be exact (up to machine precision) + auto disp_idx = system.field_store->getFieldIndex("displacement_predicted"); + FieldState final_disp = states[disp_idx]; + + FunctionalObjective>> error_obj("error", mesh_, spaces({final_disp})); + error_obj.addBodyIntegral(DependsOn<0>{}, mesh_->entireBodyName(), [=](auto /*t*/, auto X, auto U) { + auto u_val = get(U); + auto x_val = get<0>(X); + auto exact = u_exact_func(x_val); + return inner(u_val - exact, u_val - exact); + }); - std::vector> inputs; - for (size_t i = 0; i < weak_forms.size(); ++i) { - std::string wf_name = weak_forms[i]->name(); - std::vector fields_for_wk = field_store.getFields(wf_name); - inputs.push_back(fields_for_wk); - } - std::vector> params(weak_forms.size()); + double l2_error_sq = + error_obj.evaluate(TimeInfo(time + dt, dt, 0), shape_disp.get().get(), getConstFieldPointers({final_disp})); + double l2_error = std::sqrt(l2_error_sq); - return block_solve(weak_forms, index_map, field_store.getShapeDisp(), inputs, params, time_info, solver, - field_store.getBoundaryConditionManagers()); + std::cout << "Static Elasticity L2 Error (affine patch test): " << l2_error << std::endl; + EXPECT_LT(l2_error, 1e-11); // Should be machine precision for affine displacement } -TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) +TEST_F(ThermoMechanicsMeshFixture, TransientThermoMechanicsCompilation) { SMITH_MARK_FUNCTION; - FieldType> shape_disp_type("shape_displacement"); - FieldType> disp_type("displacement"); - FieldType> temperature_type("temperature"); - - std::shared_ptr d_nonlinear_solver = - buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); - double rho = 1.0; - double E = 100.0; + double E0 = 100.0; double nu = 0.25; - double c = 1.0; - double alpha = 1.0e-3; - double theta_ref = 0.0; - double k = 1.0; - GreenSaintVenantThermoelasticMaterial material{rho, E, nu, c, alpha, theta_ref, k}; - - FieldStore field_store(mesh_, 100); - - field_store.addShapeDisp(shape_disp_type); - - std::shared_ptr& disp_bc = field_store.addUnknown(disp_type); - disp_bc->setVectorBCs(mesh_->domain("left"), [](double t, smith::tensor X) { - auto bc = 0.0 * X; - bc[0] = 0.01 * t; - return bc; + double specific_heat = 1.0; + double kappa = 0.1; + GreenSaintVenantThermoelasticMaterial material{rho, E0, nu, specific_heat, 0.0, 1.0, kappa}; + + auto fast_nonlinear_opts = nonlinear_opts; + fast_nonlinear_opts.max_iterations = 5; + auto solver = buildDifferentiableNonlinearBlockSolver(fast_nonlinear_opts, linear_options, *mesh_); + + FieldType> youngs_modulus("youngs_modulus"); + auto system = buildThermoMechanicsSystem(mesh_, solver, youngs_modulus); + system.setMaterial(material, mesh_->entireBodyName()); + + system.parameter_fields[0].get()->setFromFieldFunction([=](smith::tensor) { return E0; }); + + // Add Solid Body Force (Gravity) + system.addSolidBodyForce(mesh_->entireBodyName(), [](double /*time*/, auto /*X*/, auto /*u*/, auto /*v*/, auto /*T*/, + auto /*T_dot*/, auto... /*params*/) { + tensor f{}; + f[1] = -9.81; + return f; }); - disp_bc->setFixedVectorBCs(mesh_->domain("right")); - - auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); - - std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); - temperature_bc->setFixedScalarBCs(mesh_->domain("left")); - temperature_bc->setFixedScalarBCs(mesh_->domain("right")); - auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); - - QuasiStaticFirstOrderTimeIntegrationRule disp_time_rule; - BackwardEulerFirstOrderTimeIntegrationRule temperature_time_rule; - - auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, - temperature_type, temperature_old_type); - - solid_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto disp, auto disp_old, - auto temperature, auto temperature_old) { - auto u = disp_time_rule.value(t_info, disp, disp_old); - auto v = disp_time_rule.dot(t_info, disp, disp_old); - auto T = temperature_time_rule.value(t_info, temperature, temperature_old); - GreenSaintVenantThermoelasticMaterial::State state; - auto [pk, C_v, s0, q0] = - material(t_info.dt(), state, get(u), get(v), get(T), get(T)); - return smith::tuple{smith::zero{}, pk}; + // Add Solid Traction + system.addSolidTraction("right", [](double /*time*/, auto /*X*/, auto /*n*/, auto /*u*/, auto /*v*/, auto /*T*/, + auto /*T_dot*/, auto... /*params*/) { + tensor t{}; + t[0] = 1.0; + return t; }); - auto thermal_weak_form = createWeakForm("thermal_flux", temperature_type, field_store, temperature_type, - temperature_old_type, disp_type, disp_old_type); - - thermal_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto temperature, - auto temperature_old, auto disp, auto disp_old) { - GreenSaintVenantThermoelasticMaterial::State state; - auto u = disp_time_rule.value(t_info, disp, disp_old); - auto v = disp_time_rule.dot(t_info, disp, disp_old); - auto T = temperature_time_rule.value(t_info, temperature, temperature_old); - auto T_dot = temperature_time_rule.dot(t_info, temperature, temperature_old); - auto [pk, C_v, s0, q0] = - material(t_info.dt(), state, get(u), get(v), get(T), get(T)); - auto dT_dt = get(T_dot); - return smith::tuple{C_v * dT_dt - s0, -q0}; - }); + // Add Thermal Heat Source + system.addThermalHeatSource(mesh_->entireBodyName(), + [](double /*time*/, auto /*X*/, auto /*u*/, auto /*v*/, auto /*T*/, auto /*T_dot*/, + auto... /*params*/) { return 10.0; }); - thermal_weak_form->addBodySource(smith::DependsOn<>(), mesh_->entireBodyName(), - [](auto /*t*/, auto /* x */) { return 100.0; }); + // Add Thermal Heat Flux + system.addThermalHeatFlux("left", [](double /*time*/, auto /*X*/, auto /*n*/, auto /*u*/, auto /*v*/, auto /*T*/, + auto /*T_dot*/, auto... /*params*/) { + return 5.0; // Flux into domain + }); - std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; - std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); + system.disp_bc->setVectorBCs(mesh_->domain("left"), + [](double /*t*/, smith::tensor X) { return 0.0 * X; }); + // Don't run anything, just make sure the templates build as expected. +} - // auto states = field_store.getFields(); +TEST_F(ThermoMechanicsMeshFixture, PressureBC) +{ + SMITH_MARK_FUNCTION; - EXPECT_EQ(0, 0); + double rho = 1.0; + double E0 = 100.0; + double nu = 0.25; + double specific_heat = 1.0; + double kappa = 0.1; + GreenSaintVenantThermoelasticMaterial material{rho, E0, nu, specific_heat, 0.0, 1.0, kappa}; + + auto solver = buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + FieldType> youngs_modulus("youngs_modulus"); + auto system = buildThermoMechanicsSystem(mesh_, solver, youngs_modulus); + system.setMaterial(material, mesh_->entireBodyName()); + + system.parameter_fields[0].get()->setFromFieldFunction([=](smith::tensor) { return E0; }); + + // Fixed left side + system.disp_bc->setFixedVectorBCs(mesh_->domain("left")); + system.temperature_bc->setFixedScalarBCs(mesh_->domain("left")); + system.temperature_bc->setFixedScalarBCs(mesh_->domain("right")); + + // Apply pressure on right side + double pressure_mag = 0.1; + system.addPressure("right", [pressure_mag](double, auto, auto, auto, auto, auto, auto...) { return pressure_mag; }); + + auto shape_disp = system.field_store->getShapeDisp(); + auto states = system.getStateFields(); + auto params = system.getParameterFields(); + + double dt = 0.001; + double time = 0.0; + size_t cycle = 0; + // Advance one step + TimeInfo t_info(time, dt, cycle); + auto result = system.advancer->advanceState(t_info, shape_disp, states, params); + states = result.first; + + // Check that we have some displacement + auto disp_idx = system.field_store->getFieldIndex("displacement_predicted"); + FieldState final_disp = states[disp_idx]; + + // We expect non-zero displacement + double disp_norm = final_disp.get()->Norml2(); + EXPECT_GT(disp_norm, 1e-6); } } // namespace smith diff --git a/src/smith/differentiable_numerics/thermal_system.hpp b/src/smith/differentiable_numerics/thermal_system.hpp new file mode 100644 index 0000000000..1732cf031a --- /dev/null +++ b/src/smith/differentiable_numerics/thermal_system.hpp @@ -0,0 +1,168 @@ +// 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 thermal_system.hpp + * @brief Defines the ThermalSystem struct and its factory function + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { + +/** + * @brief Container for a thermal system. + * @tparam dim Spatial dimension. + * @tparam temp_order Order of the temperature basis. + * @tparam parameter_space Finite element spaces for optional parameters. + */ +template +struct ThermalSystem : public SystemBase { + /// @brief using for ThermalWeakFormType + using ThermalWeakFormType = + TimeDiscretizedWeakForm, Parameters, parameter_space...>>; + + std::shared_ptr thermal_weak_form; ///< Thermal weak form. + std::shared_ptr temperature_bc; ///< Temperature boundary conditions. + std::shared_ptr temperature_time_rule; ///< Time integration for temperature. + + /** + * @brief Get the list of all state fields. + * @return std::vector List of state fields. + */ + std::vector getStateFields() const + { + std::vector states; + states.push_back(field_store->getField(prefix("temperature"))); + return states; + } + + /** + * @brief Create a DifferentiablePhysics object for this system. + * @param physics_name The name of the physics. + * @return std::shared_ptr The differentiable physics object. + */ + std::shared_ptr createDifferentiablePhysics(std::string physics_name) + { + return std::make_shared( + field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), + getParameterFields(), advancer, physics_name, std::vector{prefix("thermal_flux")}); + } + + /** + * @brief Set the material model for a domain. + * @tparam ThermalIntegrandType Function with signature (TimeInfo, X, T, params...) -> {residual, flux}. + * @param domain_name The name of the domain to apply the material to. + * @param integrand The thermal integrand function. + */ + template + void setThermalIntegrand(const std::string& domain_name, ThermalIntegrandType integrand) + { + auto captured_temp_rule = temperature_time_rule; + + thermal_weak_form->addBodyIntegral(domain_name, [=](auto t_info, auto X, auto temperature, auto... params) { + // Apply time integration to get current state + auto T = captured_temp_rule->value(t_info, temperature); + return integrand(t_info, X, T, params...); + }); + } + + /** + * @brief Add a body heat source to the thermal system. + * @tparam HeatSourceType Function with signature (t, X, T, params...) -> heat_source. + * @param domain_name The name of the domain to apply the source to. + * @param source_function The heat source function. + */ + template + void addBodyHeatSource(const std::string& domain_name, HeatSourceType source_function) + { + auto captured_temp_rule = temperature_time_rule; + + thermal_weak_form->addBodySource(domain_name, [=](auto t_info, auto X, auto temperature, auto... params) { + auto T = captured_temp_rule->value(t_info, temperature); + return source_function(t_info.time(), X, T, params...); + }); + } + + /** + * @brief Add a boundary heat flux to the thermal system. + * @tparam HeatFluxType Function with signature (t, X, n, T, params...) -> heat_flux. + * @param boundary_name The name of the boundary to apply the flux to. + * @param flux_function The heat flux function. + */ + template + void addBoundaryHeatFlux(const std::string& boundary_name, HeatFluxType flux_function) + { + auto captured_temp_rule = temperature_time_rule; + + thermal_weak_form->addBoundaryFlux(boundary_name, + [=](auto t_info, auto X, auto n, auto temperature, auto... params) { + auto T = captured_temp_rule->value(t_info, temperature); + return -flux_function(t_info.time(), X, n, T, params...); + }); + } +}; + +/** + * @brief Factory function to build a thermal system. + */ +template +ThermalSystem buildThermalSystem(std::shared_ptr mesh, + std::shared_ptr solver, + std::string prepend_name = "", + FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100); + + auto prefix = [&](const std::string& name) { + if (prepend_name.empty()) { + return name; + } + return prepend_name + "_" + name; + }; + + FieldType> shape_disp_type(prefix("shape_displacement")); + field_store->addShapeDisp(shape_disp_type); + + // Temperature field with quasi-static rule (1 state) + auto temperature_time_rule = std::make_shared(); + FieldType> temperature_type(prefix("temperature")); + auto temperature_bc = field_store->addIndependent(temperature_type, temperature_time_rule); + + std::vector parameter_fields; + (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); + (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + + // Thermal weak form + std::string thermal_flux_name = prefix("thermal_flux"); + auto thermal_weak_form = + std::make_shared::ThermalWeakFormType>( + thermal_flux_name, field_store->getMesh(), field_store->getField(temperature_type.name).get()->space(), + field_store->createSpaces(thermal_flux_name, temperature_type.name, temperature_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Build solver and advancer + std::vector> weak_forms{thermal_weak_form}; + auto advancer = std::make_shared(field_store, weak_forms, solver); + + return ThermalSystem{ + {field_store, solver, advancer, parameter_fields, prepend_name}, + thermal_weak_form, + temperature_bc, + temperature_time_rule}; +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/thermo_mechanics_system.hpp b/src/smith/differentiable_numerics/thermo_mechanics_system.hpp new file mode 100644 index 0000000000..13738d23a9 --- /dev/null +++ b/src/smith/differentiable_numerics/thermo_mechanics_system.hpp @@ -0,0 +1,478 @@ +// 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 thermo_mechanics_system.hpp + * @brief Defines the ThermoMechanicsSystem struct and its factory function + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { + +/** + * @brief Container for a coupled thermo-mechanical system. + * @tparam dim Spatial dimension. + * @tparam disp_order Order of the displacement basis. + * @tparam temp_order Order of the temperature basis. + * @tparam parameter_space Finite element spaces for optional parameters. + */ +template +struct ThermoMechanicsSystem : public SystemBase { + /// @brief using for SolidWeakFormType + using SolidWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, H1, H1, H1, parameter_space...>>; + + /// @brief using for ThermalWeakFormType + using ThermalWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, H1, H1, H1, parameter_space...>>; + + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + std::shared_ptr thermal_weak_form; ///< Thermal weak form. + std::shared_ptr disp_bc; ///< Displacement boundary conditions. + std::shared_ptr temperature_bc; ///< Temperature boundary conditions. + std::shared_ptr disp_time_rule; ///< Time integration for displacement. + std::shared_ptr + temperature_time_rule; ///< Time integration for temperature. + + /** + * @brief Get the list of all state fields (current and old). + * @return std::vector List of state fields. + */ + std::vector getStateFields() const + { + std::vector states; + states.push_back(field_store->getField(prefix("displacement_predicted"))); + states.push_back(field_store->getField(prefix("displacement"))); + states.push_back(field_store->getField(prefix("temperature_predicted"))); + states.push_back(field_store->getField(prefix("temperature"))); + return states; + } + + /** + * @brief Create a DifferentiablePhysics object for this system. + * @param physics_name The name of the physics. + * @return std::shared_ptr The differentiable physics object. + */ + std::shared_ptr createDifferentiablePhysics(std::string physics_name) + { + return std::make_shared( + field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), + getParameterFields(), advancer, physics_name, + std::vector{prefix("solid_force"), prefix("thermal_flux")}); + } + + /** + * @brief Set the material model for a domain, defining integrals for solid and thermal weak forms. + * @tparam MaterialType The material model type. + * @param material The material model instance. + * @param domain_name The name of the domain to apply the material to. + */ + template + void setMaterial(const MaterialType& material, const std::string& domain_name) + { + // Solid weak form: inputs are (u, u_old, temperature, temperature_old, params...) + // Manually apply time integration rules to get current state + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + solid_weak_form->addBodyIntegral(domain_name, [=](auto t_info, auto /*X*/, auto u, auto u_old, auto temperature, + auto temperature_old, auto... params) { + // Apply time integration to get current state + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto v_current = captured_disp_rule->dot(t_info, u, u_old); + auto T = captured_temp_rule->value(t_info, temperature, temperature_old); + + typename MaterialType::State state; + auto [pk, C_v, s0, q0] = material(t_info.dt(), state, get(u_current), get(v_current), + get(T), get(T), params...); + return smith::tuple{zero{}, pk}; + }); + + // Thermal weak form: inputs are (T, T_old, u, u_old, params...) + // Manually apply time integration rules to get current state + thermal_weak_form->addBodyIntegral( + domain_name, [=](auto t_info, auto /*X*/, auto T, auto T_old, auto disp, auto disp_old, auto... params) { + // Apply time integration to get current state + auto T_current = captured_temp_rule->value(t_info, T, T_old); + auto T_dot = captured_temp_rule->dot(t_info, T, T_old); + auto u = captured_disp_rule->value(t_info, disp, disp_old); + auto v = captured_disp_rule->dot(t_info, disp, disp_old); + + typename MaterialType::State state; + auto [pk, C_v, s0, q0] = material(t_info.dt(), state, get(u), get(v), + get(T_current), get(T_current), params...); + auto dT_dt = get(T_dot); + return smith::tuple{C_v * dT_dt - s0, -q0}; + }); + } + + /** + * @brief Add a body force to the solid mechanics part of the system (with DependsOn). + * @tparam active_parameters Indices of fields this force depends on. + * @tparam BodyForceType The body force function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the domain to apply the force to. + * @param force_function The force function (t, X, u, v, T, T_dot, selected params...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addSolidBodyForce(DependsOn depends_on, const std::string& domain_name, + BodyForceType force_function) + { + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + solid_weak_form->addBodySource( + depends_on, domain_name, + [=](auto t_info, auto X, auto u, auto u_old, auto temperature, auto temperature_old, auto... params) { + // Apply time integration to get current state + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto v_current = captured_disp_rule->dot(t_info, u, u_old); + auto current_T = captured_temp_rule->value(t_info, temperature, temperature_old); + auto T_dot = captured_temp_rule->dot(t_info, temperature, temperature_old); + + return force_function(t_info.time(), X, u_current, v_current, current_T, T_dot, params...); + }); + } + + /** + * @brief Add a body force to the solid mechanics part of the system. + * @tparam BodyForceType The body force function type. + * @param domain_name The name of the domain to apply the force to. + * @param force_function The force function (t, X, u, v, T, T_dot, params...). + */ + template + void addSolidBodyForce(const std::string& domain_name, BodyForceType force_function) + { + addSolidBodyForceAllParams(domain_name, force_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a surface flux (traction) to the solid mechanics part of the system (with DependsOn). + * @tparam active_parameters Indices of fields this flux depends on. + * @tparam SurfaceFluxType The surface flux function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the boundary domain to apply the flux to. + * @param flux_function The flux function (t, X, n, u, v, T, T_dot, selected params...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addSolidTraction(DependsOn depends_on, const std::string& domain_name, + SurfaceFluxType flux_function) + { + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + solid_weak_form->addBoundaryFlux( + depends_on, domain_name, + [=](auto t_info, auto X, auto n, auto u, auto u_old, auto temperature, auto temperature_old, auto... params) { + // Apply time integration to get current state + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto v_current = captured_disp_rule->dot(t_info, u, u_old); + auto current_T = captured_temp_rule->value(t_info, temperature, temperature_old); + auto T_dot = captured_temp_rule->dot(t_info, temperature, temperature_old); + + return flux_function(t_info.time(), X, n, u_current, v_current, current_T, T_dot, params...); + }); + } + + /** + * @brief Add a surface flux (traction) to the solid mechanics part of the system. + * @tparam SurfaceFluxType The surface flux function type. + * @param domain_name The name of the boundary domain to apply the flux to. + * @param flux_function The flux function (t, X, n, u, v, T, T_dot, params...). + */ + template + void addSolidTraction(const std::string& domain_name, SurfaceFluxType flux_function) + { + addSolidTractionAllParams(domain_name, flux_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a body source (heat source) to the thermal part of the system (with DependsOn). + * @tparam active_parameters Indices of fields this source depends on. + * @tparam BodySourceType The body source function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the domain to apply the source to. + * @param source_function The source function (t, X, u, v, T, T_dot, selected params...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addThermalHeatSource(DependsOn depends_on, const std::string& domain_name, + BodySourceType source_function) + { + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + thermal_weak_form->addBodySource( + depends_on, domain_name, + [=](auto t_info, auto X, auto T, auto T_old, auto disp, auto disp_old, auto... params) { + // Apply time integration to get current state + auto current_u = captured_disp_rule->value(t_info, disp, disp_old); + auto v_current = captured_disp_rule->dot(t_info, disp, disp_old); + auto T_current = captured_temp_rule->value(t_info, T, T_old); + auto T_dot = captured_temp_rule->dot(t_info, T, T_old); + + return source_function(t_info.time(), X, current_u, v_current, T_current, T_dot, params...); + }); + } + + /** + * @brief Add a body source (heat source) to the thermal part of the system. + * @tparam BodySourceType The body source function type. + * @param domain_name The name of the domain to apply the source to. + * @param source_function The source function (t, X, u, v, T, T_dot, params...). + */ + template + void addThermalHeatSource(const std::string& domain_name, BodySourceType source_function) + { + addThermalHeatSourceAllParams(domain_name, source_function, + std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a surface flux (heat flux) to the thermal part of the system (with DependsOn). + * @tparam active_parameters Indices of fields this flux depends on. + * @tparam SurfaceFluxType The surface flux function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the boundary domain to apply the flux to. + * @param flux_function The flux function (t, X, n, u, v, T, T_dot, selected params...). + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addThermalHeatFlux(DependsOn depends_on, const std::string& domain_name, + SurfaceFluxType flux_function) + { + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + thermal_weak_form->addBoundaryFlux( + depends_on, domain_name, + [=](auto t_info, auto X, auto n, auto T, auto T_old, auto disp, auto disp_old, auto... params) { + // Apply time integration to get current state + auto current_u = captured_disp_rule->value(t_info, disp, disp_old); + auto v_current = captured_disp_rule->dot(t_info, disp, disp_old); + auto T_current = captured_temp_rule->value(t_info, T, T_old); + auto T_dot = captured_temp_rule->dot(t_info, T, T_old); + + return -flux_function(t_info.time(), X, n, current_u, v_current, T_current, T_dot, params...); + }); + } + + /** + * @brief Add a surface flux (heat flux) to the thermal part of the system. + * @tparam SurfaceFluxType The surface flux function type. + * @param domain_name The name of the boundary domain to apply the flux to. + * @param flux_function The flux function (t, X, n, u, v, T, T_dot, params...). + */ + template + void addThermalHeatFlux(const std::string& domain_name, SurfaceFluxType flux_function) + { + addThermalHeatFluxAllParams(domain_name, flux_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + /** + * @brief Add a pressure boundary condition (follower force) to the solid mechanics part of the system (with + * DependsOn). + * @tparam active_parameters Indices of fields this pressure depends on. + * @tparam PressureType The pressure function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param domain_name The name of the boundary domain. + * @param pressure_function The pressure function (t, X, u, v, T, T_dot, selected params...). + * @note Pressure is applied in the current configuration: P * n_deformed. + * @note Time integration is applied to the state fields before calling the user function. + */ + template + void addPressure(DependsOn depends_on, const std::string& domain_name, + PressureType pressure_function) + { + auto captured_disp_rule = disp_time_rule; + auto captured_temp_rule = temperature_time_rule; + + solid_weak_form->addBoundaryIntegral( + depends_on, domain_name, + [=](auto t_info, auto X, auto u, auto u_old, auto temperature, auto temperature_old, auto... params) { + // Apply time integration to get current state + auto u_current = captured_disp_rule->value(t_info, u, u_old); + auto v_current = captured_disp_rule->dot(t_info, u, u_old); + auto T_current = captured_temp_rule->value(t_info, temperature, temperature_old); + auto T_dot = captured_temp_rule->dot(t_info, temperature, temperature_old); + + // Compute deformed normal and apply correction for reference configuration integration + auto x_current = X + u_current; + auto n_deformed = cross(get(x_current)); + auto n_shape_norm = norm(cross(get(X))); + + auto pressure = pressure_function(t_info.time(), get(X), u_current, v_current, T_current, T_dot, + get(params)...); + + // Return traction vector (force) + return pressure * n_deformed * (1.0 / n_shape_norm); + }); + } + + /** + * @brief Add a pressure boundary condition (follower force) to the solid mechanics part of the system. + * @tparam PressureType The pressure function type. + * @param domain_name The name of the boundary domain. + * @param pressure_function The pressure function (t, X, u, v, T, T_dot, params...). + * @note Pressure is applied in the current configuration: P * n_deformed. + */ + template + void addPressure(const std::string& domain_name, PressureType pressure_function) + { + addPressureAllParams(domain_name, pressure_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{}); + } + + private: + // Helper functions to forward non-DependsOn calls to DependsOn versions with all parameters + template + void addSolidBodyForceAllParams(const std::string& domain_name, BodyForceType force_function, + std::index_sequence) + { + addSolidBodyForce(DependsOn(Is)...>{}, domain_name, force_function); + } + + template + void addSolidTractionAllParams(const std::string& domain_name, SurfaceFluxType flux_function, + std::index_sequence) + { + addSolidTraction(DependsOn(Is)...>{}, domain_name, flux_function); + } + + template + void addPressureAllParams(const std::string& domain_name, PressureType pressure_function, std::index_sequence) + { + addPressure(DependsOn(Is)...>{}, domain_name, pressure_function); + } + + template + void addThermalHeatSourceAllParams(const std::string& domain_name, BodySourceType source_function, + std::index_sequence) + { + addThermalHeatSource(DependsOn(Is)...>{}, domain_name, source_function); + } + + template + void addThermalHeatFluxAllParams(const std::string& domain_name, SurfaceFluxType flux_function, + std::index_sequence) + { + addThermalHeatFlux(DependsOn(Is)...>{}, domain_name, flux_function); + } +}; + +/** + * @brief Factory function to build a thermo-mechanical system. + * @tparam dim Spatial dimension. + * @tparam disp_order Order of the displacement basis. + * @tparam temp_order Order of the temperature basis. + * @tparam parameter_space Finite element spaces for optional parameters. + * @param mesh The mesh. + * @param solver The differentiable block solver. + * @param prepend_name The name of the physics (used as field prefix). + * @param parameter_types Parameter field types. + * @return ThermoMechanicsSystem with all components initialized. + */ +template +ThermoMechanicsSystem buildThermoMechanicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", + FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100); + + auto prefix = [&](const std::string& name) { + if (prepend_name.empty()) { + return name; + } + return prepend_name + "_" + name; + }; + + FieldType> shape_disp_type(prefix("shape_displacement")); + field_store->addShapeDisp(shape_disp_type); + + // Displacement field with quasi-static time integration + auto disp_time_rule = std::make_shared(); + FieldType> disp_type(prefix("displacement_predicted")); + auto disp_bc = field_store->addIndependent(disp_type, disp_time_rule); + auto disp_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::VAL, prefix("displacement")); + + // Temperature field with backward Euler time integration + auto temperature_time_rule = std::make_shared(); + FieldType> temperature_type(prefix("temperature_predicted")); + auto temperature_bc = field_store->addIndependent(temperature_type, temperature_time_rule); + auto temperature_old_type = + field_store->addDependent(temperature_type, FieldStore::TimeDerivative::VAL, prefix("temperature")); + + std::vector parameter_fields; + (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); + (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + + // Solid mechanics weak form + std::string solid_force_name = prefix("solid_force"); + auto solid_weak_form = std::make_shared< + typename ThermoMechanicsSystem::SolidWeakFormType>( + solid_force_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), + field_store->createSpaces(solid_force_name, disp_type.name, disp_type, disp_old_type, temperature_type, + temperature_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Thermal weak form + std::string thermal_flux_name = prefix("thermal_flux"); + auto thermal_weak_form = std::make_shared< + typename ThermoMechanicsSystem::ThermalWeakFormType>( + thermal_flux_name, field_store->getMesh(), field_store->getField(temperature_type.name).get()->space(), + field_store->createSpaces(thermal_flux_name, temperature_type.name, temperature_type, temperature_old_type, + disp_type, disp_old_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Build solver and advancer + std::vector> weak_forms{solid_weak_form, thermal_weak_form}; + auto advancer = std::make_shared(field_store, weak_forms, solver); + + return ThermoMechanicsSystem{ + {field_store, solver, advancer, parameter_fields, prepend_name}, + solid_weak_form, + thermal_weak_form, + disp_bc, + temperature_bc, + disp_time_rule, + temperature_time_rule}; +} + +/** + * @brief Factory function to build a thermo-mechanical system (without physics name). + * @tparam dim Spatial dimension. + * @tparam disp_order Order of the displacement basis. + * @tparam temp_order Order of the temperature basis. + * @tparam parameter_space Finite element spaces for optional parameters. + * @param mesh The mesh. + * @param solver The differentiable block solver. + * @param parameter_types Parameter field types. + * @return ThermoMechanicsSystem with all components initialized. + */ +template +ThermoMechanicsSystem buildThermoMechanicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, + FieldType... parameter_types) +{ + return buildThermoMechanicsSystem(mesh, solver, "", + parameter_types...); +} + +} // namespace smith diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index 0fc42331f4..a11d423b72 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -7,8 +7,17 @@ /** * @file time_discretized_weak_form.hpp * - * @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of - * equations + * @brief Wraps FunctionalWeakForm to provide TimeInfo (time, dt, cycle) to integrands instead of just time. + * + * This class provides a thin wrapper around FunctionalWeakForm that automatically converts the time + * parameter into a TimeInfo struct containing time, timestep size (dt), and cycle number. This allows + * physics systems to access timestep information needed for time integration. + * + * Key features: + * - All integrands receive TimeInfo instead of just time + * - Systems are responsible for manually applying time integration rules + * - Supports body integrals, boundary integrals, boundary fluxes, and interior boundary integrals + * - Default behavior passes ALL input fields to integrands (can be overridden with DependsOn) */ #pragma once @@ -23,141 +32,223 @@ 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 +/** + * @brief A time-discretized weak form that provides TimeInfo to integrands. + * + * This wraps FunctionalWeakForm to pass TimeInfo (containing time, dt, and cycle) to all + * integrand functions instead of just the time value. This allows physics models to access + * timestep information for time integration. + * + * @tparam spatial_dim The spatial dimension for the problem. + * @tparam OutputSpace The output residual for the weak form (test-space). + * @tparam InputSpaces All the input FiniteElementState fields (trial-spaces). + */ template class TimeDiscretizedWeakForm> : public FunctionalWeakForm> { public: - using WeakFormT = FunctionalWeakForm>; ///< using - - /// Constructor + using Base = FunctionalWeakForm>; ///< Base class alias + + /** + * @brief Construct a time-discretized weak form. + * @param physics_name Unique name for this physics module. + * @param mesh The computational mesh. + * @param output_mfem_space The test function space (output/residual space). + * @param input_mfem_spaces Vector of trial function spaces (input field spaces). + */ 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) + const typename Base::SpacesT& input_mfem_spaces) + : Base(physics_name, mesh, output_mfem_space, input_mfem_spaces) { } - /// @overload + /** + * @brief Add a body integral with TimeInfo. + * + * The integrand receives TimeInfo (containing time, dt, cycle) instead of just time. + * The system is responsible for manually applying time integration rules inside the integrand + * to compute current state values from the raw field history. + * + * @tparam active_parameters Indices of fields this integral depends on. + * @tparam BodyIntegralType The integrand function type. + * @param depends_on Dependency specification for which input fields to pass. + * @param body_name The name of the domain. + * @param integrand Function with signature (TimeInfo, X, inputs...) -> residual. + */ 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) { + Base::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 addBodyIntegralImpl(std::string body_name, BodyForceType body_integral, - std::integer_sequence) + /** + * @brief Add a body integral with TimeInfo (defaults to all input fields). + * @tparam BodyIntegralType The integrand function type. + * @param body_name The name of the domain. + * @param integrand Function with signature (TimeInfo, X, inputs...) -> residual. + */ + template + void addBodyIntegral(std::string body_name, BodyIntegralType integrand) { - addBodyIntegral(DependsOn{}, body_name, body_integral); + constexpr int num_inputs = sizeof...(InputSpaces); + addBodyIntegralWithAllParams(body_name, integrand, std::make_integer_sequence{}); } - /// @overload - template - void addBodyIntegral(std::string body_name, BodyForceType body_integral) + /** + * @brief Add a body source (body load) with TimeInfo. + * @tparam active_parameters Indices of fields this source depends on. + * @tparam BodyLoadType The load function type. + * @param depends_on Dependency specification. + * @param body_name The name of the domain. + * @param load_function Function with signature (TimeInfo, X, inputs...) -> load vector. + */ + template + void addBodySource(DependsOn depends_on, std::string body_name, BodyLoadType load_function) { - addBodyIntegralImpl(body_name, body_integral, std::make_integer_sequence{}); + addBodyIntegral(depends_on, body_name, [load_function](auto t_info, auto X, auto... inputs) { + return smith::tuple{-load_function(t_info, get(X), get(inputs)...), smith::zero{}}; + }); } -}; -/// @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 publicly available abstract weak form is a - ///< functions of the current u, u_old, v_old, and a_old, - std::shared_ptr - final_reaction_weak_form; ///< this publicly available abstract weak form is structly a - ///< function of the current u, v, and a (no time discretization) - ///< its main purpose is to compute reaction forces after the solve is completed -}; - -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 -/// ImplicitNewmarkSecondOrderTimeIntegrationRule 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 + /// defaults to use all parameters + /// @overload + template + void addBodySource(std::string body_name, BodyLoadType load_function) + { + constexpr int num_inputs = sizeof...(InputSpaces); + addBodySourceWithAllParams(body_name, load_function, std::make_integer_sequence{}); + } - using TimeDiscretizedWeakFormT = - TimeDiscretizedWeakForm>; ///< using - using FinalReactionFormT = TimeDiscretizedWeakForm>; ///< using + /** + * @brief Add a boundary integral with TimeInfo. + * @tparam active_parameters Indices of fields this integral depends on. + * @tparam BoundaryIntegralType The boundary integrand function type. + * @param depends_on Dependency specification. + * @param boundary_name The name of the boundary. + * @param integrand Function with signature (TimeInfo, X, inputs...) -> residual. + */ + template + void addBoundaryIntegral(DependsOn depends_on, std::string boundary_name, + BoundaryIntegralType integrand) + { + const double* dt = &this->dt_; + const size_t* cycle = &this->cycle_; + Base::addBoundaryIntegral(depends_on, boundary_name, [dt, cycle, integrand](double t, auto X, auto... inputs) { + TimeInfo time_info(t, *dt, *cycle); + return integrand(time_info, X, inputs...); + }); + } - /// @brief Constructor - SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr mesh, - ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, - const mfem::ParFiniteElementSpace& output_mfem_space, - const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces) - : time_rule_(time_rule) + /// defaults to use all parameters + /// @overload + template + void addBoundaryIntegral(std::string boundary_name, BoundaryIntegralType integrand) { - time_discretized_weak_form_ = - std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces); - time_discretized_weak_form = time_discretized_weak_form_; + constexpr int num_inputs = sizeof...(InputSpaces); + addBoundaryIntegralWithAllParams(boundary_name, integrand, std::make_integer_sequence{}); + } - typename TimeDiscretizedWeakFormT::SpacesT input_mfem_spaces_trial_removed(std::next(input_mfem_spaces.begin()), - input_mfem_spaces.end()); - final_reaction_weak_form_ = - std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); - final_reaction_weak_form = final_reaction_weak_form_; + /** + * @brief Add a boundary flux with TimeInfo. + * @tparam active_parameters Indices of fields this integral depends on. + * @tparam BoundaryFluxType The flux function type. + * @param depends_on Dependency specification. + * @param boundary_name The name of the boundary. + * @param flux_function Function with signature (TimeInfo, X, n, inputs...) -> flux. + */ + template + void addBoundaryFlux(DependsOn depends_on, std::string boundary_name, + BoundaryFluxType flux_function) + { + const double* dt = &this->dt_; + const size_t* cycle = &this->cycle_; + Base::addBoundaryFlux(depends_on, boundary_name, + [dt, cycle, flux_function](double t, auto X, auto n, auto... inputs) { + TimeInfo time_info(t, *dt, *cycle); + return flux_function(time_info, X, n, inputs...); + }); } + /// defaults to use all parameters /// @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.dot(t, U, U_old, U_dot_old, U_dot_dot_old), - time_rule.ddot(t, U, U_old, U_dot_old, U_dot_dot_old), inputs...); - }); - final_reaction_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, - body_name, integrand); + template + void addBoundaryFlux(std::string boundary_name, BoundaryFluxType flux_function) + { + constexpr int num_inputs = sizeof...(InputSpaces); + addBoundaryFluxWithAllParams(boundary_name, flux_function, std::make_integer_sequence{}); } + /** + * @brief Add an interior boundary integral with TimeInfo. + * @tparam active_parameters Indices of fields this integral depends on. + * @tparam InteriorIntegralType The integrand function type. + * @param depends_on Dependency specification. + * @param interior_name The name of the interior boundary. + * @param integrand Function with signature (TimeInfo, X, inputs...) -> residual. + */ + template + void addInteriorBoundaryIntegral(DependsOn depends_on, std::string interior_name, + InteriorIntegralType integrand) + { + const double* dt = &this->dt_; + const size_t* cycle = &this->cycle_; + Base::addInteriorBoundaryIntegral(depends_on, interior_name, + [dt, cycle, integrand](double t, auto X, auto... inputs) { + TimeInfo time_info(t, *dt, *cycle); + return integrand(time_info, X, inputs...); + }); + } + + /// defaults to use all parameters /// @overload - template - void addBodyIntegral(std::string body_name, BodyForceType body_integral) + template + void addInteriorBoundaryIntegral(std::string interior_name, InteriorIntegralType integrand) { - addBodyIntegral(DependsOn<>{}, body_name, body_integral); + constexpr int num_inputs = sizeof...(InputSpaces); + addInteriorBoundaryIntegralWithAllParams(interior_name, integrand, std::make_integer_sequence{}); } 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 - final_reaction_weak_form_; ///< fully template underlying weak form (no time integration included, a function of - ///< current u, v, and a) - - ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_; ///< encodes the time integration rule + template + void addBodyIntegralWithAllParams(std::string body_name, BodyIntegralType integrand, + std::integer_sequence) + { + addBodyIntegral(DependsOn{}, body_name, integrand); + } + + template + void addBodySourceWithAllParams(std::string body_name, BodyLoadType load_function, + std::integer_sequence) + { + addBodySource(DependsOn{}, body_name, load_function); + } + + template + void addBoundaryIntegralWithAllParams(std::string boundary_name, BoundaryIntegralType integrand, + std::integer_sequence) + { + addBoundaryIntegral(DependsOn{}, boundary_name, integrand); + } + + template + void addInteriorBoundaryIntegralWithAllParams(std::string interior_name, InteriorIntegralType integrand, + std::integer_sequence) + { + addInteriorBoundaryIntegral(DependsOn{}, interior_name, integrand); + } + + template + void addBoundaryFluxWithAllParams(std::string boundary_name, BoundaryFluxType flux_function, + std::integer_sequence) + { + addBoundaryFlux(DependsOn{}, boundary_name, flux_function); + } }; } // namespace smith diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index b6b2b6d681..e3ca6b1c31 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -14,17 +14,44 @@ #pragma once #include "smith/physics/common.hpp" +#include "smith/differentiable_numerics/field_state.hpp" namespace smith { +/// @brief Abstract time integration rule for discretizing odes in time +class TimeIntegrationRule { + public: + /// @brief destructor + virtual ~TimeIntegrationRule() {} + + /// @brief update the current value of the independent variable, given the predicted value of the current independent + /// variable, followed by + virtual FieldState corrected_value(const TimeInfo& t, const std::vector& states) const = 0; + + /// @brief get the number of states required by the rule + virtual int num_args() const = 0; + + /// @brief update the current value of the independent variable's first time derivative, given the predicted value of + /// the current independent variable, followed by + virtual FieldState corrected_dot(const TimeInfo& t, const std::vector& states) const = 0; + + /// @brief update the current value of the independent variable's second time derivative, given the predicted value of + /// the current independent variable, followed by + virtual FieldState corrected_ddot(const TimeInfo& t, const std::vector& states) const = 0; +}; + /// @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 { +class BackwardEulerFirstOrderTimeIntegrationRule : public TimeIntegrationRule { + public: /// @brief Constructor BackwardEulerFirstOrderTimeIntegrationRule() {} + /// @brief get the number of states required by the rule + int num_args() const override { return 2; } + /// @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 @@ -38,6 +65,70 @@ struct BackwardEulerFirstOrderTimeIntegrationRule { { return (1.0 / t.dt()) * (field_new - field_old); } + + /// @overload + FieldState corrected_value(const TimeInfo& t, const std::vector& states) const override + { + return value(t, states[0], states[1]); + } + + /// @overload + FieldState corrected_dot(const TimeInfo& t, const std::vector& states) const override + { + return dot(t, states[0], states[1]); + } + + /// @overload + FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector& states) const override + { + SLIC_ERROR("BackwardEulerFirstOrderTimeIntegrationRule does not support second derivatives."); + return states[0]; + } +}; + +/// @brief encodes rules for time discretizing first order odes where time derivatives are zero. +/// When solving f(u, t) = 0 +/// this class provides the current discrete approximation for u as a function of u^{n+1}. +class QuasiStaticRule : public TimeIntegrationRule { + public: + /// @brief Constructor + QuasiStaticRule() {} + + /// @brief get the number of states required by the rule + int num_args() const override { return 1; } + + /// @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 + { + return field_new; + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto dot(const TimeInfo& /*t*/, const T1& /*field_new*/) const + { + return zero{}; + } + + /// @overload + FieldState corrected_value(const TimeInfo& t, const std::vector& states) const override + { + return value(t, states[0]); + } + + /// @overload + FieldState corrected_dot(const TimeInfo& /*t*/, const std::vector& states) const override + { + return zeroCopy(states[0]); + } + + /// @overload + FieldState corrected_ddot(const TimeInfo& /*t*/, const std::vector& states) const override + { + SLIC_ERROR("QuasiStaticRule does not support second derivatives."); + return states[0]; + } }; /// Alternative name for Backward Euler which makes sense when restricting what are typically second order odes, @@ -50,10 +141,14 @@ using QuasiStaticFirstOrderTimeIntegrationRule = BackwardEulerFirstOrderTimeInte /// 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 ImplicitNewmarkSecondOrderTimeIntegrationRule { +struct ImplicitNewmarkSecondOrderTimeIntegrationRule : public TimeIntegrationRule { + public: /// @brief Constructor ImplicitNewmarkSecondOrderTimeIntegrationRule() {} + /// @brief get the number of states required by the rule + int num_args() const override { return 4; } + /// @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, @@ -81,6 +176,82 @@ struct ImplicitNewmarkSecondOrderTimeIntegrationRule { auto dt = t.dt(); return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; } + + /// @overload + FieldState corrected_value(const TimeInfo& t, const std::vector& states) const override + { + return value(t, states[0], states[1], states[2], states[3]); + } + + /// @overload + FieldState corrected_dot(const TimeInfo& t, const std::vector& states) const override + { + return dot(t, states[0], states[1], states[2], states[3]); + } + + /// @overload + FieldState corrected_ddot(const TimeInfo& t, const std::vector& states) const override + { + return ddot(t, states[0], states[1], states[2], states[3]); + } +}; + +/// @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 QuasiStaticSecondOrderTimeIntegrationRule : public TimeIntegrationRule { + public: + /// @brief Constructor + QuasiStaticSecondOrderTimeIntegrationRule() {} + + /// @brief get the number of states required by the rule + int num_args() const override { return 4; } + + /// @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 dot([[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 (1.0 / t.dt()) * (field_new - field_old); + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto ddot([[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 accel_old; + } + + /// @overload + FieldState corrected_value(const TimeInfo& t, const std::vector& states) const override + { + return value(t, states[0], states[1], states[2], states[3]); + } + + /// @overload + FieldState corrected_dot(const TimeInfo& t, const std::vector& states) const override + { + return dot(t, states[0], states[1], states[2], states[3]); + } + + /// @overload + FieldState corrected_ddot(const TimeInfo& t, const std::vector& states) const override + { + return ddot(t, states[0], states[1], states[2], states[3]); + } }; } // namespace smith diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 977c4ddaf0..b6ab1df5d5 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -35,6 +35,9 @@ class NewtonSolver : public mfem::NewtonSolver { /// reconstructed smith print level mutable size_t print_level = 0; + /// Tracks if grad was monolithicized and needs deletion + mutable bool grad_monolithic = false; + public: /// constructor NewtonSolver(const NonlinearSolverOptions& nonlinear_opts) : nonlinear_options(nonlinear_opts) {} @@ -47,6 +50,12 @@ class NewtonSolver : public mfem::NewtonSolver { } #endif + /// destructor + virtual ~NewtonSolver() + { + if (grad_monolithic) delete grad; + } + /// Evaluate the residual, put in rOut and return its norm. double evaluateNorm(const mfem::Vector& x, mfem::Vector& rOut) const { @@ -65,10 +74,18 @@ class NewtonSolver : public mfem::NewtonSolver { void assembleJacobian(const mfem::Vector& x) const { SMITH_MARK_FUNCTION; + if (grad_monolithic) { + delete grad; + grad = nullptr; + grad_monolithic = false; + } grad = &oper->GetGradient(x); if (nonlinear_options.force_monolithic) { auto* grad_blocked = dynamic_cast(grad); - if (grad_blocked) grad = buildMonolithicMatrix(*grad_blocked).release(); + if (grad_blocked) { + grad = buildMonolithicMatrix(*grad_blocked).release(); + grad_monolithic = true; + } } } @@ -337,6 +354,9 @@ class TrustRegion : public mfem::NewtonSolver { /// reconstructed smith print level mutable size_t print_level = 0; + /// Tracks if grad was monolithicized and needs deletion + mutable bool grad_monolithic = false; + public: /// internal counter for hess-vecs mutable size_t num_hess_vecs = 0; @@ -358,6 +378,12 @@ class TrustRegion : public mfem::NewtonSolver { } #endif + /// destructor + virtual ~TrustRegion() + { + if (grad_monolithic) delete grad; + } + /// finds tau s.t. (z + tau*d)^2 = trSize^2 void projectToBoundaryWithCoefs(mfem::Vector& z, const mfem::Vector& d, double delta, double zz, double zd, double dd) const @@ -597,10 +623,18 @@ class TrustRegion : public mfem::NewtonSolver { { SMITH_MARK_FUNCTION; ++num_jacobian_assembles; + if (grad_monolithic) { + delete grad; + grad = nullptr; + grad_monolithic = false; + } grad = &oper->GetGradient(x); if (nonlinear_options.force_monolithic) { auto* grad_blocked = dynamic_cast(grad); - if (grad_blocked) grad = buildMonolithicMatrix(*grad_blocked).release(); + if (grad_blocked) { + grad = buildMonolithicMatrix(*grad_blocked).release(); + grad_monolithic = true; + } } } @@ -935,6 +969,21 @@ void SuperLUSolver::Mult(const mfem::Vector& input, mfem::Vector& output) const superlu_solver_.Mult(input, output); } +/** + * @brief Build a monolithic HypreParMatrix from a BlockOperator. + * + * PERFORMANCE NOTE: This function creates a NEW monolithic matrix by copying data from + * the block structure. This incurs a performance overhead: + * - Memory: Allocates new matrix storage + * - Time: Copies all block data into monolithic format + * + * This is necessary when using direct solvers (SuperLU, Strumpack) that require + * monolithic matrices. For iterative solvers, the BlockOperator can be used directly + * without this copy overhead. + * + * @param block_operator The block operator to convert. + * @return Unique pointer to the new monolithic HypreParMatrix. + */ std::unique_ptr buildMonolithicMatrix(const mfem::BlockOperator& block_operator) { int row_blocks = block_operator.NumRowBlocks(); @@ -959,7 +1008,8 @@ std::unique_ptr buildMonolithicMatrix(const mfem::BlockOpe } } - // Note that MFEM passes ownership of this matrix to the caller + // Note that MFEM passes ownership of this matrix to the caller. + // MFEM creates a new monolithic matrix (not a view), so this is a COPY operation. return std::unique_ptr(mfem::HypreParMatrixFromBlocks(hypre_blocks)); } diff --git a/src/smith/numerics/equation_solver.hpp b/src/smith/numerics/equation_solver.hpp index 33dcc42621..374b23f83e 100644 --- a/src/smith/numerics/equation_solver.hpp +++ b/src/smith/numerics/equation_solver.hpp @@ -254,13 +254,9 @@ class StrumpackSolver : public mfem::Solver { /** * @brief Function for building a monolithic parallel Hypre matrix from a block system of smaller Hypre matrices - * - * @param block_operator The block system of HypreParMatrices * @return The assembled monolithic HypreParMatrix - * - * @pre @a block_operator must have assembled HypreParMatrices for its sub-blocks */ -std::unique_ptr buildMonolithicMatrix(const mfem::BlockOperator& block_operator); +std::unique_ptr buildMonolithicMatrix(const mfem::BlockOperator&); /** * @brief Build a nonlinear solver using the nonlinear option struct diff --git a/src/smith/physics/CMakeLists.txt b/src/smith/physics/CMakeLists.txt index 19066e61a6..fca58f635c 100644 --- a/src/smith/physics/CMakeLists.txt +++ b/src/smith/physics/CMakeLists.txt @@ -38,7 +38,6 @@ set(physics_headers weak_form.hpp functional_weak_form.hpp solid_weak_form.hpp - heat_transfer_weak_form.hpp scalar_objective.hpp functional_objective.hpp constraint.hpp diff --git a/src/smith/physics/functional_objective.hpp b/src/smith/physics/functional_objective.hpp index fc9b1c8cea..a77a026fef 100644 --- a/src/smith/physics/functional_objective.hpp +++ b/src/smith/physics/functional_objective.hpp @@ -75,6 +75,21 @@ class FunctionalObjective, std::integer_ mesh_->domain(body_name)); } + /** + * @brief register a custom boundary integral calculation as part of the residual + * + * @tparam active_parameters a list of indices, describing which parameters to pass to the q-function + * @param boundary_name string specifying the boundary to integrate over + * @param qfunction a callable that returns a tuple of body-force and stress + */ + template + void addBoundaryIntegral(DependsOn, std::string boundary_name, + const FuncOfTimeSpaceAndParams& qfunction) + { + objective_->AddBoundaryIntegral(smith::Dimension{}, smith::DependsOn{}, + qfunction, mesh_->domain(boundary_name)); + } + /// @overload virtual double evaluate(TimeInfo time_info, ConstFieldPtr shape_disp, const std::vector& fields) const override diff --git a/src/smith/physics/functional_weak_form.hpp b/src/smith/physics/functional_weak_form.hpp index af7723f5b9..8044961ab8 100644 --- a/src/smith/physics/functional_weak_form.hpp +++ b/src/smith/physics/functional_weak_form.hpp @@ -61,7 +61,14 @@ class FunctionalWeakForm, axom::fmt::format("{} parameter spaces given in the template argument but {} input mfem spaces were supplied.", sizeof...(InputSpaces), input_mfem_spaces.size())); + // Validate output space + validateSpace(output_mfem_space, "output"); + + // Validate input spaces if constexpr (sizeof...(InputSpaces) > 0) { + // Validate each input space using a helper + validateInputSpaces<0>(input_mfem_spaces); + for_constexpr([&](auto i) { trial_spaces[i] = input_mfem_spaces[i]; }); for_constexpr( [&](auto i) { vector_residual_trial_spaces[i + 1] = input_mfem_spaces[i]; }); @@ -288,6 +295,7 @@ class FunctionalWeakForm, mfem::Vector residual(TimeInfo time_info, ConstFieldPtr shape_disp, const std::vector& fields, [[maybe_unused]] const std::vector& quad_fields = {}) const override { + validateFields(fields, "residual"); dt_ = time_info.dt(); cycle_ = time_info.cycle(); auto ret = (*weak_form_)(time_info.time(), *shape_disp, *fields[input_indices]...); @@ -300,6 +308,7 @@ class FunctionalWeakForm, const std::vector& jacobian_weights, [[maybe_unused]] const std::vector& quad_fields = {}) const override { + validateFields(fields, "jacobian"); dt_ = time_info.dt(); cycle_ = time_info.cycle(); @@ -338,6 +347,7 @@ class FunctionalWeakForm, [[maybe_unused]] const std::vector& v_quad_fields, DualFieldPtr jvp_reaction) const override { + validateFields(fields, "jvp"); SLIC_ERROR_IF(v_fields.size() != fields.size(), "Invalid number of field sensitivities relative to the number of fields"); @@ -363,11 +373,13 @@ class FunctionalWeakForm, DualFieldPtr vjp_shape_disp_sensitivity, const std::vector& vjp_sensitivities, [[maybe_unused]] const std::vector& vjp_quad_field_sensitivities) const override { + validateFields(fields, "vjp"); SLIC_ERROR_IF(vjp_sensitivities.size() != fields.size(), "Invalid number of field sensitivities relative to the number of fields"); dt_ = time_info.dt(); cycle_ = time_info.cycle(); + auto vecJacs = vectorJacobianFunctions(std::make_integer_sequence{}, time_info.time(), shape_disp, v_field, fields); { @@ -399,6 +411,77 @@ class FunctionalWeakForm, } protected: + /// @brief Helper to validate input spaces recursively (for constructor) + template + void validateInputSpaces(const SpacesT& input_mfem_spaces) const + { + if constexpr (I < sizeof...(InputSpaces)) { + using Space = typename std::tuple_element>::type; + validateSpace(*input_mfem_spaces[I], axom::fmt::format("input[{}]", I)); + validateInputSpaces(input_mfem_spaces); + } + } + + /// @brief Helper to validate fields passed to residual/jacobian/vjp/jvp + template + void validateFieldsRecursive(const std::vector& fields, const std::string& method_name) const + { + if constexpr (I < sizeof...(InputSpaces)) { + using Space = typename std::tuple_element>::type; + validateSpace(fields[I]->space(), + axom::fmt::format("{}(): field[{}] ('{}')", method_name, I, fields[I]->name())); + validateFieldsRecursive(fields, method_name); + } + } + + /// @brief Validate that fields vector matches the templated input spaces + void validateFields(const std::vector& fields, const std::string& method_name) const + { + SLIC_ERROR_ROOT_IF(fields.size() != sizeof...(InputSpaces), + axom::fmt::format("{}(): fields.size()={} but weak form expects {} InputSpaces", method_name, + fields.size(), sizeof...(InputSpaces))); + + if constexpr (sizeof...(InputSpaces) > 0) { + validateFieldsRecursive<0>(fields, method_name); + } + } + + /// @brief Validate that an mfem space matches the template space parameters + template + static void validateSpace(const mfem::ParFiniteElementSpace& mfem_space, const std::string& space_name) + { + const auto* fec = mfem_space.FEColl(); + + // Note: We check using the name string rather than dynamic_cast because MFEM + // has multiple FECollection types that represent H1 and L2 spaces + if constexpr (Space::family == Family::H1) { + std::string fec_name = fec->Name(); + bool is_h1 = + (fec_name.find("H1") != std::string::npos || fec_name.find("ND_") != std::string::npos || // Nedelec elements + fec_name.find("Linear") != std::string::npos); + SLIC_ERROR_ROOT_IF(!is_h1, axom::fmt::format("Space '{}': Template specifies H1 family but mfem space uses '{}'", + space_name, fec_name)); + } else if constexpr (Space::family == Family::L2) { + std::string fec_name = fec->Name(); + bool is_l2 = (fec_name.find("L2") != std::string::npos || fec_name.find("DG") != std::string::npos || + fec_name.find("Const") != std::string::npos); + SLIC_ERROR_ROOT_IF( + !is_l2, axom::fmt::format("Space '{}': Template specifies L2/DG family but mfem space uses '{}'", space_name, + fec_name)); + } + + // Validate order + SLIC_ERROR_ROOT_IF(fec->GetOrder() != Space::order, + axom::fmt::format("Space '{}': Template specifies order {} but mfem space has order {}", + space_name, Space::order, fec->GetOrder())); + + // Validate vector dimension + SLIC_ERROR_ROOT_IF( + mfem_space.GetVDim() != Space::components, + axom::fmt::format("Space '{}': Template specifies {} components but mfem space has {} components (VDim)", + space_name, Space::components, mfem_space.GetVDim())); + } + /// @brief Utility to get array of jacobian functions, one for each input field in fs template auto jacobianFunctions(std::integer_sequence, double time, ConstFieldPtr shape_disp, diff --git a/src/smith/physics/heat_transfer_weak_form.hpp b/src/smith/physics/heat_transfer_weak_form.hpp deleted file mode 100644 index 2851ab4c19..0000000000 --- a/src/smith/physics/heat_transfer_weak_form.hpp +++ /dev/null @@ -1,173 +0,0 @@ -// 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 heat_transfer_weak_form.hpp - * - * @brief Implements the WeakForm interface for heat transfer physics. - * Derives from FunctionalWeakForm. - */ - -#pragma once - -#include "smith/physics/functional_weak_form.hpp" - -namespace smith { - -template > -class HeatTransferWeakForm; - -/** - * @brief The weak form for heat transfer - * - * This uses smith::functional to compute the heat transfer residuals and tangent - * stiffness matrices. - * - * @tparam order The order of the discretization of the temperature and temperature rate - * @tparam dim The spatial dimension of the mesh - */ -template -class HeatTransferWeakForm> - : public FunctionalWeakForm, Parameters, H1, InputSpaces...>> { - public: - /// @brief typedef for underlying functional type with templates - using BaseWeakFormT = FunctionalWeakForm, Parameters, H1, InputSpaces...>>; - - // /// @brief a container holding quadrature point data of the specified type - // /// @tparam T the type of data to store at each quadrature point - // template - // using qdata_type = std::shared_ptr>; - - /// @brief temperature, temperature rate - static constexpr int NUM_STATE_VARS = 2; - - /// @brief enumeration of the required heat transfer states - enum STATE - { - TEMPERATURE, - TEMPERATURE_RATE, - NUM_STATES - }; - - /** - * @brief Construct a new HeatTransferWeakForm object - * - * @param physics_name A name for the physics module instance - * @param mesh The Smith Mesh - * @param test_space Test space - * @param parameter_fe_spaces Vector of parameters spaces - */ - HeatTransferWeakForm(std::string physics_name, std::shared_ptr mesh, - const mfem::ParFiniteElementSpace& test_space, - std::vector parameter_fe_spaces = {}) - : BaseWeakFormT(physics_name, mesh, test_space, constructAllSpaces(test_space, parameter_fe_spaces)) - { - } - - /** - * @brief Set the thermal material model for the physics module - * - * @tparam MaterialType The thermal material type - * @param body_name string name for a registered body Domain on the mesh - * @param material A material that provides a function to evaluate heat capacity and thermal flux - * @pre material must be a object that can be called with the following arguments: - * 1. `tensor x` the spatial position of the material evaluation call - * 2. `T temperature` the current temperature at the quadrature point - * 3. `tensor` the spatial gradient of the temperature at the quadrature point - * 4. `tuple{value, derivative}`, a tuple of values and derivatives for each parameter field - * specified in the `DependsOn<...>` argument. - * - * @note The actual types of these arguments passed will be `double`, `tensor` or tuples thereof - * when doing direct evaluation. When differentiating with respect to one of the inputs, its stored - * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, - * 3>`) - * - * @pre MaterialType must return a smith::tuple of volumetric heat capacity and thermal flux when operator() is called - * with the arguments listed above. - * - */ - template - void setMaterial(DependsOn, std::string body_name, const MaterialType& material) - { - ThermalMaterialFunctor material_functor(material); - BaseWeakFormT::weak_form_->AddDomainIntegral(Dimension{}, - DependsOn<0, 1, active_parameters + NUM_STATE_VARS...>{}, - std::move(material_functor), BaseWeakFormT::mesh_->domain(body_name)); - BaseWeakFormT::v_dot_weak_form_residual_->AddDomainIntegral( - Dimension{}, DependsOn<0, 1, 2, active_parameters + 1 + NUM_STATE_VARS...>{}, - [material_functor](double t, auto X, auto V, auto... params) { - auto flux = material_functor(t, X, params...); - return smith::inner(get(V), get(flux)) + - smith::inner(get(V), get(flux)); - }, - BaseWeakFormT::mesh_->domain(body_name)); - } - - /// @overload - template - void setMaterial(std::string body_name, const MaterialType& material) - { - setMaterial(DependsOn<>{}, body_name, material); - } - - protected: - /// @brief For use in the constructor, combined the correct number of state spaces (disp,velo,accel) with the vector - /// of parameters - /// @param state_space H1 temperature space - /// @param spaces parameter spaces - /// @return - std::vector constructAllSpaces( - const mfem::ParFiniteElementSpace& state_space, const std::vector& spaces) - { - std::vector all_spaces{&state_space, &state_space}; - for (auto& s : spaces) { - all_spaces.push_back(s); - } - return all_spaces; - } - - /** - * @brief Functor representing a thermal material's heat capacity and flux. - */ - template - struct ThermalMaterialFunctor { - /** - * @brief Construct a ThermalMaterialIntegrand functor with material model of type `MaterialType`. - * @param[in] material A functor representing the material model. Should be a functor, or a class/struct with - * public operator() method. Must NOT be a generic lambda, or Smith will not compile due to static asserts below. - */ - ThermalMaterialFunctor(MaterialType material) : material_(material) {} - - /// Material model - MaterialType material_; - - /** - * @brief Thermal material response call - * - * @tparam X Spatial position type - * @tparam Temperature temperature - * @tparam dT_dt temperature rate - * @tparam Params variadic parameters for call - * @param[in] x spatial position - * @param[in] temperature temperature - * @param[in] dtemp_dt temperature rate - * @param[in] params parameter pack - * @return The calculated material response (tuple of volumetric heat capacity and thermal flux) for a linear - * isotropic material - */ - template - auto operator()(double /*time*/, X x, Temperature temperature, dT_dt dtemp_dt, Params... params) const - { - // Get the value and the gradient from the input tuple - auto [u, du_dX] = temperature; - auto du_dt = get(dtemp_dt); - auto [heat_capacity, heat_flux] = material_(x, u, du_dX, params...); - return smith::tuple{heat_capacity * du_dt, -heat_flux}; - } - }; -}; - -} // namespace smith diff --git a/src/smith/physics/tests/CMakeLists.txt b/src/smith/physics/tests/CMakeLists.txt index 1834a5bd4d..066689fb9d 100644 --- a/src/smith/physics/tests/CMakeLists.txt +++ b/src/smith/physics/tests/CMakeLists.txt @@ -29,12 +29,9 @@ set(physics_serial_test_sources quasistatic_solid_adjoint.cpp finite_element_vector_set_over_domain.cpp solid_multi_material.cpp - test_solid_weak_form.cpp - test_heat_weak_form.cpp test_functional_weak_form.cpp thermomech_finite_diff.cpp thermomech_statics_patch.cpp - test_kinematic_objective.cpp ) if (TRIBOL_FOUND) diff --git a/src/smith/physics/tests/contact_solid_adjoint.cpp b/src/smith/physics/tests/contact_solid_adjoint.cpp index 39bcf3b985..e54f4eae49 100644 --- a/src/smith/physics/tests/contact_solid_adjoint.cpp +++ b/src/smith/physics/tests/contact_solid_adjoint.cpp @@ -286,8 +286,8 @@ TEST_F(ContactSensitivityFixture, ReactionShapeSensitivities) { auto solid_solver = createContactSolver(mesh, nonlinear_opts, dyn_opts, mat); auto [qoi_base, shape_sensitivity] = computeContactReactionQoiSensitivities(*solid_solver, mesh); - solid_solver->resetStates(); + FiniteElementState derivative_direction(shape_sensitivity.space(), "derivative_direction"); fillDirection(*solid_solver, derivative_direction); diff --git a/src/smith/physics/tests/test_heat_weak_form.cpp b/src/smith/physics/tests/test_heat_weak_form.cpp deleted file mode 100644 index 58dc3cf1ea..0000000000 --- a/src/smith/physics/tests/test_heat_weak_form.cpp +++ /dev/null @@ -1,205 +0,0 @@ -// Copyright 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 -#include -#include -#include - -#include "gtest/gtest.h" -#include "mpi.h" -#include "mfem.hpp" - -#include "smith/infrastructure/application_manager.hpp" -#include "smith/physics/materials/parameterized_thermal_material.hpp" -#include "smith/physics/mesh.hpp" -#include "smith/physics/state/state_manager.hpp" -#include "smith/physics/tests/physics_test_utils.hpp" -#include "smith/physics/functional_weak_form.hpp" -#include "smith/physics/heat_transfer_weak_form.hpp" -#include "smith/numerics/functional/finite_element.hpp" -#include "smith/numerics/functional/functional.hpp" -#include "smith/numerics/functional/tensor.hpp" -#include "smith/physics/common.hpp" -#include "smith/physics/field_types.hpp" -#include "smith/physics/state/finite_element_dual.hpp" -#include "smith/physics/state/finite_element_state.hpp" - -auto element_shape = mfem::Element::QUADRILATERAL; - -struct WeakFormFixture : public testing::Test { - WeakFormFixture() : time_info(time, dt) {} - - static constexpr int dim = 2; - static constexpr int order = 1; - - using ScalarSpace = smith::H1; - using ParamSpace = smith::L2; - - using ThermalMaterial = smith::heat_transfer::ParameterizedLinearIsotropicConductor; - - void SetUp() - { - MPI_Barrier(MPI_COMM_WORLD); - smith::StateManager::initialize(datastore, "heat_diffusion"); - - mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(10, 10, element_shape, true, 1.0, 1.0), - "this_mesh_name", 0, 0); - - smith::FiniteElementState temperature = smith::StateManager::newState(ScalarSpace{}, "temperature", mesh->tag()); - smith::FiniteElementState temperature_rate = - smith::StateManager::newState(ScalarSpace{}, "temperature_rate", mesh->tag()); - smith::FiniteElementState conductivity_offset = - smith::StateManager::newState(ParamSpace{}, "conductivity_offset", mesh->tag()); - - shape_disp = std::make_unique(mesh->newShapeDisplacement()); - shape_disp_dual = std::make_unique(mesh->newShapeDisplacementDual()); - - states = {temperature, temperature_rate}; - params = {conductivity_offset}; - - for (auto s : states) { - state_duals.push_back(smith::FiniteElementDual(s.space(), s.name() + "_dual")); - } - for (auto p : params) { - param_duals.push_back(smith::FiniteElementDual(p.space(), p.name() + "_dual")); - } - - state_tangents = states; - param_tangents = params; - - std::string physics_name = "heat"; - - auto heat_transfer_weak_form = std::make_shared( - physics_name, mesh, states[HeatWeakFormT::TEMPERATURE].space(), getSpaces(params)); - - ThermalMaterial mat(1.0, 1.0, 1.0); - heat_transfer_weak_form->setMaterial(smith::DependsOn<0>{}, mesh->entireBodyName(), mat); - - std::string source_name = "source"; - mesh->addDomainOfBoundaryElements(source_name, smith::by_attr({1, 2})); - - heat_transfer_weak_form->addBoundaryFlux(source_name, [](auto /* t */, auto /* x */, auto /* n */) { return 1.0; }); - - for (auto& s : state_tangents) { - pseudoRand(s); - } - for (auto& p : param_tangents) { - pseudoRand(p); - } - - // used to test that vjp acts via +=, add initial value to shape displacement dual - state_duals[HeatWeakFormT::TEMPERATURE] = 1.0; - - states[HeatWeakFormT::TEMPERATURE].setFromFieldFunction( - [](smith::tensor x) { return 0.1 * std::pow(std::pow(x[0], 2.0) + std::pow(x[1], 2.0), 0.5); }); - states[HeatWeakFormT::TEMPERATURE_RATE].setFromFieldFunction([](smith::tensor x) { - return 0.01 * std::pow(std::pow(x[0], 2.0) + std::pow(0.5 * x[1], 2.0), 0.5); - }); - params[0] = 0.5; - - weak_form = heat_transfer_weak_form; - } - - using HeatWeakFormT = smith::HeatTransferWeakForm>; - - const double time = 0.0; - const double dt = 1.0; - smith::TimeInfo time_info; - - axom::sidre::DataStore datastore; - std::shared_ptr mesh; - std::shared_ptr weak_form; - - std::unique_ptr shape_disp; - std::unique_ptr shape_disp_dual; - - std::vector states; - std::vector params; - - std::vector state_duals; - std::vector param_duals; - - std::vector state_tangents; - std::vector param_tangents; -}; - -TEST_F(WeakFormFixture, VjpConsistency) -{ - auto input_fields = getConstFieldPointers(states, params); - - smith::FiniteElementDual res_vector(states[HeatWeakFormT::TEMPERATURE].space(), "residual"); - res_vector = weak_form->residual(time_info, shape_disp.get(), input_fields); - ASSERT_NE(0.0, res_vector.Norml2()); - - auto jacobian_weights = [&](size_t i) { - std::vector tangents(input_fields.size()); - tangents[i] = 1.0; - return tangents; - }; - - // test vjp - smith::FiniteElementState v(res_vector.space(), "v"); - pseudoRand(v); - auto field_vjps = getFieldPointers(state_duals, param_duals); - - weak_form->vjp(time_info, shape_disp.get(), input_fields, {}, &v, shape_disp_dual.get(), field_vjps, {}); - - for (size_t i = 0; i < input_fields.size(); ++i) { - smith::FiniteElementState vjp = *input_fields[i]; - vjp = 0.0; - auto J = weak_form->jacobian(time_info, shape_disp.get(), input_fields, jacobian_weights(i)); - J->MultTranspose(v, vjp); - if (i == HeatWeakFormT::TEMPERATURE) vjp += 1.0; // make sure vjp uses += - EXPECT_NEAR(vjp.Norml2(), field_vjps[i]->Norml2(), 1e-12); - } -} - -TEST_F(WeakFormFixture, JvpConsistency) -{ - auto input_fields = getConstFieldPointers(states, params); - - smith::FiniteElementDual res_vector(states[HeatWeakFormT::TEMPERATURE].space(), "residual"); - res_vector = weak_form->residual(time_info, shape_disp.get(), input_fields); - ASSERT_NE(0.0, res_vector.Norml2()); - - auto jacobian_weights = [&](size_t i) { - std::vector tangents(input_fields.size()); - tangents[i] = 1.0; - return tangents; - }; - - auto selectStates = [&](size_t i) { - auto field_tangents = getConstFieldPointers(state_tangents, param_tangents); - for (size_t j = 0; j < field_tangents.size(); ++j) { - if (i != j) { - field_tangents[j] = nullptr; - } - } - return field_tangents; - }; - - smith::FiniteElementDual jvp_slow(states[HeatWeakFormT::TEMPERATURE].space(), "jvp_slow"); - smith::FiniteElementDual jvp(states[HeatWeakFormT::TEMPERATURE].space(), "jvp"); - jvp = 4.0; - - auto field_tangents = getFieldPointers(state_tangents, param_tangents); - - for (size_t i = 0; i < input_fields.size(); ++i) { - auto J = weak_form->jacobian(time_info, shape_disp.get(), input_fields, jacobian_weights(i)); - J->Mult(*field_tangents[i], jvp_slow); - weak_form->jvp(time_info, shape_disp.get(), input_fields, {}, nullptr, selectStates(i), {}, &jvp); - EXPECT_NEAR(jvp_slow.Norml2(), jvp.Norml2(), 1e-12); - } -} - -int main(int argc, char* argv[]) -{ - ::testing::InitGoogleTest(&argc, argv); - smith::ApplicationManager applicationManager(argc, argv); - return RUN_ALL_TESTS(); -}