diff --git a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp index c6161392a..e038af13e 100644 --- a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp +++ b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp @@ -35,7 +35,7 @@ namespace smith { template auto buildSolidMechanics(std::shared_ptr mesh, std::shared_ptr d_solid_nonlinear_solver, - smith::SecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, std::string physics_name, const std::vector& param_names = {}) { auto graph = std::make_shared(num_checkpoints); diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp index f06563dcb..009ef1396 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -266,6 +266,133 @@ std::vector LinearDifferentiableBlockSolver return u_duals; } +NonlinearDifferentiableBlockSolver::NonlinearDifferentiableBlockSolver(std::unique_ptr s) + : nonlinear_solver_(std::move(s)) +{ +} + +void NonlinearDifferentiableBlockSolver::completeSetup(const std::vector&) +{ + // eventually may need something like: initializeSolver(&nonlinear_solver_->preconditioner(), u); +} + +std::vector NonlinearDifferentiableBlockSolver::solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residual_funcs, + std::function>(const std::vector&)> jacobian_funcs) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_guesses.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_guesses[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_u = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_u->GetBlock(row_i) = *u_guesses[static_cast(row_i)]; + } + + auto block_r = std::make_unique(block_offsets); + + auto residual_op_ = std::make_unique( + block_u->Size(), + [&residual_funcs, num_rows, &u_guesses, &block_r](const mfem::Vector& u_, mfem::Vector& r_) { + const mfem::BlockVector* u = dynamic_cast(&u_); + SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a blocl vector"); + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); + } + auto residuals = residual_funcs(u_guesses); + SLIC_ERROR_IF(!block_r, "Invalid r cast in block differentiable solver to a block vector"); + for (int row_i = 0; row_i < num_rows; ++row_i) { + auto r = residuals[static_cast(row_i)]; + block_r->GetBlock(row_i) = r; + } + r_ = *block_r; + }, + [this, &block_offsets, &u_guesses, jacobian_funcs, num_rows](const mfem::Vector& u_) -> mfem::Operator& { + const mfem::BlockVector* u = dynamic_cast(&u_); + SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a block vector"); + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); + } + block_jac_ = std::make_unique(block_offsets); + matrix_of_jacs_ = jacobian_funcs(u_guesses); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + auto& J = matrix_of_jacs_[static_cast(i)][static_cast(j)]; + if (J) { + block_jac_->SetBlock(i, j, J.get()); + } + } + } + return *block_jac_; + }); + nonlinear_solver_->setOperator(*residual_op_); + nonlinear_solver_->solve(*block_u); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] = block_u->GetBlock(row_i); + } + + return u_guesses; +} + +std::vector NonlinearDifferentiableBlockSolver::solveAdjoint( + const std::vector& u_bars, std::vector>& jacobian_transposed) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_bars.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + std::vector u_duals(static_cast(num_rows)); + for (int row_i = 0; row_i < num_rows; ++row_i) { + u_duals[static_cast(row_i)] = std::make_shared( + u_bars[static_cast(row_i)]->space(), "u_dual_" + std::to_string(row_i)); + } + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_bars[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_ds = std::make_unique(block_offsets); + *block_ds = 0.0; + + auto block_r = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_r->GetBlock(row_i) = *u_bars[static_cast(row_i)]; + } + + auto block_jac = std::make_unique(block_offsets); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + block_jac->SetBlock(i, j, jacobian_transposed[static_cast(i)][static_cast(j)].get()); + } + } + + auto& linear_solver = nonlinear_solver_->linearSolver(); + linear_solver.SetOperator(*block_jac); + linear_solver.Mult(*block_r, *block_ds); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_duals[static_cast(row_i)] = block_ds->GetBlock(row_i); + } + + return u_duals; +} + std::shared_ptr buildDifferentiableLinearSolver(LinearSolverOptions linear_opts, const smith::Mesh& mesh) { @@ -273,11 +400,19 @@ std::shared_ptr buildDifferentiableLinearSolver(Line return std::make_shared(std::move(linear_solver), std::move(precond)); } -std::shared_ptr buildDifferentiableNonlinearSolver( - smith::NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) +std::shared_ptr buildDifferentiableNonlinearSolver(NonlinearSolverOptions nonlinear_opts, + LinearSolverOptions linear_opts, + const smith::Mesh& mesh) +{ + auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); + return std::make_shared(std::move(solid_solver)); +} + +std::shared_ptr buildDifferentiableNonlinearBlockSolver( + NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) { - auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); - return std::make_shared(std::move(solid_solver)); + auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); + return std::make_shared(std::move(solid_solver)); } } // namespace smith diff --git a/src/smith/differentiable_numerics/differentiable_solver.hpp b/src/smith/differentiable_numerics/differentiable_solver.hpp index bf6b738aa..2276a205f 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.hpp +++ b/src/smith/differentiable_numerics/differentiable_solver.hpp @@ -19,6 +19,7 @@ namespace mfem { class Solver; class Vector; class HypreParMatrix; +class BlockOperator; } // namespace mfem namespace smith { @@ -175,6 +176,36 @@ class LinearDifferentiableBlockSolver : public DifferentiableBlockSolver { mutable std::unique_ptr mfem_preconditioner; ///< stored mfem block preconditioner }; +/// @brief Implementation of the DifferentiableBlockSolver interface for the special case of nonlinear solves with +/// linear adjoint solves +class NonlinearDifferentiableBlockSolver : public DifferentiableBlockSolver { + public: + /// @brief Construct from a linear solver and linear block precondition which may be used by the linear solver + NonlinearDifferentiableBlockSolver(std::unique_ptr s); + + /// @overload + void completeSetup(const std::vector& us) override; + + /// @overload + std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residuals, + std::function>(const std::vector&)> jacobians) const override; + + /// @overload + std::vector solveAdjoint(const std::vector& u_bars, + std::vector>& jacobian_transposed) const override; + + mutable std::unique_ptr + block_jac_; ///< Need to hold an instance of a block operator to work with the mfem solver interface + mutable std::vector> + matrix_of_jacs_; ///< Holding vectors of block matrices to that do not going out of scope before the mfem solver + ///< is done with using them in the block_jac_ + + mutable std::unique_ptr + nonlinear_solver_; ///< the nonlinear equation solver used for the forward pass +}; + /// @brief Create a differentiable linear solver /// @param linear_opts linear options struct /// @param mesh mesh @@ -189,4 +220,11 @@ std::shared_ptr buildDifferentiableNonlinearSolve LinearSolverOptions linear_opts, const smith::Mesh& mesh); +/// @brief Create a differentiable nonlinear block solver +/// @param nonlinear_opts nonlinear options struct +/// @param linear_opts linear options struct +/// @param mesh mesh +std::shared_ptr buildDifferentiableNonlinearBlockSolver( + NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh); + } // namespace smith diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp index c148350df..a83819176 100644 --- a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -72,74 +72,56 @@ class DirichletBoundaryConditions { /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. /// @param domain All dofs in this domain have boundary conditions applied to it. - /// @param component component direction to apply boundary condition to if the underlying field is a vector-field. /// @param applied_displacement applied_displacement is a functor which takes time, and a /// smith::tensor corresponding to the spatial coordinate. The functor must return a double. For /// example: [](double t, smith::tensor X) { return 1.0; } template - void setVectorBCs(const Domain& domain, int component, AppliedDisplacementFunction applied_displacement) + void setScalarBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) { - const int field_dim = space_.GetVDim(); - SLIC_ERROR_IF(component >= field_dim, - axom::fmt::format("Trying to set boundary conditions on a field with dim {}, using component {}", - field_dim, component)); auto mfem_coefficient_function = [applied_displacement](const mfem::Vector& X_mfem, double t) { auto X = make_tensor([&X_mfem](int k) { return X_mfem[k]; }); return applied_displacement(t, X); }; auto dof_list = domain.dof_list(&space_); - // scalar ldofs -> vector ldofs - space_.DofsToVDofs(component, dof_list); + space_.DofsToVDofs(static_cast(0), dof_list); auto component_disp_bdr_coef_ = std::make_shared(mfem_coefficient_function); - bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, component); - } - - /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. - /// @param domain All dofs in this domain have boundary conditions applied to it. - /// @param applied_displacement applied_displacement is a functor which takes time, and a - /// smith::tensor corresponding to the spatial coordinate. The functor must return a double. For - /// example: [](double t, smith::tensor X) { return 1.0; } - template - void setScalarBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) - { - setScalarBCs(domain, 0, applied_displacement); + bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, 0); } /// @brief Constrain the dofs of a scalar field over a domain template void setFixedScalarBCs(const Domain& domain) { - setVectorBCs(domain, [](auto, auto) { return 0.0; }); + setScalarBCs(domain, [](auto, auto) { return 0.0; }); } - /// @brief Constrain the dofs of a scalar field over a domain - template - void setFixedVectorBCs(const Domain& domain, int component) + /// @brief Constrain the vector dofs over a domain corresponding to a subset of the vector components + template + void setFixedVectorBCs(const Domain& domain, std::vector components) { - setScalarBCs(domain, component, [](auto, auto) { return 0.0; }); + setVectorBCs(domain, components, [](auto, auto) { return smith::tensor{}; }); } - /// @brief Constrain the vector dofs over a domain corresponding to a subset of the vector components - template - void setFixedVectorBCs(const Domain& domain, std::vector components) + /// @brief Constrain the dofs of a scalar field over a domain + template + void setFixedVectorBCs(const Domain& domain, int component) { - setVectorBCs(domain, components, [](auto, auto) { return smith::tensor{}; }); + std::vector components{component}; + setVectorBCs(domain, components); } /// @brief Constrain all the vector dofs over a domain - template + template void setFixedVectorBCs(const Domain& domain) { - const int field_dim = space_.GetVDim(); - SLIC_ERROR_IF(field_dim != spatial_dim, - "Vector boundary conditions current only work if they match the spatial dimension"); + SLIC_ERROR_IF(field_dim != space_.GetVDim(), "Vector boundary condition field_dim does not match the fields vdim"); std::vector components(static_cast(field_dim)); for (int component = 0; component < field_dim; ++component) { components[static_cast(component)] = component; } - setFixedVectorBCs(domain, components); + setFixedVectorBCs(domain, components); } /// @brief Return the smith BoundaryConditionManager diff --git a/src/smith/differentiable_numerics/field_state.hpp b/src/smith/differentiable_numerics/field_state.hpp index e9528faf8..4f2322c30 100644 --- a/src/smith/differentiable_numerics/field_state.hpp +++ b/src/smith/differentiable_numerics/field_state.hpp @@ -212,22 +212,30 @@ inline std::vector spaces(const std::vector< }; /// @brief Get a vector of FieldPtr or DualFieldPtr from a vector of FieldState -inline std::vector getFieldPointers(std::vector& states) +inline std::vector getFieldPointers(std::vector& states, + std::vector params = {}) { std::vector pointers; for (auto& t : states) { pointers.push_back(t.get().get()); } + for (auto& p : params) { + pointers.push_back(p.get().get()); + } return pointers; } /// @brief Get a vector of ConstFieldPtr or ConstDualFieldPtr from a vector of FieldState -inline std::vector getConstFieldPointers(const std::vector& states) +inline std::vector getConstFieldPointers(const std::vector& states, + const std::vector& params = {}) { std::vector pointers; for (auto& t : states) { pointers.push_back(t.get().get()); } + for (auto& p : params) { + pointers.push_back(p.get().get()); + } return pointers; } diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 2f31498fc..37777a173 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -272,7 +272,7 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers) + const std::vector& bc_managers) { SMITH_MARK_FUNCTION; size_t num_rows_ = residual_evals.size(); @@ -305,7 +305,6 @@ std::vector block_solve(const std::vector& residual_evals } } allFields.push_back(shape_disp); - struct ZeroDualVectors { std::vector operator()(const std::vector& fs) { @@ -319,15 +318,14 @@ std::vector block_solve(const std::vector& residual_evals FieldVecState sol = shape_disp.create_state, std::vector>(allFields, ZeroDualVectors()); - sol.set_eval([=](const gretl::UpstreamStates& upstreams, gretl::DownstreamState& downstream) { SMITH_MARK_BEGIN("solve forward"); const size_t num_rows = num_state_inputs.size(); std::vector> input_fields(num_rows); - SLIC_ERROR_IF(num_rows != num_param_inputs.size(), "row count for params and columns are inconsistent"); + SLIC_ERROR_IF(num_rows != num_param_inputs.size(), "row count for params and states are inconsistent"); - // The order of inputs in upstreams seems to be: - // states of residual 0 -> states of residual 1 -> params of residual 0 -> params of residual 1 + // The order of inputs in upstreams is: + // states of residual 0, states of residual 1, ... , params of residual 0, params of residual 1, ... size_t field_count = 0; for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t state_i = 0; state_i < num_state_inputs[row_i]; ++state_i) { @@ -342,12 +340,23 @@ 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) { - diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][block_indices[row_i][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)"); + diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][prime_unknown_row_i]); + } + + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + FEFieldPtr primal_field_row_i = diagonal_fields[row_i]; + applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); } for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + 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]; + } } } @@ -363,13 +372,11 @@ std::vector block_solve(const std::vector& residual_evals *primal_field_row_i = *unknowns[row_i]; applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); } - for (size_t row_i = 0; row_i < num_rows; ++row_i) { residuals[row_i] = residual_evals[row_i]->residual(time_info, shape_disp_ptr.get(), getConstFieldPointers(input_fields[row_i])); residuals[row_i].SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); } - return residuals; }; @@ -389,29 +396,37 @@ std::vector block_solve(const std::vector& residual_evals std::vector tangent_weights(row_field_inputs.size(), 0.0); for (size_t col_j = 0; col_j < num_rows; ++col_j) { size_t field_index_to_diff = block_indices[row_i][col_j]; - tangent_weights[field_index_to_diff] = 1.0; - auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), - getConstFieldPointers(row_field_inputs), tangent_weights); - jacobians[row_i].emplace_back(std::move(jac_ij)); - tangent_weights[field_index_to_diff] = 0.0; + if (field_index_to_diff != invalid_block_index) { + tangent_weights[field_index_to_diff] = 1.0; + auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), + getConstFieldPointers(row_field_inputs), tangent_weights); + jacobians[row_i].emplace_back(std::move(jac_ij)); + tangent_weights[field_index_to_diff] = 0.0; + } else { + jacobians[row_i].emplace_back(nullptr); + } } } // Apply BCs to the block system for (size_t row_i = 0; row_i < num_rows; ++row_i) { - mfem::HypreParMatrix* Jii = - jacobians[row_i][row_i]->EliminateRowsCols(bc_managers[row_i]->allEssentialTrueDofs()); - delete Jii; + if (jacobians[row_i][row_i]) { + jacobians[row_i][row_i]->EliminateBC(bc_managers[row_i]->allEssentialTrueDofs(), + mfem::Operator::DiagonalPolicy::DIAG_ONE); + } for (size_t col_j = 0; col_j < num_rows; ++col_j) { if (col_j != row_i) { - jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); - mfem::HypreParMatrix* Jji = - jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); - delete Jji; + if (jacobians[row_i][col_j]) { + jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); + } + if (jacobians[col_j][row_i]) { + mfem::HypreParMatrix* Jji = + jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jji; + } } } } - return jacobians; }; @@ -558,7 +573,11 @@ std::vector block_solve(const std::vector& residual_evals for (size_t i = 0; i < num_rows_; ++i) { FieldState s = gretl::create_state( zero_dual_from_state(), - [i](const std::vector& sols) { return std::make_shared(*sols[i]); }, + [i](const std::vector& sols) { + auto state_copy = std::make_shared(sols[i]->space(), sols[i]->name()); + *state_copy = *sols[i]; + return state_copy; + }, [i](const std::vector&, const FEFieldPtr&, std::vector& sols_, const FEDualPtr& output_) { *sols_[i] += *output_; }, sol); diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index 3529a3006..eebfbdd40 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -24,6 +24,9 @@ class DifferentiableBlockSolver; class BoundaryConditionManager; class DirichletBoundaryConditions; +/// @brief magic number for represending a field which is not an argument of the weak form. +static constexpr size_t invalid_block_index = std::numeric_limits::max() - 1; + /// @brief Solve a nonlinear system of equations as defined by the weak form, assuming that the field indexed by /// unknown_index is the unknown field /// @param residual_eval The weak form which defines the equations to be solved @@ -61,6 +64,6 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers); + const std::vector& bc_managers); } // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp index 7a9b888e4..67e2b66a9 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -16,7 +16,7 @@ namespace smith { SolidMechanicsStateAdvancer::SolidMechanicsStateAdvancer( std::shared_ptr solver, std::shared_ptr vector_bcs, std::shared_ptr solid_dynamic_weak_forms, - smith::SecondOrderTimeIntegrationRule time_rule) + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule) : solver_(solver), vector_bcs_(vector_bcs), solid_dynamic_weak_forms_(solid_dynamic_weak_forms), @@ -31,25 +31,23 @@ std::vector SolidMechanicsStateAdvancer::advanceState(const TimeInfo { std::vector states_old = states_old_; if (time_info.cycle() == 0) { - states_old[ACCELERATION] = solve(*solid_dynamic_weak_forms_->quasi_static_weak_form, shape_disp, states_old_, + states_old[ACCELERATION] = solve(*solid_dynamic_weak_forms_->final_reaction_weak_form, shape_disp, states_old_, params, time_info, *solver_, *vector_bcs_, ACCELERATION); } - TimeInfo final_time_info = time_info.endTimeInfo(); - std::vector solid_inputs{states_old[DISPLACEMENT], states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]}; auto displacement = solve(*solid_dynamic_weak_forms_->time_discretized_weak_form, shape_disp, solid_inputs, params, - final_time_info, *solver_, *vector_bcs_); + time_info, *solver_, *vector_bcs_); std::vector states = states_old; states[DISPLACEMENT] = displacement; - states[VELOCITY] = time_rule_.derivative(final_time_info, displacement, states_old[DISPLACEMENT], - states_old[VELOCITY], states_old[ACCELERATION]); - states[ACCELERATION] = time_rule_.second_derivative(final_time_info, displacement, states_old[DISPLACEMENT], - states_old[VELOCITY], states_old[ACCELERATION]); + states[VELOCITY] = + time_rule_.dot(time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]); + states[ACCELERATION] = time_rule_.ddot(time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], + states_old[ACCELERATION]); return states; } @@ -61,7 +59,7 @@ std::vector SolidMechanicsStateAdvancer::computeReactions(const T { std::vector solid_inputs{states[DISPLACEMENT], states[VELOCITY], states[ACCELERATION]}; solid_inputs.insert(solid_inputs.end(), params.begin(), params.end()); - return {evaluateWeakForm(solid_dynamic_weak_forms_->quasi_static_weak_form, time_info, shape_disp, solid_inputs, + return {evaluateWeakForm(solid_dynamic_weak_forms_->final_reaction_weak_form, time_info, shape_disp, solid_inputs, states[DISPLACEMENT])}; } diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp index b81c514d8..121872491 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp @@ -40,7 +40,7 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { SolidMechanicsStateAdvancer(std::shared_ptr solid_solver, std::shared_ptr vector_bcs, std::shared_ptr solid_dynamic_weak_forms, - SecondOrderTimeIntegrationRule time_rule); + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule); /// State enum for indexing convenience enum STATE @@ -70,7 +70,7 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { /// application you will get back: shape_disp, states, params, time, and solid_mechanics_weak_form template static auto buildWeakFormAndStates(const std::shared_ptr& mesh, const std::shared_ptr& graph, - SecondOrderTimeIntegrationRule time_rule, std::string physics_name, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, std::string physics_name, const std::vector& param_names, double initial_time = 0.0) { auto shape_disp = createFieldState(*graph, ShapeDispSpace{}, physics_name + "_shape_displacement", mesh->tag()); @@ -111,8 +111,9 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { std::shared_ptr solid_dynamic_weak_forms_; ///< Solid mechanics time discretized weak forms, user must setup the appropriate ///< integrals. Has both the time discretized and the undiscretized weak forms. - SecondOrderTimeIntegrationRule time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, - ///< given the current predicted u and the previous u, u_dot, u_dot_dot + ImplicitNewmarkSecondOrderTimeIntegrationRule + time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, + ///< given the current predicted u and the previous u, u_dot, u_dot_dot }; } // namespace smith diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index dc5925a79..06681ec85 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -10,6 +10,7 @@ set(differentiable_numerics_test_source test_field_state.cpp test_explicit_dynamics.cpp test_solid_mechanics_state_advancer.cpp + test_thermo_mechanics.cpp ) smith_add_tests( SOURCES ${differentiable_numerics_test_source} diff --git a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp index a49e2a737..e22e123fb 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp @@ -90,7 +90,7 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule; auto [physics, solid_weak_form, bcs] = buildSolidMechanics( @@ -181,20 +181,16 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) EXPECT_NEAR(0.0, u_err, 1e-14); } -TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) +auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) { - SMITH_MARK_FUNCTION; - - std::string physics_name = "solid"; - + size_t num_checkpoints = 200; std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); - auto [physics, solid_weak_form, bcs] = buildSolidMechanics( - mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); + mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule(), num_checkpoints, + physics_name, {"bulk", "shear"}); bcs->setFixedVectorBCs(mesh->domain("right")); bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { @@ -221,7 +217,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) auto shape_disp = physics->getShapeDispFieldState(); auto params = physics->getFieldParams(); - auto initial_states = physics->getInitialFieldStates(); + auto states = physics->getInitialFieldStates(); params[0].get()->setFromFieldFunction([=](smith::tensor) { double scaling = 1.0; @@ -235,6 +231,15 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) physics->resetStates(); + return std::make_tuple(physics, shape_disp, states, params); +} + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) +{ + SMITH_MARK_FUNCTION; + std::string physics_name = "solid"; + auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); + auto pv_writer = smith::createParaviewWriter(*mesh, physics->getFieldStatesAndParamStates(), physics_name); pv_writer.write(0, physics->time(), physics->getFieldStatesAndParamStates()); for (size_t m = 0; m < num_steps_; ++m) { @@ -311,56 +316,8 @@ void adjointBackward(std::shared_ptr physics, smith::FiniteElementD TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) { SMITH_MARK_FUNCTION; - std::string physics_name = "solid"; - - std::shared_ptr d_solid_nonlinear_solver = - buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - - smith::SecondOrderTimeIntegrationRule time_rule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); - - auto [physics, solid_weak_form, bcs] = - buildSolidMechanics( - mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); - - bcs->setFixedVectorBCs(mesh->domain("right")); - bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { - auto bc = 0.0 * X; - bc[0] = 0.01 * t; - bc[1] = -0.05 * t; - return bc; - }); - - double E = 100.0; - double nu = 0.25; - auto K = E / (3.0 * (1.0 - 2 * nu)); - auto G = E / (2.0 * (1.0 + nu)); - using MaterialType = solid_mechanics::ParameterizedNeoHookeanSolid; - MaterialType material{.density = 10.0, .K0 = K, .G0 = G}; - - solid_weak_form->addBodyIntegral( - smith::DependsOn<0, 1>{}, mesh->entireBodyName(), - [material](const auto& /*time_info*/, auto /*X*/, auto u, auto /*v*/, auto a, auto bulk, auto shear) { - MaterialType::State state; - auto pk_stress = material(state, get(u), bulk, shear); - return smith::tuple{get(a) * material.density, pk_stress}; - }); - - auto shape_disp = physics->getShapeDispFieldState(); - auto params = physics->getFieldParams(); - auto states = physics->getInitialFieldStates(); - - params[0].get()->setFromFieldFunction([=](smith::tensor) { - double scaling = 1.0; - return scaling * material.K0; - }); - - params[1].get()->setFromFieldFunction([=](smith::tensor) { - double scaling = 1.0; - return scaling * material.G0; - }); - - physics->resetStates(); + auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); double qoi = integrateForward(physics, num_steps_, dt_); SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp new file mode 100644 index 000000000..9ab89d1b0 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -0,0 +1,432 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Smith Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "gtest/gtest.h" + +#include "smith/smith_config.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/solver_config.hpp" + +#include "smith/physics/state/state_manager.hpp" + +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/differentiable_numerics/differentiable_test_utils.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" + +namespace smith { + +/** + * @brief Compute Green's strain from the displacement gradient + */ +template +auto greenStrain(const tensor& grad_u) +{ + return 0.5 * (grad_u + transpose(grad_u) + dot(transpose(grad_u), grad_u)); +} + +/// @brief Green-Saint Venant isotropic thermoelastic model +struct GreenSaintVenantThermoelasticMaterial { + double density; ///< density + double E; ///< Young's modulus + double nu; ///< Poisson's ratio + double C_v; ///< volumetric heat capacity + double alpha; ///< thermal expansion coefficient + double theta_ref; ///< datum temperature for thermal expansion + double kappa; ///< thermal conductivity + + using State = Empty; + + /** + * @brief Evaluate constitutive variables for thermomechanics + * + * @tparam T1 Type of the displacement gradient components (number-like) + * @tparam T2 Type of the temperature (number-like) + * @tparam T3 Type of the temperature gradient components (number-like) + * + * @param[in] grad_u Displacement gradient + * @param[in] theta Temperature + * @param[in] grad_theta Temperature gradient + * @param[in,out] state State variables for this material + * + * @return[out] tuple of constitutive outputs. Contains the + * First Piola stress, the volumetric heat capacity in the reference + * configuration, the heat generated per unit volume during the time + * step (units of energy), and the referential heat flux (units of + * energy per unit time and per unit area). + */ + template + auto operator()(double, State&, const tensor& grad_u, const tensor& grad_v, T3 theta, + const tensor& grad_theta) const + { + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + const auto Eg = greenStrain(grad_u); + const auto trEg = tr(Eg); + + // stress + static constexpr auto I = Identity(); + const auto S = 2.0 * G * dev(Eg) + K * (trEg - dim * alpha * (theta - theta_ref)) * I; + auto F = grad_u + I; + const auto Piola = dot(F, S); + + // internal heat power + auto greenStrainRate = + 0.5 * (grad_v + transpose(grad_v) + dot(transpose(grad_v), grad_u) + dot(transpose(grad_u), grad_v)); + const auto s0 = -dim * K * alpha * (theta + 273.1) * tr(greenStrainRate); + + // heat flux + const auto q0 = -kappa * grad_theta; + + return smith::tuple{Piola, C_v, s0, q0}; + } +}; + +smith::LinearSolverOptions linear_options{.linear_solver = smith::LinearSolver::Strumpack, + .relative_tol = 1e-8, + .absolute_tol = 1e-8, + .max_iterations = 200, + .print_level = 0}; + +smith::NonlinearSolverOptions nonlinear_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, + .relative_tol = 1.9e-6, + .absolute_tol = 1.0e-10, + .max_iterations = 500, + .max_line_search_iterations = 50, + .print_level = 2}; + +static constexpr int dim = 3; +static constexpr int order = 1; + +struct SolidMechanicsMeshFixture : public testing::Test { + double length = 1.0; + double width = 0.04; + int num_elements_x = 12; + int num_elements_y = 2; + int num_elements_z = 2; + double elem_size = length / num_elements_x; + + void SetUp() + { + smith::StateManager::initialize(datastore_, "solid"); + auto mfem_shape = mfem::Element::QUADRILATERAL; + mesh_ = std::make_shared( + mfem::Mesh::MakeCartesian3D(num_elements_x, num_elements_y, num_elements_z, mfem_shape, length, width, width), + "mesh", 0, 0); + mesh_->addDomainOfBoundaryElements("left", smith::by_attr(3)); + mesh_->addDomainOfBoundaryElements("right", smith::by_attr(5)); + } + + static constexpr double total_simulation_time_ = 1.1; + static constexpr size_t num_steps_ = 4; + static constexpr double dt_ = total_simulation_time_ / num_steps_; + + axom::sidre::DataStore datastore_; + std::shared_ptr mesh_; +}; + +template +struct FieldType { + FieldType(std::string n, int unknown_index_ = -1) : name(n), unknown_index(unknown_index_) {} + std::string name; + int unknown_index; +}; + +struct FieldStore { + FieldStore(std::shared_ptr mesh, size_t storage_size = 50) + : mesh_(mesh), data_store_(std::make_shared(storage_size)) + { + } + + template + void addShapeDisp(FieldType type) + { + shape_disp_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); + } + + template + std::shared_ptr& addUnknown(FieldType& type) + { + type.unknown_index = static_cast(num_unknowns_); + to_fields_index_[type.name] = fields_.size(); + to_unknown_index_[type.name] = num_unknowns_; + FieldState new_field = smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag()); + fields_.push_back(new_field); + ++num_unknowns_; + boundary_conditions_.push_back( + std::make_shared(mesh_->mfemParMesh(), new_field.get()->space())); + SLIC_ERROR_IF(num_unknowns_ != boundary_conditions_.size(), + "Inconcistency between num unknowns and boundary condition size"); + return boundary_conditions_.back(); + } + + template + auto addDerived(FieldType, std::string name) + { + to_fields_index_[name] = fields_.size(); + fields_.push_back(smith::createFieldState(*data_store_, Space{}, name, mesh_->tag())); + return FieldType(name); + } + + void addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index) + { + FieldLabel argument_name_and_index{.field_name = argument_name, .field_index_in_residual = argument_index}; + if (weak_form_name_to_unknown_name_index_.count(weak_form_name)) { + weak_form_name_to_unknown_name_index_.at(weak_form_name).push_back(argument_name_and_index); + } else { + weak_form_name_to_unknown_name_index_[weak_form_name] = std::vector{argument_name_and_index}; + } + } + + void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) + { + size_t field_index = to_fields_index_.at(argument_name); + if (weak_form_name_to_field_indices_.count(weak_form_name)) { + weak_form_name_to_field_indices_.at(weak_form_name).push_back(field_index); + } else { + weak_form_name_to_field_indices_[weak_form_name] = std::vector{field_index}; + } + SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_indices_.at(weak_form_name).size(), + "Invalid order for adding weak form arguments."); + } + + void printMap() + { + for (auto& keyval : weak_form_name_to_unknown_name_index_) { + std::cout << "for residual: " << keyval.first << " "; + for (auto& name_index : keyval.second) { + std::cout << "arg " << name_index.field_name << " at " << name_index.field_index_in_residual << ", "; + } + std::cout << std::endl; + } + } + + std::vector> indexMap(const std::vector& residual_names) const + { + std::vector> block_indices(residual_names.size()); + + for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { + std::vector& res_indices = block_indices[res_i]; + res_indices = std::vector(num_unknowns_, invalid_block_index); + const std::string& res_name = residual_names[res_i]; + const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); + + for (const auto& field_name_and_arg_index : arg_info) { + const std::string field_name = field_name_and_arg_index.field_name; + size_t unknown_index = to_unknown_index_.at(field_name); + SLIC_ASSERT(unknown_index < num_unknowns_); + res_indices[unknown_index] = field_name_and_arg_index.field_index_in_residual; + } + } + + return block_indices; + } + + std::vector getBoundaryConditionManagers() const + { + std::vector bcs; + for (auto& bc : boundary_conditions_) { + bcs.push_back(&bc->getBoundaryConditionManager()); + } + return bcs; + } + + size_t getFieldIndex(const std::string& field_name) const { return to_fields_index_.at(field_name); } + + const FieldState& getField(const std::string& field_name) const + { + size_t field_index = getFieldIndex(field_name); + return fields_[field_index]; + } + + void setField(const std::string& field_name, FieldState updated_field) + { + size_t field_index = getFieldIndex(field_name); + fields_[field_index] = updated_field; + } + + const FieldState& getShapeDisp() const { return shape_disp_[0]; } + + const std::vector& getAllFields() const { return fields_; } + + std::vector getFields(const std::string& weak_form_name) const + { + auto unknown_field_indices = weak_form_name_to_field_indices_.at(weak_form_name); + std::vector fields_for_residual; + for (auto& i : unknown_field_indices) { + fields_for_residual.push_back(fields_[i]); + } + return fields_for_residual; + } + + const std::shared_ptr& getMesh() const { return mesh_; } + + private: + std::shared_ptr mesh_; + std::shared_ptr data_store_; + + std::vector shape_disp_; + std::vector fields_; + std::map to_fields_index_; + + size_t num_unknowns_ = 0; + std::map to_unknown_index_; + std::vector> boundary_conditions_; + + struct FieldLabel { + std::string field_name; + size_t field_index_in_residual; + }; + + std::map> weak_form_name_to_unknown_name_index_; + + std::map> weak_form_name_to_field_indices_; +}; + +template +void createSpaces(const std::string& weak_form_name, FieldStore& field_store, + std::vector& spaces, size_t arg_num, FirstType type, + Types... types) +{ + SLIC_ERROR_IF(spaces.size() != arg_num, "Error creating spaces recursively"); + spaces.push_back(&field_store.getField(type.name).get()->space()); + field_store.addWeakFormArg(weak_form_name, type.name, arg_num); + if (type.unknown_index >= 0) { + field_store.addWeakFormUnknownArg(weak_form_name, type.name, arg_num); + } + if constexpr (sizeof...(types) > 0) { + createSpaces(weak_form_name, field_store, spaces, arg_num + 1, types...); + } +} + +template +auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, + FieldType... field_types) +{ + const mfem::ParFiniteElementSpace& test_space = field_store.getField(test_type.name).get()->space(); + std::vector input_spaces; + createSpaces(name, field_store, input_spaces, 0, field_types...); + return std::make_shared>>( + name, field_store.getMesh(), test_space, input_spaces); +} + +std::vector solve(const std::vector& weak_forms, const FieldStore& field_store, + const DifferentiableBlockSolver* solver, const TimeInfo& time_info) +{ + std::vector weak_form_names; + for (const auto& wf : weak_forms) { + weak_form_names.push_back(wf->name()); + } + std::vector> index_map = field_store.indexMap(weak_form_names); + + std::vector> inputs; + for (size_t i = 0; i < weak_forms.size(); ++i) { + std::string wf_name = weak_forms[i]->name(); + std::vector fields_for_wk = field_store.getFields(wf_name); + inputs.push_back(fields_for_wk); + } + std::vector> params(weak_forms.size()); + + return block_solve(weak_forms, index_map, field_store.getShapeDisp(), inputs, params, time_info, solver, + field_store.getBoundaryConditionManagers()); +} + +TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) +{ + SMITH_MARK_FUNCTION; + + FieldType> shape_disp_type("shape_displacement"); + FieldType> disp_type("displacement"); + FieldType> temperature_type("temperature"); + + std::shared_ptr d_nonlinear_solver = + buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + + double rho = 1.0; + double E = 100.0; + double nu = 0.25; + double c = 1.0; + double alpha = 1.0e-3; + double theta_ref = 0.0; + double k = 1.0; + GreenSaintVenantThermoelasticMaterial material{rho, E, nu, c, alpha, theta_ref, k}; + + FieldStore field_store(mesh_, 100); + + field_store.addShapeDisp(shape_disp_type); + + std::shared_ptr& disp_bc = field_store.addUnknown(disp_type); + disp_bc->setVectorBCs(mesh_->domain("left"), [](double t, smith::tensor X) { + auto bc = 0.0 * X; + bc[0] = 0.01 * t; + return bc; + }); + disp_bc->setFixedVectorBCs(mesh_->domain("right")); + + auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); + + std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); + temperature_bc->setFixedScalarBCs(mesh_->domain("left")); + temperature_bc->setFixedScalarBCs(mesh_->domain("right")); + + auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); + + QuasiStaticFirstOrderTimeIntegrationRule disp_time_rule; + BackwardEulerFirstOrderTimeIntegrationRule temperature_time_rule; + + auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, + temperature_type, temperature_old_type); + + solid_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto disp, auto disp_old, + auto temperature, auto temperature_old) { + auto u = disp_time_rule.value(t_info, disp, disp_old); + auto v = disp_time_rule.dot(t_info, disp, disp_old); + auto T = temperature_time_rule.value(t_info, temperature, temperature_old); + GreenSaintVenantThermoelasticMaterial::State state; + auto [pk, C_v, s0, q0] = + material(t_info.dt(), state, get(u), get(v), get(T), get(T)); + return smith::tuple{smith::zero{}, pk}; + }); + + auto thermal_weak_form = createWeakForm("thermal_flux", temperature_type, field_store, temperature_type, + temperature_old_type, disp_type, disp_old_type); + + thermal_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto temperature, + auto temperature_old, auto disp, auto disp_old) { + GreenSaintVenantThermoelasticMaterial::State state; + auto u = disp_time_rule.value(t_info, disp, disp_old); + auto v = disp_time_rule.dot(t_info, disp, disp_old); + auto T = temperature_time_rule.value(t_info, temperature, temperature_old); + auto T_dot = temperature_time_rule.dot(t_info, temperature, temperature_old); + auto [pk, C_v, s0, q0] = + material(t_info.dt(), state, get(u), get(v), get(T), get(T)); + auto dT_dt = get(T_dot); + return smith::tuple{C_v * dT_dt - s0, -q0}; + }); + + thermal_weak_form->addBodySource(smith::DependsOn<>(), mesh_->entireBodyName(), + [](auto /*t*/, auto /* x */) { return 100.0; }); + + std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; + std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); + + // auto states = field_store.getFields(); + + EXPECT_EQ(0, 0); +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index bde1c31fe..a98c6d307 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -55,11 +55,19 @@ class TimeDiscretizedWeakForm + void addBodyIntegralImpl(std::string body_name, BodyForceType body_integral, + std::integer_sequence) + { + addBodyIntegral(DependsOn{}, body_name, body_integral); + } + /// @overload template void addBodyIntegral(std::string body_name, BodyForceType body_integral) { - addBodyIntegral(DependsOn<>{}, body_name, body_integral); + addBodyIntegralImpl(body_name, body_integral, std::make_integer_sequence{}); } }; @@ -67,10 +75,12 @@ class TimeDiscretizedWeakForm time_discretized_weak_form; ///< this publically available abstract weak form is a + std::shared_ptr time_discretized_weak_form; ///< this publicaly available abstract weak form is a ///< functions of the current u, u_old, v_old, and a_old, - std::shared_ptr quasi_static_weak_form; ///< this publically available abstract weak form is structly a - ///< function of the current u, v, and a (no time discretization) + std::shared_ptr + final_reaction_weak_form; ///< this publicaly available abstract weak form is structly a + ///< function of the current u, v, and a (no time discretization) + ///< its main purpose is to compute reaction forces after the solve is completed }; template > @@ -78,8 +88,8 @@ class SecondOrderTimeDiscretizedWeakForm; /// @brief Useful for time-discretized PDEs of second order (involves for first and second derivatives of time). Users /// write q-functions in terns of u, u_dot, u_dot_dot, and the weak form is transformed by the -/// SecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, u_dot_dot_old, with u as -/// the distinct unknown for the time discretized system. +/// ImplicitNewmarkSecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, +/// u_dot_dot_old, with u as the distinct unknown for the time discretized system. /// @tparam spatial_dim Spatial dimension, 2 or 3. /// @tparam OutputSpace The space corresponding to the output residual for the weak form (test-space). /// @tparam TrialInputSpace The space corresponding to the predicted solution u, i.e., the trial solution, the unique @@ -93,13 +103,12 @@ class SecondOrderTimeDiscretizedWeakForm>; ///< using - using QuasiStaticWeakFormT = - TimeDiscretizedWeakForm>; ///< using + TimeDiscretizedWeakForm>; ///< using + using FinalReactionFormT = TimeDiscretizedWeakForm>; ///< using /// @brief Constructor SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr mesh, - SecondOrderTimeIntegrationRule time_rule, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, const mfem::ParFiniteElementSpace& output_mfem_space, const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces) : time_rule_(time_rule) @@ -110,9 +119,9 @@ class SecondOrderTimeDiscretizedWeakForm(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); - quasi_static_weak_form = quasi_static_weak_form_; + final_reaction_weak_form_ = + std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); + final_reaction_weak_form = final_reaction_weak_form_; } /// @overload @@ -126,11 +135,11 @@ class SecondOrderTimeDiscretizedWeakFormaddBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, body_name, - integrand); + final_reaction_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, + body_name, integrand); } /// @overload @@ -144,11 +153,11 @@ class SecondOrderTimeDiscretizedWeakForm time_discretized_weak_form_; ///< fully templated time discretized weak form (with time integration rule injected ///< into the q-function) - std::shared_ptr - quasi_static_weak_form_; ///< fully template underlying weak form (no time integration included, a function of - ///< current u, v, and a) + std::shared_ptr + final_reaction_weak_form_; ///< fully template underlying weak form (no time integration included, a function of + ///< current u, v, and a) - SecondOrderTimeIntegrationRule time_rule_; ///< encodes the time integration rule + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_; ///< encodes the time integration rule }; } // namespace smith diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index dc4e0cd41..79a78a67e 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -34,27 +34,25 @@ struct BackwardEulerFirstOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto derivative(const TimeInfo& t, const T1& field_new, const T2& field_old) const + SMITH_HOST_DEVICE auto dot(const TimeInfo& t, const T1& field_new, const T2& field_old) const { return (1.0 / t.dt()) * (field_new - field_old); } }; -/// @brief Options for second order time integration methods -enum class SecondOrderTimeIntegrationMethod -{ - IMPLICIT_NEWMARK, /// implicit newmark discretization - QUASI_STATIC /// quasi-static, specifies current field, velocity is central difference (for quasi-static artificial - /// viscosity), and acceleration is lagged for cases where it is set to some fixed value in time. -}; +/// Altenative name for Backward Euler which makes sense when restricting what are typically second order odes, +/// for example transient solid mechanics, to the quasi-static approximation. It happens that the implementation is +/// identical to backward-Euler applied to first order systems as we want to be able to capture current velocity +/// dependencies. +using QuasiStaticFirstOrderTimeIntegrationRule = BackwardEulerFirstOrderTimeIntegrationRule; /// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). /// When solving f(u, u_dot, u_dot_dot, t) = 0 /// this class provides the current discrete approximation for u, u_dot, and u_dot_dot as a function of /// (u^{n+1},u^n,u_dot^n,u_dot_dot^n). -struct SecondOrderTimeIntegrationRule { +struct ImplicitNewmarkSecondOrderTimeIntegrationRule { /// @brief Constructor - SecondOrderTimeIntegrationRule(SecondOrderTimeIntegrationMethod method) : method_(method) {} + ImplicitNewmarkSecondOrderTimeIntegrationRule() {} /// @brief evaluate value of the ode state as used by the integration rule template @@ -67,24 +65,58 @@ struct SecondOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, - [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, - [[maybe_unused]] const T4& accel_old) const + SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const { return (2.0 / t.dt()) * (field_new - field_old) - velo_old; } /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto second_derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, - [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, - [[maybe_unused]] const T4& accel_old) const + SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const { auto dt = t.dt(); return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; } +}; + +/// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). +/// When solving f(u, u_dot, u_dot_dot, t) = 0 +/// this class provides the current discrete approximation for u, u_dot, and u_dot_dot as a function of +/// (u^{n+1},u^n,u_dot^n,u_dot_dot^n). +struct QuasiStaticSecondOrderTimeIntegrationRule { + /// @brief Constructor + QuasiStaticSecondOrderTimeIntegrationRule() {} - SecondOrderTimeIntegrationMethod method_; ///< method specifying time integration rule to inject into the q-function. + /// @brief evaluate value of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto value([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return field_new; + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return (1.0 / t.dt()) * (field_new - field_old); + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return accel_old; + } }; } // namespace smith diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 6c24646dc..977c4ddaf 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -1019,20 +1019,19 @@ void StrumpackSolver::SetOperator(const mfem::Operator& op) #endif -std::unique_ptr buildNonlinearSolver(const NonlinearSolverOptions& nonlinear_opts, +std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions& linear_opts, mfem::Solver& prec, MPI_Comm comm) { std::unique_ptr nonlinear_solver; if (nonlinear_opts.nonlin_solver == NonlinearSolver::Newton) { - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "Newton's method does not support nonzero min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "Newton's method does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm, nonlinear_opts); - // nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::LBFGS) { - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "LBFGS does not support nonzero min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "LBFGS does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::NewtonLineSearch) { nonlinear_solver = std::make_unique(comm, nonlinear_opts); @@ -1052,9 +1051,8 @@ std::unique_ptr buildNonlinearSolver(const NonlinearSolverOp // KINSOL else { #ifdef SMITH_USE_SUNDIALS - - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "kinsol solvers do not support min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "kinsol solvers do not support min_iterations"); int kinsol_strat = KIN_NONE; diff --git a/src/smith/numerics/equation_solver.hpp b/src/smith/numerics/equation_solver.hpp index 83bc4673f..33dcc4262 100644 --- a/src/smith/numerics/equation_solver.hpp +++ b/src/smith/numerics/equation_solver.hpp @@ -271,7 +271,7 @@ std::unique_ptr buildMonolithicMatrix(const mfem::BlockOpe * @param comm The MPI communicator for the supplied nonlinear operators and HypreParVectors * @return The constructed nonlinear solver */ -std::unique_ptr buildNonlinearSolver(const NonlinearSolverOptions& nonlinear_opts, +std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions& linear_opts, mfem::Solver& preconditioner, MPI_Comm comm = MPI_COMM_WORLD); diff --git a/src/smith/numerics/solver_config.hpp b/src/smith/numerics/solver_config.hpp index 3ab439046..56890a4b3 100644 --- a/src/smith/numerics/solver_config.hpp +++ b/src/smith/numerics/solver_config.hpp @@ -431,7 +431,7 @@ struct NonlinearSolverOptions { int max_iterations = 20; /// Maximum line search cutbacks - int max_line_search_iterations = 0; + int max_line_search_iterations = 8; /// Debug print level int print_level = 0; diff --git a/src/smith/physics/common.hpp b/src/smith/physics/common.hpp index d65f78b6d..cde66b33d 100644 --- a/src/smith/physics/common.hpp +++ b/src/smith/physics/common.hpp @@ -22,8 +22,11 @@ struct TimeInfo { { } - /// @brief accessor for time - double time() const { return time_.first; } + /// @brief accessor for time previous time + double time_old() const { return time_.first; } + + /// @brief accessor for the current time + double time() const { return time_.first + dt_.first; } /// @brief accessor for dt double dt() const { return dt_.first; } @@ -31,9 +34,7 @@ struct TimeInfo { /// @brief accessor for cycle size_t cycle() const { return cycle_; } - /// @brief return the time info corresponding to the end of this cycle (time + dt, dt, cycle) - TimeInfo endTimeInfo() const { return TimeInfo(time() + dt(), dt(), cycle()); } - + private: std::pair time_; ///< time and its dual std::pair dt_; ///< timestep and its dual size_t cycle_; ///< cycle, step, iteration count diff --git a/src/smith/physics/solid_mechanics.hpp b/src/smith/physics/solid_mechanics.hpp index 988b1ccee..6ff77a66c 100644 --- a/src/smith/physics/solid_mechanics.hpp +++ b/src/smith/physics/solid_mechanics.hpp @@ -1647,9 +1647,16 @@ class SolidMechanics, std::integer_se } if (use_warm_start_) { - // Update external forcing + // We want to solve using the forcing from the actual end of step, but, there are cases where the initial + // configuration is already out of equilibrium. In that case, we really want to consider the change in the + // forcing. For most runs, the previous state will be solved for equilibrium, so evaluating with displacement_, + // at time_ would give ~zero for unconstrained nodes. auto r = (*residual_)(time_ + dt, shapeDisplacement(), displacement_, acceleration_, *parameters_[parameter_indices].state...); + auto r_previous = (*residual_)(time_, shapeDisplacement(), displacement_, acceleration_, + *parameters_[parameter_indices].state...); + + r -= r_previous; // use the most recently evaluated Jacobian auto [_, drdu] = (*residual_)(time_, shapeDisplacement(), differentiate_wrt(displacement_), acceleration_,