Development supporting p-multigrid

include/base/dof_map.h

@@ -1070,6 +1070,33 @@ class DofMap : public DofMapBase,
    */
   void create_dof_constraints (const MeshBase &, Real time=0);
 
+#ifdef LIBMESH_ENABLE_DIRICHLET
+  /**
+   * Computes the value each \p DirichletBoundary in \p dirichlets prescribes for
+   * every degree of freedom it reaches, by the same local per-entity projection
+   * create_dof_constraints() uses, and stores those values in \p values keyed on
+   * global degree of freedom index.
+   *
+   * This constrains nothing: the DofMap is left exactly as it was found, so a
+   * caller that wants the projected coefficients themselves need not add a
+   * boundary, sweep, and remove it again.  A projected value is a coefficient in
+   * whatever basis is current, which makes it meaningful on a modal basis as well
+   * as an interpolatory one.
+   *
+   * A degree of freedom this DofMap already constrains is omitted, exactly as
+   * create_dof_constraints() leaves such a degree of freedom to the constraint
+   * that already holds it.
+   *
+   * Only local degrees of freedom are computed, as in create_dof_constraints().
+   *
+   * A time is specified for use with time-dependent Dirichlet functions.
+   */
+  void compute_dirichlet_values (const DirichletBoundaries & dirichlets,
+                                 const MeshBase & mesh,
+                                 Real time,
+                                 DofConstraintValueMap & values) const;
+#endif // LIBMESH_ENABLE_DIRICHLET
+
   /**
    * Gathers constraint equation dependencies from other processors
    */

include/base/libmesh_logging.h

@@ -66,27 +66,35 @@ struct PerfItem
   PerfItem(const char * label,
            const char * header,
            bool enabled=true,
-           PerfLog * my_perflog=&perflog) :
+           PerfLog * my_perflog=&perflog,
+           bool nvtx_range=true) :
     _label(label),
     _header(header),
     _enabled(enabled),
-    _perflog(*my_perflog)
+    _perflog(*my_perflog),
+    _nvtx_range(nvtx_range)
   {
     if (_enabled)
-      _perflog.fast_push(label, header);
+      _perflog.fast_push(label, header, _nvtx_range);
   }
 
   ~PerfItem()
   {
     if (_enabled)
-      _perflog.fast_pop(_label, _header);
+      _perflog.fast_pop(_label, _header, _nvtx_range);
   }
 
 private:
   const char * _label;
   const char * _header;
   bool _enabled;
   PerfLog & _perflog;
+
+  /**
+   * Whether to open an NVTX range alongside the performance log event; see
+   * PerfLog::fast_push().
+   */
+  bool _nvtx_range;
 };
 
 
@@ -106,6 +114,10 @@ struct PerfItem
 #  define LOG_SCOPE(a,b)   libMesh::PerfItem TOKENPASTE2(perf_item_, __LINE__)(a,b);
 #  define LOG_SCOPE_IF(a,b,enabled)   libMesh::PerfItem TOKENPASTE2(perf_item_, __LINE__)(a,b,enabled);
 #  define LOG_SCOPE_WITH(a,b,logger)   libMesh::PerfItem TOKENPASTE2(perf_item_, __LINE__)(a,b,true,&logger);
+// Times the event in the performance log while leaving it out of any NVTX
+// capture.  Use for events that fire once per mesh entity or more often; see
+// PerfLog::fast_push().
+#  define LOG_SCOPE_NO_NVTX(a,b)   libMesh::PerfItem TOKENPASTE2(perf_item_, __LINE__)(a,b,true,&libMesh::perflog,false);
 
 #else
 
@@ -116,6 +128,7 @@ struct PerfItem
 #  define LOG_SCOPE(a,b)   {}
 #  define LOG_SCOPE_IF(a,b,enabled) {}
 #  define LOG_SCOPE_WITH(a,b,logger) {}
+#  define LOG_SCOPE_NO_NVTX(a,b) {}
 
 #endif
 

include/fe/fe.h

@@ -26,6 +26,7 @@
 #include "libmesh/libmesh.h"
 
 // C++ includes
+#include <cmath>
 #include <cstddef>
 
 namespace libMesh
@@ -1558,6 +1559,36 @@ OutputShape fe_fdm_deriv(const ElemType type,
                             const Elem *, const unsigned int,
                             const Point &));
 
+/**
+ * The scaling that gives the \p i'th one-dimensional HIERARCHIC bubble function unit \f$H^1\f$
+ * seminorm on the reference interval, for \p i greater than one.
+ *
+ * The bubbles are \f$\xi^i - 1\f$ for even \p i and \f$\xi^i - \xi\f$ for odd \p i, up to this
+ * scaling. Integrating the square of their derivatives over \f$[-1,1]\f$ gives a seminorm of
+ * \f$\sqrt{2/(2i-1)}/(i-1)!\f$ in the even case and \f$(i-1)\sqrt{2/(2i-1)}/i!\f$ in the odd one, so
+ * the reciprocal of the seminorm cancels a factorial and what remains grows only as \f$\sqrt{i}\f$.
+ *
+ * The normalization matters because a shape function carrying \f$1/i!\f$ instead shrinks factorially
+ * with its order, and an operator assembled from such a basis inherits that spread on its diagonal.
+ * At order eight in two dimensions the smallest diagonal entry falls below the roundoff of the
+ * largest, which leaves the discretization numerically singular in double precision however it is
+ * solved.
+ *
+ * Only the bubbles are scaled. The two vertex functions are interpolatory, so their coefficients are
+ * values of the finite element solution at the vertices, and scaling them would change what a nodal
+ * boundary condition or a nodal output of a HIERARCHIC variable means.
+ */
+inline Real fe_hierarchic_bubble_scaling(const unsigned int i)
+{
+  libmesh_assert_greater(i, 1);
+
+  // An even bubble differentiates to xi^(i-1)/(i-1)!, whose square integrates to 2/(2i-1) over the
+  // interval. The linear term an odd bubble carries turns the i^2 of that calculation into (i-1)^2.
+  const Real denominator = (i % 2) ? Real(i) - 1. : Real(i);
+
+  return std::sqrt((2. * Real(i) - 1.) / 2.) / denominator;
+}
+
 
 template <typename OutputShape>
 OutputShape

include/geom/elem.h

@@ -1159,6 +1159,34 @@ class Elem : public ReferenceCountedObject<Elem>,
    */
   bool positive_face_orientation(const unsigned int i) const;
 
+  /**
+   * \returns The orientation of edge \p i, which is 1 when the edge is positively
+   * oriented and 0 when it is not.
+   *
+   * A basis whose shape functions depend on the order of an edge's vertices, such as a
+   * hierarchic basis of order two or above, consumes this index to select the shape
+   * functions of the degrees of freedom that the edge owns.
+   */
+  unsigned int edge_orientation(const unsigned int i) const;
+
+  /**
+   * \returns The orientation of face \p i, an index over the symmetries of the face.
+   *
+   * The orientation of a face is the ordering that the positions of the face's vertices
+   * induce on them: which vertex is the lexicographically least, together with the
+   * direction in which the remaining vertices run from it. This index carries both, as
+   * twice the position that the least vertex takes in the face's node map plus one when
+   * the face is positively oriented, so it runs over the six symmetries of a triangular
+   * face or the eight of a quadrilateral face.
+   *
+   * A basis whose shape functions depend on the order of a face's vertices, such as a
+   * hierarchic basis of order three or above, consumes this index to select the shape
+   * functions of the degrees of freedom that the face owns. The single bit that
+   * \p positive_face_orientation carries fixes the direction alone, which leaves the
+   * least vertex of a face with more than three vertices open.
+   */
+  unsigned int face_orientation(const unsigned int i) const;
+
   /**
    * \returns \p true iff, for an edge \p e on side \p s, the node map for
    * side \p s is such that the first vertex (i.e. zeroth node) of \p e is

include/numerics/petsc_matrix_shell_matrix.h

@@ -58,6 +58,27 @@ class PetscMatrixShellMatrix : public PetscMatrixBase<T>
 
   virtual bool require_sparsity_pattern() const override { return false; }
 
+  virtual void zero() override;
+  virtual std::unique_ptr<SparseMatrix<T>> zero_clone() const override;
+  virtual std::unique_ptr<SparseMatrix<T>> clone() const override;
+  virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override;
+  virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override;
+  virtual void add_matrix(const DenseMatrix<T> & dm,
+                          const std::vector<numeric_index_type> & rows,
+                          const std::vector<numeric_index_type> & cols) override;
+  virtual void add_matrix(const DenseMatrix<T> & dm,
+                          const std::vector<numeric_index_type> & dof_indices) override;
+  virtual void add(const T a, const SparseMatrix<T> & X) override;
+  virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override;
+  virtual Real l1_norm() const override;
+  virtual Real linfty_norm() const override;
+  virtual void print_personal(std::ostream & os = libMesh::out) const override;
+  virtual void get_diagonal(NumericVector<T> & dest) const override;
+  virtual void get_transpose(SparseMatrix<T> & dest) const override;
+  virtual void get_row(numeric_index_type i,
+                       std::vector<numeric_index_type> & indices,
+                       std::vector<T> & values) const override;
+
 private:
   // Make this private because we mark as initialized after we've done our initialization, and we
   // don't want derived classes to mistakenly register their data as initialized (or not)

include/numerics/petsc_mffd_matrix.h

@@ -48,7 +48,15 @@ class PetscMFFDMatrix : public PetscMatrixBase<T>
 
   explicit PetscMFFDMatrix(const Parallel::Communicator & comm_in);
 
-  PetscMFFDMatrix & operator=(Mat m);
+  /**
+   * Adopt an existing, externally-owned Mat, without destroying it when this
+   * object goes out of scope. \p set_context controls whether we attach a
+   * context pointer to \p m allowing \p get_context() to recover this object
+   * from the Mat later; skip this when this wrapper is short-lived (e.g. a
+   * function-local variable) so we don't leave a dangling context on \p m
+   * after we're destroyed.
+   */
+  void assign(Mat m, bool set_context);
 
   virtual void init(const numeric_index_type,
                     const numeric_index_type,
@@ -101,11 +109,14 @@ PetscMFFDMatrix<T>::PetscMFFDMatrix(const Parallel::Communicator & comm_in)
 }
 
 template <typename T>
-PetscMFFDMatrix<T> &
-PetscMFFDMatrix<T>::operator=(Mat m)
+void
+PetscMFFDMatrix<T>::assign(Mat m, bool set_context)
 {
   this->_mat = m;
-  return *this;
+  this->_is_initialized = true;
+  this->_destroy_mat_on_exit = false;
+  if (set_context)
+    this->set_context();
 }
 
 template <typename T>

include/solvers/nonlinear_solver.h

@@ -102,6 +102,21 @@ class NonlinearSolver : public ReferenceCountedObject<NonlinearSolver<T>>,
                                                const double,      // Stopping tolerance
                                                const unsigned int) = 0; // N. Iterations
 
+  /**
+   * Solves the nonlinear system using \p jac_in as the actual Jacobian operator and \p pre_in as
+   * the preconditioning matrix -- which may be the same object (the common case) or genuinely
+   * distinct (e.g. a matrix-free operator paired with an assembled preconditioning matrix).
+   */
+  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> & /* jac_in */,
+                                               SparseMatrix<T> & /* pre_in */,
+                                               NumericVector<T> & /* x_in */,
+                                               NumericVector<T> & /* r_in */,
+                                               const double /* tol */,
+                                               const unsigned int /* m_its */)
+  {
+    libmesh_not_implemented();
+  }
+
   /**
    * Prints a useful message about why the latest nonlinear solve
    * con(di)verged.

include/solvers/petsc_nonlinear_solver.h

@@ -107,11 +107,23 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
   SNES snes(const char * name = nullptr);
 
   /**
-   * Call the Petsc solver.  It calls the method below, using the
-   * same matrix for the system and preconditioner matrices.
+   * Call the Petsc solver, using the same matrix for the system and preconditioner matrices. Calls
+   * the two-matrix overload below with \p pre_in for both.
    */
   virtual std::pair<unsigned int, Real>
-  solve (SparseMatrix<T> &,                     // System Jacobian Matrix
+  solve (SparseMatrix<T> & pre_in,               // System Preconditioning Matrix
+         NumericVector<T> &,                    // Solution vector
+         NumericVector<T> &,                    // Residual vector
+         const double,                         // Stopping tolerance
+         const unsigned int) override; // N. Iterations
+
+  /**
+   * Call the Petsc solver, using \p jac_in as the actual SNES Jacobian operator (Amat) and
+   * \p pre_in as the preconditioning matrix (Pmat).
+   */
+  virtual std::pair<unsigned int, Real>
+  solve (SparseMatrix<T> & jac_in,               // Jacobian operator matrix (Amat)
+         SparseMatrix<T> & pre_in,               // Preconditioning matrix (Pmat)
          NumericVector<T> &,                    // Solution vector
          NumericVector<T> &,                    // Residual vector
          const double,                         // Stopping tolerance
@@ -171,6 +183,27 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    */
   bool snes_mf_reuse_base() const { return _snesmf_reuse_base; }
 
+  /**
+   * Set whether evaluating the residual twice at the same solution yields bitwise identical
+   * results. The debug-mode check that reusing the nonlinear residual as the matrix-free base
+   * vector was safe compares the recomputed base residual against the stored one bitwise, so it
+   * says something only when assembly sums each degree of freedom's contributions in a fixed
+   * order. MPI-parallel and threaded assembly do not, and are already excluded from that check;
+   * an application that performs its own parallel assembly, for instance accumulating with
+   * device atomics whose completion order is not reproducible, has to report that itself.
+   *
+   * Setting this to false disables only that check. It does not change what is computed, and in
+   * particular does not stop the base vector from being reused; where the difference between two
+   * evaluations is at the level of round-off, reusing it remains sound.
+   */
+  void set_reproducible_residual(bool state) { _reproducible_residual = state; }
+
+  /**
+   * @return Whether evaluating the residual twice at the same solution yields bitwise identical
+   * results
+   */
+  bool reproducible_residual() const { return _reproducible_residual; }
+
   /**
    * Set whether we are computing the base vector for matrix-free finite-differencing
    */
@@ -268,6 +301,12 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
    */
   bool _snesmf_reuse_base;
 
+  /**
+   * True if two residual evaluations at the same solution agree bitwise, which is what makes the
+   * debug-mode base vector reuse check meaningful. See set_reproducible_residual()
+   */
+  bool _reproducible_residual;
+
   void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
                             void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
                             MatNullSpace *);
@@ -291,9 +330,6 @@ class PetscNonlinearSolver : public NonlinearSolver<T>
   PetscDMWrapper _dm_wrapper;
 #endif
 
-  /// Wrapper for matrix-free finite-difference Jacobians
-  PetscMFFDMatrix<Number> _mffd_jac;
-
 private:
   friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
   friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);

include/systems/nonlinear_implicit_system.h

@@ -257,10 +257,20 @@ class NonlinearImplicitSystem : public ImplicitSystem
   virtual void reinit () override;
 
   /**
-   * Assembles & solves the nonlinear system R(x) = 0.
+   * Assembles & solves the nonlinear system R(x) = 0. If a matrix has been registered via
+   * set_operator_matrix(), it is used as the actual Jacobian operator (Amat) while \p matrix
+   * remains the preconditioning matrix (Pmat); otherwise \p matrix is used for both, as usual.
    */
   virtual void solve () override;
 
+  /**
+   * Set a matrix to use as the actual Jacobian operator (Amat) on the next solve(), distinct from
+   * \p matrix, which continues to be used as the preconditioning matrix (Pmat). Pass nullptr (the
+   * default) to restore the ordinary behavior of using \p matrix for both. Only solve() itself
+   * reads this value back, so there is no public getter.
+   */
+  void set_operator_matrix(SparseMatrix<Number> * mat) { _operator_matrix = mat; }
+
   /**
    * \returns An integer corresponding to the upper iteration count
    * limit and a Real corresponding to the convergence tolerance to
@@ -335,6 +345,12 @@ class NonlinearImplicitSystem : public ImplicitSystem
    * The final residual for the nonlinear system R(x)
    */
   Real _final_nonlinear_residual;
+
+  /**
+   * An optional matrix to use as the actual Jacobian operator (Amat), distinct from \p matrix
+   * (used as the preconditioning matrix, Pmat). See set_operator_matrix().
+   */
+  SparseMatrix<Number> * _operator_matrix;
 };
 
 } // namespace libMesh