Merge pull request #1198 from LLNL/atkins12/contact-bug

Add cylinder example and contact force monitoring

examples/CMakeLists.txt

@@ -40,3 +40,13 @@ install(
     DESTINATION
         examples/serac/contact
 )
+
+# Add the buckling examples
+add_subdirectory(buckling)
+
+install(
+    FILES
+        buckling/cylinder.cpp
+    DESTINATION
+        examples/serac/buckling
+)

examples/buckling/CMakeLists.txt

@@ -0,0 +1,13 @@
+# Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+# other Serac Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+if(TRIBOL_FOUND)
+    blt_add_executable( NAME        buckling_cylinder
+                        SOURCES     cylinder.cpp
+                        OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
+                        DEPENDS_ON  serac_physics serac_mesh
+                        )
+endif()

examples/buckling/cylinder.cpp

@@ -0,0 +1,211 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+// other Serac Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/**
+ * @file cylinder.cpp
+ *
+ * @brief A buckling cylinder under compression, run with or without contact
+ *
+ * @note Run with mortar contact and PETSc preconditioners:
+ * @code{.sh}
+ * ./build/examples/buckling_cylinder --contact-type 1 --preconditioner 6 \
+ *     -options_file examples/buckling/cylinder_petsc_options.yml
+ * @endcode
+ * @note Run with penalty contact and HYPRE BoomerAMG preconditioner
+ * @code{.sh}
+ * ./build/examples/buckling_cylinder --penalty 1e3
+ * @endcode
+ * @note Run without contact:
+ * @code{.sh}
+ * ./build/examples/buckling_cylinder --no-contact
+ * @endcode
+ */
+
+#include <set>
+#include <string>
+
+#include "axom/slic/core/SimpleLogger.hpp"
+#include "axom/inlet.hpp"
+
+#include "mfem.hpp"
+
+#include "serac/physics/solid_mechanics_contact.hpp"
+#include "serac/infrastructure/terminator.hpp"
+#include "serac/mesh/mesh_utils.hpp"
+#include "serac/physics/state/state_manager.hpp"
+#include "serac/physics/materials/parameterized_solid_material.hpp"
+#include "serac/serac_config.hpp"
+
+using namespace serac;
+
+std::function<std::string(const std::string&)> petscPCTypeValidator = [](const std::string& in) -> std::string {
+  return std::to_string(static_cast<int>(mfem_ext::stringToPetscPCType(in)));
+};
+
+/**
+ * @brief Run buckling cylinder example
+ *
+ * @note Based on doi:10.1016/j.cma.2014.08.012
+ */
+int main(int argc, char* argv[])
+{
+  constexpr int dim = 3;
+  constexpr int p   = 1;
+
+  // Command line arguments
+  // Mesh options
+  int    serial_refinement   = 0;
+  int    parallel_refinement = 0;
+  double dt                  = 0.1;
+
+  // Solver options
+  NonlinearSolverOptions nonlinear_options     = solid_mechanics::default_nonlinear_options;
+  LinearSolverOptions    linear_options        = solid_mechanics::default_linear_options;
+  nonlinear_options.nonlin_solver              = serac::NonlinearSolver::TrustRegion;
+  nonlinear_options.relative_tol               = 1e-6;
+  nonlinear_options.absolute_tol               = 1e-10;
+  nonlinear_options.min_iterations             = 1;
+  nonlinear_options.max_iterations             = 500;
+  nonlinear_options.max_line_search_iterations = 20;
+  nonlinear_options.print_level                = 1;
+#ifdef SERAC_USE_PETSC
+  linear_options.linear_solver  = serac::LinearSolver::GMRES;
+  linear_options.preconditioner = serac::Preconditioner::HypreAMG;
+  linear_options.relative_tol   = 1e-8;
+  linear_options.absolute_tol   = 1e-16;
+  linear_options.max_iterations = 2000;
+#endif
+
+  // Contact specific options
+  double penalty      = 1e3;
+  bool   use_contact  = true;
+  auto   contact_type = serac::ContactEnforcement::Penalty;
+
+  // Initialize Serac and all of the external libraried
+  serac::initialize(argc, argv);
+
+  // Handle command line arguments
+  axom::CLI::App app{"Hollow cylinder buckling example"};
+  // Mesh options
+  app.add_option("--serial-refinement", serial_refinement, "Serial refinement steps", true);
+  app.add_option("--parallel-refinement", parallel_refinement, "Parallel refinement steps", true);
+  // Solver options
+  app.add_option("--nonlinear-solver", nonlinear_options.nonlin_solver, "Nonlinear solver", true);
+  app.add_option("--linear-solver", linear_options.linear_solver, "Linear solver", true);
+  app.add_option("--preconditioner", linear_options.preconditioner, "Preconditioner", true);
+  app.add_option("--petsc-pc-type", linear_options.petsc_preconditioner, "Petsc preconditioner", true)
+      ->transform(
+          [](const std::string& in) -> std::string {
+            return std::to_string(static_cast<int>(mfem_ext::stringToPetscPCType(in)));
+          },
+          "Convert string to PetscPCType", "PetscPCTypeTransform");
+  app.add_option("--dt", dt, "Size of pseudo-time step pre-contact", true);
+  // Contact options
+  app.add_flag("--contact,!--no-contact", use_contact, "Use contact for the inner faces of the cylinder");
+  app.add_option("--contact-type", contact_type,
+                 "Type of contact enforcement, 0 for penalty or 1 for Lagrange multipliers", true)
+      ->needs("--contact");
+  app.add_option("--penalty", penalty, "Penalty for contact", true)->needs("--contact");
+
+  // Need to allow extra arguments for PETSc support
+  app.set_help_flag("--help");
+  app.allow_extras()->parse(argc, argv);
+
+  nonlinear_options.force_monolithic = linear_options.preconditioner != Preconditioner::Petsc;
+
+  // Create DataStore
+  std::string            name     = use_contact ? "buckling_cylinder_contact" : "buckling_cylinder";
+  std::string            mesh_tag = "mesh";
+  axom::sidre::DataStore datastore;
+  serac::StateManager::initialize(datastore, name + "_data");
+
+  // Create and refine mesh
+  std::string filename = SERAC_REPO_DIR "/data/meshes/hollow-cylinder.mesh";
+  auto        mesh     = serac::buildMeshFromFile(filename);
+  auto        pmesh    = mesh::refineAndDistribute(std::move(mesh), serial_refinement, parallel_refinement);
+  serac::StateManager::setMesh(std::move(pmesh), mesh_tag);
+
+  // Surface attributes for boundary conditions
+  std::set<int> xneg{2};
+  std::set<int> xpos{3};
+  std::set<int> bottom{1};
+  std::set<int> top{4};
+
+  // Create solver, either with or without contact
+  std::unique_ptr<SolidMechanics<p, dim>> solid_solver;
+  if (use_contact) {
+    auto solid_contact_solver = std::make_unique<serac::SolidMechanicsContact<p, dim>>(
+        nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options,
+        serac::GeometricNonlinearities::On, name, mesh_tag);
+
+    // Add the contact interaction
+    serac::ContactOptions contact_options{.method      = serac::ContactMethod::SingleMortar,
+                                          .enforcement = contact_type,
+                                          .type        = serac::ContactType::Frictionless,
+                                          .penalty     = penalty};
+    auto                  contact_interaction_id = 0;
+    solid_contact_solver->addContactInteraction(contact_interaction_id, xpos, xneg, contact_options);
+    solid_solver = std::move(solid_contact_solver);
+  } else {
+    solid_solver = std::make_unique<serac::SolidMechanics<p, dim>>(nonlinear_options, linear_options,
+                                                                   serac::solid_mechanics::default_quasistatic_options,
+                                                                   serac::GeometricNonlinearities::On, name, mesh_tag);
+    auto domain  = serac::Domain::ofBoundaryElements(
+         StateManager::mesh(mesh_tag), [&](std::vector<vec3>, int attr) { return xpos.find(attr) != xpos.end(); });
+    solid_solver->setPressure([&](auto&, double t) { return 0.01 * t; }, domain);
+  }
+
+  // Define a Neo-Hookean material
+  auto                        lambda = 1.0;
+  auto                        G      = 0.1;
+  solid_mechanics::NeoHookean mat{.density = 1.0, .K = (3 * lambda + 2 * G) / 3, .G = G};
+
+  solid_solver->setMaterial(mat);
+
+  // Set up essential boundary conditions
+  // Bottom of cylinder is fixed
+  auto clamp = [](const mfem::Vector&, mfem::Vector& u) {
+    u.SetSize(dim);
+    u = 0.0;
+  };
+  solid_solver->setDisplacementBCs(bottom, clamp);
+
+  // Top of cylinder has prescribed displacement of magnitude in x-z direction
+  auto compress = [&](const mfem::Vector&, double t, mfem::Vector& u) {
+    u.SetSize(dim);
+    u    = 0.0;
+    u(0) = u(2) = -1.5 / std::sqrt(2.0) * t;
+  };
+  solid_solver->setDisplacementBCs(top, compress);
+
+  // Finalize the data structures
+  solid_solver->completeSetup();
+
+  // Save initial state
+  std::string paraview_name = name + "_paraview";
+  solid_solver->outputStateToDisk(paraview_name);
+
+  // Perform the quasi-static solve
+  SLIC_INFO_ROOT(axom::fmt::format("Running hollow cylinder bucking example with {} displacement dofs",
+                                   solid_solver->displacement().GlobalSize()));
+  SLIC_INFO_ROOT("Starting pseudo-timestepping.");
+  serac::logger::flush();
+  while (solid_solver->time() < 1.0 && std::abs(solid_solver->time() - 1) > DBL_EPSILON) {
+    SLIC_INFO_ROOT(axom::fmt::format("time = {} (out of 1.0)", solid_solver->time()));
+    serac::logger::flush();
+
+    // Refine dt as contact starts
+    auto next_dt = solid_solver->time() < 0.65 ? dt : dt * 0.1;
+    solid_solver->advanceTimestep(next_dt);
+
+    // Output the sidre-based plot files
+    solid_solver->outputStateToDisk(paraview_name);
+  }
+  SLIC_INFO_ROOT(axom::fmt::format("final time = {}", solid_solver->time()));
+
+  // Exit without error
+  serac::exitGracefully();
+}

examples/buckling/cylinder_petsc_options.yml

@@ -0,0 +1,21 @@
+ksp:
+  converged_reason:
+  type: fgmres
+
+snes:
+  converged_reason:
+
+nestpc:
+  type: fieldsplit
+  fieldsplit:
+    detect_saddle_point:
+    type: schur
+    schur_fact_type: diag
+    schur_precondition: a11
+    # type: multiplicative
+nestfieldsplit:
+  0:
+    ksp_max_it: 10
+  1:
+    pc_type: jacobi
+    # pc_type: jacobi

src/serac/numerics/equation_solver.cpp

@@ -57,6 +57,10 @@ class NewtonSolver : public mfem::NewtonSolver {
   {
     SERAC_MARK_FUNCTION;
     grad = &oper->GetGradient(x);
+    if (nonlinear_options.force_monolithic) {
+      auto* grad_blocked = dynamic_cast<mfem::BlockOperator*>(grad);
+      if (grad_blocked) grad = buildMonolithicMatrix(*grad_blocked).release();
+    }
   }
 
   /// set the preconditioner for the linear solver
@@ -428,6 +432,10 @@ class TrustRegion : public mfem::NewtonSolver {
   {
     SERAC_MARK_FUNCTION;
     grad = &oper->GetGradient(x);
+    if (nonlinear_options.force_monolithic) {
+      auto* grad_blocked = dynamic_cast<mfem::BlockOperator*>(grad);
+      if (grad_blocked) grad = buildMonolithicMatrix(*grad_blocked).release();
+    }
   }
 
   /// evaluate the nonlinear residual
@@ -685,14 +693,6 @@ void SuperLUSolver::Mult(const mfem::Vector& input, mfem::Vector& output) const
   superlu_solver_.Mult(input, output);
 }
 
-/**
- * @brief Function for building a monolithic parallel Hypre matrix from a block system of smaller Hypre matrices
- *
- * @param block_operator The block system of HypreParMatrices
- * @return The assembled monolithic HypreParMatrix
- *
- * @pre @a block_operator must have assembled HypreParMatrices for its sub-blocks
- */
 std::unique_ptr<mfem::HypreParMatrix> buildMonolithicMatrix(const mfem::BlockOperator& block_operator)
 {
   int row_blocks = block_operator.NumRowBlocks();

src/serac/numerics/equation_solver.hpp

@@ -250,6 +250,16 @@ class StrumpackSolver : public mfem::Solver {
 
 #endif
 
+/**
+ * @brief Function for building a monolithic parallel Hypre matrix from a block system of smaller Hypre matrices
+ *
+ * @param block_operator The block system of HypreParMatrices
+ * @return The assembled monolithic HypreParMatrix
+ *
+ * @pre @a block_operator must have assembled HypreParMatrices for its sub-blocks
+ */
+std::unique_ptr<mfem::HypreParMatrix> buildMonolithicMatrix(const mfem::BlockOperator& block_operator);
+
 /**
  * @brief Build a nonlinear solver using the nonlinear option struct
  *

src/serac/numerics/solver_config.hpp

@@ -392,6 +392,9 @@ struct NonlinearSolverOptions {
 
   /// Debug print level
   int print_level = 0;
+
+  /// Should the gradient be converted to a monolithic matrix
+  bool force_monolithic = false;
 };
 // _nonlinear_options_end
 

src/serac/physics/contact/contact_data.cpp

@@ -43,6 +43,9 @@ void ContactData::addContactInteraction(int interaction_id, const std::set<int>&
 
 void ContactData::update(int cycle, double time, double& dt)
 {
+  cycle_ = cycle;
+  time_  = time;
+  dt_    = dt;
   // This updates the redecomposed surface mesh based on the current displacement, then transfers field quantities to
   // the updated mesh.
   tribol::updateMfemParallelDecomposition();
@@ -236,14 +239,13 @@ void ContactData::residualFunction(const mfem::Vector& u, mfem::Vector& r)
   mfem::Vector r_blk(r, 0, disp_size);
   mfem::Vector g_blk(r, disp_size, numPressureDofs());
 
-  double dt = 1.0;
   setDisplacements(u_blk);
   // we need to call update first to update gaps
-  update(1, 1.0, dt);
+  update(cycle_, time_, dt_);
   // with updated gaps, we can update pressure for contact interactions with penalty enforcement
   setPressures(p_blk);
   // call update again with the right pressures
-  update(1, 1.0, dt);
+  update(cycle_, time_, dt_);
 
   r_blk += forces();
   // calling mergedGaps() with true will zero out gap on inactive dofs (so the residual converges and the linearized

src/serac/physics/contact/contact_data.hpp

@@ -265,6 +265,10 @@ class ContactData {
    * @note This is mutable so it can be updated when pressures/gaps/Jacobians are retrieved.
    */
   mutable mfem::Array<HYPRE_BigInt> global_pressure_dof_offsets_;
+
+  int    cycle_{0};
+  double time_{0.0};
+  double dt_{1.0};
 };
 
 }  // namespace serac