Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Robust gradual incrementing of boundary conditions. #1252

Open
wants to merge 9 commits into
base: develop
Choose a base branch
from
Open

Robust gradual incrementing of boundary conditions. #1252

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a36e9dd
fix parallel issue.
tupek2 Oct 17, 2024
70e5a3a
Add option to symmetrize the linear system from the physics side.
tupek2 Oct 24, 2024
0d965d2
Add catches to cut back load step internally whenever the nonlinear s…
tupek2 Oct 24, 2024
8910405
Cleanup output, function for computing matrix norms.
tupek2 Oct 25, 2024
3d71efa
revert a cmake change.
tupek2 Oct 28, 2024
c955801
Tidy up the more robust boundary conditions a bit.
tupek2 Oct 28, 2024
5fa7df0
Apply style updates
Oct 28, 2024
ff8e419
Merge branch 'develop' into tupek/reaction_sensitivities_with_robust_bcs
tupek2 Oct 28, 2024
29c922e
Merge branch 'develop' into tupek/reaction_sensitivities_with_robust_bcs
tupek2 Oct 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 39 additions & 19 deletions src/serac/numerics/equation_solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,16 @@

namespace serac {

class SolveException : public std::exception {
public:
SolveException(const std::string& message) : msg(message) {}

const char* what() const noexcept override { return msg.c_str(); }

private:
std::string msg;
};

/// Newton solver with a 2-way line-search. Reverts to regular Newton if max_line_search_iterations is set to 0.
class NewtonSolver : public mfem::NewtonSolver {
protected:
Expand Down Expand Up @@ -107,8 +117,11 @@ class NewtonSolver : public mfem::NewtonSolver {
}

if (norm != norm) {
mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
mfem::out << "Newton: No convergence!\n";
if (print_options.first_and_last) {
mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
mfem::out << "Newton: No convergence!\n";
}
throw SolveException("Initial residual for Newton iteration is undefined/nan.");
return;
}

Expand Down Expand Up @@ -195,9 +208,14 @@ class NewtonSolver : public mfem::NewtonSolver {
if (print_options.summary || (!converged && print_options.warnings) || print_options.first_and_last) {
mfem::out << "Newton: Number of iterations: " << final_iter << '\n' << " ||r|| = " << final_norm << '\n';
}

if (!converged && (print_options.summary || print_options.warnings)) {
mfem::out << "Newton: No convergence!\n";
}

if (!converged) {
throw SolveException("Newton algorithm failed to converge.");
}
}
};

Expand Down Expand Up @@ -522,14 +540,17 @@ class TrustRegion : public mfem::NewtonSolver {
mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm);
mfem::out << ", x_incr = " << std::setw(13) << d.Norml2();
} else {
mfem::out << ", norm goal = " << std::setw(13) << norm_goal << "\n";
mfem::out << ", norm goal = " << std::setw(13) << norm_goal;
}
mfem::out << '\n';
mfem::out << " :: ";
}

if (norm != norm) {
mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
mfem::out << "Newton: No convergence!\n";
if (print_options.first_and_last) {
mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
mfem::out << "Newton: No convergence!\n";
}
throw SolveException("Initial residual for trust-region iteration is undefined/nan.");
return;
}

Expand Down Expand Up @@ -595,7 +616,7 @@ class TrustRegion : public mfem::NewtonSolver {

doglegStep(trResults.cauchy_point, trResults.z, tr_size, d);

static constexpr double roundOffTol = 0.0; // 1e-14;
static constexpr double roundOffTol = 1e-15;

hess_vec_func(d, Hd);
double dHd = Dot(d, Hd);
Expand All @@ -606,15 +627,9 @@ class TrustRegion : public mfem::NewtonSolver {
double realObjective = std::numeric_limits<double>::max();
double normPred = std::numeric_limits<double>::max();
try {
normPred = computeResidual(x_pred, r_pred);
double obj1 = 0.5 * (Dot(r, d) + Dot(r_pred, d)) - roundOffTol;

// midstep work estimate
// add(X, 0.5, d, x_mid);
// computeResidual(x_mid, r_mid);
// double obj2 = Dot(r_mid, d);

realObjective = obj1; // 0.5 * obj1 + 0.5 * obj2;
normPred = computeResidual(x_pred, r_pred);
double obj1 = 0.5 * (Dot(r, d) + Dot(r_pred, d)) - roundOffTol;
realObjective = obj1;
} catch (const std::exception&) {
realObjective = std::numeric_limits<double>::max();
normPred = std::numeric_limits<double>::max();
Expand All @@ -627,8 +642,8 @@ class TrustRegion : public mfem::NewtonSolver {
happyAboutTrSize = true;
if (print_options.iterations) {
printTrustRegionInfo(realObjective, modelObjective, trResults.cg_iterations_count, tr_size, true);
trResults.cg_iterations_count =
0; // zero this output so it doesn't look like the linesearch is doing cg iterations
// zero this output so it doesn't look like the linesearch is doing cg iterations
trResults.cg_iterations_count = 0;
}
break;
}
Expand All @@ -641,7 +656,7 @@ class TrustRegion : public mfem::NewtonSolver {
if (print_options.iterations || print_options.warnings) {
mfem::out << "Found a positive model objective increase. Debug if you see this.\n";
}
rho = realImprove / -modelImprove;
rho = 0.0; // realImprove / -modelImprove;
}

if (!(rho >= settings.eta2)) { // write it this way to handle NaNs
Expand Down Expand Up @@ -675,9 +690,14 @@ class TrustRegion : public mfem::NewtonSolver {
final_iter = it;
final_norm = norm;

if (!converged) {
throw SolveException("Trust-region algorithm failed to converge.");
}

if (print_options.summary || (!converged && print_options.warnings) || print_options.first_and_last) {
mfem::out << "Newton: Number of iterations: " << final_iter << '\n' << " ||r|| = " << final_norm << '\n';
}

if (!converged && (print_options.summary || print_options.warnings)) {
mfem::out << "Newton: No convergence!\n";
}
Expand Down
6 changes: 5 additions & 1 deletion src/serac/physics/solid_mechanics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_,
mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_,
mfem::HypreParVector& adjoint_essential, BoundaryConditionManager& bcs_,
mfem::Solver& lin_solver)
mfem::Solver& lin_solver, bool symmetrize)
{
// there are hard-coded here for now
static constexpr double beta = 0.25;
Expand All @@ -36,6 +36,10 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
auto J_ = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(1.0, *m_mat, fac3 * dt_n * dt_n, *k_mat));
auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());

if (symmetrize) {
J_T = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
}

// recall that temperature_adjoint_load_vector and d_temperature_dt_adjoint_load_vector were already multiplied by
// -1 above

Expand Down
105 changes: 87 additions & 18 deletions src/serac/physics/solid_mechanics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,17 @@

namespace serac {

static constexpr bool debug_matrix_symmetry = false;

inline double matrixNorm(std::unique_ptr<mfem::HypreParMatrix>& K)
{
mfem::HypreParMatrix* H = K.get();
hypre_ParCSRMatrix* Hhypre = static_cast<hypre_ParCSRMatrix*>(*H);
double Hfronorm;
hypre_ParCSRMatrixNormFro(Hhypre, &Hfronorm);
return Hfronorm;
}

namespace solid_mechanics {

namespace detail {
Expand All @@ -39,7 +50,7 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_,
mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_,
mfem::HypreParVector& adjoint_essential, BoundaryConditionManager& bcs_,
mfem::Solver& lin_solver);
mfem::Solver& lin_solver, bool symmetrize);
} // namespace detail

/**
Expand Down Expand Up @@ -875,7 +886,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
}

auto flux = dot(stress, transpose(inv(dx_dX))) * det(dx_dX);

return serac::tuple{material_.density * d2u_dt2, flux};
}
};
Expand Down Expand Up @@ -1181,8 +1191,25 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
*parameters_[parameter_indices].state...);
J_.reset();
J_ = assemble(drdu);

if constexpr (debug_matrix_symmetry) {
auto Jt = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
Jt = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, -0.5, *Jt));
double n1 = matrixNorm(J_);
double n2 = matrixNorm(Jt);
if (mpi_rank_ == 0) {
std::cout << "matrix norm of J = " << n1 << ", matrix norm of skew(J) = " << n2 << std::endl;
}
}

if (symmetrize) {
auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
J_ = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
}

J_e_.reset();
J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_);

return *J_;
});
}
Expand Down Expand Up @@ -1250,6 +1277,12 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
// J = M + c0 * K
J_.reset();
J_.reset(mfem::Add(1.0, *m_mat, c0_, *k_mat));

if (symmetrize) {
auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
J_ = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
}

J_e_.reset();
J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_);

Expand Down Expand Up @@ -1302,7 +1335,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
}

if (is_quasistatic_) {
quasiStaticSolve(dt);
quasiStaticSolve(dt, 0.0, 1.0);
} else {
// The current ode interface tracks 2 times, one internally which we have a handle to via time_,
// and one here via the step interface.
Expand Down Expand Up @@ -1427,6 +1460,9 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
J_ = assemble(drdu);

auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
if (symmetrize) {
J_T = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
}

J_e_.reset();
J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T);
Expand Down Expand Up @@ -1467,7 +1503,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
solid_mechanics::detail::adjoint_integrate(
dt_n_to_np1, dt_np1_to_np2, m_mat.get(), k_mat.get(), displacement_adjoint_load_, velocity_adjoint_load_,
acceleration_adjoint_load_, adjoint_displacement_, implicit_sensitivity_displacement_start_of_step_,
implicit_sensitivity_velocity_start_of_step_, reactions_adjoint_bcs_, bcs_, lin_solver);
implicit_sensitivity_velocity_start_of_step_, reactions_adjoint_bcs_, bcs_, lin_solver, symmetrize);
}

time_end_step_ = time_;
Expand Down Expand Up @@ -1626,6 +1662,9 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
/// @brief A flag denoting whether to compute geometric nonlinearities in the residual
GeometricNonlinearities geom_nonlin_;

/// @brief A flag denoting whether to symmetrize linear operators passed to the solvers
bool symmetrize = false; // true;

/// @brief A flag denoting whether to compute the warm start for improved robustness
bool use_warm_start_;

Expand All @@ -1647,16 +1686,39 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
}...};

/// @brief Solve the Quasi-static Newton system
virtual void quasiStaticSolve(double dt)
virtual void quasiStaticSolve(double dt, double a, double b, int level = 0)
{
// warm start must be called prior to the time update so that the previous Jacobians can be used consistently
// throughout.
warmStartDisplacement(dt);
time_ += dt;

// this method is essentially equivalent to the 1-liner
// u += dot(inv(J), dot(J_elim[:, dofs], (U(t + dt) - u)[dofs]));
nonlin_solver_->solve(displacement_);
if (level >= 6) {
if (mpi_rank_ == 0)
std::cout << "Too many boundary condition cutbacks, accepting solution even though there may be issues. Try "
"increasing the number of load steps "
<< std::endl;
return;
}

bool solver_success = true;
try {
// warm start must be called prior to the time update so that the previous Jacobians can be used consistently
// throughout.
if (a == 0.0) time_ += dt;
warmStartDisplacement(dt, b);
nonlin_solver_->solve(displacement_);
} catch (const std::exception& e) {
if (mpi_rank_ == 0) mfem::out << "caught nonlinear solver exception: " << e.what() << std::endl;
displacement_ -= du_;
solver_success = false;
quasiStaticSolve(dt, 1.0, 0.5 * b, level + 1);
quasiStaticSolve(dt, 1.0, 1.0, level + 1);
}

if (solver_success) {
if (b == 1.0) {
if (mpi_rank_ == 0)
mfem::out << "SolidMechanics solve succeeded for time " << time_ << " dt = " << dt << std::endl;
} else {
if (mpi_rank_ == 0) mfem::out << "substep solve succeeded for time " << time_ << " dt = " << dt << std::endl;
}
}
}

/**
Expand Down Expand Up @@ -1754,14 +1816,14 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
nonlinear solvers will ensure positive definiteness at equilibrium. *Once any parameter is changed, it is no longer
certain to be positive definite, which will cause issues for many types linear solvers.
*/
void warmStartDisplacement(double dt)
void warmStartDisplacement(double dt, double displacement_scale_factor)
{
SERAC_MARK_FUNCTION;

du_ = 0.0;
for (auto& bc : bcs_.essentials()) {
// apply the future boundary conditions, but use the most recent Jacobians stiffness.
bc.setDofs(du_, time_ + dt);
bc.setDofs(du_, time_);
}

auto& constrained_dofs = bcs_.allEssentialTrueDofs();
Expand All @@ -1770,16 +1832,24 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
du_[j] -= displacement_(j);
}

du_ *= displacement_scale_factor;

if (use_warm_start_) {
// Update the linearized Jacobian matrix
auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_,
auto r = (*residual_)(time_, shape_displacement_, displacement_, acceleration_,
*parameters_[parameter_indices].state...);

// use the most recently evaluated Jacobian
auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
auto [_, drdu] = (*residual_)(time_ - dt, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
*parameters_[parameter_indices].previous_state...);
J_.reset();
J_ = assemble(drdu);

if (symmetrize) {
auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
J_ = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
}

J_e_.reset();
J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_);

Expand All @@ -1794,7 +1864,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
auto& lin_solver = nonlin_solver_->linearSolver();

lin_solver.SetOperator(*J_);

lin_solver.Mult(r, du_);
}

Expand Down
9 changes: 6 additions & 3 deletions src/serac/physics/tests/dynamic_solid_adjoint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ std::unique_ptr<SolidMechanics<p, dim>> createNonlinearSolidMechanicsSolver(
.preconditioner = Preconditioner::HypreJacobi,
.relative_tol = 1.0e-9,
.absolute_tol = 1.0e-16,
.max_iterations = 2000,
.max_iterations = 200,
.print_level = 0};

bool checkpoint_to_disk = true;
Expand Down Expand Up @@ -258,8 +258,11 @@ struct SolidMechanicsSensitivityFixture : public ::testing::Test {
axom::sidre::DataStore dataStore;
mfem::ParMesh* mesh;

NonlinearSolverOptions nonlinear_opts{
.nonlin_solver = NonlinearSolver::TrustRegion, .relative_tol = 1.0e-15, .absolute_tol = 1.0e-14};
NonlinearSolverOptions nonlinear_opts{.nonlin_solver = NonlinearSolver::TrustRegion,
.relative_tol = 1.0e-14,
.absolute_tol = 5.0e-15,
.max_iterations = 50,
.print_level = 0};

bool dispBc = true;
TimesteppingOptions dyn_opts{.timestepper = TimestepMethod::Newmark,
Expand Down