--- /dev/null
+//-----------------------------------------------------------
+//
+// 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();
+ }
+}
--- /dev/null
+
+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.
+--------------------------------------------------------
+
projects / dealii.git / commitdiff