diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 950392b794..4feef4fe35 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -31,6 +31,7 @@ The Smith project release numbers follow [Semantic Versioning](http://semver.org - Added axisymmetric solid mechanics materials and loads for 2D `(r, z)` meshes. - Added host-configs for CZMatrix and RZVector machines - Added CI pipelines for CZMatrix machine +- Added differentiable plastic system with internal variables represented by discontinuous galerkin nodal fields ### Removed diff --git a/src/smith/differentiable_numerics/CMakeLists.txt b/src/smith/differentiable_numerics/CMakeLists.txt index 71ade68cc9..2daffb1c5c 100644 --- a/src/smith/differentiable_numerics/CMakeLists.txt +++ b/src/smith/differentiable_numerics/CMakeLists.txt @@ -47,6 +47,7 @@ set(differentiable_numerics_headers combined_system.hpp system_base.hpp differentiable_test_utils.hpp + plastic_mechanics_system.hpp ) set(differentiable_numerics_depends diff --git a/src/smith/differentiable_numerics/differentiable_physics.cpp b/src/smith/differentiable_numerics/differentiable_physics.cpp index 4928db1a2b..58353be536 100644 --- a/src/smith/differentiable_numerics/differentiable_physics.cpp +++ b/src/smith/differentiable_numerics/differentiable_physics.cpp @@ -67,9 +67,8 @@ DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr mesh, std::sh reaction_names_.reserve(reaction_infos_.size()); for (size_t i = 0; i < reaction_infos_.size(); ++i) { - SLIC_ERROR_IF( - reaction_infos_[i].space == nullptr, - axom::fmt::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_)); + SLIC_ERROR_IF(reaction_infos_[i].space == nullptr, + std::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_)); reaction_names_.push_back(reaction_infos_[i].name); reaction_name_to_reaction_index_[reaction_infos_[i].name] = i; } @@ -135,7 +134,7 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction if (reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end()) { std::string available; for (auto& n : reaction_names_) available += n + " "; - SLIC_ERROR(axom::fmt::format( + SLIC_ERROR(std::format( "Could not find reaction named {0} in mesh with tag \"{1}\" to get. Available reactions (size {2}): {3}", reaction_name, mesh_->tag(), reaction_names_.size(), available)); } @@ -153,11 +152,10 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction FiniteElementDual DifferentiablePhysics::loadCheckpointedDual(const std::string& reaction_name, int cycle) { - SLIC_ERROR_IF( - cycle != cycle_, - axom::fmt::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for " - "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}", - cycle, cycle_)); + SLIC_ERROR_IF(cycle != cycle_, + std::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for " + "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}", + cycle, cycle_)); return dual(reaction_name); } @@ -218,9 +216,8 @@ void DifferentiablePhysics::setAdjointLoad( for (auto string_reaction_pair : string_to_reaction) { std::string field_name = string_reaction_pair.first; const smith::FiniteElementDual& reaction = string_reaction_pair.second; - SLIC_ERROR_IF( - state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), - axom::fmt::format("Could not find reaction named {0} in mesh with tag {1}", field_name, mesh_->tag())); + SLIC_ERROR_IF(state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), + std::format("Could not find reaction named {0} in mesh with tag {1}", field_name, mesh_->tag())); size_t state_index = state_name_to_field_index_.at(field_name); *field_states_[state_index].get_dual() += reaction; } @@ -232,10 +229,9 @@ void DifferentiablePhysics::setDualAdjointBcs( for (auto string_bc_pair : string_to_bc) { std::string reaction_name = string_bc_pair.first; const smith::FiniteElementState& reaction_adjoint_load = string_bc_pair.second; - SLIC_ERROR_IF( - reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(), - axom::fmt::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}", - reaction_name, mesh_->tag())); + SLIC_ERROR_IF(reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(), + std::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}", + reaction_name, mesh_->tag())); size_t reaction_index = reaction_name_to_reaction_index_.at(reaction_name); *reaction_states_[reaction_index].get_dual() += reaction_adjoint_load; } diff --git a/src/smith/differentiable_numerics/differentiable_test_utils.hpp b/src/smith/differentiable_numerics/differentiable_test_utils.hpp index 4ecf7e79d2..f93ec2d7a3 100644 --- a/src/smith/differentiable_numerics/differentiable_test_utils.hpp +++ b/src/smith/differentiable_numerics/differentiable_test_utils.hpp @@ -79,6 +79,7 @@ inline double checkGradWrt(const gretl::State& objective, smith::FieldSt std::vector grad_errors; auto [grad, grad_fd] = checkGradients(objective, input, *dual_vec, objectiveBase, graph, eps); + if (printmore) std::cout << "grad = " << grad << "\ngrad fd = " << grad_fd << std::endl; grad_errors.push_back(std::abs(grad - grad_fd)); for (size_t step = 1; step < num_fd_steps; ++step) { diff --git a/src/smith/differentiable_numerics/field_store.cpp b/src/smith/differentiable_numerics/field_store.cpp index 20d0869ade..04ff05cf4a 100644 --- a/src/smith/differentiable_numerics/field_store.cpp +++ b/src/smith/differentiable_numerics/field_store.cpp @@ -31,16 +31,6 @@ std::shared_ptr FieldStore::addBoundaryConditions(F return std::make_shared(mesh_->mfemParMesh(), field->space()); } -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 @@ -55,71 +45,43 @@ void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument void FieldStore::printMap() { - for (auto& keyval : weak_form_name_to_unknown_name_index_) { + for (const auto& keyval : weak_form_name_to_field_names_) { 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 << ", "; + for (size_t argument_index = 0; argument_index < keyval.second.size(); ++argument_index) { + std::cout << "arg " << keyval.second[argument_index] << " at " << argument_index << ", "; } std::cout << std::endl; } } -std::vector> FieldStore::indexMap(const std::vector& residual_names) const +BlockArgumentMap FieldStore::indexMap(const std::vector& residual_names) const { - // Build a local column space: each residual in the subsystem contributes one local column, - // corresponding to its "self" diagonal unknown. The self-unknown is preferably the residual's - // reaction (test) field if that field appears in the unknown-arg list for this weak form; - // otherwise fall back on the first unknown argument (handles cases like the cycle-zero - // acceleration solve, where the reaction field is a dependent/history field). std::map global_state_to_local_col; for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { const std::string& res_name = residual_names[res_i]; - size_t global_state_idx = invalid_block_index; - - std::string reaction_name; - for (const auto& kv : weak_form_to_test_field_) { - if (kv.first == res_name) { - reaction_name = kv.second; - break; - } - } - - // Check if the reaction field is one of the registered unknown args for this weak form. - bool reaction_is_unknown = false; - if (!reaction_name.empty() && weak_form_name_to_unknown_name_index_.count(res_name)) { - for (const auto& label : weak_form_name_to_unknown_name_index_.at(res_name)) { - if (label.field_name == reaction_name) { - reaction_is_unknown = true; - break; - } - } - } - - if (reaction_is_unknown) { - global_state_idx = to_states_index_.at(reaction_name); - } else { - const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); - SLIC_ERROR_IF(arg_info.empty(), - "Weak form '" << res_name << "' has no unknown arguments; cannot build index map."); - global_state_idx = to_states_index_.at(arg_info.front().field_name); - } - global_state_to_local_col[global_state_idx] = res_i; + size_t global_state_idx = to_states_index_.at(getWeakFormFieldNames(res_name).unknown); + auto [owner, inserted] = global_state_to_local_col.emplace(global_state_idx, res_i); + SLIC_ERROR_IF(!inserted, std::format("Weak forms '{}' and '{}' both own state index {}", + residual_names[owner->second], res_name, global_state_idx)); } - std::vector> block_indices(residual_names.size()); + BlockArgumentMap block_indices(residual_names.size(), + std::vector(residual_names.size(), BlockArgumentIndices{})); 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(residual_names.size(), 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); + const auto& field_names = weak_form_name_to_field_names_.at(res_name); + size_t state_slot = 0; + for (const auto& field_name : field_names) { + if (!to_states_index_.count(field_name)) { + continue; + } - for (const auto& field_name_and_arg_index : arg_info) { - size_t global_state_index = to_states_index_.at(field_name_and_arg_index.field_name); - auto it = global_state_to_local_col.find(global_state_index); - if (it != global_state_to_local_col.end()) { - res_indices[it->second] = field_name_and_arg_index.field_index_in_residual; + size_t global_state_index = to_states_index_.at(field_name); + auto column = global_state_to_local_col.find(global_state_index); + if (column != global_state_to_local_col.end()) { + block_indices[res_i][column->second].push_back(state_slot); } - // else: field belongs to a different subsystem; treat as fixed input here. + ++state_slot; } } @@ -132,7 +94,7 @@ std::vector FieldStore::getBoundaryConditionMan std::vector field_names; field_names.reserve(weak_form_names.size()); for (const auto& wf_name : weak_form_names) { - field_names.push_back(getWeakFormReaction(wf_name)); + field_names.push_back(getWeakFormFieldNames(wf_name).unknown); } return getBoundaryConditionManagersForFields(field_names); } @@ -277,9 +239,9 @@ 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 addReaction()?", - weak_form_name)); + SLIC_ERROR_ROOT_IF( + weak_form_name_to_field_names_.count(weak_form_name) == 0, + std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name)); auto field_names = weak_form_name_to_field_names_.at(weak_form_name); std::vector fields_for_residual; @@ -287,7 +249,7 @@ std::vector FieldStore::getStates(const std::string& weak_form_name) // 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)); + std::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 @@ -303,9 +265,9 @@ std::vector FieldStore::getStatesFromVectors(const std::string& weak 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 addReaction()?", - weak_form_name)); + SLIC_ERROR_ROOT_IF( + weak_form_name_to_field_names_.count(weak_form_name) == 0, + std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name)); auto field_names = weak_form_name_to_field_names_.at(weak_form_name); std::vector fields_for_residual; @@ -313,21 +275,21 @@ std::vector FieldStore::getStatesFromVectors(const std::string& weak // 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)); + SLIC_ERROR_ROOT_IF( + idx >= state_fields.size(), + std::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)); + std::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)); + SLIC_ERROR_ROOT(std::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore", name, + weak_form_name)); } } return fields_for_residual; @@ -359,26 +321,27 @@ void FieldStore::markWeakFormInternal(const std::string& weak_form_name) internal_weak_forms_.insert(weak_form_name); } -void FieldStore::addWeakFormReaction(std::string weak_form_name, std::string field_name) +void FieldStore::addWeakFormFieldNames(std::string weak_form_name, UnknownAndTestFieldNames field_names) { - for (auto& kv : weak_form_to_test_field_) { + for (auto& kv : weak_form_field_names_) { if (kv.first == weak_form_name) { - kv.second = field_name; + kv.second = std::move(field_names); return; } } - weak_form_to_test_field_.push_back({weak_form_name, field_name}); + weak_form_field_names_.push_back({std::move(weak_form_name), std::move(field_names)}); } -std::string FieldStore::getWeakFormReaction(const std::string& weak_form_name) const +const UnknownAndTestFieldNames& FieldStore::getWeakFormFieldNames(const std::string& weak_form_name) const { - for (const auto& kv : weak_form_to_test_field_) { + for (const auto& kv : weak_form_field_names_) { if (kv.first == weak_form_name) { return kv.second; } } - SLIC_ERROR("Reaction field not found for weak form " << weak_form_name); - return ""; + SLIC_ERROR("Field names not found for weak form " << weak_form_name); + static const UnknownAndTestFieldNames empty; + return empty; } const std::vector& FieldStore::getParameterFields() const { return params_; } @@ -405,12 +368,12 @@ std::vector FieldStore::getOutputFieldStates() const std::vector FieldStore::getReactionInfos() const { std::vector infos; - for (const auto& kv : weak_form_to_test_field_) { + for (const auto& kv : weak_form_field_names_) { const std::string& weak_form_name = kv.first; if (internal_weak_forms_.count(weak_form_name)) { continue; } - const std::string& field_name = kv.second; + const std::string& field_name = kv.second.test; infos.push_back({weak_form_name, &getField(field_name).get()->space()}); } return infos; diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp index d85b3a6338..8dcd9edada 100644 --- a/src/smith/differentiable_numerics/field_store.hpp +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -7,6 +7,7 @@ #pragma once #include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "smith/differentiable_numerics/time_integration_rule.hpp" #include "smith/physics/functional_weak_form.hpp" #include "smith/physics/mesh.hpp" @@ -30,16 +31,14 @@ struct ReactionInfo { const mfem::ParFiniteElementSpace* space = nullptr; ///< The finite element space of the dual field. }; +/// Names identifying a weak form's test field and solver-owned unknown. +struct UnknownAndTestFieldNames { + std::string unknown; ///< Field owned and solved by this residual row. + std::string test; ///< Field defining residual test space and reaction output. +}; + /** - * @brief Representation of a field type with a name and a flag indicating whether it is an - * active Jacobian unknown in the current weak-form context. - * - * @c is_unknown is intentionally a per-instance flag, not a global property. The same field - * may be a Jacobian variable in one weak form (e.g. displacement in the main solid solve) and - * a fixed input in another (e.g. displacement in the stress projection). Callers control this - * by passing the same @c FieldType object (set by @c addIndependent) or a plain copy with the - * default @c is_unknown = false. - * + * @brief Representation of a named field type. * @tparam Space The finite element space type. * @tparam Time The time integration type (unused by default). */ @@ -47,14 +46,9 @@ template struct FieldType { using space_type = Space; ///< The finite element space type. - /** - * @brief Construct a new FieldType object. - * @param n Name of the field. - * @param is_unknown_ Whether this field is a Jacobian unknown in the current context. - */ - FieldType(std::string n, bool is_unknown_ = false) : name(n), is_unknown(is_unknown_) {} - std::string name; ///< Name of the field. - bool is_unknown = false; ///< True if this field is a Jacobian variable in the current weak-form context. + /// Construct a field type with the given name. + FieldType(std::string n) : name(n) {} + std::string name; ///< Name of the field. }; /** @@ -115,9 +109,9 @@ struct FieldStore { // Verify the space matches by checking vdim (== Space::components). auto& existing = params_[to_params_index_.at(type.name)]; SLIC_ERROR_ROOT_IF(existing.get()->space().GetVDim() != Space::components, - axom::fmt::format("Parameter '{}' re-registered with a different space " - "(existing vdim={}, new vdim={})", - type.name, existing.get()->space().GetVDim(), Space::components)); + std::format("Parameter '{}' re-registered with a different space " + "(existing vdim={}, new vdim={})", + type.name, existing.get()->space().GetVDim(), Space::components)); return; } to_params_index_[type.name] = params_.size(); @@ -127,13 +121,10 @@ struct FieldStore { /** * @brief Add an independent field (a solver unknown) to the store. * - * Registers the field as an unknown by setting @c type.is_unknown = true, so the same - * @c FieldType object can later be passed to @c createSpaces to mark this argument - * as an active Jacobian variable. Also creates a boundary-condition slot keyed by field - * name that callers can populate after this call returns. + * Also creates a boundary-condition slot keyed by field name that callers can populate. * * @tparam Space The finite element space type. - * @param type The field type specification; @c type.is_unknown is set to @c true on return. + * @param type The field type specification. * @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. @@ -143,7 +134,6 @@ struct FieldStore { std::shared_ptr time_rule) { type.name = prefix(type.name); - type.is_unknown = true; to_states_index_[type.name] = states_.size(); FieldState new_field = smith::createFieldState(*graph_, Space{}, type.name, mesh_->tag()); states_.push_back(new_field); @@ -170,12 +160,8 @@ struct FieldStore { * 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 and @c is_unknown is @c false. 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. To make it the Jacobian variable for a specific weak form (e.g. - * the cycle-zero acceleration solve), copy the returned object and set @c is_unknown = true - * before passing to @c createSpaces. + * Returns a descriptor for the newly registered field so callers can pass it directly to + * @c createSpaces when assembling the weak-form argument list. * * @tparam Space The finite element space type (must match the independent field). * @param independent_field The @c FieldType of the independent (predicted) field. @@ -220,14 +206,6 @@ struct FieldStore { 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. @@ -236,17 +214,8 @@ struct FieldStore { */ 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); + /// Register test-space and solver-owned field names for a weak form. + void addWeakFormFieldNames(std::string weak_form_name, UnknownAndTestFieldNames field_names); /** * @brief Mark a weak form as internal so it is excluded from getReactionInfos(). @@ -256,56 +225,57 @@ struct FieldStore { */ void markWeakFormInternal(const std::string& weak_form_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; + /// Get explicit test-space and solver-owned field names for a weak form. + const UnknownAndTestFieldNames& getWeakFormFieldNames(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. If @c type.is_unknown is @c true, the field is also registered as an active - * Jacobian unknown for this weak form. + * 1. Registers explicit test-space and solver-owned field names. + * 2. Registers every @c FieldType in @p types as an ordered input argument. * 3. Returns the ordered vector of finite element spaces. * - * A field may have @c is_unknown = true in one weak form and @c false in another (e.g. - * displacement is a Jacobian variable in the main solid solve but a fixed input in the stress - * projection). Callers control this by passing the @c FieldType returned from @c addIndependent - * (has @c is_unknown = true) or a plain copy constructed from the field name (has @c is_unknown - * = false). Similarly, a dependent field can be made the Jacobian variable for a specific weak - * form (e.g. acceleration in the cycle-zero solve) by copying the returned @c FieldType and - * setting @c is_unknown = true. - * * @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 field_names Names of the test-space and solver-owned fields. * @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, + UnknownAndTestFieldNames field_names, FieldTypes... types) { - addWeakFormReaction(weak_form_name, reaction_field_name); + addWeakFormFieldNames(weak_form_name, field_names); std::vector spaces; size_t arg_num = 0; + bool unknown_found = false; auto register_field = [&](auto type) { spaces.push_back(&getField(type.name).get()->space()); addWeakFormArg(weak_form_name, type.name, arg_num); - if (type.is_unknown) { - addWeakFormUnknownArg(weak_form_name, type.name, arg_num); - } + unknown_found = unknown_found || type.name == field_names.unknown; ++arg_num; }; (register_field(types), ...); + SLIC_ERROR_IF(!unknown_found, "Unknown field '" << field_names.unknown << "' is not an argument of weak form '" + << weak_form_name << "'"); return spaces; } + /** + * @brief Register input fields when test and unknown field names match. + * @param weak_form_name Name of the weak form being constructed. + * @param field_name Shared test-space and solver-owned field name. + * @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& field_name, FieldTypes... types) + { + return createSpaces(weak_form_name, {.unknown = field_name, .test = field_name}, types...); + } + /** * @brief Mapping between primary and history/derivative fields for time integration. */ @@ -332,18 +302,18 @@ struct FieldStore { /** * @brief Generate an index map for the residuals. * @param residual_names Names of the residuals. - * @return std::vector> The index map. + * @return Slot-list matrix mapping residual rows and solved columns. */ - std::vector> indexMap(const std::vector& residual_names) const; + BlockArgumentMap indexMap(const std::vector& residual_names) const; /** * @brief Get the boundary condition managers for the given weak forms, one per residual row. * - * For each weak form in @p weak_form_names the reaction (test) field name is looked up. The + * For each weak form in @p weak_form_names the solver-owned field name is looked up. The * returned manager is selected by consulting the registered @c TimeIntegrationMapping s: - * - reaction = primary or history slot -> value-level BC manager - * - reaction = second-derivative (ddot) slot -> second-derivative BC manager - * - reaction has its own DBC entry not tied to a mapping -> that DBC's value manager + * - unknown = primary or history slot -> value-level BC manager + * - unknown = second-derivative (ddot) slot -> second-derivative BC manager + * - unknown has its own DBC entry not tied to a mapping -> that DBC's value manager * - otherwise -> @c nullptr (solver skips null entries) * * The second-derivative manager is rebuilt on each call, so late value-BC additions are reflected. @@ -484,19 +454,11 @@ struct FieldStore { /// Boundary conditions keyed by primary unknown field name. std::map> 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::vector> weak_form_to_test_field_; + std::vector> weak_form_field_names_; std::set internal_weak_forms_; ///< weak forms excluded from getReactionInfos() (subsystem-internal) std::vector, TimeIntegrationMapping>> time_integration_rules_; @@ -509,13 +471,15 @@ struct 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) +template +auto createWeakForm(std::string name, FieldType test_type, FieldType unknown_type, + FieldStore& field_store, FieldType... field_types) { - return std::make_shared>>( + return std::make_shared< + FunctionalWeakForm>>( name, field_store.getMesh(), field_store.getField(test_type.name).get()->space(), - field_store.createSpaces(name, test_type.name, field_types...)); + field_store.createSpaces(name, {.unknown = unknown_type.name, .test = test_type.name}, unknown_type, + field_types...)); } /** diff --git a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp index 7979726a51..57f4d75f04 100644 --- a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp @@ -25,8 +25,8 @@ MultiphysicsTimeIntegrator::MultiphysicsTimeIntegrator(std::shared_ptrweak_forms.size(); ++i) { const std::string wf_name = system_->weak_forms[i]->name(); - const std::string reaction_name = system_->field_store->getWeakFormReaction(wf_name); - main_unknown_name_to_local_idx_[reaction_name] = i; + const std::string unknown_name = system_->field_store->getWeakFormFieldNames(wf_name).unknown; + main_unknown_name_to_local_idx_[unknown_name] = i; } } @@ -75,7 +75,7 @@ std::pair, std::vector> MultiphysicsTimeI for (size_t i = 0; i < cz_sys->weak_forms.size(); ++i) { const std::string result_field_name = cz_sys->solve_result_field_names.empty() - ? system_->field_store->getWeakFormReaction(cz_sys->weak_forms[i]->name()) + ? system_->field_store->getWeakFormFieldNames(cz_sys->weak_forms[i]->name()).unknown : cz_sys->solve_result_field_names[i]; size_t result_field_state_idx = system_->field_store->getFieldIndex(result_field_name); current_states[result_field_state_idx] = cycle_zero_unknowns[i]; @@ -122,8 +122,8 @@ std::pair, std::vector> MultiphysicsTimeI for (const auto& ps : post_solve_systems_) { auto ps_unknowns = ps->solve(time_info); for (size_t i = 0; i < ps->weak_forms.size(); ++i) { - const std::string reaction_name = ps->field_store->getWeakFormReaction(ps->weak_forms[i]->name()); - size_t u_idx = ps->field_store->getFieldIndex(reaction_name); + const std::string unknown_name = ps->field_store->getWeakFormFieldNames(ps->weak_forms[i]->name()).unknown; + size_t u_idx = ps->field_store->getFieldIndex(unknown_name); ps->field_store->setField(u_idx, ps_unknowns[i]); } } diff --git a/src/smith/differentiable_numerics/nonlinear_block_solver.hpp b/src/smith/differentiable_numerics/nonlinear_block_solver.hpp index d75b2b1f0a..fd8dab0095 100644 --- a/src/smith/differentiable_numerics/nonlinear_block_solver.hpp +++ b/src/smith/differentiable_numerics/nonlinear_block_solver.hpp @@ -98,6 +98,9 @@ class NonlinearBlockSolverBase { /// @brief Set an inner-solve tolerance multiplier, e.g. for staggered solves. virtual void setInnerToleranceMultiplier(double multiplier) = 0; + /// @brief Return the configured solver print level. + virtual int printLevel() const { return 0; } + protected: mutable bool is_setup_ = false; ///< Records if this block solver has its preconditioner initialized. mutable NonlinearConvergenceContext convergence_context_ = {}; ///< Solver-owned convergence state for one solve. @@ -140,6 +143,9 @@ class NonlinearBlockSolver : public NonlinearBlockSolverBase { /// @brief Set the inner tolerance multiplier. void setInnerToleranceMultiplier(double multiplier) override { inner_tol_multiplier_ = multiplier; } + /// @brief Return the configured solver print level. + int printLevel() const override { return retained_nonlinear_options_ ? retained_nonlinear_options_->print_level : 0; } + /// @brief Build a fresh solver instance from retained config. std::shared_ptr cloneFresh() const; diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index e600f6cec6..42244f7d20 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -26,9 +26,8 @@ void applyBoundaryConditions(double time, const smith::BoundaryConditionManager* } } -std::vector block_solve(const std::vector& residual_evals, - const std::vector> block_indices, const FieldState& shape_disp, - const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers) @@ -50,18 +49,20 @@ std::vector block_solve(const std::vector& residual_evals // Validate block_indices bounds against states sizes for (size_t row = 0; row < num_rows_; ++row) { + std::vector slot_owners(states[row].size(), invalid_block_index); for (size_t col = 0; col < num_rows_; ++col) { - size_t idx = block_indices[row][col]; - if (idx != invalid_block_index) { + for (size_t idx : block_indices[row][col]) { 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())); + std::format("block_indices[{}][{}] contains {} outside states[{}] size {}", row, col, idx, row, + states[row].size())); + SLIC_ERROR_IF( + slot_owners[idx] != invalid_block_index, + std::format("State slot {} in row {} belongs to columns {} and {}", idx, row, slot_owners[idx], col)); + slot_owners[idx] = col; } } - // 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)); + SLIC_ERROR_IF(block_indices[row][row].empty(), + std::format("block_indices[{}][{}] diagonal entry must not be empty", row, row)); } std::vector num_state_inputs; @@ -115,9 +116,7 @@ std::vector block_solve(const std::vector& residual_evals std::vector diagonal_fields(num_rows); for (size_t row_i = 0; row_i < num_rows; ++row_i) { - size_t prime_unknown_row_i = block_indices[row_i][row_i]; - SLIC_ERROR_IF(prime_unknown_row_i == invalid_block_index, - "The primary unknown field (field index for block_index[n][n], must not be invalid)"); + size_t prime_unknown_row_i = block_indices[row_i][row_i].front(); diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][prime_unknown_row_i]); } @@ -128,9 +127,8 @@ std::vector block_solve(const std::vector& residual_evals for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - size_t prime_unknown_ij = block_indices[row_i][col_j]; - if (prime_unknown_ij != invalid_block_index) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + for (size_t slot : block_indices[row_i][col_j]) { + input_fields[row_i][slot] = diagonal_fields[col_j]; } } } @@ -172,13 +170,17 @@ std::vector block_solve(const std::vector& residual_evals std::vector row_field_inputs = input_fields[row_i]; std::vector tangent_weights(row_field_inputs.size(), 0.0); for (size_t col_j = 0; col_j < num_rows; ++col_j) { - size_t field_index_to_diff = block_indices[row_i][col_j]; - if (field_index_to_diff != invalid_block_index) { - tangent_weights[field_index_to_diff] = 1.0; + const auto& field_indices_to_diff = block_indices[row_i][col_j]; + if (!field_indices_to_diff.empty()) { + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 1.0; + } auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), getConstFieldPointers(row_field_inputs), tangent_weights); jacobians[row_i].emplace_back(std::move(jac_ij)); - tangent_weights[field_index_to_diff] = 0.0; + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 0.0; + } } else { jacobians[row_i].emplace_back(nullptr); } @@ -244,14 +246,15 @@ std::vector block_solve(const std::vector& residual_evals // make a copy so we don't accidentally override the original copy std::vector diagonal_fields(num_rows); for (size_t row_i = 0; row_i < num_rows; ++row_i) { - diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][block_indices[row_i][row_i]]); + diagonal_fields[row_i] = + std::make_shared(*input_fields[row_i][block_indices[row_i][row_i].front()]); *diagonal_fields[row_i] = *s[row_i]; } for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - if (block_indices[row_i][col_j] != invalid_block_index) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + for (size_t slot : block_indices[row_i][col_j]) { + input_fields[row_i][slot] = diagonal_fields[col_j]; } } } @@ -272,13 +275,17 @@ std::vector block_solve(const std::vector& residual_evals std::vector row_field_inputs = input_fields[row_i]; std::vector tangent_weights(row_field_inputs.size(), 0.0); for (size_t col_j = 0; col_j < num_rows; ++col_j) { - size_t field_index_to_diff = block_indices[row_i][col_j]; - if (field_index_to_diff != invalid_block_index) { - tangent_weights[field_index_to_diff] = 1.0; + const auto& field_indices_to_diff = block_indices[row_i][col_j]; + if (!field_indices_to_diff.empty()) { + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 1.0; + } auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), getConstFieldPointers(row_field_inputs), tangent_weights); jacobians[row_i].emplace_back(std::move(jac_ij)); - tangent_weights[field_index_to_diff] = 0.0; + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 0.0; + } } else { jacobians[row_i].emplace_back(nullptr); } @@ -351,8 +358,8 @@ std::vector block_solve(const std::vector& residual_evals // No sensitivity needed for primal fields for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - if (block_indices[row_i][col_j] != invalid_block_index) { - field_sensitivities[row_i][block_indices[row_i][col_j]] = nullptr; + for (size_t slot : block_indices[row_i][col_j]) { + field_sensitivities[row_i][slot] = nullptr; } } } diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index c083bd018b..1825bb7ab5 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -28,15 +28,21 @@ class DirichletBoundaryConditions; /// @brief magic number for representing a field which is not an argument of the weak form. static constexpr size_t invalid_block_index = std::numeric_limits::max() - 1; +/// Local weak-form argument slots associated with one solved field. +using BlockArgumentIndices = std::vector; + +/// Matrix mapping residual rows and solved columns to local argument slots. +using BlockArgumentMap = std::vector>; + /// @brief Solve a block nonlinear system of equations as defined by the vector of weak form /// @param residual_evals Vector of weak forms which define the equations to be solved -/// @param block_indices Matrix of index arguments specifying where in each WeakForm the unknown fields are passed in. +/// @param block_indices Matrix specifying every argument slot where each solved field is passed. /// Example: for a 2 weak-form system, with weak-forms, r1, r2 /// r1(a,b,c) /// r2(b,d,e,a) // with unknowns (with respect to the solver) being a, and b. -// r1 has unknowns a,b in the ‘slots’ 0, 1 -// r2 has unknowns a,b, in the ‘slots’ 3,0 +// r1 has unknowns a,b in slots {0}, {1} +// r2 has unknowns a,b in slots {3}, {0} /// @param shape_disp The mesh-morphed shape displacement /// @param states The time varying states as inputs to the weak form /// @param params The fixed field parameters as inputs to the weak form @@ -44,9 +50,8 @@ static constexpr size_t invalid_block_index = std::numeric_limits::max() /// @param solver The nonlinear block solver used to solve the system of equations /// @param bc_managers Holds information about which degrees of freedom (DOFS) /// @return Vector of field solutions satisfying the weak forms -std::vector block_solve(const std::vector& residual_evals, - const std::vector> block_indices, const FieldState& shape_disp, - const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers); @@ -67,8 +72,9 @@ inline FieldState solve(const WeakForm& residual_eval, const FieldState& shape_d const DirichletBoundaryConditions& bcs, size_t unknown_state_index = 0) { std::vector bc_managers{&bcs.getBoundaryConditionManager()}; - auto solutions = block_solve({const_cast(&residual_eval)}, {{unknown_state_index}}, shape_disp, {states}, - {params}, time_info, &solver, bc_managers); + BlockArgumentMap block_indices{{BlockArgumentIndices{unknown_state_index}}}; + auto solutions = block_solve({const_cast(&residual_eval)}, block_indices, shape_disp, {states}, {params}, + time_info, &solver, bc_managers); return solutions[0]; } diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp new file mode 100644 index 0000000000..22b2830d08 --- /dev/null +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -0,0 +1,275 @@ +// 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 plastic_mechanics_system.hpp + * @brief Defines the PlasticMechanicsSystem struct and its factory function. + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/nonlinear_block_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/differentiable_physics.hpp" +#include "smith/physics/functional_weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { +/** + * @brief Container for a plastic solid mechanics system. + * @tparam dim Spatial dimension. + * @tparam disp_order Order of the displacement basis. + * @tparam parameter_space Finite element spaces for optional parameters. + * + * This system relies on a fixed-point iteration to solve plasticity, equivalent to a traditional + * inner-outer Newton method. + * n in [0, N] denotes the timestep and k in [0, L] denotes the fixed-point iteration count. + * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> + * u^{n+1,k+1}. Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, + * epsilon_p^{n+1,k+1}]. + */ +template +struct PlasticMechanicsSystem : public SystemBase { + using SystemBase::SystemBase; + // Primary weak forms to solve for displacement + /// @brief using for SolidWeakFormType with inputs [u^{n+1, k+1}, Fp^{n+1, k}, epsilon_p^{n+1, k}] + using SolidWeakFormType = FunctionalWeakForm< + dim, H1, + Parameters, L2, L2, parameter_space...>>; + + /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, + /// u^{n+1, k+1}] + using PlasticDeformWeakFormType = + FunctionalWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; + + /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, + /// u^{n+1, k+1}] + using PlasticStrainWeakFormType = + FunctionalWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; + + // Primary weak forms + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + + // Internal variable weak forms to update in staggered solve + std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. + std::shared_ptr plastic_strain_weak_form; ///< Plastic strain weak form. + + // Primary variable bcs + std::shared_ptr disp_bc; ///< Displacement boundary conditions. + + // Internal variable bcs + std::shared_ptr plastic_deform_bc; ///< Plastic deformation gradient conditions. + std::shared_ptr plastic_strain_bc; ///< Plastic strain conditions + + std::shared_ptr quasistatic_time_rule; ///< Quasistatic time integration rule + std::shared_ptr + backward_euler_time_rule; ///< Backward euler time integration rule + + /** + * @brief Transform tensor variable stored as [dim * dim, 1] to [dim, dim] + */ + template + SMITH_HOST_DEVICE tensor recoverTensor(const tensor& f_state) const + { + tensor f_state_tensor{}; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + f_state_tensor(i, j) = f_state(i * dim + j); + } + } + + return f_state_tensor; + } + + /** + * @brief Transform tensor variable stored as [dim, dim] to [dim * dim, 1] + */ + template + SMITH_HOST_DEVICE tensor flattenTensor(const tensor& f_state) const + { + tensor f_state_flattened{}; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + f_state_flattened(i * dim + j) = f_state(i, j); + } + } + + return f_state_flattened; + } + + /** + * @brief Get the list of physical state fields for visualization. + * @return std::vector List of physical fields suitable for output. + */ + std::vector getOutputFieldStates() const + { + return {field_store->getField(field_store->prefix("displacement")), + field_store->getField(field_store->prefix("plastic_defgrad")), + field_store->getField(field_store->prefix("plastic_strain"))}; + } + + /** + * @brief Get information about reaction fields for this system. + * @return List of ReactionInfo structures. + */ + std::vector getReactionInfos() const + { + return {{field_store->prefix("solid_force"), + &field_store->getField(field_store->prefix("displacement")).get()->space()}}; + } + + /** + * @brief Set the material model for a domain, defining integrals for weak form. + * @tparam MaterialType The material model type. + * @param domain_name The name of the domain to apply the material to. + * @param mat The material model instance. + */ + template + void setMaterial(const std::string& domain_name, const MaterialType& mat) + { + solid_weak_form->addBodyIntegral( + domain_name, [&](auto t_info, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { + auto du_dX = get(u); + auto Fp_tensor = this->recoverTensor(get(Fp)); + auto dt = t_info.dt(); + + auto P = mat.firstPiolaStress(dt, du_dX, Fp_tensor, get(epsilon_p), params...); + + return tuple{zero{}, P}; + }); + } + + /** + * @brief Set the plasticity update model for a domain, defining integrals for weak form. + * @tparam MaterialType The material model type. + * @param domain_name The name of the domain to apply the material to. + * @param mat The material model instance. + */ + template + void setPlasticity(const std::string& domain_name, const MaterialType& mat) + { + auto captured_strain_rule = backward_euler_time_rule; + + plastic_deform_weak_form->addBodyIntegral( + domain_name, [=, this](auto t_info, auto /* X */, auto Fp, auto Fp_old, auto epsilon_p, auto epsilon_p_old, + auto u, auto... params) { + auto du_dX = get(u); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto epsilon_dot = captured_strain_rule->dot(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, get(epsilon_dot), du_dX, params...); + auto Fp_predict = this->flattenTensor(Fp_predict_tensor); + + return tuple(get(Fp) - Fp_predict, zero{}); + }); + + plastic_strain_weak_form->addBodyIntegral(domain_name, [=, this](auto t_info, auto /* X */, auto epsilon_p, + auto epsilon_p_old, auto /* Fp */, auto Fp_old, + auto u, auto... params) { + auto du_dX = get(u); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + auto epsilon_dot_predict = + mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX, params...); + return tuple(get(epsilon_dot) - epsilon_dot_predict, zero{}); + }); + } +}; + +/** + * @brief Factory function to build a plastic mechanics system with L2 internal variables. + */ +template +std::shared_ptr> buildPlasticMechanicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", + FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100, prepend_name); + + FieldType> shape_disp_type("shape_displacement"); + field_store->addShapeDisp(shape_disp_type); + + auto quasistatic_time_rule = std::make_shared(); + auto backward_euler_time_rule = std::make_shared(); + + // Displacement with quasi-static rule + FieldType> disp_type("displacement"); + auto disp_bc = field_store->addIndependent(disp_type, quasistatic_time_rule); + + // State variable fields + FieldType> plastic_defgrad_type("plastic_defgrad"); + auto plastic_defgrad_bc = field_store->addIndependent(plastic_defgrad_type, backward_euler_time_rule); + auto plastic_defgrad_old_type = + field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, "plastic_defgrad_old"); + + FieldType> plastic_strain_type("plastic_strain"); + auto plastic_strain_bc = field_store->addIndependent(plastic_strain_type, backward_euler_time_rule); + auto plastic_strain_old_type = + field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, "plastic_strain_old"); + + // Parameters + auto register_parameter = [&](auto& parameter_type) { + parameter_type.name = "param_" + parameter_type.name; + field_store->addParameter(parameter_type); + }; + + if constexpr (sizeof...(parameter_space) > 0) { + (register_parameter(parameter_types), ...); + } + + using SystemType = PlasticMechanicsSystem; + + // Main weak forms for mechanics + std::string solid_res_name = field_store->prefix("solid_residual"); + auto solid_weak_form = std::make_shared( + 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, plastic_defgrad_type, plastic_strain_type, + parameter_types...)); + + std::string plastic_defgrad_res_name = field_store->prefix("plastic_defgrad_residual"); + auto plastic_defgrad_weak_form = std::make_shared( + plastic_defgrad_res_name, field_store->getMesh(), field_store->getField(plastic_defgrad_type.name).get()->space(), + field_store->createSpaces(plastic_defgrad_res_name, plastic_defgrad_type.name, plastic_defgrad_type, + plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, + parameter_types...)); + + std::string plastic_strain_res_name = field_store->prefix("plastic_strain_residual"); + auto plastic_strain_weak_form = std::make_shared( + plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), + field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, + plastic_strain_old_type, plastic_defgrad_type, plastic_defgrad_old_type, disp_type, + parameter_types...)); + + // Build solver and advancer + std::vector> weak_forms{solid_weak_form, plastic_defgrad_weak_form, + plastic_strain_weak_form}; + + auto sys = std::make_shared(field_store, solver, weak_forms); + + sys->solid_weak_form = solid_weak_form; + sys->plastic_deform_weak_form = plastic_defgrad_weak_form; + sys->plastic_strain_weak_form = plastic_strain_weak_form; + sys->disp_bc = disp_bc; + sys->plastic_deform_bc = plastic_defgrad_bc; + sys->plastic_strain_bc = plastic_strain_bc; + sys->quasistatic_time_rule = quasistatic_time_rule; + sys->backward_euler_time_rule = backward_euler_time_rule; + + return sys; +} +} // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_system.hpp b/src/smith/differentiable_numerics/solid_mechanics_system.hpp index cc3a885984..9e9a603239 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_system.hpp @@ -307,7 +307,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons auto disp_time_rule_ptr = std::make_shared(); FieldType> shape_disp_type(field_store->prefix("shape_displacement")); - FieldType> disp_type(field_store->prefix("displacement_solve_state"), true); + FieldType> disp_type(field_store->prefix("displacement_solve_state")); FieldType> disp_old_type(field_store->prefix("displacement")); FieldType> velo_old_type(field_store->prefix("velocity")); FieldType> accel_old_type(field_store->prefix("acceleration")); @@ -351,7 +351,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons } if (has_stress_output) { - FieldType> stress_type(field_store->prefix("stress"), true); + FieldType> stress_type(field_store->prefix("stress")); FieldType> disp_as_input(disp_type.name); std::string stress_name = field_store->prefix("stress_projection"); sys->stress_weak_form = detail::buildWeakFormWithCoupling( diff --git a/src/smith/differentiable_numerics/state_variable_system.hpp b/src/smith/differentiable_numerics/state_variable_system.hpp index 25c237d2de..752dd1d768 100644 --- a/src/smith/differentiable_numerics/state_variable_system.hpp +++ b/src/smith/differentiable_numerics/state_variable_system.hpp @@ -160,7 +160,7 @@ auto buildInternalVariableSystemImpl(std::shared_ptr field_store, co { auto internal_variable_time_rule = std::make_shared(); - FieldType state_type(field_store->prefix("state_solve_state"), true); + FieldType state_type(field_store->prefix("state_solve_state")); FieldType state_old_type(field_store->prefix("state")); auto internal_variable_bc = field_store->getBoundaryConditions(state_type.name); diff --git a/src/smith/differentiable_numerics/system_base.cpp b/src/smith/differentiable_numerics/system_base.cpp index e13589eb90..7908039ce2 100644 --- a/src/smith/differentiable_numerics/system_base.cpp +++ b/src/smith/differentiable_numerics/system_base.cpp @@ -8,6 +8,8 @@ #include "smith/differentiable_numerics/reaction.hpp" #include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include + namespace smith { std::vector SystemBase::solve(const TimeInfo& time_info) const @@ -16,44 +18,58 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const for (const auto& wf : weak_forms) { weak_form_names.push_back(wf->name()); } - std::vector> index_map = field_store->indexMap(weak_form_names); + BlockArgumentMap index_map; + if (solve_result_field_names.empty()) { + index_map = field_store->indexMap(weak_form_names); + } std::vector> inputs; + std::vector> input_field_indices; if (!solve_input_field_names.empty()) { SLIC_ERROR_IF(solve_input_field_names.size() != weak_forms.size(), "solve_input_field_names size must match weak_forms size"); } for (size_t i = 0; i < weak_forms.size(); ++i) { std::vector fields_for_wk; + std::vector indices_for_wk; if (solve_input_field_names.empty()) { fields_for_wk = field_store->getStates(weak_forms[i]->name()); } else { for (const auto& field_name : solve_input_field_names[i]) { fields_for_wk.push_back(field_store->getField(field_name)); + indices_for_wk.push_back(field_store->getFieldIndex(field_name)); } } inputs.push_back(fields_for_wk); + input_field_indices.push_back(indices_for_wk); } std::vector bc_field_names; if (!solve_result_field_names.empty()) { SLIC_ERROR_IF(solve_result_field_names.size() != weak_forms.size(), "solve_result_field_names size must match weak_forms size"); + SLIC_ERROR_IF(solve_input_field_names.empty(), "solve_input_field_names must accompany solve_result_field_names"); bc_field_names = solve_result_field_names; + index_map = BlockArgumentMap(weak_forms.size(), std::vector(weak_forms.size())); + std::set result_field_indices; + for (const auto& field_name : solve_result_field_names) { + size_t field_index = field_store->getFieldIndex(field_name); + SLIC_ERROR_IF(!result_field_indices.insert(field_index).second, + "Each solve result field must have exactly one owning row"); + } for (size_t row = 0; row < weak_forms.size(); ++row) { for (size_t col = 0; col < weak_forms.size(); ++col) { - index_map[row][col] = invalid_block_index; - for (size_t arg = 0; arg < inputs[row].size(); ++arg) { - if (inputs[row][arg].get()->name() == solve_result_field_names[col]) { - index_map[row][col] = arg; - break; + size_t result_field_index = field_store->getFieldIndex(solve_result_field_names[col]); + for (size_t arg = 0; arg < input_field_indices[row].size(); ++arg) { + if (input_field_indices[row][arg] == result_field_index) { + index_map[row][col].push_back(arg); } } } - SLIC_ERROR_IF(index_map[row][row] == invalid_block_index, "Requested solve result field '" - << solve_result_field_names[row] - << "' is not an argument of weak form '" - << weak_form_names[row] << "'"); + SLIC_ERROR_IF(index_map[row][row].empty(), "Requested solve result field '" + << solve_result_field_names[row] + << "' is not an argument of weak form '" << weak_form_names[row] + << "'"); } } @@ -78,7 +94,7 @@ std::vector SystemBase::computeReactions(const TimeInfo& time_inf auto params = field_store->getParameterFields(); 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()); + std::string test_field_name = field_store->getWeakFormFieldNames(wf->name()).test; 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, field_store->getShapeDisp(), wf_fields, test_field)); diff --git a/src/smith/differentiable_numerics/system_base.hpp b/src/smith/differentiable_numerics/system_base.hpp index 4a4d9e7ec1..fde91bd22b 100644 --- a/src/smith/differentiable_numerics/system_base.hpp +++ b/src/smith/differentiable_numerics/system_base.hpp @@ -61,7 +61,7 @@ struct SystemBase { ///< independently; cycle-zero solves do not couple across subsystems. std::vector> post_solve_systems; ///< Optional systems solved after main state update. std::vector - solve_result_field_names; ///< Optional per-weak-form fields to solve/update instead of reaction fields. + solve_result_field_names; ///< Optional per-weak-form fields overriding registered solver-owned fields. std::vector> solve_input_field_names; ///< Optional per-weak-form input field ordering used during solve. diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index 5d234a8c20..284505110e 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -13,10 +13,7 @@ #include #include -#include -#include -#include -#include +#include "axom/slic.hpp" namespace smith { @@ -37,7 +34,7 @@ void SystemSolver::addSubsystemSolver(const std::vector& block_indices, { SLIC_ERROR_IF(!solver, "SystemSolver stage solver must be non-null"); SLIC_ERROR_IF(relaxation_factor <= 0.0 || relaxation_factor > 1.0, - axom::fmt::format("Stage relaxation_factor {} must be in (0, 1]", relaxation_factor)); + std::format("Stage relaxation_factor {} must be in (0, 1]", relaxation_factor)); stages_.push_back(Stage{block_indices, std::move(solver), relaxation_factor}); } @@ -55,8 +52,8 @@ void SystemSolver::appendStagesWithBlockMapping(const SystemSolver& subsystem_so remapped_block_indices.reserve(stage.block_indices.size()); for (size_t local_block_index : stage.block_indices) { SLIC_ERROR_IF(local_block_index >= global_block_indices.size(), - axom::fmt::format("Local block index {} exceeds subsystem size {}", local_block_index, - global_block_indices.size())); + std::format("Local block index {} exceeds subsystem size {}", local_block_index, + global_block_indices.size())); remapped_block_indices.push_back(global_block_indices[local_block_index]); } } @@ -65,8 +62,7 @@ void SystemSolver::appendStagesWithBlockMapping(const SystemSolver& subsystem_so } std::vector SystemSolver::solve(const std::vector& residual_evals, - const std::vector>& block_indices, - const FieldState& shape_disp, + const BlockArgumentMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, @@ -75,6 +71,21 @@ std::vector SystemSolver::solve(const std::vector& residu SLIC_ERROR_IF(stages_.empty(), "SystemSolver has no stages defined."); size_t num_residuals = residual_evals.size(); + SLIC_ERROR_IF(states.size() != num_residuals, "State rows must match residual count"); + SLIC_ERROR_IF(block_indices.size() != num_residuals, "Block index rows must match residual count"); + for (size_t row = 0; row < num_residuals; ++row) { + SLIC_ERROR_IF(block_indices[row].size() != num_residuals, + std::format("Block index columns for row {} must match residual count", row)); + for (size_t col = 0; col < num_residuals; ++col) { + for (size_t slot : block_indices[row][col]) { + SLIC_ERROR_IF(slot >= states[row].size(), + std::format("Block slot {} for row {}, column {} exceeds state count {}", slot, row, col, + states[row].size())); + } + } + SLIC_ERROR_IF(block_indices[row][row].empty(), + std::format("Block index diagonal for row {} must not be empty", row)); + } std::vector active_stages = stages_; for (auto& stage : active_stages) { if (stage.block_indices.empty()) { @@ -83,7 +94,7 @@ std::vector SystemSolver::solve(const std::vector& residu } for (size_t block_index : stage.block_indices) { SLIC_ERROR_IF(block_index >= num_residuals, - axom::fmt::format("Stage block index {} exceeds residual count {}", block_index, num_residuals)); + std::format("Stage block index {} exceeds residual count {}", block_index, num_residuals)); } } // Set the inner tolerance factor based on the number of stages. For single-stage @@ -95,6 +106,10 @@ std::vector SystemSolver::solve(const std::vector& residu stage.solver->setInnerToleranceMultiplier(inner_tol_factor); } + const bool print_staggered_progress = std::all_of(active_stages.begin(), active_stages.end(), [](const auto& stage) { + return stage.solver && stage.solver->printLevel() > 0; + }); + // Reset each stage solver's convergence tracking (e.g. initial residual norm for rel-tol) for (const auto& stage : active_stages) { stage.solver->resetConvergenceState(); @@ -104,16 +119,6 @@ std::vector SystemSolver::solve(const std::vector& residu // Working copy of states, updated in-place as stages solve std::vector> current_states = states; - // Pre-compute name -> (row, slot) routing so the propagation loop avoids O(N*M) string compares - // on every staggered iteration. Field-name identity within current_states is invariant across - // the iteration loop: only values are replaced, never the underlying name. - std::unordered_map>> field_routing; - for (size_t r = 0; r < num_residuals; ++r) { - for (size_t slot = 0; slot < current_states[r].size(); ++slot) { - field_routing[current_states[r][slot].get()->name()].emplace_back(r, slot); - } - } - // Helper lambda to assemble input pointers, evaluate residual, and zero essential BCs auto eval_residual_and_zero_bcs = [&](size_t global_row) { std::vector input_ptrs; @@ -142,13 +147,23 @@ std::vector SystemSolver::solve(const std::vector& residu } for (int iter = 0; iter < max_staggered_iterations_; ++iter) { + if (print_staggered_progress) { + SLIC_INFO_ROOT("Solving staggered iteration " << iter); + smith::logger::flush(); + } + // --- Run each stage --- for (size_t stage_idx = 0; stage_idx < active_stages.size(); ++stage_idx) { + if (print_staggered_progress) { + SLIC_INFO_ROOT("Solving stage " << stage_idx); + smith::logger::flush(); + } + const auto& stage = active_stages[stage_idx]; size_t num_stage_blocks = stage.block_indices.size(); std::vector stage_residuals; - std::vector> stage_block_indices; + BlockArgumentMap stage_block_indices; std::vector> stage_states; std::vector> stage_params; std::vector stage_bc_managers; @@ -160,7 +175,7 @@ std::vector SystemSolver::solve(const std::vector& residu stage_states.push_back(current_states[global_row]); stage_params.push_back(params[global_row]); - std::vector row_indices(num_stage_blocks, invalid_block_index); + std::vector row_indices(num_stage_blocks); for (size_t col_idx = 0; col_idx < num_stage_blocks; ++col_idx) { size_t global_col = stage.block_indices[col_idx]; row_indices[col_idx] = block_indices[global_row][global_col]; @@ -172,23 +187,20 @@ std::vector SystemSolver::solve(const std::vector& residu block_solve(stage_residuals, stage_block_indices, shape_disp, stage_states, stage_params, time_info, stage.solver.get(), stage_bc_managers); - // Propagate updated fields to every residual input that references the solved field. - // Match by field name (looked up via the pre-computed routing map): coupling fields appear - // as fixed inputs in other rows and therefore do not have a valid unknown-block entry there. + // Propagate updated fields to every residual input sharing the solved FieldStore index. // Apply relaxation: x_new = omega * x_solved + (1 - omega) * x_k. for (size_t i = 0; i < num_stage_blocks; ++i) { size_t global_col = stage.block_indices[i]; FieldState new_state = stage_solutions[i]; if (stage.relaxation_factor != 1.0) { - FieldState old_state = current_states[global_col][block_indices[global_col][global_col]]; + FieldState old_state = current_states[global_col][block_indices[global_col][global_col].front()]; new_state = weighted_average(new_state, old_state, stage.relaxation_factor); } - auto it = field_routing.find(new_state.get()->name()); - if (it != field_routing.end()) { - for (const auto& [r, slot] : it->second) { - current_states[r][slot] = new_state; + for (size_t row = 0; row < num_residuals; ++row) { + for (size_t slot : block_indices[row][global_col]) { + current_states[row][slot] = new_state; } } } @@ -213,7 +225,7 @@ std::vector SystemSolver::solve(const std::vector& residu } } if (all_converged) { - SLIC_INFO_ROOT(axom::fmt::format("Staggered iteration converged after {} iteration(s)", iter + 1)); + SLIC_INFO_ROOT(std::format("Staggered iteration converged after {} iteration(s)", iter + 1)); break; } } @@ -223,7 +235,7 @@ std::vector SystemSolver::solve(const std::vector& residu std::vector final_solutions; final_solutions.reserve(num_residuals); for (size_t r = 0; r < num_residuals; ++r) { - size_t s_idx = block_indices[r][r]; + size_t s_idx = block_indices[r][r].front(); final_solutions.push_back(current_states[r][s_idx]); } diff --git a/src/smith/differentiable_numerics/system_solver.hpp b/src/smith/differentiable_numerics/system_solver.hpp index 71351f1c53..1d8157391c 100644 --- a/src/smith/differentiable_numerics/system_solver.hpp +++ b/src/smith/differentiable_numerics/system_solver.hpp @@ -10,6 +10,7 @@ #include #include #include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "smith/numerics/solver_config.hpp" #include "smith/physics/common.hpp" @@ -59,16 +60,15 @@ class SystemSolver { /// @brief Solves the multiphysics system using staggered iterations. /// @param residual_evals Vector of WeakForm evaluations for each block. - /// @param block_indices Block indices for each residual evaluation. + /// @param block_indices Argument slot lists for each residual and solved field. /// @param shape_disp Current shape displacement. /// @param states Nested vector of field states. /// @param params Nested vector of parameters. /// @param time_info Current time information. /// @param bc_managers Managers for boundary conditions. /// @return Updated field states. - std::vector solve(const std::vector& residual_evals, - const std::vector>& block_indices, const FieldState& shape_disp, - const std::vector>& states, + std::vector solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const std::vector& bc_managers) const; diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index d9c676c94d..f74653b13f 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -17,6 +17,7 @@ set(differentiable_numerics_test_source test_thermal_static.cpp test_solid_static_with_internal_vars.cpp test_multiphysics_time_integrator.cpp + test_differentiable_plasticity.cpp test_thermo_mechanics_with_internal_vars.cpp test_mixed_poisson.cpp test_amg_vdim_block_preconditioner.cpp diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp new file mode 100644 index 0000000000..d90cb74b3b --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -0,0 +1,401 @@ +// 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/mesh_utils/mesh_utils.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/solver_config.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/solid_material.hpp" + +#include "smith/differentiable_numerics/nonlinear_block_solver.hpp" +#include "smith/differentiable_numerics/system_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/plastic_mechanics_system.hpp" + +namespace smith { + +LinearSolverOptions lin_opts{.linear_solver = LinearSolver::CG, + .preconditioner = Preconditioner::HypreAMG, + .relative_tol = 1e-15, + .absolute_tol = 1e-15, + .max_iterations = 200, + .print_level = 0}; + +NonlinearSolverOptions nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1e-9, + .absolute_tol = 1e-10, + .max_iterations = 25, + .print_level = 2}; + +/// @brief Differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening +template +class DifferentiableJ2SmallStrain { + public: + DifferentiableJ2SmallStrain(HardeningType hardeningModel, double youngsModulus, double poissonsRatio, + double hardeningModulus, double density) + : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density) + { + } + + /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material + * state */ + template + SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, + const tensor& Fp, const T3& /* epsilon_p */) const + { + constexpr auto I = Identity(); + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + + return s + p * I; + } + + /** @brief calculate the plastic deformation gradient */ + template + SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, const T2& epsilon_dot, + const tensor& du_dX) const + { + using std::sqrt; + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto Np = 1.5 * eta / q; + + return Fp_old + epsilon_dot * dt * Np; + } + + /** @brief calculate the plastic strain */ + template + SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, + const tensor& Fp_old, const tensor& du_dX) const + { + using std::sqrt; + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto epsilon_dot_predict = q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + auto non_negativity = epsilon_dot_predict >= 0.0; + return epsilon_dot_predict * non_negativity; + } + + private: + HardeningType hardening; ///< Flow stress hardening model + double E; ///< Young's modulus + double nu; ///< Poisson's ratio + double Hk; ///< Kinematic hardening modulus + double rho; ///< Mass density +}; + +TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) +{ + MPI_Barrier(MPI_COMM_WORLD); + + NonlinearSolverOptions primal_nonlin_opts = nonlin_opts; + LinearSolverOptions primal_lin_opts = lin_opts; + + NonlinearSolverOptions state_nonlin_opts = nonlin_opts; + state_nonlin_opts.relative_tol = 1e-7; + state_nonlin_opts.absolute_tol = 1e-8; + LinearSolverOptions state_lin_opts = lin_opts; + state_lin_opts.preconditioner = Preconditioner::HypreJacobi; + state_lin_opts.relative_tol = 1e-10; + state_lin_opts.absolute_tol = 1e-10; + + int serial_refinement = 0; + int parallel_refinement = 0; + + static constexpr int dim = 3; + static constexpr int order = 2; + + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "plasticity_small_strain"); + + std::string filename = SMITH_REPO_DIR "/data/meshes/beam-hex.mesh"; + const std::string meshtag = "mesh"; + auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); + + auto staggered_coupled_solver = std::make_shared(50); + auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); + auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + + staggered_coupled_solver->addSubsystemSolver({0}, primal_block_solver); + staggered_coupled_solver->addSubsystemSolver({2}, strain_block_solver); + staggered_coupled_solver->addSubsystemSolver({1}, defgrad_block_solver); + + auto plastic_mechanics_system = buildPlasticMechanicsSystem(mesh, staggered_coupled_solver); + + using Hardening = solid_mechanics::LinearHardening; + Hardening hardening{.sigma_y = 40.0, .Hi = 50.0, .eta = 0.0}; + DifferentiableJ2SmallStrain mat(hardening, 1e+4, 0.25, 5.0, 1.0); + + plastic_mechanics_system->setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system->setPlasticity(mesh->entireBodyName(), mat); + + // prescribe zero displacement at the supported end of the beam, + mesh->addDomainOfBoundaryElements("support", by_attr(1)); + plastic_mechanics_system->disp_bc->setFixedVectorBCs(mesh->domain("support")); + + // apply a displacement along z to the the tip of the beam + mesh->addDomainOfBoundaryElements("tip", by_attr(2)); + auto translated_in_z = [](double t, tensor) { + tensor u{}; + u[2] = t * (t - 1); + return u; + }; + plastic_mechanics_system->disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + + double dt = 0.1; + double time = 0.0; + auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); + auto states = plastic_mechanics_system->field_store->getStateFields(); + auto params = plastic_mechanics_system->field_store->getParameterFields(); + const double initial_plastic_defgrad_l2 = norm( + *plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")) + .get(), + 2.0); + const double initial_plastic_strain_l2 = norm( + *plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_strain")) + .get(), + 2.0); + + auto advancer = makeAdvancer(plastic_mechanics_system); + + auto pv_writer = createParaviewWriter(*mesh, states, "J2_small_strain"); + pv_writer.write(0, 0.0, states); + + std::vector reactions; + for (size_t step = 0; step < 10; ++step) { + std::tie(states, reactions) = advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); + time += dt; + pv_writer.write(step + 1, time, states); + } + + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); + auto plastic_defgrad = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")); + auto plastic_strain = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_strain")); + double final_disp_l2 = norm(*displacement.get(), 2.0); + EXPECT_GT(final_disp_l2, 0.0); + EXPECT_GT(norm(*plastic_defgrad.get(), 2.0), initial_plastic_defgrad_l2 + 1.0e-10); + EXPECT_GT(norm(*plastic_strain.get(), 2.0), initial_plastic_strain_l2 + 1.0e-10); +} + +/// @brief Parameterized differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening +template +class ParameterizedDifferentiableJ2SmallStrain { + public: + ParameterizedDifferentiableJ2SmallStrain(HardeningType hardeningModel, double baseYoungsModulus, double poissonsRatio, + double baseHardeningModulus, double density) + : hardening(hardeningModel), E0(baseYoungsModulus), nu(poissonsRatio), Hk0(baseHardeningModulus), rho(density) + { + } + + /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material + * state */ + template + SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, + const tensor& Fp, const T3& /* epsilon_p */, const T4& deltaE, + const T5& /* deltaHk */) const + { + constexpr auto I = Identity(); + auto E = E0 * (1.0 + get(deltaE)); + auto K = E / (3.0 * (1.0 - 2.0 * nu)); + auto G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + + return s + p * I; + } + + /** @brief calculate the plastic deformation gradient */ + template + SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, const T2& epsilon_dot, + const tensor& du_dX, const T4& deltaE, const T5& deltaHk) const + { + using std::sqrt; + auto E = E0 * (1.0 + get(deltaE)); + auto Hk = Hk0 * (1.0 + get(deltaHk)); + auto G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto Np = 1.5 * eta / q; + + return Fp_old + epsilon_dot * dt * Np; + } + + /** @brief calculate the plastic strain */ + template + SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, + const tensor& Fp_old, const tensor& du_dX, + const T5& deltaE, const T6& deltaHk) const + { + using std::sqrt; + auto E = E0 * (1.0 + get(deltaE)); + auto Hk = Hk0 * (1.0 + get(deltaHk)); + auto G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto epsilon_dot_predict = q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + auto non_negativity = epsilon_dot_predict >= 0.0; + return epsilon_dot_predict * non_negativity; + } + + private: + HardeningType hardening; ///< Flow stress hardening model + double E0; ///< Base Young's modulus + double nu; ///< Poisson's ratio + double Hk0; ///< Base kinematic hardening modulus + double rho; ///< Mass density +}; + +TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) +{ + MPI_Barrier(MPI_COMM_WORLD); + + std::string physics_name = "differentiable_plasticity"; + + NonlinearSolverOptions primal_nonlin_opts = nonlin_opts; + primal_nonlin_opts.print_level = 0; + LinearSolverOptions primal_lin_opts = lin_opts; + + NonlinearSolverOptions state_nonlin_opts = nonlin_opts; + state_nonlin_opts.relative_tol = 1e-7; + state_nonlin_opts.absolute_tol = 1e-8; + state_nonlin_opts.print_level = 0; + LinearSolverOptions state_lin_opts = lin_opts; + state_lin_opts.preconditioner = Preconditioner::HypreJacobi; + state_lin_opts.relative_tol = 1e-10; + state_lin_opts.absolute_tol = 1e-10; + + int serial_refinement = 0; + int parallel_refinement = 0; + + static constexpr int dim = 3; + static constexpr int order = 1; + + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "plasticity_small_strain"); + + std::string filename = SMITH_REPO_DIR "/data/meshes/beam-hex.mesh"; + const std::string meshtag = "mesh"; + auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); + + auto staggered_coupled_solver = std::make_shared(100); + auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); + auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + + staggered_coupled_solver->addSubsystemSolver({0}, primal_block_solver); + staggered_coupled_solver->addSubsystemSolver({2}, strain_block_solver); + staggered_coupled_solver->addSubsystemSolver({1}, defgrad_block_solver); + + auto plastic_mechanics_system = buildPlasticMechanicsSystem( + mesh, staggered_coupled_solver, physics_name, FieldType>("deltaE"), FieldType>("deltaHk")); + + using Hardening = solid_mechanics::LinearHardening; + Hardening hardening{.sigma_y = 20.0, .Hi = 50.0, .eta = 0.0}; + ParameterizedDifferentiableJ2SmallStrain mat(hardening, 1.5e+4, 0.25, 500.0, 1.0); + + plastic_mechanics_system->setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system->setPlasticity(mesh->entireBodyName(), mat); + + // prescribe zero displacement at the supported end of the beam, + mesh->addDomainOfBoundaryElements("support", by_attr(1)); + plastic_mechanics_system->disp_bc->setFixedVectorBCs(mesh->domain("support")); + + // apply a displacement along z to the the tip of the beam + mesh->addDomainOfBoundaryElements("tip", by_attr(2)); + auto translated_in_z = [](double t, tensor) { + tensor u{}; + u[2] = t; + return u; + }; + plastic_mechanics_system->disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + + auto advancer = makeAdvancer(plastic_mechanics_system); + std::unique_ptr differentiable_plasticity = + makeDifferentiablePhysics(plastic_mechanics_system, advancer, physics_name); + + auto pv_writer = createParaviewWriter(*mesh, differentiable_plasticity->getFieldStatesAndParamStates(), physics_name, + smith::ParaviewWriter::Options{.write_duals = false}); + pv_writer.write(0, differentiable_plasticity->time(), differentiable_plasticity->getFieldStatesAndParamStates()); + + double dt = 0.1; + for (size_t m = 0; m < 2; ++m) { + differentiable_plasticity->advanceTimestep(dt); + pv_writer.write(m + 1, differentiable_plasticity->time(), + differentiable_plasticity->getFieldStatesAndParamStates()); + } + + auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); + auto params = plastic_mechanics_system->field_store->getParameterFields(); + + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); + auto disp_squared = 0.5 * innerProduct(displacement, displacement); + gretl::set_as_objective(disp_squared); + + EXPECT_GT(checkGradWrt(disp_squared, shape_disp, 2e-4, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[0], 1e-2, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[1], 1e-1, 4, true), 0.95); + + std::vector sensitivity_fields{shape_disp, params[0], params[1]}; + auto sens_writer = createParaviewWriter(*mesh, sensitivity_fields, physics_name + "_final_sensitivities"); + + auto& graph = disp_squared.data_store(); + graph.reset(); + gretl::set_as_objective(disp_squared); + graph.back_prop(); + + sens_writer.write(0, differentiable_plasticity->time(), sensitivity_fields); +} + +} // 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_mixed_poisson.cpp b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp index 9c21313fd3..e9f8e5b61c 100644 --- a/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp +++ b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp @@ -216,7 +216,7 @@ TEST_P(BlockPreconditionerTest, BlockSolve) auto nonlinear_block_solver = smith::buildNonlinearBlockSolver(nonlin_opts, linear_options, *mesh); - auto sols = block_solve({&con_form, &bal_form}, {{0, 1}, {0, 1}}, shape_disp, {con_arguments, bal_arguments}, + auto sols = block_solve({&con_form, &bal_form}, {{{0}, {1}}, {{0}, {1}}}, shape_disp, {con_arguments, bal_arguments}, {flux_params, potential_params}, smith::TimeInfo(time.get(), dt.get(), cycle), nonlinear_block_solver.get(), {flux_bc_manager.get(), potential_bc_manager.get()}); diff --git a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp index c2cb034247..2225ef95aa 100644 --- a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp @@ -121,6 +121,24 @@ class ConstantNonlinearBlockSolver : public NonlinearBlockSolverBase { mutable int solve_calls_ = 0; }; +class ResidualRecordingNonlinearBlockSolver : public NoOpNonlinearBlockSolver { + public: + std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residual, + std::function>(const std::vector&)> jacobian) const override + { + auto residuals = residual(u_guesses); + residual_norm_ = residuals.front().Norml2(); + return NoOpNonlinearBlockSolver::solve(u_guesses, std::move(residual), std::move(jacobian)); + } + + double residualNorm() const { return residual_norm_; } + + private: + mutable double residual_norm_ = 0.0; +}; + class NeedsInitialSolveRule : public QuasiStaticRule { public: bool requiresInitialAccelerationSolve() const override { return true; } @@ -173,7 +191,7 @@ bool allNodesOnBoundary(const std::vector& nodes, double x_target) } // namespace -TEST(MultiphysicsTimeIntegrator, CycleZeroUsesBcsForReactionFieldNotUnknownZero) +TEST(MultiphysicsTimeIntegrator, CycleZeroUsesOwnedFieldBoundaryConditions) { axom::sidre::DataStore datastore; StateManager::initialize(datastore, "multiphysics_time_integrator"); @@ -195,9 +213,6 @@ TEST(MultiphysicsTimeIntegrator, CycleZeroUsesBcsForReactionFieldNotUnknownZero) FieldType> displacement_type("displacement"); auto displacement_bc = field_store->addIndependent(displacement_type, quasi_static); - ASSERT_TRUE(temperature_type.is_unknown); - ASSERT_TRUE(displacement_type.is_unknown); - temperature_bc->setScalarBCs<2>(mesh->domain("left"), [](double, tensor) { return 0.0; }); displacement_bc->setScalarBCs<2>(mesh->domain("right"), [](double, tensor) { return 1.0; }); @@ -514,6 +529,125 @@ TEST(SystemSolver, AppendsStagesWithBlockMappingForCombinedSubsystems) StateManager::reset(); } +TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) +{ + constexpr int dim = 2; + constexpr int order = 1; + using ShapeSpace = H1; + using FieldBTestSpace = H1; + using FieldBUnknownSpace = H1; + using FieldAFirstOccurrenceSpace = H1; + using FieldASecondOccurrenceSpace = H1; + using FieldAHistorySpace = H1; + using FieldBParameters = + Parameters; + using FieldBWeakForm = FunctionalWeakForm; + + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "system_solver_repeated_dependency_routing"); + + auto mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL, true, 1.0, 1.0), + "system_solver_repeated_dependency_routing_mesh"); + auto field_store = std::make_shared(mesh, 20); + FieldType shape_disp_type("shape_displacement"); + field_store->addShapeDisp(shape_disp_type); + + auto static_rule = std::make_shared(); + auto backward_euler_rule = std::make_shared(); + FieldType field_a_type("field_a"); + FieldType field_b_type("field_b"); + field_store->addIndependent(field_a_type, backward_euler_rule); + auto field_a_old_type = field_store->addDependent(field_a_type, FieldStore::TimeDerivative::VAL, "field_a_old"); + field_store->addIndependent(field_b_type, static_rule); + + auto field_a_weak_form = buildScalarDiffusionWeakForm("field_a_residual", mesh, field_store, field_a_type); + auto field_b_weak_form = std::make_shared( + "field_b_residual", mesh, field_store->getField(field_b_type.name).get()->space(), + field_store->createSpaces("field_b_residual", field_b_type.name, field_b_type, field_a_type, field_a_type, + field_a_old_type)); + field_b_weak_form->addBodySource(mesh->entireBodyName(), + [](auto, auto, auto, auto field_a_first, auto field_a_second, auto field_a_old) { + return field_a_first + field_a_second - 2.0 * field_a_old; + }); + + *field_store->getField(field_a_old_type.name).get() = *field_store->getField(field_a_type.name).get(); + + const std::vector residual_names{field_a_weak_form->name(), field_b_weak_form->name()}; + const auto block_indices = field_store->indexMap(residual_names); + EXPECT_EQ(block_indices[0][0], BlockArgumentIndices({0})); + EXPECT_TRUE(block_indices[0][1].empty()); + EXPECT_EQ(block_indices[1][0], BlockArgumentIndices({1, 2})); + EXPECT_EQ(block_indices[1][1], BlockArgumentIndices({0})); + + // Force field_a to a known nonzero solution without testing nonlinear convergence. Relaxation produces a + // predictable update, isolating whether SystemSolver propagates every repeated current-field slot while leaving + // history unchanged. + constexpr double field_a_solved_value = 3.0; + constexpr double field_a_relaxation_factor = 0.5; + constexpr int max_staggered_iterations = 1; + constexpr bool exact_staggered_steps = true; + auto field_a_solver = std::make_shared(field_a_solved_value); + auto field_b_solver = std::make_shared(); + SystemSolver solver(max_staggered_iterations, exact_staggered_steps); + solver.addSubsystemSolver({0}, field_a_solver, field_a_relaxation_factor); + solver.addSubsystemSolver({1}, field_b_solver); + + const std::vector residuals{field_a_weak_form.get(), field_b_weak_form.get()}; + const std::vector> states{field_store->getStates(field_a_weak_form->name()), + field_store->getStates(field_b_weak_form->name())}; + const std::vector> params(residuals.size()); + const auto bc_managers = field_store->getBoundaryConditionManagers(residual_names); + const TimeInfo time_info(0.0, 1.0, 0); + auto solutions = + solver.solve(residuals, block_indices, field_store->getShapeDisp(), states, params, time_info, bc_managers); + + // field_a starts at zero, ConstantNonlinearBlockSolver returns 3.0, and 0.5 relaxation produces 1.5. + EXPECT_NEAR((*solutions[0].get())[0], field_a_relaxation_factor * field_a_solved_value, 1.0e-12); + + // Expected field_b inputs contain relaxed field_a in both current-value slots and unchanged field_a_old history. + auto expected_field_b_states = states[1]; + expected_field_b_states[1] = solutions[0]; + expected_field_b_states[2] = solutions[0]; + auto expected_residual = field_b_weak_form->residual(time_info, field_store->getShapeDisp().get().get(), + getConstFieldPointers(expected_field_b_states)); + + EXPECT_NEAR(field_b_solver->residualNorm(), expected_residual.Norml2(), 1.0e-12); + EXPECT_GT(expected_residual.Norml2(), 1.0e-8); + + StateManager::reset(); +} + +TEST(FieldStore, UsesExplicitUnknownFieldInsteadOfTestField) +{ + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "field_store_explicit_unknown_field"); + + auto mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL, true, 1.0, 1.0), + "field_store_explicit_unknown_field_mesh"); + auto field_store = std::make_shared(mesh, 10); + FieldType> test_type("test_field"); + FieldType> unknown_type("unknown_field"); + auto static_rule = std::make_shared(); + field_store->addIndependent(test_type, static_rule); + auto unknown_bcs = field_store->addIndependent(unknown_type, static_rule); + + // Most weak forms use one field for both solver-owned unknown and test space. This case separates those roles: + // unknown_field owns the solve and boundary conditions, while test_field defines residual and reaction spaces. + field_store->createSpaces("residual", {.unknown = unknown_type.name, .test = test_type.name}, test_type, + unknown_type); + + const auto block_indices = field_store->indexMap({"residual"}); + EXPECT_EQ(field_store->getWeakFormFieldNames("residual").test, test_type.name); + EXPECT_EQ(field_store->getWeakFormFieldNames("residual").unknown, unknown_type.name); + // Unlike the matching-name case, the diagonal unknown occupies argument slot 1 instead of test-field slot 0. + EXPECT_EQ(block_indices[0][0], BlockArgumentIndices({1})); + // Boundary conditions follow the solver-owned unknown; reaction output follows the test space. + EXPECT_EQ(field_store->getBoundaryConditionManagers({"residual"})[0], &unknown_bcs->getBoundaryConditionManager()); + EXPECT_EQ(field_store->getReactionInfos()[0].space, &field_store->getField(test_type.name).get()->space()); + + StateManager::reset(); +} + } // namespace smith int main(int argc, char* argv[]) diff --git a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp index af2b9429d2..b52a69b763 100644 --- a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp +++ b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp @@ -303,7 +303,7 @@ TEST_P(BlockPreconditionerTest, BlockSolve) auto nonlinear_block_solver = smith::buildNonlinearBlockSolver(nonlin_opts, linear_options, *mesh); - auto sols = block_solve({&T1_form, &T2_form}, {{0, 1}, {0, 1}}, shape_disp, {T1_arguments, T2_arguments}, + auto sols = block_solve({&T1_form, &T2_form}, {{{0}, {1}}, {{0}, {1}}}, shape_disp, {T1_arguments, T2_arguments}, {T1_params, T2_params}, smith::TimeInfo(time.get(), dt.get(), cycle), nonlinear_block_solver.get(), {T1_bc_manager.get(), T2_bc_manager.get()}); diff --git a/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp index 333e8d39c0..4f0bf3ca80 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp @@ -418,15 +418,15 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) 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(std::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())); + SLIC_INFO_ROOT( + std::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())); + SLIC_INFO_ROOT( + std::format("{} {} {}", params[p].get()->name(), params[p].get()->Norml2(), params[p].get_dual()->Norml2())); } } @@ -473,7 +473,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) createSolidMechanicsBasePhysics(physics_name, mesh, solid_nonlinear_opts, solid_linear_options); double qoi = integrateForward(*physics, num_steps_, dt_, physics_name + "_reactions"); - SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); + SLIC_INFO_ROOT(std::format("{}", qoi)); // Check that reaction forces are zero away from Dirichlet DOFs auto reactions = physics->getReactionStates(); @@ -491,14 +491,13 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) 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(std::format("{} {}", name_and_state_sensitivity.first, name_and_state_sensitivity.second.Norml2())); } - SLIC_INFO_ROOT(axom::fmt::format("{} {}", shape_sensitivity.name(), shape_sensitivity.Norml2())); + SLIC_INFO_ROOT(std::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())); + SLIC_INFO_ROOT(std::format("{} {}", parameter_sensitivities[p].name(), parameter_sensitivities[p].Norml2())); } } @@ -557,13 +556,13 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) // 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)); + SLIC_INFO_ROOT(std::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)); + SLIC_INFO_ROOT(std::format("Parameter {} sensitivity difference: {:.6e}", p, param_diff)); EXPECT_LT(param_diff, tol); } @@ -572,7 +571,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) 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)); + SLIC_INFO_ROOT(std::format("Comparing sensitivity for {}: {} vs {}", suffix, nameG, nameB)); // Find Gretl dual const FiniteElementDual* dualG = nullptr; @@ -585,7 +584,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) 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)); + SLIC_INFO_ROOT(std::format("Initial state {} sensitivity difference: {:.6e}", suffix, state_diff)); EXPECT_LT(state_diff, tol); } } diff --git a/src/smith/differentiable_numerics/thermal_system.hpp b/src/smith/differentiable_numerics/thermal_system.hpp index 1a7f54baee..f82ca2ef30 100644 --- a/src/smith/differentiable_numerics/thermal_system.hpp +++ b/src/smith/differentiable_numerics/thermal_system.hpp @@ -214,7 +214,7 @@ auto buildThermalSystemImpl(std::shared_ptr field_store, const Coupl auto temperature_time_rule_ptr = std::make_shared(); FieldType> shape_disp_type(field_store->prefix("shape_displacement")); - FieldType> temperature_type(field_store->prefix("temperature_solve_state"), true); + FieldType> temperature_type(field_store->prefix("temperature_solve_state")); FieldType> temperature_old_type(field_store->prefix("temperature")); auto temperature_bc = field_store->getBoundaryConditions(temperature_type.name); diff --git a/src/smith/numerics/block_preconditioner.cpp b/src/smith/numerics/block_preconditioner.cpp index 38478a85cf..a0a543f499 100644 --- a/src/smith/numerics/block_preconditioner.cpp +++ b/src/smith/numerics/block_preconditioner.cpp @@ -7,7 +7,6 @@ #include "mfem.hpp" #include "axom/slic/core/SimpleLogger.hpp" -#include "axom/fmt.hpp" namespace smith { @@ -64,10 +63,9 @@ void BlockDiagonalPreconditioner::SetOperator(const mfem::Operator& jacobian) block_jacobian_ = dynamic_cast(&jacobian); MFEM_VERIFY(block_jacobian_, "Jacobian must be a BlockOperator"); - SLIC_ERROR_ROOT_IF( - block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, - axom::fmt::format("BlockDiagonalPreconditioner solver count ({}) must match block operator size ({}x{})", - num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + SLIC_ERROR_ROOT_IF(block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, + std::format("BlockDiagonalPreconditioner solver count ({}) must match block operator size ({}x{})", + num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); solver_diag_ = std::make_unique(block_offsets_); @@ -219,8 +217,8 @@ void BlockTriangularPreconditioner::SetOperator(const mfem::Operator& jacobian) SLIC_ERROR_ROOT_IF( block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, - axom::fmt::format("BlockTriangularPreconditioner solver count ({}) must match block operator size ({}x{})", - num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + std::format("BlockTriangularPreconditioner solver count ({}) must match block operator size ({}x{})", num_blocks_, + block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); @@ -406,11 +404,10 @@ void BlockSchurPreconditioner::SetOperator(const mfem::Operator& jacobian) MFEM_VERIFY(block_jacobian_, "Jacobian must be a BlockOperator"); SLIC_ERROR_ROOT_IF(block_jacobian_->NumRowBlocks() != 2 || block_jacobian_->NumColBlocks() != 2, - axom::fmt::format("BlockSchurPreconditioner requires a 2x2 block operator, got {}x{}", - block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); - SLIC_ERROR_ROOT_IF( - mfem_solvers_.size() != 2, - axom::fmt::format("BlockSchurPreconditioner requires exactly 2 solvers, got {}", mfem_solvers_.size())); + std::format("BlockSchurPreconditioner requires a 2x2 block operator, got {}x{}", + block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + SLIC_ERROR_ROOT_IF(mfem_solvers_.size() != 2, + std::format("BlockSchurPreconditioner requires exactly 2 solvers, got {}", mfem_solvers_.size())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); solver_diag_ = std::make_unique(block_offsets_); diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index d7adc577e6..89a26c9a0f 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -123,9 +123,9 @@ bool monolithicizeOperatorIfNeeded(const LinearSolverOptions& linear_options, mf gradient_operator = buildMonolithicMatrix(*block_gradient).release(); SLIC_DEBUG_ROOT( - axom::fmt::format("Automatically monolithicizing block Jacobian for linear solver {} with " - "preconditioner {}", - linearName(linear_options.linear_solver), preconditionerName(linear_options.preconditioner))); + std::format("Automatically monolithicizing block Jacobian for linear solver {} with " + "preconditioner {}", + linearName(linear_options.linear_solver), preconditionerName(linear_options.preconditioner))); return true; } diff --git a/src/smith/numerics/nonlinear_convergence.cpp b/src/smith/numerics/nonlinear_convergence.cpp index 454e6c8118..cd036849d3 100644 --- a/src/smith/numerics/nonlinear_convergence.cpp +++ b/src/smith/numerics/nonlinear_convergence.cpp @@ -9,7 +9,6 @@ #include #include -#include "axom/fmt.hpp" #include "axom/slic.hpp" namespace smith { diff --git a/src/smith/physics/functional_weak_form.hpp b/src/smith/physics/functional_weak_form.hpp index 521f7ef9f1..2eb7a29dba 100644 --- a/src/smith/physics/functional_weak_form.hpp +++ b/src/smith/physics/functional_weak_form.hpp @@ -333,6 +333,9 @@ class FunctionalWeakForm, [[maybe_unused]] const std::vector& quad_fields = {}) const override { validateFields(fields, "jacobian"); + SLIC_ERROR_IF(jacobian_weights.size() != fields.size(), + std::format("jacobian(): jacobian_weights.size()={} but fields.size()={}", jacobian_weights.size(), + fields.size())); SetCurrentTimeInfoRAII clear_current_time_info_on_exit(current_time_info_, time_info); std::unique_ptr J; @@ -474,7 +477,7 @@ class FunctionalWeakForm, { if constexpr (I < sizeof...(InputSpaces)) { using Space = typename std::tuple_element>::type; - validateSpace(*input_mfem_spaces[I], axom::fmt::format("input[{}]", I)); + validateSpace(*input_mfem_spaces[I], std::format("input[{}]", I)); validateInputSpaces(input_mfem_spaces); } } @@ -486,7 +489,7 @@ class FunctionalWeakForm, 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())); + std::format("{}(): field[{}] ('{}')", method_name, I, fields[I]->name())); validateFieldsRecursive(fields, method_name); } } @@ -495,8 +498,8 @@ class FunctionalWeakForm, 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))); + std::format("{}(): fields.size()={} but weak form expects {} InputSpaces", method_name, + fields.size(), sizeof...(InputSpaces))); if constexpr (sizeof...(InputSpaces) > 0) { validateFieldsRecursive<0>(fields, method_name); @@ -516,27 +519,26 @@ class FunctionalWeakForm, 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)); + SLIC_ERROR_ROOT_IF(!is_h1, std::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)); + SLIC_ERROR_ROOT_IF(!is_l2, std::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())); + std::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())); + std::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 diff --git a/src/smith/physics/tests/test_functional_weak_form.cpp b/src/smith/physics/tests/test_functional_weak_form.cpp index 2112439811..d0900eda85 100644 --- a/src/smith/physics/tests/test_functional_weak_form.cpp +++ b/src/smith/physics/tests/test_functional_weak_form.cpp @@ -104,6 +104,8 @@ struct WeakFormFixture : public testing::Test { [](auto /*t_info*/, auto /*x*/, auto n, auto... /*args*/) { return 1.0 * n; }); f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto /*x*/, auto u, auto... /*args*/) { return u; }); + f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto /*x*/, auto /*u*/, auto velocity, + auto /*density*/) { return 0.25 * velocity; }); f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto x, auto... /*args*/) { return 0.5 * x; }); @@ -251,6 +253,30 @@ TEST_F(WeakFormFixture, JvpConsistency) } } +TEST_F(WeakFormFixture, MultipleJacobianWeightsSumArgumentDerivatives) +{ + auto input_fields = getConstFieldPointers(states, params); + const std::vector displacement_weights = {1.0, 0.0, 0.0}; + const std::vector velocity_weights = {0.0, 1.0, 0.0}; + const std::vector combined_weights = {1.0, 1.0, 0.0}; + + auto displacement_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, displacement_weights); + auto velocity_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, velocity_weights); + auto combined_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, combined_weights); + + smith::FiniteElementDual expected(states[DISP].space(), "expected"); + smith::FiniteElementDual velocity_contribution(states[DISP].space(), "velocity_contribution"); + smith::FiniteElementDual actual(states[DISP].space(), "actual"); + displacement_jacobian->Mult(state_tangents[DISP], expected); + velocity_jacobian->Mult(state_tangents[DISP], velocity_contribution); + expected += velocity_contribution; + combined_jacobian->Mult(state_tangents[DISP], actual); + + EXPECT_GT(velocity_contribution.Norml2(), 1.0e-12); + actual -= expected; + EXPECT_NEAR(actual.Norml2(), 0.0, 1.0e-12); +} + TEST_F(WeakFormFixture, ForwardsOriginalTimeInfoToIntegrands) { using TrialSpace = VectorSpace;