From 9cdca8aa39c501079eb85c01f78c541643d937c7 Mon Sep 17 00:00:00 2001 From: Agent Style Date: Tue, 10 Sep 2024 20:44:16 -0700 Subject: [PATCH] Apply style updates --- src/serac/physics/solid_mechanics.cpp | 7 ++-- src/serac/physics/solid_mechanics.hpp | 12 +++--- src/serac/physics/solid_mechanics_contact.hpp | 19 ++++----- .../physics/tests/contact_solid_adjoint.cpp | 39 +++++++++---------- 4 files changed, 36 insertions(+), 41 deletions(-) diff --git a/src/serac/physics/solid_mechanics.cpp b/src/serac/physics/solid_mechanics.cpp index 9001d7b4e..4d49ec49b 100644 --- a/src/serac/physics/solid_mechanics.cpp +++ b/src/serac/physics/solid_mechanics.cpp @@ -15,10 +15,9 @@ namespace detail { void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat, mfem::HypreParMatrix* k_mat, mfem::HypreParVector& disp_adjoint_load_vector, mfem::HypreParVector& velo_adjoint_load_vector, mfem::HypreParVector& accel_adjoint_load_vector, mfem::HypreParVector& adjoint_displacement_, - mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_, - mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_, - BoundaryConditionManager& bcs_, - mfem::Solver& lin_solver) + mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_, + mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_, + BoundaryConditionManager& bcs_, mfem::Solver& lin_solver) { // there are hard-coded here for now static constexpr double beta = 0.25; diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 7b17404d5..7825be899 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -36,10 +36,9 @@ namespace detail { void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat, mfem::HypreParMatrix* k_mat, mfem::HypreParVector& disp_adjoint_load_vector, mfem::HypreParVector& velo_adjoint_load_vector, mfem::HypreParVector& accel_adjoint_load_vector, mfem::HypreParVector& adjoint_displacement_, - mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_, - mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_, - BoundaryConditionManager& bcs_, - mfem::Solver& lin_solver); + mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_, + mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_, + BoundaryConditionManager& bcs_, mfem::Solver& lin_solver); } // namespace detail /** @@ -1193,7 +1192,7 @@ class SolidMechanics, std::integer_se add(1.0, u_, c0_, d2u_dt2, predicted_displacement_); // K := dR/du - auto K = serac::get((*residual_)(time_, shape_displacement_, + auto K = serac::get((*residual_)(time_, shape_displacement_, differentiate_wrt(predicted_displacement_), d2u_dt2, *parameters_[parameter_indices].state...)); std::unique_ptr k_mat(assemble(K)); @@ -1537,7 +1536,6 @@ class SolidMechanics, std::integer_se /// serac::Functional that is used to calculate the residual and its derivatives std::unique_ptr> residual_; - /// mfem::Operator that calculates the residual after applying essential boundary conditions std::unique_ptr residual_with_bcs_; @@ -1624,7 +1622,7 @@ class SolidMechanics, std::integer_se adjoint_essential = 0.0; auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].state...); + *parameters_[parameter_indices].state...); auto jacobian = assemble(drdu); auto J_T = std::unique_ptr(jacobian->Transpose()); diff --git a/src/serac/physics/solid_mechanics_contact.hpp b/src/serac/physics/solid_mechanics_contact.hpp index 1e2a8f5b5..5b483d91f 100644 --- a/src/serac/physics/solid_mechanics_contact.hpp +++ b/src/serac/physics/solid_mechanics_contact.hpp @@ -118,7 +118,7 @@ class SolidMechanicsContact, // tracking strategy // See https://github.com/mfem/mfem/issues/3531 mfem::Vector r_blk(r, 0, displacement_.Size()); - r_blk = res; + r_blk = res; mfem::Vector uPlusShapeDisp = u; uPlusShapeDisp += shape_displacement_; contact_.residualFunction(uPlusShapeDisp, r); @@ -264,19 +264,20 @@ class SolidMechanicsContact, /// @brief Solve the Quasi-static Newton system void quasiStaticAdjointSolve(double /*dt*/) override { - SLIC_ERROR_ROOT_IF(contact_.haveLagrangeMultipliers(), "Lagrange multiplier contact does not currently support sensitivities/adjoints."); + SLIC_ERROR_ROOT_IF(contact_.haveLagrangeMultipliers(), + "Lagrange multiplier contact does not currently support sensitivities/adjoints."); // By default, use a homogeneous essential boundary condition mfem::HypreParVector adjoint_essential(displacement_adjoint_load_); adjoint_essential = 0.0; - + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].state...); + *parameters_[parameter_indices].state...); auto jacobian = assemble(drdu); mfem::Vector uPlusShapeDisp = displacement_; uPlusShapeDisp += shape_displacement_; - auto block_J = contact_.jacobianFunction(uPlusShapeDisp, jacobian.release()); + auto block_J = contact_.jacobianFunction(uPlusShapeDisp, jacobian.release()); block_J->owns_blocks = false; jacobian = std::unique_ptr(static_cast(&block_J->GetBlock(0, 0))); auto J_T = std::unique_ptr(jacobian->Transpose()); @@ -301,10 +302,10 @@ class SolidMechanicsContact, mfem::Vector uPlusShapeDisp = displacement_; uPlusShapeDisp += shape_displacement_; - auto block_J = contact_.jacobianFunction(uPlusShapeDisp, drdshape_mat.release()); + auto block_J = contact_.jacobianFunction(uPlusShapeDisp, drdshape_mat.release()); block_J->owns_blocks = false; drdshape_mat = std::unique_ptr(static_cast(&block_J->GetBlock(0, 0))); - + drdshape_mat->MultTranspose(adjoint_displacement_, *shape_displacement_sensitivity_); return *shape_displacement_sensitivity_; @@ -322,10 +323,10 @@ class SolidMechanicsContact, using BasePhysics::states_; using BasePhysics::time_; using SolidMechanicsBase::acceleration_; + using SolidMechanicsBase::adjoint_displacement_; using SolidMechanicsBase::d_residual_d_; using SolidMechanicsBase::DERIVATIVE; using SolidMechanicsBase::displacement_; - using SolidMechanicsBase::adjoint_displacement_; using SolidMechanicsBase::displacement_adjoint_load_; using SolidMechanicsBase::du_; using SolidMechanicsBase::J_; @@ -334,8 +335,8 @@ class SolidMechanicsContact, using SolidMechanicsBase::ode_time_point_; using SolidMechanicsBase::residual_; using SolidMechanicsBase::residual_with_bcs_; - using SolidMechanicsBase::warmStartDisplacement; using SolidMechanicsBase::time_end_step_; + using SolidMechanicsBase::warmStartDisplacement; /// Pointer to the Jacobian operator (J_ if no Lagrange multiplier contact, J_constraint_ otherwise) mfem::Operator* J_operator_; diff --git a/src/serac/physics/tests/contact_solid_adjoint.cpp b/src/serac/physics/tests/contact_solid_adjoint.cpp index e819cc7a0..0b6f09097 100644 --- a/src/serac/physics/tests/contact_solid_adjoint.cpp +++ b/src/serac/physics/tests/contact_solid_adjoint.cpp @@ -57,28 +57,30 @@ void computeStepAdjointLoad(const FiniteElementState& displacement, FiniteElemen using SolidMechT = serac::SolidMechanicsContact; -std::unique_ptr createContactSolver( - const NonlinearSolverOptions& nonlinear_opts, const TimesteppingOptions& dyn_opts, const SolidMaterial& mat) +std::unique_ptr createContactSolver(const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, const SolidMaterial& mat) { static int iter = 0; - auto solid = - std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, dyn_opts, - geoNonlinear, physics_prefix + std::to_string(iter++), mesh_tag); + auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, dyn_opts, + geoNonlinear, physics_prefix + std::to_string(iter++), mesh_tag); solid->setMaterial(mat); solid->setDisplacementBCs({2}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { disp = 0.0; }); - solid->setDisplacementBCs({4}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { disp = 0.0; disp[1] = -0.1; }); + solid->setDisplacementBCs({4}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { + disp = 0.0; + disp[1] = -0.1; + }); - auto contact_type = serac::ContactEnforcement::Penalty; + auto contact_type = serac::ContactEnforcement::Penalty; double element_length = 1.0; - double penalty = 105.1 * mat.K / element_length; + double penalty = 105.1 * mat.K / element_length; serac::ContactOptions contact_options{.method = serac::ContactMethod::SingleMortar, .enforcement = contact_type, .type = serac::ContactType::Frictionless, .penalty = penalty}; - auto contact_interaction_id = 0; + auto contact_interaction_id = 0; solid->addContactInteraction(contact_interaction_id, {3}, {5}, contact_options); solid->completeSetup(); @@ -96,8 +98,7 @@ double computeSolidMechanicsQoi(BasePhysics& solid_solver, const TimeSteppingInf return computeStepQoi(solid_solver.state("displacement")); } -auto computeContactQoiSensitivities( - BasePhysics& solid_solver, const TimeSteppingInfo& ts_info) +auto computeContactQoiSensitivities(BasePhysics& solid_solver, const TimeSteppingInfo& ts_info) { EXPECT_EQ(0, solid_solver.cycle()); @@ -117,7 +118,6 @@ auto computeContactQoiSensitivities( return std::make_tuple(qoi, shape_sensitivity); } - double computeSolidMechanicsQoiAdjustingShape(SolidMechanics& solid_solver, const TimeSteppingInfo& ts_info, const FiniteElementState& shape_derivative_direction, double pertubation) { @@ -129,7 +129,6 @@ double computeSolidMechanicsQoiAdjustingShape(SolidMechanics& solid_solv return computeSolidMechanicsQoi(solid_solver, ts_info); } - struct ContactSensitivityFixture : public ::testing::Test { void SetUp() override { @@ -138,9 +137,9 @@ struct ContactSensitivityFixture : public ::testing::Test { std::string filename = std::string(SERAC_REPO_DIR) + "/data/meshes/contact_two_blocks.g"; mesh = &StateManager::setMesh(mesh::refineAndDistribute(buildMeshFromFile(filename), 0), mesh_tag); - mat.density = 1.0; - mat.K = 1.0; - mat.G = 0.1; + mat.density = 1.0; + mat.K = 1.0; + mat.G = 0.1; } void fillDirection(FiniteElementState& direction) const @@ -157,7 +156,7 @@ struct ContactSensitivityFixture : public ::testing::Test { NonlinearSolverOptions nonlinear_opts{.relative_tol = 1.0e-10, .absolute_tol = 1.0e-12}; bool dispBc = true; - TimesteppingOptions dyn_opts{.timestepper = TimestepMethod::QuasiStatic}; + TimesteppingOptions dyn_opts{.timestepper = TimestepMethod::QuasiStatic}; SolidMaterial mat; TimeSteppingInfo tsInfo; @@ -165,7 +164,6 @@ struct ContactSensitivityFixture : public ::testing::Test { static constexpr double eps = 2e-7; }; - TEST_F(ContactSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSameObjectiveAndGradient) { auto solid_solver = createContactSolver(nonlinear_opts, dyn_opts, mat); @@ -184,7 +182,6 @@ TEST_F(ContactSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSameObjec EXPECT_EQ(directional_deriv1, directional_deriv2); } - TEST_F(ContactSensitivityFixture, QuasiStaticShapeSensitivities) { auto solid_solver = createContactSolver(nonlinear_opts, dyn_opts, mat); @@ -196,10 +193,10 @@ TEST_F(ContactSensitivityFixture, QuasiStaticShapeSensitivities) double qoi_plus = computeSolidMechanicsQoiAdjustingShape(*solid_solver, tsInfo, derivative_direction, eps); - double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); + double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); double directional_deriv_fd = (qoi_plus - qoi_base) / eps; // std::cout << "qoi, derivs = " << qoi_base << " " << directional_deriv << " " << directional_deriv_fd << std::endl; - EXPECT_NEAR(directional_deriv, directional_deriv_fd, 0.2); // These are very large tolerances + EXPECT_NEAR(directional_deriv, directional_deriv_fd, 0.2); // These are very large tolerances } } // namespace serac