Options for restricting solves to unconstrained DoFs only

examples/adjoints/adjoints_ex3/adjoints_ex3.C

@@ -790,6 +790,17 @@ int main (int argc, char ** argv)
 
   equation_systems.init ();
 
+  if (infile.search("--restrict_solves_to_unconstrained"))
+    {
+      libMesh::out << "Restricting solves to unconstrained DoFs" << std::endl;
+      system.restrict_solves_to_unconstrained(true);
+      libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(),
+                              &system.get_dof_map());
+      libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained());
+    }
+  else
+    libMesh::out << "Not restricting solves to unconstrained DoFs" << std::endl;
+
   // Print information about the mesh and system to the screen.
   mesh.print_info();
   equation_systems.print_info();
@@ -830,9 +841,17 @@ int main (int argc, char ** argv)
                    << std::endl
                    << std::endl;
 
+      libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(),
+                              &system.get_dof_map());
+      libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained());
+
       // Solve the forward system
       system.solve();
 
+      libmesh_assert_equal_to(system.get_linear_solver()->restrictions_to_unconstrained(),
+                              &system.get_dof_map());
+      libmesh_assert(system.get_time_solver().diff_solver()->get_restrict_solves_to_unconstrained());
+
       // Write the output
       write_output(equation_systems, 0, a_step, "primal", param);
 

examples/adjoints/adjoints_ex3/run.sh

@@ -8,7 +8,7 @@ example_name=adjoints_ex3
 
 example_dir=examples/adjoints/$example_name
 
-run_example "$example_name" $options
+run_example "$example_name" --restrict_solves_to_unconstrained $options
 
 # This example needs better solvers on larger meshes + processor counts
 #benchmark_example 1 "$example_name" "max_adaptivesteps=18 -pc_type bjacobi -sub_pc_factor_levels 4"

include/base/dof_map.h

@@ -1251,6 +1251,20 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   void constrain_nothing (std::vector<dof_id_type> & dofs) const;
 
+  /**
+   * Constrains the numeric vector \p v, which represents a solution defined on
+   * the mesh.  This may need to be used after a linear solve, if your linear
+   * solver's solutions do not satisfy your DoF constraints to a tight enough
+   * tolerance.
+   *
+   * If \p homogeneous == true, heterogeneous constraints are enforced
+   * as if they were homogeneous.  This might be appropriate for e.g. a
+   * vector representing a difference between two
+   * heterogeneously-constrained solutions.
+   */
+  void enforce_constraints_exactly (NumericVector<Number> & v,
+                                    bool homogeneous = false) const;
+
   /**
    * Constrains the numeric vector \p v, which represents a solution defined on
    * the mesh.  This may need to be used after a linear solve, if your linear
@@ -1263,6 +1277,8 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    * as if they were homogeneous.  This might be appropriate for e.g. a
    * vector representing a difference between two
    * heterogeneously-constrained solutions.
+   *
+   * Will be deprecated eventually; use the non-System& overload.
    */
   void enforce_constraints_exactly (const System & system,
                                     NumericVector<Number> * v = nullptr,

include/solvers/diff_solver.h

@@ -139,6 +139,15 @@ class DiffSolver : public ReferenceCountedObject<DiffSolver>,
    */
   sys_type & system () { return _system; }
 
+  /**
+   * After calling this method with \p true, all successive solves
+   * will be restricted to unconstrained dofs, with constraints
+   * applied separately as necessary.
+   */
+  virtual void restrict_solves_to_unconstrained(bool restricting) = 0;
+
+  virtual bool get_restrict_solves_to_unconstrained() = 0;
+
   /**
    * Each linear solver step should exit after \p max_linear_iterations
    * is exceeded.

include/solvers/linear_solver.h

@@ -40,6 +40,7 @@ template <typename T> class SparseMatrix;
 template <typename T> class NumericVector;
 template <typename T> class ShellMatrix;
 template <typename T> class Preconditioner;
+class DofMap;
 class System;
 class SolverConfiguration;
 enum SolverPackage : int;
@@ -152,6 +153,26 @@ class LinearSolver : public ReferenceCountedObject<LinearSolver<T>>,
   virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
                                   const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
 
+  /**
+   * After calling this method, all successive solves will be
+   * restricted to unconstrained dofs.  This mode can be disabled by
+   * calling this method with \p dof_map being a \p nullptr.
+   *
+   * This does a pure submatrix \ subvector solve, with no correction
+   * of the rhs to account for ignored columns (so the application of
+   * constraints to the matrix should be done symmetrically
+   * beforehand) and no post-solution application of ignored rows (so
+   * homogeneous or heterogeneous primal or heterogeneous adjoint
+   * constraints should be applied manually afterward)
+   */
+  virtual void restrict_solve_to_unconstrained (const DofMap * dof_map);
+
+  /**
+   * Returns the current restriction (or lack thereof) state with
+   * regard to constrained/unconstrained dofs.
+   */
+  virtual const DofMap * restrictions_to_unconstrained();
+
   /**
    * This function calls the solver \p _solver_type preconditioned
    * with the \p _preconditioner_type preconditioner.

include/solvers/newton_solver.h

@@ -88,6 +88,15 @@ class NewtonSolver : public DiffSolver
     return *_linear_solver;
   }
 
+  /**
+   * After calling this method with \p true, all successive solves
+   * will be restricted to unconstrained dofs, with constraints
+   * applied separately as necessary.
+   */
+  virtual void restrict_solves_to_unconstrained(bool restricting) override;
+
+  virtual bool get_restrict_solves_to_unconstrained() override;
+
   /**
    * If this is set to true, the solver is forced to test the residual
    * after each Newton step, and to reduce the length of its steps

include/solvers/petsc_diff_solver.h

@@ -86,13 +86,40 @@ class PetscDiffSolver : public DiffSolver
    */
   void clear ();
 
+  /**
+   * After calling this method with \p true, all successive solves
+   * will be restricted to unconstrained dofs, with constraints
+   * applied separately as necessary.
+   */
+  virtual void restrict_solves_to_unconstrained(bool restricting) override;
+
+  virtual bool get_restrict_solves_to_unconstrained() override;
+
   /**
    * This method performs a solve.  What occurs in
    * this method will depend on the PETSc SNES settings.
    * See the PETSc documentation for more details.
    */
   virtual unsigned int solve () override;
 
+  /*
+   * Scatter connecting full vectors with unconstrained subvectors
+   */
+  WrappedPetsc<VecScatter> scatter;
+
+  /*
+   * Submatrix of Jacobian with unconstrained-unconstrained coupling
+   */
+  WrappedPetsc<Mat> submat;
+
+  bool submat_created;
+
+  /**
+   * PETSc index set containing unconstrained dofs on which to solve
+   * (\p nullptr means solve on all dofs).
+   */
+  WrappedPetsc<IS> _unconstrained_dofs_is;
+
 protected:
 
   /**
@@ -109,6 +136,11 @@ class PetscDiffSolver : public DiffSolver
 #endif
 #endif
 
+  /**
+   * Are we restricting solves to unconstrained DoFs?
+   */
+  bool _restrict_to_unconstrained;
+
 private:
 
   /**

include/solvers/petsc_linear_solver.h

@@ -147,6 +147,26 @@ class PetscLinearSolver : public LinearSolver<T>
   virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
                                   const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override;
 
+  /**
+   * After calling this method, all successive solves will be
+   * restricted to unconstrained dofs.  This mode can be disabled by
+   * calling this method with \p dof_map being a \p nullptr.
+   *
+   * This does a pure submatrix \ subvector solve, with no correction
+   * of the rhs to account for ignored columns (so the application of
+   * constraints to the matrix should be done symmetrically
+   * beforehand) and no post-solution application of ignored rows (so
+   * homogeneous or heterogeneous primal or heterogeneous adjoint
+   * constraints should be applied manually afterward)
+   */
+  virtual void restrict_solve_to_unconstrained (const DofMap * dof_map) override;
+
+  /**
+   * Returns the current restriction (or lack thereof) state with
+   * regard to constrained/unconstrained dofs.
+   */
+  virtual const DofMap * restrictions_to_unconstrained() override;
+
   /**
    * Call the Petsc solver.  It calls the method below, using the
    * same matrix for the system and preconditioner matrices.
@@ -344,9 +364,21 @@ class PetscLinearSolver : public LinearSolver<T>
   WrappedPetsc<IS> _restrict_solve_to_is_complement;
 
   /**
-   * \returns The local size of \p _restrict_solve_to_is.
+   * DoFMap to use if we're restricting solves to unconstrained DoFs.
+   * (\p nullptr means solve on constrained DoFs too)
+   */
+  const DofMap * _restrict_to_unconstrained_dofmap;
+
+  /**
+   * PETSc index set containing unconstrained dofs on which to solve
+   * (\p nullptr means solve on all dofs).
+   */
+  WrappedPetsc<IS> _unconstrained_dofs_is;
+
+  /**
+   * \returns The local size of \p is.
    */
-  PetscInt restrict_solve_to_is_local_size() const;
+  PetscInt local_is_size(IS & is) const;
 
   /**
    * Creates \p _restrict_solve_to_is_complement to contain all

include/systems/diff_system.h

@@ -304,6 +304,13 @@ class DifferentiableSystem : public ImplicitSystem,
    */
   virtual void side_postprocess (DiffContext &) {}
 
+  /**
+   * After calling this method with \p true, all successive solves
+   * will be restricted to unconstrained dofs, with constraints
+   * applied separately as necessary.
+   */
+  virtual void restrict_solves_to_unconstrained(bool restricting) override;
+
   /**
    * For a given second order (in time) variable var, this method will return
    * the index to the corresponding "dot" variable. For FirstOrderUnsteadySolver

include/systems/implicit_system.h

@@ -308,6 +308,13 @@ class ImplicitSystem : public ExplicitSystem
    */
   SparseMatrix<Number> & get_system_matrix();
 
+  /**
+   * After calling this method with \p true, all successive solves
+   * will be restricted to unconstrained dofs, with constraints
+   * applied separately as necessary.
+   */
+  virtual void restrict_solves_to_unconstrained(bool restricting);
+
   /**
    * The system matrix.  Implicit systems are characterized by
    * the need to solve the linear system Ax=b.  This is the

src/base/dof_map_constraints.C

@@ -2895,8 +2895,7 @@ void DofMap::constrain_nothing (std::vector<dof_id_type> & dofs) const
 
 
 
-void DofMap::enforce_constraints_exactly (const System & system,
-                                          NumericVector<Number> * v,
+void DofMap::enforce_constraints_exactly (NumericVector<Number> & v,
                                           bool homogeneous) const
 {
   parallel_object_only();
@@ -2906,53 +2905,49 @@ void DofMap::enforce_constraints_exactly (const System & system,
 
   LOG_SCOPE("enforce_constraints_exactly()","DofMap");
 
-  if (!v)
-    v = system.solution.get();
-
   NumericVector<Number> * v_local  = nullptr; // will be initialized below
   NumericVector<Number> * v_global = nullptr; // will be initialized below
   std::unique_ptr<NumericVector<Number>> v_built;
-  if (v->type() == SERIAL)
+  if (v.type() == SERIAL)
     {
       v_built = NumericVector<Number>::build(this->comm());
       v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
       v_built->close();
 
       for (dof_id_type i=v_built->first_local_index();
            i<v_built->last_local_index(); i++)
-        v_built->set(i, (*v)(i));
+        v_built->set(i, v(i));
       v_built->close();
       v_global = v_built.get();
 
-      v_local = v;
+      v_local = &v;
       libmesh_assert (v_local->closed());
     }
-  else if (v->type() == PARALLEL)
+  else if (v.type() == PARALLEL)
     {
       v_built = NumericVector<Number>::build(this->comm());
-      v_built->init (v->size(), v->local_size(),
+      v_built->init (v.size(), v.local_size(),
                      this->get_send_list(), true,
                      GHOSTED);
-      v->localize(*v_built, this->get_send_list());
+      v.localize(*v_built, this->get_send_list());
       v_built->close();
       v_local = v_built.get();
 
-      v_global = v;
+      v_global = &v;
     }
-  else if (v->type() == GHOSTED)
+  else if (v.type() == GHOSTED)
     {
-      v_local = v;
-      v_global = v;
+      v_local = &v;
+      v_global = &v;
     }
-  else // unknown v->type()
-    libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v->type()));
+  else // unknown v.type()
+    libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v.type()));
 
   // We should never hit these asserts because we should error-out in
   // else clause above.  Just to be sure we don't try to use v_local
   // and v_global uninitialized...
   libmesh_assert(v_local);
   libmesh_assert(v_global);
-  libmesh_assert_equal_to (this, &(system.get_dof_map()));
 
   for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
     {
@@ -2975,14 +2970,26 @@ void DofMap::enforce_constraints_exactly (const System & system,
 
   // If the old vector was serial, we probably need to send our values
   // to other processors
-  if (v->type() == SERIAL)
+  if (v.type() == SERIAL)
     {
 #ifndef NDEBUG
       v_global->close();
 #endif
-      v_global->localize (*v);
+      v_global->localize (v);
     }
-  v->close();
+  v.close();
+}
+
+void DofMap::enforce_constraints_exactly (const System & system,
+                                          NumericVector<Number> * v,
+                                          bool homogeneous) const
+{
+  libmesh_assert_equal_to (this, &(system.get_dof_map()));
+
+  if (!v)
+    v = system.solution.get();
+
+  this->enforce_constraints_exactly(*v, homogeneous);
 }
 
 void DofMap::enforce_constraints_on_residual (const NonlinearImplicitSystem & system,