diff --git a/src/serac/numerics/odes.cpp b/src/serac/numerics/odes.cpp index 8a63f99aa..9768a6433 100644 --- a/src/serac/numerics/odes.cpp +++ b/src/serac/numerics/odes.cpp @@ -21,6 +21,8 @@ SecondOrderODE::SecondOrderODE(int n, State&& state, const EquationSolver& solve void SecondOrderODE::SetTimestepper(const serac::TimestepMethod timestepper) { + timestepper_ = timestepper; + switch (timestepper) { case serac::TimestepMethod::Newmark: second_order_ode_solver_ = std::make_unique(); @@ -258,6 +260,8 @@ FirstOrderODE::FirstOrderODE(int n, FirstOrderODE::State&& state, const Equation void FirstOrderODE::SetTimestepper(const serac::TimestepMethod timestepper) { + timestepper_ = timestepper; + switch (timestepper) { case serac::TimestepMethod::BackwardEuler: ode_solver_ = std::make_unique(); diff --git a/src/serac/numerics/odes.hpp b/src/serac/numerics/odes.hpp index 33a61b480..a0ba72ce2 100644 --- a/src/serac/numerics/odes.hpp +++ b/src/serac/numerics/odes.hpp @@ -163,6 +163,13 @@ class SecondOrderODE : public mfem::SecondOrderTimeDependentOperator { */ const State& GetState() { return state_; } + /** + * @brief Query the timestep method for the ode solver + * + * @return The timestep method used by the underlying ode solver + */ + TimestepMethod GetTimestepper() { return timestepper_; } + private: /** * @brief Internal implementation used for mfem::SOTDO::Mult and mfem::SOTDO::ImplicitSolve @@ -214,6 +221,8 @@ class SecondOrderODE : public mfem::SecondOrderTimeDependentOperator { mutable mfem::Vector U_plus_; mutable mfem::Vector dU_dt_; mutable mfem::Vector d2U_dt2_; + + serac::TimestepMethod timestepper_; }; /** @@ -338,6 +347,13 @@ class FirstOrderODE : public mfem::TimeDependentOperator { } } + /** + * @brief Query the timestep method for the ode solver + * + * @return The timestep method used by the underlying ode solver + */ + TimestepMethod GetTimestepper() { return timestepper_; } + private: /** * @brief Internal implementation used for mfem::TDO::Mult and mfem::TDO::ImplicitSolve\ @@ -378,6 +394,8 @@ class FirstOrderODE : public mfem::TimeDependentOperator { mutable mfem::Vector U_; mutable mfem::Vector U_plus_; mutable mfem::Vector dU_dt_; + + TimestepMethod timestepper_; }; } // namespace serac::mfem_ext diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index ead9b44e5..71b077f8a 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -115,7 +115,7 @@ class BasePhysics { * @param parameter_name The name of the parameter to generate * * @note The user is responsible for managing the lifetime of this object. It is required - * to exist whenever advanceTimestep, solveAdjoint, or computeSensitivity is called. + * to exist whenever advanceTimestep, reverseAdjointTimestep, or computeTimestepSensitivity is called. * * @note The finite element space for this object is generated from the parameter * discretization space (e.g. L2, H1) and the computational mesh given in the physics module constructor. @@ -190,9 +190,9 @@ class BasePhysics { * * @return The sensitivity with respect to the parameter * - * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method. + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method. */ - virtual FiniteElementDual& computeSensitivity(size_t /* parameter_index */) + virtual FiniteElementDual& computeTimestepSensitivity(size_t /* parameter_index */) { SLIC_ERROR_ROOT(axom::fmt::format("Parameter sensitivities not enabled in physics module {}", name_)); return *parameters_[0].sensitivity; @@ -204,14 +204,28 @@ class BasePhysics { * * @return The sensitivity with respect to the shape displacement * - * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method. + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method. */ - virtual FiniteElementDual& computeShapeSensitivity() + virtual FiniteElementDual& computeTimestepShapeSensitivity() { SLIC_ERROR_ROOT(axom::fmt::format("Shape sensitivities not enabled in physics module {}", name_)); return shape_displacement_sensitivity_; } + /** + * @brief Compute the implicit sensitivity of the quantity of interest with respect to the initial condition fields + * + * @return Fields states corresponding to the sensitivities with respect to the initial condition fields + * + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method as many times as + * the forward advance is called. + */ + virtual const std::unordered_map computeInitialConditionSensitivity() + { + SLIC_ERROR_ROOT(axom::fmt::format("Initial condition sensitivities not enabled in physics module {}", name_)); + return {}; + } + /** * @brief Advance the state variables according to the chosen time integrator * @@ -227,7 +241,7 @@ class BasePhysics { * * @return The computed adjoint finite element states */ - virtual const std::unordered_map solveAdjoint( + virtual const std::unordered_map reverseAdjointTimestep( std::unordered_map /* adjoint_loads */, std::unordered_map /* adjoint_with_essential_boundary */ = {}) { diff --git a/src/serac/physics/heat_transfer.hpp b/src/serac/physics/heat_transfer.hpp index 107bd5b3d..f41ab05e4 100644 --- a/src/serac/physics/heat_transfer.hpp +++ b/src/serac/physics/heat_transfer.hpp @@ -102,7 +102,7 @@ class HeatTransfer, std::integer_sequ * @param[in] pmesh The mesh to conduct the simulation on, if different than the default mesh */ HeatTransfer(const NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions lin_opts, - const serac::TimesteppingOptions timestepping_opts, const std::string& name = "", + const TimesteppingOptions timestepping_opts, const std::string& name = "", mfem::ParMesh* pmesh = nullptr) : HeatTransfer(std::make_unique(nonlinear_opts, lin_opts, StateManager::mesh(StateManager::collectionID(pmesh)).GetComm()), @@ -126,14 +126,16 @@ class HeatTransfer, std::integer_sequ FiniteElementState::Options{ .order = order, .vector_dim = 1, .name = detail::addPrefix(name, "temperature")}, sidre_datacoll_id_)), + temperature_rate_(temperature_), adjoint_temperature_(StateManager::newState( FiniteElementState::Options{ .order = order, .vector_dim = 1, .name = detail::addPrefix(name, "adjoint_temperature")}, sidre_datacoll_id_)), + implicit_sensitivity_temperature_start_of_step_(adjoint_temperature_.space(), "total_deriv_wrt_temperature"), residual_with_bcs_(temperature_.space().TrueVSize()), nonlin_solver_(std::move(solver)), ode_(temperature_.space().TrueVSize(), - {.time = ode_time_point_, .u = u_, .dt = dt_, .du_dt = previous_, .previous_dt = previous_dt_}, + {.time = ode_time_point_, .u = u_, .dt = dt_, .du_dt = u_rate_start_of_step_, .previous_dt = previous_dt_}, *nonlin_solver_, bcs_) { SLIC_ERROR_ROOT_IF( @@ -189,15 +191,17 @@ class HeatTransfer, std::integer_sequ u_.SetSize(true_size); u_predicted_.SetSize(true_size); - previous_.SetSize(true_size); - previous_ = 0.0; + u_rate_start_of_step_.SetSize(true_size); + u_rate_start_of_step_ = 0.0; zero_.SetSize(true_size); zero_ = 0.0; - shape_displacement_ = 0.0; - temperature_ = 0.0; - adjoint_temperature_ = 0.0; + shape_displacement_ = 0.0; + temperature_ = 0.0; + temperature_rate_ = 0.0; + adjoint_temperature_ = 0.0; + implicit_sensitivity_temperature_start_of_step_ = 0.0; } /** @@ -556,7 +560,6 @@ class HeatTransfer, std::integer_sequ J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); return *J_; }); - } else { residual_with_bcs_ = mfem_ext::StdFunctionOperator( temperature_.space().TrueVSize(), @@ -601,29 +604,61 @@ class HeatTransfer, std::integer_sequ * @pre The adjoint load maps are expected to contain a single entry named "temperature" * @note If the essential boundary dual is not specified, homogeneous essential boundary conditions are applied to * the adjoint system + * @note We provide a quick derivation for the discrete adjoint equations and notation used here for backward + * Euler. There are two equations satisfied at each step of the forward solve using backward Euler \n 1). \f$r^n(u^n, + * v^n; d) = 0,\f$ \n where \f$u^n\f$ is the end-step primal value, \f$d\f$ are design parameters, and \f$v^n\f$ is + * the central difference velocity, satisfying \n 2). \f$\Delta t \, v^n = u^n-u^{n-1}.\f$ \n We are interesting in + * the implicit sensitivity of a qoi (quantity of interest), \f$\sum_{n=1}^N \pi^n(u^n, v^n; d)\f$, while maintaining + * the above constraints. We construct a Lagrangian that adds 'zero' to this qoi using the free multipliers + * \f$\lambda^n\f$ (which we call the adjoint temperature) and + * \f$\mu^n\f$ (which can eventually be intepreted as the implicit sensitivity of the qoi with respect to the + * start-of-step primal value \f$u^{n-1}\f$) with \f[ \mathcal{L} := \sum_{n=1}^N \pi(u^n, v^n; d) + \lambda^n \cdot + * r^n(u^n, v^n; d) + \mu^n \cdot (\Delta t \, v^n - u^n + u^{n-1}).\f] We are interesting in the total derivative + * \f[\frac{d\mathcal{L}}{dp} = \sum_{n=1}^N \pi^n_{,u} u^n_{,d} + \pi^n_{,v} v^n_{,d} + \pi^n_{,d} + \lambda^n \cdot + * \left( r^n_{,u} u^n_{,d} + r^n_{,v} v^n_{,d} + r^n_{,d} \right) + \mu^n \cdot (\Delta t \, v^n_{,d} - u^n_{,d} + + * u^{n-1}_{,d}). \f] We are free to choose \f$\lambda^n\f$ and \f$\mu^n\f$ in any convenient way, and the way we + * choose is by grouping terms involving \f$u^n_{,d}\f$ and \f$v^n_{,d}\f$, and setting the multipliers such that + * those terms cancel out in the final expression of the qoi sensitivity. In particular, by choosing \f[ \lambda^n = - + * \left[ r_{,u}^n + \frac{r_{,v}^n}{\Delta t} \right]^{-T} \left( \pi^n_{,u} + \frac{\pi^n_{,v}}{\Delta t} + + * \mu^{n+1} \right) \f] and \f[ \mu^n = -\frac{1}{\Delta t}\left( \pi^n_{,v} + \lambda^n \cdot r^n_{,v} \right), \f] + * we find + * \f[ \frac{d\mathcal{L}}{dp} = \sum_{n=1}^N \left( \pi^n_{,d} + \lambda^n \cdot r^n_{,d} \right) + \mu^1 \cdot + * u^{0}_{,d}, \f] where the multiplier/adjoint equations are solved backward in time starting at \f$n=N\f$, + * \f$\mu^{N+1} = 0\f$, and \f$u^{0}_{,d}\f$ is the sensitivity of the initial primal variable with respect to the + * parameters.\n We call the quantities \f$\pi^n_{,u}\f$ and \f$\pi^n_{,v}\f$ the temperature and temperature-rate + * adjoint loads, respectively. * * @param adjoint_loads An unordered map containing finite element duals representing the RHS of the adjoint equations - * indexed by their name + * indexed by their name. It must have a load named "temperature", which is the sensitivity of the qoi over that + * timestep to the end-of-step temperature. It may optionally have a load named "temperature_rate", which is the + * sensitivity of the qoi over that timestep with respect to the end-of_step temperature rate of change. * @param adjoint_with_essential_boundary An unordered map containing finite element states representing the * non-homogeneous essential boundary condition data for the adjoint problem indexed by their name * @return An unordered map of the adjoint solutions indexed by their name. It has a single entry named - * "adjoint_temperature" + * "adjoint_temperature". */ - const std::unordered_map solveAdjoint( + const std::unordered_map reverseAdjointTimestep( std::unordered_map adjoint_loads, std::unordered_map adjoint_with_essential_boundary = {}) override { - SLIC_ERROR_ROOT_IF(adjoint_loads.size() != 1, - "Adjoint load container is not the expected size of 1 in the heat transfer module."); + SLIC_ERROR_ROOT_IF(adjoint_loads.size() == 0, + "Adjoint load container size must be greater than 0 in the heat transfer module."); - auto temp_adjoint_load = adjoint_loads.find("temperature"); + auto temp_adjoint_load = adjoint_loads.find("temperature"); + auto temp_rate_adjoint_load = adjoint_loads.find("temperature_rate"); // does not need to be specified SLIC_ERROR_ROOT_IF(temp_adjoint_load == adjoint_loads.end(), "Adjoint load for \"temperature\" not found."); - mfem::HypreParVector adjoint_load_vector(temp_adjoint_load->second); - + mfem::HypreParVector temperature_adjoint_load_vector(temp_adjoint_load->second); // Add the sign correction to move the term to the RHS - adjoint_load_vector *= -1.0; + temperature_adjoint_load_vector *= -1.0; + + mfem::HypreParVector temperature_rate_adjoint_load_vector(temperature_adjoint_load_vector); + temperature_rate_adjoint_load_vector = 0.0; + if (temp_rate_adjoint_load != adjoint_loads.end()) { + temperature_rate_adjoint_load_vector = temp_rate_adjoint_load->second; + temperature_rate_adjoint_load_vector *= -1.0; + } auto& lin_solver = nonlin_solver_->linearSolver(); @@ -631,11 +666,6 @@ class HeatTransfer, std::integer_sequ mfem::HypreParVector adjoint_essential(temp_adjoint_load->second); adjoint_essential = 0.0; - auto [r, drdu] = (*residual_)(differentiate_wrt(temperature_), zero_, shape_displacement_, - *parameters_[parameter_indices].state...); - auto jacobian = assemble(drdu); - auto J_T = std::unique_ptr(jacobian->Transpose()); - // If we have a non-homogeneous essential boundary condition, extract it from the given state auto essential_adjoint_temp = adjoint_with_essential_boundary.find("temperature"); @@ -649,34 +679,113 @@ class HeatTransfer, std::integer_sequ "boundary condition named \"temperature\"."); } + if (is_quasistatic_) { + // We store the previous timestep's temperature as the current temperature for use in the lambdas computing the + // sensitivities. + + auto [r, drdu] = (*residual_)(differentiate_wrt(temperature_), zero_, shape_displacement_, + *parameters_[parameter_indices].state...); + auto jacobian = assemble(drdu); + auto J_T = std::unique_ptr(jacobian->Transpose()); + + for (const auto& bc : bcs_.essentials()) { + bc.apply(*J_T, temperature_adjoint_load_vector, adjoint_essential); + } + + lin_solver.SetOperator(*J_T); + lin_solver.Mult(temperature_adjoint_load_vector, adjoint_temperature_); + + // Reset the equation solver to use the full nonlinear residual operator + nonlin_solver_->setOperator(residual_with_bcs_); + + return {{"adjoint_temperature", adjoint_temperature_}}; + } + + SLIC_ERROR_ROOT_IF(ode_.GetTimestepper() != TimestepMethod::BackwardEuler, + "Only backward Euler implemented for transient adjoint heat conduction."); + + SLIC_ERROR_ROOT_IF(cycle_ == 0, + "Maximum number of adjoint timesteps exceeded! The number of adjoint timesteps must equal the " + "number of forward timesteps"); + + // Load the temperature from the previous cycle from disk + serac::FiniteElementState temperature_n_minus_1(temperature_); + StateManager::loadCheckpointedStates(cycle_, {temperature_}); + StateManager::loadCheckpointedStates(cycle_ - 1, {temperature_n_minus_1}); + + temperature_rate_ = temperature_; + temperature_rate_.Add(-1.0, temperature_n_minus_1); + temperature_rate_ /= dt_; + + // K := dR/du + auto K = serac::get((*residual_)(differentiate_wrt(temperature_), temperature_rate_, + shape_displacement_, *parameters_[parameter_indices].state...)); + std::unique_ptr k_mat(assemble(K)); + + // M := dR/du_dot + auto M = serac::get((*residual_)(temperature_, differentiate_wrt(temperature_rate_), + shape_displacement_, *parameters_[parameter_indices].state...)); + std::unique_ptr m_mat(assemble(M)); + + J_.reset(mfem::Add(1.0, *m_mat, dt_, *k_mat)); + auto J_T = std::unique_ptr(J_->Transpose()); + + // recall that temperature_adjoint_load_vector and d_temperature_dt_adjoint_load_vector were already multiplied by + // -1 above + mfem::HypreParVector modified_RHS(temperature_adjoint_load_vector); + modified_RHS *= dt_; + modified_RHS.Add(1.0, temperature_rate_adjoint_load_vector); + modified_RHS.Add(-dt_, implicit_sensitivity_temperature_start_of_step_); + for (const auto& bc : bcs_.essentials()) { - bc.apply(*J_T, adjoint_load_vector, adjoint_essential); + bc.apply(*J_T, modified_RHS, adjoint_essential); } lin_solver.SetOperator(*J_T); - lin_solver.Mult(adjoint_load_vector, adjoint_temperature_); + lin_solver.Mult(modified_RHS, adjoint_temperature_); + + // This multiply is technically on M transposed. However, this matrix should be symmetric unless + // the thermal capacity is a function of the temperature rate of change, which is thermodynamically + // impossible, and fortunately not possible with our material interface. + // Not doing the transpose here to avoid doing unnecessary work. + m_mat->Mult(adjoint_temperature_, implicit_sensitivity_temperature_start_of_step_); + implicit_sensitivity_temperature_start_of_step_ *= -1.0 / dt_; + implicit_sensitivity_temperature_start_of_step_.Add( + 1.0 / dt_, temperature_rate_adjoint_load_vector); // already multiplied by -1 // Reset the equation solver to use the full nonlinear residual operator nonlin_solver_->setOperator(residual_with_bcs_); + time_ -= dt_; + cycle_--; + return {{"adjoint_temperature", adjoint_temperature_}}; } + /** + * @brief Get a temperature which has been previously checkpointed at the give cycle + * @param cycle The previous timestep where the temperature solution is requested + * @return The temperature FiniteElementState at a previous cycle + */ + FiniteElementState loadCheckpointedTemperature(int cycle) const + { + FiniteElementState previous_temperature(temperature_); + StateManager::loadCheckpointedStates(cycle, {previous_temperature}); + return previous_temperature; + } + /** * @brief Compute the implicit sensitivity of the quantity of interest used in defining the load for the adjoint - * problem with respect to the parameter field + * problem with respect to the parameter field for the last computed adjoint timestep * * @tparam parameter_field The index of the parameter to take a derivative with respect to * @return The sensitivity with respect to the parameter * - * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method. + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method. */ - FiniteElementDual& computeSensitivity(size_t parameter_field) override + FiniteElementDual& computeTimestepSensitivity(size_t parameter_field) override { - SLIC_ASSERT_MSG(parameter_field < sizeof...(parameter_indices), - axom::fmt::format("Invalid parameter index '{}' reqested for sensitivity.")); - - auto drdparam = serac::get<1>(d_residual_d_[parameter_field]()); + auto drdparam = serac::get(d_residual_d_[parameter_field]()); auto drdparam_mat = assemble(drdparam); drdparam_mat->MultTranspose(adjoint_temperature_, *parameters_[parameter_field].sensitivity); @@ -686,17 +795,16 @@ class HeatTransfer, std::integer_sequ /** * @brief Compute the implicit sensitivity of the quantity of interest used in defining the load for the adjoint - * problem with respect to the shape displacement field + * problem with respect to the shape displacement field for the last computed adjoint timestep * * @return The sensitivity with respect to the shape displacement * - * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method. + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method. */ - FiniteElementDual& computeShapeSensitivity() override + FiniteElementDual& computeTimestepShapeSensitivity() override { - auto drdshape = serac::get((*residual_)(DifferentiateWRT{}, temperature_, zero_, + auto drdshape = serac::get((*residual_)(DifferentiateWRT{}, temperature_, temperature_rate_, shape_displacement_, *parameters_[parameter_indices].state...)); - auto drdshape_mat = assemble(drdshape); drdshape_mat->MultTranspose(adjoint_temperature_, shape_displacement_sensitivity_); @@ -704,6 +812,18 @@ class HeatTransfer, std::integer_sequ return shape_displacement_sensitivity_; } + /** + * @brief Compute the implicit sensitivity of the quantity of interest with respect to the initial temperature + * + * @return The sensitivity with respect to the initial temperature + * + * @pre `reverseAdjointTimestep` must be called as many times as the forward solver was advanced before this is called + */ + const std::unordered_map computeInitialConditionSensitivity() override + { + return {{"temperature", implicit_sensitivity_temperature_start_of_step_}}; + } + /// Destroy the Thermal Solver object virtual ~HeatTransfer() = default; @@ -720,9 +840,15 @@ class HeatTransfer, std::integer_sequ /// The temperature finite element state serac::FiniteElementState temperature_; - /// The adjoint temperature finite element state + /// Rate of change in temperature at the current adjoint timestep + FiniteElementState temperature_rate_; + + /// The adjoint temperature finite element states, the multiplier on the residual for a given timestep serac::FiniteElementState adjoint_temperature_; + /// The total/implicit sensitivity of the qoi with respect to the start of the previous timestep's temperature + serac::FiniteElementDual implicit_sensitivity_temperature_start_of_step_; + /// serac::Functional that is used to calculate the residual and its derivatives std::unique_ptr> residual_; @@ -771,16 +897,16 @@ class HeatTransfer, std::integer_sequ mfem::Vector u_predicted_; /// Previous value of du_dt used to prime the pump for the nonlinear solver - mfem::Vector previous_; + mfem::Vector u_rate_start_of_step_; /// @brief Array functions computing the derivative of the residual with respect to each given parameter /// @note This is needed so the user can ask for a specific sensitivity at runtime as opposed to it being a /// template parameter. - std::array{}, temperature_, zero_, shape_displacement_, - *parameters_[parameter_indices].state...))()>, + std::array{}, temperature_, temperature_rate_, + shape_displacement_, *parameters_[parameter_indices].state...))()>, sizeof...(parameter_indices)> d_residual_d_ = {[&]() { - return (*residual_)(DifferentiateWRT{}, temperature_, zero_, + return (*residual_)(DifferentiateWRT{}, temperature_, temperature_rate_, shape_displacement_, *parameters_[parameter_indices].state...); }...}; }; diff --git a/src/serac/physics/materials/thermal_material.hpp b/src/serac/physics/materials/thermal_material.hpp index 2b0de9e93..33920aad0 100644 --- a/src/serac/physics/materials/thermal_material.hpp +++ b/src/serac/physics/materials/thermal_material.hpp @@ -65,6 +65,65 @@ struct LinearIsotropicConductor { double conductivity_; }; +/// Nonlinear isotropic thermal conduction material model +struct IsotropicConductorWithLinearConductivityVsTemperature { + /** + * @brief Construct a Isotropic Conductor with Conductivity linear with Temparture object + * + * @param density Density of the material (mass/volume) + * @param specific_heat_capacity Specific heat capacity of the material (energy / (mass * temp)) + * @param reference_conductivity Reference thermal conductivity of the material at temp = 0 (power / (length * temp)) + * @param d_conductivity_d_temperature Slope for the thermal conductivity as a function of temperature + */ + IsotropicConductorWithLinearConductivityVsTemperature(double density = 1.0, double specific_heat_capacity = 1.0, + double reference_conductivity = 1.0, + double d_conductivity_d_temperature = 0.0) + : density_(density), + specific_heat_capacity_(specific_heat_capacity), + reference_conductivity_(reference_conductivity), + d_conductivity_d_temperature_(d_conductivity_d_temperature) + { + SLIC_ERROR_ROOT_IF(density_ < 0.0, "Density must be positive in the linear isotropic conductor material model."); + + SLIC_ERROR_ROOT_IF(specific_heat_capacity_ < 0.0, + "Specific heat capacity must be positive in the linear isotropic conductor material model."); + } + + /** + * @brief Material response call for a linear isotropic material with linear conductivity vs temperature + * + * @tparam T1 Spatial position type + * @tparam T2 Temperature type + * @tparam T3 Temperature gradient type + * @param[in] temperature Temperature + * @param[in] temperature_gradient Temperature gradient + * @return The calculated material response (tuple of volumetric heat capacity and thermal flux) for a linear + * isotropic material + */ + template + SERAC_HOST_DEVICE auto operator()(const T1& /* x */, const T2& temperature, const T3& temperature_gradient) const + { + const auto currentConductivity = reference_conductivity_ + d_conductivity_d_temperature_ * temperature; + SLIC_ERROR_ROOT_IF( + serac::get_value(currentConductivity) < 0.0, + "Conductivity in the IsotropicConductorWithLinearConductivityVsTemperature model has gone negative."); + return serac::tuple{density_ * specific_heat_capacity_, -1.0 * currentConductivity * temperature_gradient}; + } + +private: + /// Density + double density_; + + /// Specific heat capacity + double specific_heat_capacity_; + + /// Reference isotropic thermal conductivity + double reference_conductivity_; + + /// Slope of nonlinear thermal conductivity dependence on temperature + double d_conductivity_d_temperature_; +}; + /** * @brief Linear anisotropic thermal material model * diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 3c74adf4a..25eba2410 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1077,7 +1077,7 @@ class SolidMechanics, std::integer_se * @return An unordered map of the adjoint solutions indexed by their name. It has a single entry named * "adjoint_displacement" */ - const std::unordered_map solveAdjoint( + const std::unordered_map reverseAdjointTimestep( std::unordered_map adjoint_loads, std::unordered_map adjoint_with_essential_boundary = {}) override { @@ -1135,9 +1135,9 @@ class SolidMechanics, std::integer_se * @param parameter_field The index of the parameter to take a derivative with respect to * @return The sensitivity with respect to the parameter * - * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method. + * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method. */ - FiniteElementDual& computeSensitivity(size_t parameter_field) override + FiniteElementDual& computeTimestepSensitivity(size_t parameter_field) override { SLIC_ASSERT_MSG(parameter_field < sizeof...(parameter_indices), axom::fmt::format("Invalid parameter index '{}' requested for sensitivity.")); @@ -1157,7 +1157,7 @@ class SolidMechanics, std::integer_se * * @return The sensitivity with respect to the shape displacement */ - FiniteElementDual& computeShapeSensitivity() override + FiniteElementDual& computeTimestepShapeSensitivity() override { auto drdshape = serac::get((*residual_)(DifferentiateWRT<2>{}, displacement_, zero_, shape_displacement_, *parameters_[parameter_indices].state...)); diff --git a/src/serac/physics/state/finite_element_vector.cpp b/src/serac/physics/state/finite_element_vector.cpp index f3b423afc..fb753983f 100644 --- a/src/serac/physics/state/finite_element_vector.cpp +++ b/src/serac/physics/state/finite_element_vector.cpp @@ -150,4 +150,37 @@ double min(const FiniteElementVector& fe_vector) return global_min; } +/** + * @brief Check if two finite element spaces are the same + * + * @param left + * @param right + * @return Bool which is true if the spaces are the same, otherwise false + */ +bool sameFiniteElementSpace(const mfem::FiniteElementSpace& left, const mfem::FiniteElementSpace& right) +{ + bool sameMesh = (left.GetMesh() == right.GetMesh()); + bool equivalentFEColl = strcmp(left.FEColl()->Name(), right.FEColl()->Name()) == 0; + bool sameVectorDimension = (left.GetVDim() == right.GetVDim()); + bool sameOrdering = (left.GetOrdering() == right.GetOrdering()); + return sameMesh && equivalentFEColl && sameVectorDimension && sameOrdering; +} + +double innerProduct(const FiniteElementVector& v1, const FiniteElementVector& v2) +{ + SLIC_ERROR_IF( + v1.Size() != v2.Size(), + axom::fmt::format("Finite element vector of size '{}' can not inner product with another vector of size '{}'", + v1.Size(), v2.Size())); + SLIC_ERROR_IF(v1.comm() != v2.comm(), + "Cannot compute inner products between finite element vectors with different mpi communicators"); + SLIC_ERROR_IF(!sameFiniteElementSpace(v1.space(), v2.space()), + "Currently cannot compute inner products between finite element vectors with different mfem spaces"); + + double global_ip; + double local_ip = mfem::InnerProduct(v1, v2); + MPI_Allreduce(&local_ip, &global_ip, 1, MPI_DOUBLE, MPI_SUM, v1.comm()); + return global_ip; +} + } // namespace serac diff --git a/src/serac/physics/state/finite_element_vector.hpp b/src/serac/physics/state/finite_element_vector.hpp index be625fd7f..a3d6acc5b 100644 --- a/src/serac/physics/state/finite_element_vector.hpp +++ b/src/serac/physics/state/finite_element_vector.hpp @@ -229,4 +229,13 @@ double max(const FiniteElementVector& fe_vector); */ double min(const FiniteElementVector& fe_vector); +/** + * @brief Find the inner prodcut between two finite element vectors across all dofs + * + * @param vec1 The first vector + * @param vec2 The second vector + * @return The inner prodcut between finite element vectors + */ +double innerProduct(const FiniteElementVector& vec1, const FiniteElementVector& vec2); + } // namespace serac diff --git a/src/serac/physics/state/state_manager.cpp b/src/serac/physics/state/state_manager.cpp index b1cffa733..9a84b76d6 100644 --- a/src/serac/physics/state/state_manager.cpp +++ b/src/serac/physics/state/state_manager.cpp @@ -92,6 +92,28 @@ double StateManager::newDataCollection(const std::string& name, const std::optio return datacoll.GetTime(); } +void StateManager::loadCheckpointedStates(int cycle_to_load, + std::vector> states_to_load) +{ + std::string mesh_name = collectionID(&states_to_load.begin()->get().mesh()); + + std::string coll_name = mesh_name + "_datacoll"; + + axom::sidre::MFEMSidreDataCollection previous_datacoll(coll_name); + + previous_datacoll.SetComm(states_to_load.begin()->get().mesh().GetComm()); + previous_datacoll.SetPrefixPath(output_dir_); + previous_datacoll.Load(cycle_to_load); + + for (auto state : states_to_load) { + SLIC_ERROR_ROOT_IF(collectionID(&state.get().mesh()) != mesh_name, + "Loading FiniteElementStates from two different meshes at one time is not allowed."); + mfem::ParGridFunction* datacoll_owned_grid_function = previous_datacoll.GetParField(state.get().name()); + + state.get().setFromGridFunction(*datacoll_owned_grid_function); + } +} + void StateManager::initialize(axom::sidre::DataStore& ds, const std::string& output_directory) { // If the global object has already been initialized, clear it out diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index 3955770b7..c36d77f1f 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -193,6 +193,15 @@ class StateManager { */ static FiniteElementState& shapeDisplacement(const std::string& mesh_tag = default_mesh_name_); + /** + * @brief loads the finite element states from a previously checkpointed cycle + * + * @param cycle_to_load + * @param states_to_load + */ + static void loadCheckpointedStates(int cycle_to_load, + std::vector> states_to_load); + /** * @brief Get the shape displacement sensitivity finite element dual * diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index 092a052d3..984eed3c1 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -30,6 +30,7 @@ set(solver_tests solid_shape.cpp thermal_shape.cpp thermal_mechanics.cpp + dynamic_thermal_adjoint.cpp ) serac_add_tests( SOURCES ${solver_tests} diff --git a/src/serac/physics/tests/dynamic_thermal_adjoint.cpp b/src/serac/physics/tests/dynamic_thermal_adjoint.cpp new file mode 100644 index 000000000..178a09856 --- /dev/null +++ b/src/serac/physics/tests/dynamic_thermal_adjoint.cpp @@ -0,0 +1,315 @@ +// Copyright (c) 2019-2023, Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include "serac/physics/heat_transfer.hpp" + +#include +#include +#include + +#include "axom/slic/core/SimpleLogger.hpp" +#include +#include "mfem.hpp" + +#include "serac/mesh/mesh_utils.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/materials/thermal_material.hpp" +#include "serac/physics/materials/parameterized_thermal_material.hpp" +#include "serac/serac_config.hpp" + +namespace serac { + +constexpr int dim = 2; +constexpr int p = 1; + +const std::string thermal_prefix = "thermal"; +const std::string parametrized_thermal_prefix = "thermal_with_param"; + +struct TimeSteppingInfo { + double totalTime = 0.6; + int num_timesteps = 4; +}; + +double computeStepQoi(const FiniteElementState& temperature, double dt) +{ + return 0.5 * dt * innerProduct(temperature, temperature); +} + +void computeStepAdjointLoad(const FiniteElementState& temperature, FiniteElementDual& d_qoi_d_temperature, double dt) +{ + d_qoi_d_temperature = temperature; + d_qoi_d_temperature *= dt; +} + +std::unique_ptr> create_nonlinear_heat_transfer( + axom::sidre::DataStore& /*data_store*/, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, + const heat_transfer::IsotropicConductorWithLinearConductivityVsTemperature& mat) +{ + // eventually figure out how to clear out sidre state + // auto saveMesh = std::make_unique(StateManager::mesh()); + // StateManager::reset(); + // static int iter = 0; + // StateManager::initialize(data_store, "thermal_dynamic_solve"+std::to_string(iter++)); + // std::string filename = std::string(SERAC_REPO_DIR) + "/data/meshes/star.mesh"; + // mfem::ParMesh* mesh = StateManager::setMesh(std::move(saveMesh)); + static int iter = 0; + auto thermal = std::make_unique>(nonlinear_opts, heat_transfer::direct_linear_options, dyn_opts, + thermal_prefix + std::to_string(iter++)); + thermal->setMaterial(mat); + thermal->setTemperature([](const mfem::Vector&, double) { return 0.0; }); + thermal->setTemperatureBCs({1}, [](const mfem::Vector&, double) { return 0.0; }); + thermal->setSource([](auto /* X */, auto /* time */, auto /* u */, auto /* du_dx */) { return 1.0; }); + thermal->completeSetup(); + return thermal; +} + +using ParametrizedHeatTransferT = HeatTransfer>, std::integer_sequence>; + +std::unique_ptr create_parameterized_heat_transfer( + axom::sidre::DataStore& /*data_store*/, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, const heat_transfer::ParameterizedLinearIsotropicConductor& mat) +{ + // eventually figure out how to clear out sidre state + // auto saveMesh = std::make_unique(StateManager::mesh()); + // StateManager::reset(); + // static int iter = 0; + // StateManager::initialize(data_store, "thermal_dynamic_solve"+std::to_string(iter++)); + // std::string filename = std::string(SERAC_REPO_DIR) + "/data/meshes/star.mesh"; + // mfem::ParMesh* mesh = StateManager::setMesh(std::move(saveMesh)); + static int iter = 0; + auto thermal = + std::make_unique(nonlinear_opts, heat_transfer::direct_linear_options, dyn_opts, + parametrized_thermal_prefix + std::to_string(iter++)); + auto user_defined_conductivity_ptr = thermal->generateParameter("pcond", 0); + *user_defined_conductivity_ptr = 1.1; + thermal->setMaterial(DependsOn<0>{}, mat); + thermal->setTemperature([](const mfem::Vector&, double) { return 0.0; }); + thermal->setTemperatureBCs({1}, [](const mfem::Vector&, double) { return 0.0; }); + thermal->setSource([](auto /* X */, auto /* time */, auto /* u */, auto /* du_dx */) { return 1.0; }); + thermal->completeSetup(); + return thermal; // std::make_pair(thermal, user_defined_conductivity); +} + +double computeThermalQoiAdjustingInitalTemperature( + axom::sidre::DataStore& data_store, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, + const heat_transfer::IsotropicConductorWithLinearConductivityVsTemperature& mat, const TimeSteppingInfo& ts_info, + const FiniteElementState& init_temp_derivative_direction, double pertubation) +{ + auto thermal = create_nonlinear_heat_transfer(data_store, nonlinear_opts, dyn_opts, mat); + + auto& temperature = thermal->temperature(); + SLIC_ASSERT_MSG(temperature.Size() == init_temp_derivative_direction.Size(), + "Shape displacement and intended derivative direction FiniteElementState sizes do not agree."); + + temperature.Add(pertubation, init_temp_derivative_direction); + + double qoi = 0.0; + thermal->outputState(); + for (int i = 0; i < ts_info.num_timesteps; ++i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + thermal->advanceTimestep(dt); + thermal->outputState(); + qoi += computeStepQoi(thermal->temperature(), dt); + } + return qoi; +} + +double computeThermalQoiAdjustingShape(axom::sidre::DataStore& data_store, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, + const heat_transfer::IsotropicConductorWithLinearConductivityVsTemperature& mat, + const TimeSteppingInfo& ts_info, + const FiniteElementState& shape_derivative_direction, double pertubation) +{ + auto thermal = create_nonlinear_heat_transfer(data_store, nonlinear_opts, dyn_opts, mat); + + auto& shapeDisp = thermal->shapeDisplacement(); + SLIC_ASSERT_MSG(shapeDisp.Size() == shape_derivative_direction.Size(), + "Shape displacement and intended derivative direction FiniteElementState sizes do not agree."); + + shapeDisp.Add(pertubation, shape_derivative_direction); + + double qoi = 0.0; + thermal->outputState(); + for (int i = 0; i < ts_info.num_timesteps; ++i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + thermal->advanceTimestep(dt); + thermal->outputState(); + qoi += computeStepQoi(thermal->temperature(), dt); + } + return qoi; +} + +std::tuple computeThermalQoiAndInitialTemperatureAndShapeSensitivity( + axom::sidre::DataStore& data_store, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, + const heat_transfer::IsotropicConductorWithLinearConductivityVsTemperature& mat, const TimeSteppingInfo& ts_info) +{ + auto thermal = create_nonlinear_heat_transfer(data_store, nonlinear_opts, dyn_opts, mat); + + double qoi = 0.0; + thermal->outputState(); + for (int i = 0; i < ts_info.num_timesteps; ++i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + thermal->advanceTimestep(dt); + thermal->outputState(); + qoi += computeStepQoi(thermal->temperature(), dt); + } + + FiniteElementDual initial_temperature_sensitivity(thermal->temperature().space(), "init_temp_sensitivity"); + initial_temperature_sensitivity = 0.0; + FiniteElementDual shape_sensitivity(thermal->shapeDisplacement().space(), "shape_sensitivity"); + shape_sensitivity = 0.0; + + FiniteElementDual adjoint_load(thermal->temperature().space(), "adjoint_load"); + + for (int i = ts_info.num_timesteps; i > 0; --i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + FiniteElementState temperature_end_of_step = thermal->loadCheckpointedTemperature(thermal->cycle()); + computeStepAdjointLoad(temperature_end_of_step, adjoint_load, dt); + thermal->reverseAdjointTimestep({{"temperature", adjoint_load}}); + shape_sensitivity += thermal->computeTimestepShapeSensitivity(); + } + + EXPECT_EQ(0, thermal->cycle()); // we are back to the start + auto initialConditionSensitivities = thermal->computeInitialConditionSensitivity(); + auto initialTemperatureSensitivityIter = initialConditionSensitivities.find("temperature"); + SLIC_ASSERT_MSG(initialTemperatureSensitivityIter != initialConditionSensitivities.end(), + "Could not find temperature in the computed initial condition sensitivities."); + initial_temperature_sensitivity += initialTemperatureSensitivityIter->second; + + return std::make_tuple(qoi, initial_temperature_sensitivity, shape_sensitivity); +} + +std::tuple computeThermalConductivitySensitivity( + axom::sidre::DataStore& data_store, const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, const heat_transfer::ParameterizedLinearIsotropicConductor& mat, + const TimeSteppingInfo& ts_info) +{ + auto thermal = create_parameterized_heat_transfer(data_store, nonlinear_opts, dyn_opts, mat); + + double qoi = 0.0; + thermal->outputState(); + for (int i = 0; i < ts_info.num_timesteps; ++i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + thermal->advanceTimestep(dt); + thermal->outputState(); + qoi += computeStepQoi(thermal->temperature(), dt); + } + + FiniteElementDual initial_temperature_sensitivity(thermal->temperature().space(), "init_temp_sensitivity"); + initial_temperature_sensitivity = 0.0; + FiniteElementDual shape_sensitivity(thermal->shapeDisplacement().space(), "shape_sensitivity"); + shape_sensitivity = 0.0; + + FiniteElementDual adjoint_load(thermal->temperature().space(), "adjoint_load"); + + for (int i = ts_info.num_timesteps; i > 0; --i) { + double dt = ts_info.totalTime / ts_info.num_timesteps; + FiniteElementState temperature_end_of_step = thermal->loadCheckpointedTemperature(thermal->cycle()); + computeStepAdjointLoad(temperature_end_of_step, adjoint_load, dt); + thermal->reverseAdjointTimestep({{"temperature", adjoint_load}}); + shape_sensitivity += thermal->computeTimestepShapeSensitivity(); + } + + EXPECT_EQ(0, thermal->cycle()); // we are back to the start + auto initialConditionSensitivities = thermal->computeInitialConditionSensitivity(); + auto initialTemperatureSensitivityIter = initialConditionSensitivities.find("temperature"); + SLIC_ASSERT_MSG(initialTemperatureSensitivityIter != initialConditionSensitivities.end(), + "Could not find temperature in the computed initial condition sensitivities."); + initial_temperature_sensitivity += initialTemperatureSensitivityIter->second; + + return std::make_tuple(qoi, initial_temperature_sensitivity, shape_sensitivity); +} + +struct HeatTransferSensitivityFixture : public ::testing::Test { + void SetUp() override + { + MPI_Barrier(MPI_COMM_WORLD); + StateManager::initialize(dataStore, "thermal_dynamic_solve"); + std::string filename = std::string(SERAC_REPO_DIR) + "/data/meshes/star.mesh"; + mesh = StateManager::setMesh(mesh::refineAndDistribute(buildMeshFromFile(filename), 0)); + } + + void fillDirection(FiniteElementState& direction) const { direction = 1.1; } + + // Create DataStore + axom::sidre::DataStore dataStore; + mfem::ParMesh* mesh; + + // Solver options + NonlinearSolverOptions nonlinear_opts{.relative_tol = 5.0e-13, .absolute_tol = 5.0e-13}; + + TimesteppingOptions dyn_opts{.timestepper = TimestepMethod::BackwardEuler, + .enforcement_method = DirichletEnforcementMethod::DirectControl}; + + heat_transfer::IsotropicConductorWithLinearConductivityVsTemperature mat{1.0, 1.0, 1.0, 2.0}; + heat_transfer::ParameterizedLinearIsotropicConductor parameterizedMat; + + TimeSteppingInfo tsInfo{.totalTime = 0.5, .num_timesteps = 4}; +}; + +TEST_F(HeatTransferSensitivityFixture, InitialTemperatureSensitivities) +{ + auto [qoi_base, temperature_sensitivity, _] = + computeThermalQoiAndInitialTemperatureAndShapeSensitivity(dataStore, nonlinear_opts, dyn_opts, mat, tsInfo); + + FiniteElementState derivative_direction(temperature_sensitivity.space(), "derivative_direction"); + fillDirection(derivative_direction); + + const double eps = 1e-7; + double qoi_plus = computeThermalQoiAdjustingInitalTemperature(dataStore, nonlinear_opts, dyn_opts, mat, tsInfo, + derivative_direction, eps); + double directional_deriv = innerProduct(derivative_direction, temperature_sensitivity); + EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); +} + +TEST_F(HeatTransferSensitivityFixture, ShapeSensitivities) +{ + auto [qoi_base, _, shape_sensitivity] = + computeThermalQoiAndInitialTemperatureAndShapeSensitivity(dataStore, nonlinear_opts, dyn_opts, mat, tsInfo); + + FiniteElementState derivative_direction(shape_sensitivity.space(), "derivative_direction"); + fillDirection(derivative_direction); + + const double eps = 1e-7; + + double qoi_plus = + computeThermalQoiAdjustingShape(dataStore, nonlinear_opts, dyn_opts, mat, tsInfo, derivative_direction, eps); + double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); + EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); +} + +TEST_F(HeatTransferSensitivityFixture, DISABLED_ConductivityParameterSensitivities) +{ + auto [qoi_base, _, shape_sensitivity] = + computeThermalConductivitySensitivity(dataStore, nonlinear_opts, dyn_opts, parameterizedMat, tsInfo); + + FiniteElementState derivative_direction(shape_sensitivity.space(), "derivative_direction"); + fillDirection(derivative_direction); + + const double eps = 1e-7; + + double qoi_plus = + computeThermalQoiAdjustingShape(dataStore, nonlinear_opts, dyn_opts, mat, tsInfo, derivative_direction, eps); + double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); + EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); +} + +} // namespace serac + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + + axom::slic::SimpleLogger logger; + + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + + return result; +} diff --git a/src/serac/physics/tests/parameterized_thermal.cpp b/src/serac/physics/tests/parameterized_thermal.cpp index 4c6bd18f8..02b516aef 100644 --- a/src/serac/physics/tests/parameterized_thermal.cpp +++ b/src/serac/physics/tests/parameterized_thermal.cpp @@ -42,7 +42,7 @@ TEST(Thermal, ParameterizedMaterial) // Define a boundary attribute set std::set ess_bdr = {1}; - // Construct and initialized the user-defined conductivity to be used as a differentiable parameter in + // Construct and initialize the user-defined conductivity to be used as a differentiable parameter in // the thermal conduction physics module. FiniteElementState user_defined_conductivity( StateManager::newState(FiniteElementState::Options{.order = 1, .name = "parameterized_conductivity"})); @@ -106,10 +106,10 @@ TEST(Thermal, ParameterizedMaterial) adjoint_load = 1.0; // Solve the adjoint problem - thermal_solver.solveAdjoint({{"temperature", adjoint_load}}); + thermal_solver.reverseAdjointTimestep({{"temperature", adjoint_load}}); // Compute the sensitivity (d QOI/ d state * d state/d parameter) given the current adjoint solution - auto& sensitivity = thermal_solver.computeSensitivity(conductivity_parameter_index); + auto& sensitivity = thermal_solver.computeTimestepSensitivity(conductivity_parameter_index); EXPECT_NEAR(1.7890782925134845, mfem::ParNormlp(sensitivity, 2, MPI_COMM_WORLD), 1.0e-6); } diff --git a/src/serac/physics/tests/parameterized_thermomechanics_example.cpp b/src/serac/physics/tests/parameterized_thermomechanics_example.cpp index 9de87400b..72f22d2d4 100644 --- a/src/serac/physics/tests/parameterized_thermomechanics_example.cpp +++ b/src/serac/physics/tests/parameterized_thermomechanics_example.cpp @@ -184,9 +184,9 @@ TEST(Thermomechanics, ParameterizedMaterial) check_gradient(qoi, simulation.displacement()); - simulation.solveAdjoint({{"displacement", adjoint_load}}); + simulation.reverseAdjointTimestep({{"displacement", adjoint_load}}); - auto& dqoi_dalpha = simulation.computeSensitivity(1); + auto& dqoi_dalpha = simulation.computeTimestepSensitivity(1); double epsilon = 1.0e-5; auto dalpha = alpha.CreateCompatibleVector(); diff --git a/src/serac/physics/tests/solid.cpp b/src/serac/physics/tests/solid.cpp index 5b8394950..5b5599881 100644 --- a/src/serac/physics/tests/solid.cpp +++ b/src/serac/physics/tests/solid.cpp @@ -371,8 +371,8 @@ void functional_parameterized_solid_test(double expected_disp_norm) // are not used, but running them as part of this test // checks the index-translation part of the derivative // kernels is working - solid_solver.computeSensitivity(0); - solid_solver.computeSensitivity(1); + solid_solver.computeTimestepSensitivity(0); + solid_solver.computeTimestepSensitivity(1); // Output the sidre-based plot files solid_solver.outputState(); diff --git a/src/serac/physics/tests/solid_finite_diff.cpp b/src/serac/physics/tests/solid_finite_diff.cpp index 15b9e8456..cb6825225 100644 --- a/src/serac/physics/tests/solid_finite_diff.cpp +++ b/src/serac/physics/tests/solid_finite_diff.cpp @@ -116,10 +116,10 @@ TEST(SolidMechanics, FiniteDifferenceParameter) adjoint_load = *assembled_vector; // Solve the adjoint problem - solid_solver.solveAdjoint({{"displacement", adjoint_load}}); + solid_solver.reverseAdjointTimestep({{"displacement", adjoint_load}}); // Compute the sensitivity (d QOI/ d state * d state/d parameter) given the current adjoint solution - [[maybe_unused]] auto& sensitivity = solid_solver.computeSensitivity(bulk_parameter_index); + [[maybe_unused]] auto& sensitivity = solid_solver.computeTimestepSensitivity(bulk_parameter_index); // Perform finite difference on each bulk modulus value // to check if computed qoi sensitivity is consistent @@ -279,10 +279,10 @@ void finite_difference_shape_test(LoadingType load) adjoint_load = *assembled_vector; // Solve the adjoint problem - solid_solver.solveAdjoint({{"displacement", adjoint_load}}); + solid_solver.reverseAdjointTimestep({{"displacement", adjoint_load}}); // Compute the sensitivity (d QOI/ d state * d state/d parameter) given the current adjoint solution - [[maybe_unused]] auto& sensitivity = solid_solver.computeShapeSensitivity(); + [[maybe_unused]] auto& sensitivity = solid_solver.computeTimestepShapeSensitivity(); // Perform finite difference on each shape velocity value // to check if computed qoi sensitivity is consistent diff --git a/src/serac/physics/tests/thermal_finite_diff.cpp b/src/serac/physics/tests/thermal_finite_diff.cpp index bf00eb802..5fdfd9180 100644 --- a/src/serac/physics/tests/thermal_finite_diff.cpp +++ b/src/serac/physics/tests/thermal_finite_diff.cpp @@ -101,10 +101,10 @@ TEST(Thermal, FiniteDifference) adjoint_load = *assembled_vector; // Solve the adjoint problem - thermal_solver.solveAdjoint({{"temperature", adjoint_load}}); + thermal_solver.reverseAdjointTimestep({{"temperature", adjoint_load}}); // Compute the sensitivity (d QOI/ d state * d state/d parameter) given the current adjoint solution - [[maybe_unused]] auto& sensitivity = thermal_solver.computeSensitivity(conductivity_parameter_index); + [[maybe_unused]] auto& sensitivity = thermal_solver.computeTimestepSensitivity(conductivity_parameter_index); // Perform finite difference on each conduction value // to check if computed qoi sensitivity is consistent @@ -216,10 +216,10 @@ TEST(HeatTransfer, FiniteDifferenceShape) adjoint_load = *assembled_vector; // Solve the adjoint problem - thermal_solver.solveAdjoint({{"temperature", adjoint_load}}); + thermal_solver.reverseAdjointTimestep({{"temperature", adjoint_load}}); // Compute the sensitivity (d QOI/ d state * d state/d parameter) given the current adjoint solution - [[maybe_unused]] auto& sensitivity = thermal_solver.computeShapeSensitivity(); + [[maybe_unused]] auto& sensitivity = thermal_solver.computeTimestepShapeSensitivity(); // Perform finite difference on each shape displacement value // to check if computed qoi sensitivity is consistent