Skip to content

Commit

Permalink
Apply style updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@tupek2
Agent Style authored and tupek2 committed Sep 13, 2024
1 parent 730bbc6 commit 9cdca8a
Show file tree
Hide file tree
Showing 4 changed files with 36 additions and 41 deletions.
7 changes: 3 additions & 4 deletions src/serac/physics/solid_mechanics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
12 changes: 5 additions & 7 deletions src/serac/physics/solid_mechanics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

/**
Expand Down Expand Up @@ -1193,7 +1192,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
add(1.0, u_, c0_, d2u_dt2, predicted_displacement_);

// K := dR/du
auto K = serac::get<DERIVATIVE>((*residual_)(time_, shape_displacement_,
auto K = serac::get<DERIVATIVE>((*residual_)(time_, shape_displacement_,
differentiate_wrt(predicted_displacement_), d2u_dt2,
*parameters_[parameter_indices].state...));
std::unique_ptr<mfem::HypreParMatrix> k_mat(assemble(K));
Expand Down Expand Up @@ -1537,7 +1536,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
/// serac::Functional that is used to calculate the residual and its derivatives
std::unique_ptr<ShapeAwareFunctional<shape_trial, test(trial, trial, parameter_space...)>> residual_;


/// mfem::Operator that calculates the residual after applying essential boundary conditions
std::unique_ptr<mfem_ext::StdFunctionOperator> residual_with_bcs_;

Expand Down Expand Up @@ -1624,7 +1622,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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<mfem::HypreParMatrix>(jacobian->Transpose());

Expand Down
19 changes: 10 additions & 9 deletions src/serac/physics/solid_mechanics_contact.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
// 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);
Expand Down Expand Up @@ -264,19 +264,20 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
/// @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<mfem::HypreParMatrix>(static_cast<mfem::HypreParMatrix*>(&block_J->GetBlock(0, 0)));
auto J_T = std::unique_ptr<mfem::HypreParMatrix>(jacobian->Transpose());
Expand All @@ -301,10 +302,10 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,

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<mfem::HypreParMatrix>(static_cast<mfem::HypreParMatrix*>(&block_J->GetBlock(0, 0)));

drdshape_mat->MultTranspose(adjoint_displacement_, *shape_displacement_sensitivity_);

return *shape_displacement_sensitivity_;
Expand All @@ -322,10 +323,10 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
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_;
Expand All @@ -334,8 +335,8 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
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_;
Expand Down
39 changes: 18 additions & 21 deletions src/serac/physics/tests/contact_solid_adjoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,28 +57,30 @@ void computeStepAdjointLoad(const FiniteElementState& displacement, FiniteElemen

using SolidMechT = serac::SolidMechanicsContact<p, dim>;

std::unique_ptr<SolidMechT> createContactSolver(
const NonlinearSolverOptions& nonlinear_opts, const TimesteppingOptions& dyn_opts, const SolidMaterial& mat)
std::unique_ptr<SolidMechT> createContactSolver(const NonlinearSolverOptions& nonlinear_opts,
const TimesteppingOptions& dyn_opts, const SolidMaterial& mat)
{
static int iter = 0;

auto solid =
std::make_unique<SolidMechT>(nonlinear_opts, solid_mechanics::direct_linear_options, dyn_opts,
geoNonlinear, physics_prefix + std::to_string(iter++), mesh_tag);
auto solid = std::make_unique<SolidMechT>(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();
Expand All @@ -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());

Expand All @@ -117,7 +118,6 @@ auto computeContactQoiSensitivities(
return std::make_tuple(qoi, shape_sensitivity);
}


double computeSolidMechanicsQoiAdjustingShape(SolidMechanics<p, dim>& solid_solver, const TimeSteppingInfo& ts_info,
const FiniteElementState& shape_derivative_direction, double pertubation)
{
Expand All @@ -129,7 +129,6 @@ double computeSolidMechanicsQoiAdjustingShape(SolidMechanics<p, dim>& solid_solv
return computeSolidMechanicsQoi(solid_solver, ts_info);
}


struct ContactSensitivityFixture : public ::testing::Test {
void SetUp() override
{
Expand All @@ -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
Expand All @@ -157,15 +156,14 @@ 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;

static constexpr double eps = 2e-7;
};


TEST_F(ContactSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSameObjectiveAndGradient)
{
auto solid_solver = createContactSolver(nonlinear_opts, dyn_opts, mat);
Expand All @@ -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);
Expand All @@ -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
Expand Down

0 comments on commit 9cdca8a

Please sign in to comment.