Merge branch 'develop' into tupek/bugfix/warm-start-solves-linear-pro…

…blems

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/functional/dof_numbering.hpp

@@ -277,15 +277,18 @@ struct GradientAssemblyLookupTables {
     };
   };
 
+  /// dummy default ctor to enable deferred initialization
+  GradientAssemblyLookupTables() : initialized{false} {};
+
   /**
    * @param block_test_dofs object containing information about dofs for the test space
    * @param block_trial_dofs object containing information about dofs for the trial space
    *
    * @brief create lookup tables describing which degrees of freedom
    * correspond to each domain/boundary element
    */
-  GradientAssemblyLookupTables(const serac::BlockElementRestriction& block_test_dofs,
-                               const serac::BlockElementRestriction& block_trial_dofs)
+  void init(const serac::BlockElementRestriction& block_test_dofs,
+            const serac::BlockElementRestriction& block_trial_dofs)
   {
     // we start by having each element and boundary element emit the (i,j) entry that it
     // touches in the global "stiffness matrix", and also keep track of some metadata about
@@ -342,6 +345,8 @@ struct GradientAssemblyLookupTables {
     }
 
     row_ptr.back() = static_cast<int>(nnz);
+
+    initialized = true;
   }
 
   /**
@@ -368,6 +373,9 @@ struct GradientAssemblyLookupTables {
    * corresponding to the (i,j) entry
    */
   std::unordered_map<Entry, uint32_t, Entry::Hasher> nz_LUT;
+
+  /// @brief specifies if the table has already been initialized or not
+  bool initialized;
 };
 
 }  // namespace serac

src/serac/numerics/functional/functional.hpp

@@ -538,7 +538,6 @@ class Functional<test(trials...), exec> {
     Gradient(Functional<test(trials...), exec>& f, uint32_t which = 0)
         : mfem::Operator(f.test_space_->GetTrueVSize(), f.trial_space_[which]->GetTrueVSize()),
           form_(f),
-          lookup_tables(f.G_test_[Domain::Type::Elements], f.G_trial_[Domain::Type::Elements][which]),
           which_argument(which),
           test_space_(f.test_space_),
           trial_space_(f.trial_space_[which]),
@@ -576,6 +575,11 @@ class Functional<test(trials...), exec> {
 
       constexpr bool col_ind_is_sorted = true;
 
+      if (!lookup_tables.initialized) {
+        lookup_tables.init(form_.G_test_[Domain::Type::Elements],
+                           form_.G_trial_[Domain::Type::Elements][which_argument]);
+      }
+
       double* values = new double[lookup_tables.nnz]{};
 
       std::map<mfem::Geometry::Type, ExecArray<double, 3, exec>> element_gradients[Domain::num_types];

src/serac/numerics/petsc_solvers.cpp

@@ -431,15 +431,12 @@ void PetscGAMGSolver::SetupNearNullSpace()
     mfem::VectorFunctionCoefficient coeff_coords(sdim, func_coords);
     mfem::ParGridFunction           gf_coords(fespace_coords);
     gf_coords.ProjectCoefficient(coeff_coords);
-    int                   num_nodes   = fespace_->GetNDofs();
     mfem::HypreParVector* hvec_coords = gf_coords.ParallelProject();
     auto data_coords = const_cast<PetscScalar*>(mfem::Read(hvec_coords->GetMemory(), hvec_coords->Size(), false));
-    PetscCallAbort(GetComm(), PCSetCoordinates(*this, sdim, num_nodes, data_coords));
 
     Vec pvec_coords;
-    PetscCallAbort(GetComm(), VecCreateMPIWithArray(GetComm(), sdim, hvec_coords->Size(), hvec_coords->GlobalSize(),
-                                                    data_coords, &pvec_coords));
-    VecViewFromOptions(pvec_coords, NULL, "-view_coords");
+    PetscCallAbort(GetComm(), VecCreateMPIWithArray(GetComm(), sdim, hvec_coords->Size(),
+                                                    fespace_coords->GlobalTrueVSize(), data_coords, &pvec_coords));
     PetscCallAbort(GetComm(), MatNullSpaceCreateRigidBody(pvec_coords, &nnsp));
     PetscCallAbort(GetComm(), MatSetNearNullSpace(pA, nnsp));
     PetscCallAbort(GetComm(), MatNullSpaceDestroy(&nnsp));