From 05957f979fd069deda653e880bf2327a5028f4b3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 25 May 2023 14:06:35 -0600 Subject: [PATCH] Add a test that shows that PETSc TS cannot recover from errors. --- tests/petsc/petsc_ts_04.cc | 275 +++++++++++++++++++++++++++++++++ tests/petsc/petsc_ts_04.output | 200 ++++++++++++++++++++++++ 2 files changed, 475 insertions(+) create mode 100644 tests/petsc/petsc_ts_04.cc create mode 100644 tests/petsc/petsc_ts_04.output diff --git a/tests/petsc/petsc_ts_04.cc b/tests/petsc/petsc_ts_04.cc new file mode 100644 index 0000000000..0d1d2ad155 --- /dev/null +++ b/tests/petsc/petsc_ts_04.cc @@ -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 + +#include +#include +#include + +#include + +#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; +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 index 0000000000..26db1c105a --- /dev/null +++ b/tests/petsc/petsc_ts_04.output @@ -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. +-------------------------------------------------------- + -- 2.39.5