Tutorials for SNES and TS

doc/doxygen/tutorial/tutorial.h.in

@@ -684,6 +684,12 @@
  *       <br/> Keywords: FEInterfaceValues, NonMatching::FEImmersedSurfaceValues
  *       </td></tr>
  *
+ *   <tr valign="top">
+ *       <td>step-86</td>
+ *       <td> The nonlinear minimal surface equation revisited.
+ *       Interfacing with PETSc' SNES nonlinear solver.
+ *       </td></tr>
+ *
  * </table>
  *
  *
@@ -804,7 +810,8 @@
  *       step-18,
  *       step-40,
  *       step-50,
- *       step-55
+ *       step-55,
+ *       step-86
  *     </td>
  *   </tr>
  *
@@ -917,7 +924,9 @@
  *       step-42,
  *       step-43,
  *       step-57,
- *       step-70
+ *       step-70,
+ *       step-77,
+ *       step-86
  *     </td>
  *   </tr>
  *
@@ -1081,7 +1090,8 @@
  *       step-58,
  *       step-62,
  *       step-81,
- *       step-82
+ *       step-82,
+ *       step-86
  *     </td>
  *   </tr>
  *
@@ -1121,7 +1131,8 @@
  *       step-55,
  *       step-59,
  *       step-66,
- *       step-75
+ *       step-75,
+ *       step-86
  *     </td>
  *   </tr>
  *
@@ -1157,7 +1168,9 @@
  *       step-41,
  *       step-42,
  *       step-44,
- *       step-57
+ *       step-57,
+ *       step-77,
+ *       step-86
  *     </td>
  *   </tr>
  *
@@ -1217,7 +1230,9 @@
  *     <td> Minimal surface equation
  *     </td>
  *     <td>
- *       step-15
+ *       step-15,
+ *       step-77,
+ *       step-86
  *     </td>
  *   </tr>
  *

examples/step-67/doc/results.dox

@@ -3,7 +3,7 @@
 <h3>Program output</h3>
 
 Running the program with the default settings on a machine with 40 processes
-produces the following output:
+produces the following output if PETSc is not installed:
 @code
 Running with 40 MPI processes
 Vectorization over 8 doubles = 512 bits (AVX512)
@@ -703,7 +703,8 @@ efficiency in terms of step size per stage before they violate the CFL
 condition. An interesting extension would be to compare the low-storage
 variants proposed here with standard Runge--Kutta integrators or to use vector
 operations that are run separate from the mass matrix operation and compare
-performance.
+performance. This could be done using the PETSc TimeStepper solver, as already
+implemented in this tutorial.
 
 
 <h4>More advanced numerical flux functions and skew-symmetric formulations</h4>

examples/step-67/doc/tooltip

@@ -1 +1,2 @@
 An explicit time integrator for the Euler equations with matrix-free implementation.
+The tutorial takes advantage of PETSc's ODE integrators if available.

examples/step-67/step-67.cc

@@ -15,10 +15,13 @@
 
  *
  * Author: Martin Kronbichler, 2020
+ 2022
  */
 
 // The include files are similar to the previous matrix-free tutorial programs
-// step-37, step-48, and step-59
+// step-37, step-48, and step-59. This tutorial will use the PETSc ODE solver
+// interface if available. We thus include the relevant headers.
+
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/logstream.h>
@@ -45,6 +48,11 @@
 
 #include <deal.II/numerics/data_out.h>
 
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/lac/petsc_ts.h>
+#  include <deal.II/lac/petsc_vector.h>
+#endif
+
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -372,7 +380,6 @@ namespace Euler_DG
   };
 
 
-
   // @sect3{Implementation of point-wise operations of the Euler equations}
 
   // In the following functions, we implement the various problem-specific
@@ -2248,6 +2255,82 @@ namespace Euler_DG
     // mostly done by the TimerOutput::print_wall_time_statistics() function.
     unsigned int timestep_number = 0;
 
+#ifdef DEAL_II_WITH_PETSC
+    // Here we detail how to use PETSc TS explicit solvers of Runge Kutta type.
+    // We use two small helper functions to copy to and from
+    // instances of `LinearAlgebra::distributed::Vector`, since the
+    // current implementation does not directly support this vector class.
+    using VectorType = PETScWrappers::MPI::Vector;
+
+    auto copy_from_petsc = [](const VectorType &                          src,
+                              LinearAlgebra::distributed::Vector<Number> &dst) {
+      auto lr = src.local_range();
+      for (unsigned int i = lr.first; i < lr.second; i++)
+        dst[i] = src[i];
+    };
+
+    auto copy_to_petsc =
+      [](const LinearAlgebra::distributed::Vector<Number> &src,
+         VectorType &                                      dst) {
+        auto lr = dst.local_range();
+        for (unsigned int i = lr.first; i < lr.second; i++)
+          dst[i] = src[i];
+        dst.compress(VectorOperation::insert);
+      };
+
+    // We construct our solver with parameters like initial time step and final
+    // time, together with the solver type for Runge-Kutta. We also decide to
+    // use adaptive time stepping based on <a
+    // href="https://dl.acm.org/doi/10.1145/641876.641877"> digital signal
+    // processing</a>, and set adaptive tolerances to -1 to use default
+    // tolerances from PETSc
+    PETScWrappers::TimeStepperData tsdata;
+    tsdata.final_time         = final_time;
+    tsdata.initial_step_size  = time_step;
+    tsdata.ts_type            = "rk";
+    tsdata.ts_adapt_type      = "dsp";
+    tsdata.absolute_tolerance = -1.0;
+    tsdata.relative_tolerance = -1.0;
+
+    PETScWrappers::TimeStepper<VectorType> petsc_time_stepper(tsdata);
+
+    // PETScWrappers::TimeStepper can solve ODE/DAE problems in explicit form
+    // i.e., $\dot{u} = G(u,t)$, or in fully implicit form $F(\dot{u},u,t) = 0$.
+    // This function implements the explicit interface for $G$.
+    petsc_time_stepper.explicit_function =
+      [&](const double t, const VectorType &y, VectorType &f) {
+        time = t;
+        copy_from_petsc(y, rk_register_1);
+        euler_operator.apply(t, rk_register_1, rk_register_2);
+        copy_to_petsc(rk_register_2, f);
+        return 0;
+      };
+
+    // We can also attach a monitoring callback to inspect the solution.
+    petsc_time_stepper.monitor = [&](const double      t,
+                                     const VectorType &y,
+                                     const unsigned int) {
+      time      = t;
+      time_step = petsc_time_stepper.get_time_step();
+      copy_from_petsc(y, solution);
+
+      static int prev_tick = 0;
+      int        curr_tick = static_cast<int>(t / output_tick);
+      if (curr_tick != prev_tick)
+        output_results(static_cast<unsigned int>(std::round(t / output_tick)));
+      prev_tick = curr_tick;
+      return 0;
+    };
+
+    // We finally set the initial condition and integrate the ODE.
+    VectorType psolution(solution.locally_owned_elements(),
+                         solution.get_mpi_communicator());
+    copy_to_petsc(solution, psolution);
+    petsc_time_stepper.solve(psolution);
+    copy_from_petsc(psolution, solution);
+
+#else
+
     while (time < final_time - 1e-12)
       {
         ++timestep_number;
@@ -2275,7 +2358,7 @@ namespace Euler_DG
           output_results(
             static_cast<unsigned int>(std::round(time / output_tick)));
       }
-
+#endif
     timer.print_wall_time_statistics(MPI_COMM_WORLD);
     pcout << std::endl;
   }

examples/step-86/CMakeLists.txt

@@ -0,0 +1,49 @@
+##
+#  CMake script for the step-86 tutorial program:
+##
+
+# Set the name of the project and target:
+set(TARGET "step-86")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  file(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  file(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  set(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+set(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+cmake_minimum_required(VERSION 3.3.0)
+
+find_package(deal.II 9.5.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+if(NOT ${deal.II_FOUND})
+  message(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+endif()
+
+if(NOT DEAL_II_WITH_PETSC) # keep in one line
+  message(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_PETSC = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+    DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
+This conflicts with the requirements."
+    )
+endif()
+
+deal_ii_initialize_cached_variables()
+project(${TARGET})
+deal_ii_invoke_autopilot()

examples/step-86/doc/builds-on

@@ -0,0 +1,2 @@
+step-15
+step-77

examples/step-86/doc/intro.dox

@@ -0,0 +1,89 @@
+<br>
+
+<i>
+This program was contributed by Stefano Zampini, King Abdullah University of Science and Technology.
+
+</i>
+<br>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+The step-15 program solved the following, nonlinear equation
+describing the minimal surface problem:
+@f{align*}{
+    -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad
+    \qquad &&\textrm{in} ~ \Omega
+    \\
+    u&=g \qquad\qquad &&\textrm{on} ~ \partial \Omega.
+@f}
+step-15 uses a Newton method, and
+Newton's method works by repeatedly solving a *linearized* problem for
+an update $\delta u_k$ -- called the "search direction" --, computing a
+"step length"
+$\alpha_k$, and then combining them to compute the new
+guess for the solution via
+@f{align*}{
+    u_{k+1} = u_k + \alpha_k \, \delta u_k.
+@f}
+
+Here, we use PETSc to solve the nonlinear system of equations with the SNES
+package, with a very similar interface as discussed for KINSOL in step-77.
+
+<h3> How deal.II interfaces with SNES </h3>
+
+SNES, like many similar packages, works in a pretty abstract way. At
+its core, it sees a nonlinear problem of the form
+@f{align*}{
+    F(U) = 0
+@f}
+and constructs a sequence of iterates $U_k$ which, in general, are
+vectors of the same length as the vector returned by the function
+$F$. To do this, there are a few things it needs from the user:
+- A way to evaluate, for a given vector $U$, the function $F(U)$. This
+  function is generally called the "residual" operation because the
+  goal is of course to find a point $U^\ast$ for which $F(U^\ast)=0$;
+  if $F(U)$ returns a nonzero vector, then this is the
+  <a href="https://en.wikipedia.org/wiki/Residual_(numerical_analysis)">"residual"</a>
+  (i.e., the "rest", or whatever is "left over"). The function
+  that will do this is in essence the same as the computation of
+  the right hand side vector in step-15, but with an important difference:
+  there, the right hand side denoted the *negative* of the residual,
+  so we have to switch a sign.
+- A way to compute the matrix $J_k$ if that is necessary in the
+  current iteration. This operation will generally be called the
+  "jacobian" operation.
+
+The deal.II interface also offers an alternative approach, more
+in line with the deal.II philosophy:
+- A way to compute the matrix $J_k$ if that is necessary in the
+  current iteration, along with possibly a preconditioner or other
+  data structures. This operation will generally be called the
+  "setup_jacobian" operation.
+- A way to solve a linear system $\tilde J_k x = b$ with whatever
+  matrix $\tilde J_k$ was last computed. This operation will generally
+  be called the "solve_with_jacobian" operation.
+
+These operations must be provided in the form of
+[std::function](https://en.cppreference.com/w/cpp/utility/functional/function)
+objects that take the appropriate set of arguments and that generally
+return an integer that indicates success (a zero return value) or
+failure (a nonzero return value). Specifically, the object we will
+access is PETScWrappers::NonlinearSolver, and its methods:
+PETScWrappers::NonlinearSolver::residual,
+PETScWrappers::NonlinearSolver::setup_jacobian, and
+PETScWrappers::NonlinearSolver::solve_with_jacobian.
+In our implementation, we will use
+[lambda functions](https://en.cppreference.com/w/cpp/language/lambda)
+to implement these "callbacks" that in turn can call member functions;
+SNES will then call these callbacks whenever its internal algorithms
+think it is useful.
+
+
+<h3> Details of the implementation </h3>
+
+The majority of the code of this tutorial program is as in step-15 and
+in step-77, and we will not comment on it, except for one thing.
+We don't need to explicitly factorize the Jacobian matrix ourselves as
+it is done in step-77, since PETSc will automatically do it for
+us when needed.

examples/step-86/doc/kind

@@ -0,0 +1 @@
+techniques