]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a test that shows that PETSc TS cannot recover from errors.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 25 May 2023 20:06:35 +0000 (14:06 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 25 May 2023 21:06:23 +0000 (15:06 -0600)
tests/petsc/petsc_ts_04.cc [new file with mode: 0644]
tests/petsc/petsc_ts_04.output [new file with mode: 0644]

diff --git a/tests/petsc/petsc_ts_04.cc b/tests/petsc/petsc_ts_04.cc
new file mode 100644 (file)
index 0000000..0d1d2ad
--- /dev/null
@@ -0,0 +1,275 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal.II library is free software; you can use it, redistribute
+//    it, and/or modify it under the terms of the GNU Lesser General
+//    Public License as published by the Free Software Foundation; either
+//    version 2.1 of the License, or (at your option) any later version.
+//    The full text of the license can be found in the file LICENSE.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_matrix_base.h>
+#include <deal.II/lac/petsc_ts.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+/**
+ * Like the _03 test, but throw an error whenever we find that the
+ * time step for evaluating the residual is too large.
+ *
+ * Solves the equations of exponential decay:
+ *
+ * u' = -k u
+ * u (t0) = 1
+ *
+ * using the PETScWrappers::TimeStepper class
+ * that interfaces PETSc TS ODE solver object.
+ * The ODE can be an EXPLICIT first order ode:
+ *
+ * y[0]' = - k   y[0]     -> y' = G(y,t)
+ *
+ * or specified in IMPLICIT form:
+ *
+ * y[0]' +  k  y[1]  = 0  -> F(y',y,t) = 0
+ *
+ * The exact solution is, in either case,
+ *
+ * y[0](t) = exp(-kt)
+ *
+ * We use the same class to test both formulations.
+ * The class also supports a third approach where
+ * control for linear systems generation and
+ * solution is completely on users. In this case
+ * users are in charge of solving for the
+ * implicit Jacobian.
+ * The model used to wrap PETSc's TS is the same
+ * used for the nonlinear solver SNES. Check
+ * petsc_snes_00.cc for additional information.
+ *
+ */
+using VectorType  = PETScWrappers::MPI::Vector;
+using MatrixType  = PETScWrappers::MatrixBase;
+using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
+using real_type   = TimeStepper::real_type;
+
+class HarmonicOscillator
+{
+public:
+  HarmonicOscillator(real_type                                      _kappa,
+                     const typename PETScWrappers::TimeStepperData &data,
+                     bool                                           setjac,
+                     bool                                           implicit,
+                     bool                                           user)
+    : time_stepper(data)
+    , kappa(_kappa)
+    , last_eval_time(0)
+  {
+    // In this case we use the implicit form
+    if (implicit)
+      {
+        time_stepper.implicit_function = [&](const real_type   t,
+                                             const VectorType &y,
+                                             const VectorType &y_dot,
+                                             VectorType &      res) -> void {
+          deallog << "Evaluating implicit function at t=" << t << std::endl;
+
+          if (t > last_eval_time + 0.1)
+            {
+              deallog << "Time step too large: last_eval_time="
+                      << last_eval_time << ", t=" << t << std::endl;
+              throw RecoverableUserCallbackError();
+            }
+
+
+          res(0) = y_dot(0) + kappa * y(0);
+          res.compress(VectorOperation::insert);
+
+          last_eval_time = t;
+        };
+
+        // We either have the possibility of using PETSc standard
+        // functionalities for Jacobian setup and solve, or take full control of
+        // the Jacobian solutions. The latter aligns more with deal.II tutorials
+        // For our testing class, 'user' discriminate between those two
+        // approaches.
+        if (!user && setjac)
+          {
+            // Callback for Jacobian evaluation
+            // This callback populates the P matrix with
+            //     shift * dF/dydot + dF/dy
+            // A and P are the same matrix if not otherwise specified
+            // PETSc can compute a Jacobian by finite-differencing the residual
+            // evaluation, and it can be selected at command line.
+            time_stepper.implicit_jacobian = [&](const real_type   t,
+                                                 const VectorType &y,
+                                                 const VectorType &y_dot,
+                                                 const real_type   shift,
+                                                 MatrixType &      A,
+                                                 MatrixType &      P) -> void {
+              deallog << "Evaluating implicit Jacobian at t=" << t << std::endl;
+              P.set(0, 0, shift + kappa);
+              P.compress(VectorOperation::insert);
+            };
+          }
+        else if (user)
+          {
+            // In this code block, users are completely in charge of
+            // setting up the Jacobian system and solve for it.
+            // In this example we only store the solver shift
+            // during setup.
+            time_stepper.setup_jacobian = [&](const real_type   t,
+                                              const VectorType &y,
+                                              const VectorType &y_dot,
+                                              const real_type   shift) -> void {
+              deallog << "Setting up Jacobian at t=" << t << std::endl;
+              myshift = shift;
+            };
+
+            // In the solve phase we se the stored shift to solve
+            // for the implicit Jacobian system
+            time_stepper.solve_with_jacobian = [&](const VectorType &src,
+                                                   VectorType &dst) -> void {
+              deallog << "Solving with Jacobian" << std::endl;
+              dst(0) = src(0) / (myshift + kappa);
+              dst.compress(VectorOperation::insert);
+            };
+          }
+      }
+    else
+      { // Here we instead use the explicit form
+        // This is the only function one would populate in case an explicit
+        // solver is used.
+        time_stepper.explicit_function =
+          [&](const real_type t, const VectorType &y, VectorType &res) -> void {
+          deallog << "Evaluating explicit function at t=" << t << std::endl;
+
+          if (t > last_eval_time + 0.1)
+            {
+              deallog << "Time step too large: last_eval_time="
+                      << last_eval_time << ", t=" << t << std::endl;
+              throw RecoverableUserCallbackError();
+            }
+
+          res(0) = -kappa * y(0);
+          res.compress(VectorOperation::insert);
+
+          last_eval_time = t;
+        };
+
+        // The explicit Jacobian callback is not needed in case
+        // an explicit solver is used. We need it when moving to
+        // an implicit solver like in the parameter file used for
+        // this test
+        if (setjac)
+          {
+            time_stepper.explicit_jacobian = [&](const real_type   t,
+                                                 const VectorType &y,
+                                                 MatrixType &      A,
+                                                 MatrixType &      P) -> void {
+              deallog << "Evaluating explicit Jacobian at t=" << t << std::endl;
+              P.set(0, 0, -kappa);
+              P.compress(VectorOperation::insert);
+            };
+          }
+      }
+
+    // Monitoring routine. Here we print diagnostic for the exact
+    // solution to the log file.
+    time_stepper.monitor = [&](const real_type   t,
+                               const VectorType &sol,
+                               const unsigned int /*step_number*/) -> void {
+      deallog << "Intermediate output:" << std::endl;
+      deallog << "  t =" << t << std::endl;
+      deallog << "  y =" << sol[0] << "  (exact: " << std::exp(-kappa * t)
+              << ')' << std::endl;
+    };
+  }
+
+  void
+  run()
+  {
+    // Set initial conditions
+    // This test uses a nonzero initial time
+    // We can access the time value with the
+    // get_time method
+    auto t = time_stepper.get_time();
+
+    VectorType y(MPI_COMM_SELF, 1, 1);
+    y[0] = 1;
+    y.compress(VectorOperation::insert);
+
+    // Integrate the ODE.
+    try
+      {
+        auto nt = time_stepper.solve(y);
+        deallog << "# Number of steps taken: " << nt << std::endl;
+      }
+    catch (const std::exception &exc)
+      {
+        deallog << "Time stepper aborted with an exception:" << std::endl
+                << exc.what() << std::endl;
+      }
+  }
+
+private:
+  TimeStepper time_stepper;
+  real_type   kappa;   // Defines the oscillator
+  real_type   myshift; // Used by the user solve
+  double      last_eval_time;
+};
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  PETScWrappers::TimeStepperData data;
+  ParameterHandler               prm;
+
+  data.add_parameters(prm);
+  deallog << "# Default Parameters" << std::endl;
+  prm.print_parameters(deallog.get_file_stream(), ParameterHandler::ShortText);
+
+  std::ifstream ifile(SOURCE_DIR "/petsc_ts_03_in.prm");
+  prm.parse_input(ifile);
+
+  deallog << "# Testing Parameters" << std::endl;
+  prm.print_parameters(deallog.get_file_stream(), ParameterHandler::ShortText);
+
+  for (int setjaci = 0; setjaci < 2; setjaci++)
+    {
+      bool setjac = setjaci ? true : false;
+
+      {
+        deallog << "# Test explicit interface (J " << setjac << ")"
+                << std::endl;
+        HarmonicOscillator ode_expl(1.0, data, setjac, false, false);
+        ode_expl.run();
+      }
+
+      {
+        deallog << "# Test implicit interface (J " << setjac << ")"
+                << std::endl;
+        HarmonicOscillator ode_impl(1.0, data, setjac, true, false);
+        ode_impl.run();
+      }
+    }
+
+  {
+    deallog << "# Test user interface" << std::endl;
+    HarmonicOscillator ode_user(1.0, data, true, true, true);
+    ode_user.run();
+  }
+}
diff --git a/tests/petsc/petsc_ts_04.output b/tests/petsc/petsc_ts_04.output
new file mode 100644 (file)
index 0000000..26db1c1
--- /dev/null
@@ -0,0 +1,200 @@
+
+DEAL::# Default Parameters
+subsection Error control
+  set absolute error tolerance = -1
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = -1
+  set minimum step size        = -1
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 0
+  set initial step size       = 0
+  set initial time            = 0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = 
+end
+DEAL::# Testing Parameters
+subsection Error control
+  set absolute error tolerance = 0.0001
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = 1e+7
+  set minimum step size        = 1e-7
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 2.0
+  set initial step size       = 0.2
+  set initial time            = 0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = bdf
+end
+DEAL::# Test explicit interface (J 0)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in line <0> of file <> in function
+    
+The violated condition was: 
+    
+Additional information: 
+    A user call-back function encountered a recoverable error, but the
+    underlying library that called the call-back did not manage to recover
+    from the error and aborted its operation.
+    
+    See the glossary entry on user call-back functions for more
+    information.
+--------------------------------------------------------
+
+DEAL::# Test implicit interface (J 0)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in line <0> of file <> in function
+    
+The violated condition was: 
+    
+Additional information: 
+    A user call-back function encountered a recoverable error, but the
+    underlying library that called the call-back did not manage to recover
+    from the error and aborted its operation.
+    
+    See the glossary entry on user call-back functions for more
+    information.
+--------------------------------------------------------
+
+DEAL::# Test explicit interface (J 1)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit Jacobian at t=0.100000
+DEAL::Evaluating explicit function at t=0.100000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Evaluating explicit Jacobian at t=0.200000
+DEAL::Evaluating explicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating explicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in line <0> of file <> in function
+    
+The violated condition was: 
+    
+Additional information: 
+    A user call-back function encountered a recoverable error, but the
+    underlying library that called the call-back did not manage to recover
+    from the error and aborted its operation.
+    
+    See the glossary entry on user call-back functions for more
+    information.
+--------------------------------------------------------
+
+DEAL::# Test implicit interface (J 1)
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit Jacobian at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit Jacobian at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in line <0> of file <> in function
+    
+The violated condition was: 
+    
+Additional information: 
+    A user call-back function encountered a recoverable error, but the
+    underlying library that called the call-back did not manage to recover
+    from the error and aborted its operation.
+    
+    See the glossary entry on user call-back functions for more
+    information.
+--------------------------------------------------------
+
+DEAL::# Test user interface
+DEAL::Intermediate output:
+DEAL::  t =0.00000
+DEAL::  y =1.00000  (exact: 1.00000)
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Setting up Jacobian at t=0.100000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.100000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Setting up Jacobian at t=0.200000
+DEAL::Solving with Jacobian
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Evaluating implicit function at t=0.200000
+DEAL::Intermediate output:
+DEAL::  t =0.200000
+DEAL::  y =0.823864  (exact: 0.818731)
+DEAL::Evaluating implicit function at t=0.400000
+DEAL::Time step too large: last_eval_time=0.200000, t=0.400000
+DEAL::Time stepper aborted with an exception:
+DEAL::
+--------------------------------------------------------
+An error occurred in line <0> of file <> in function
+    
+The violated condition was: 
+    
+Additional information: 
+    A user call-back function encountered a recoverable error, but the
+    underlying library that called the call-back did not manage to recover
+    from the error and aborted its operation.
+    
+    See the glossary entry on user call-back functions for more
+    information.
+--------------------------------------------------------
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.