<< "Please provide an implementation for the function \""
<< arg1 << "\"");
+ /**
+ * This exception is used if some user function returns nonzero exit codes.
+ */
+ DeclException2(ExcFunctionNonzeroReturn,
+ std::string,
+ int,
+ << "Function \"" << arg1 << "\" returned " << arg2);
+
/**
* This exception is used if some object is found uninitialized.
*/
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+//---------------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#ifndef dealii_petsc_ts_h
+#define dealii_petsc_ts_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+# include <deal.II/base/mpi.h>
+# include <deal.II/base/parameter_handler.h>
+# include <deal.II/base/smartpointer.h>
+
+# include <deal.II/lac/petsc_matrix_base.h>
+# include <deal.II/lac/petsc_precondition.h>
+# include <deal.II/lac/petsc_vector_base.h>
+
+# include <petscts.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace PETScWrappers
+{
+ /**
+ * Additional parameters that can be passed to the TimeStepper class.
+ */
+ class TimeStepperData
+ {
+ public:
+ /**
+ * Type that holds real-valued numbers.
+ *
+ * Used to represent time and norms tolerances.
+ */
+ using real_type = PetscReal;
+
+ /**
+ * Initialization parameters for TimeStepper.
+ *
+ * Running parameters:
+ *
+ * @param opts_prefix The string indicating the options prefix for command line customization.
+ * @param tstype The string indicating the PETSc solver type.
+ * @param initial_time Initial simulation time.
+ * @param final_time Final simulation time.
+ * @param initial_step_size Initial step size.
+ * @param max_steps Maximum number of steps allowed.
+ * @param match_step Whether or not to exactly stop at final time or step over it.
+ *
+ * Error parameters:
+ *
+ * @param tsadapttype The string indicating the PETSc time step adaptor type.
+ * @param minimum_step_size Minimum step size allowed.
+ * @param maximum_step_size Maximum step size allowed.
+ * @param absolute_tolerance Absolute error tolerance.
+ * @param relative_tolerance Relative error tolerance.
+ * @param ignore_algebraic_lte Ignore algebraic terms for
+ * error computations
+ *
+ * Note that one between `final_time` or `max_steps` must
+ * be specified by the user, otherwise PETSc will complain.
+ * Adaptive time stepping is disabled by default.
+ * Negative values indicate using PETSc's default.
+ *
+ * Note that all parameters values specified here can be overriden by
+ * command line choices.
+ */
+ TimeStepperData(
+ // Running parameters
+ const std::string &opts_prefix = "",
+ const std::string &tstype = "",
+ const real_type initial_time = 0.0,
+ const real_type final_time = 0.0,
+ const real_type initial_step_size = 0.0,
+ const int max_steps = -1,
+ const bool match_step = false,
+ // Error parameters
+ const std::string &tsadapttype = "none",
+ const real_type minimum_step_size = -1.0,
+ const real_type maximum_step_size = -1.0,
+ const real_type absolute_tolerance = -1.0,
+ const real_type relative_tolerance = -1.0,
+ const bool ignore_algebraic_lte = true)
+ : opts_prefix(opts_prefix)
+ , tstype(tstype)
+ , initial_time(initial_time)
+ , final_time(final_time)
+ , initial_step_size(initial_step_size)
+ , max_steps(max_steps)
+ , match_step(match_step)
+ , tsadapttype(tsadapttype)
+ , minimum_step_size(minimum_step_size)
+ , maximum_step_size(maximum_step_size)
+ , absolute_tolerance(absolute_tolerance)
+ , relative_tolerance(relative_tolerance)
+ , ignore_algebraic_lte(ignore_algebraic_lte)
+ {}
+
+ void
+ add_parameters(ParameterHandler &prm)
+ {
+ prm.enter_subsection("Running parameters");
+ prm.add_parameter(
+ "options prefix",
+ opts_prefix,
+ "The string indicating the options prefix for command line customization.");
+ prm.add_parameter("solver type",
+ tstype,
+ "The string indicating the PETSc TS type.");
+ prm.add_parameter("initial time",
+ initial_time,
+ "The value for the initial time.");
+ prm.add_parameter("final time",
+ final_time,
+ "The value for the final time.");
+ prm.add_parameter("initial step size",
+ initial_step_size,
+ "The value for the initial time step.");
+ prm.add_parameter("maximum number of steps",
+ max_steps,
+ "Maximum number of time steps allowed.");
+ prm.add_parameter(
+ "match final time",
+ match_step,
+ "Whether or not to exactly stop at final time or step over it.");
+ prm.leave_subsection();
+
+ prm.enter_subsection("Error control");
+ prm.add_parameter("adaptor type",
+ tsadapttype,
+ "The string for the TSAdapt type.");
+ prm.add_parameter("minimum step size",
+ minimum_step_size,
+ "Minimum time step size allowed.");
+ prm.add_parameter("maximum step size",
+ maximum_step_size,
+ "Maximum time step size allowed.");
+ prm.add_parameter("absolute error tolerance",
+ absolute_tolerance,
+ "Absolute error tolerance.");
+ prm.add_parameter("relative error tolerance",
+ relative_tolerance,
+ "Absolute error tolerance.");
+ prm.add_parameter(
+ "ignore algebraic lte",
+ ignore_algebraic_lte,
+ "Indicate whether or not to suppress algebraic variables "
+ "in the local truncation error test.");
+ prm.leave_subsection();
+ }
+
+ /**
+ * Options database prefix.
+ */
+ std::string opts_prefix;
+
+ /**
+ * PETSc solver type.
+ */
+ std::string tstype;
+
+ /**
+ * Initial time for the DAE.
+ */
+ real_type initial_time;
+
+ /**
+ * Final time.
+ */
+ real_type final_time;
+
+ /**
+ * Initial step size.
+ * Non-positive values are ignored.
+ */
+ real_type initial_step_size;
+
+ /**
+ * Maximum number of steps to be taken.
+ * Negative values are ignored.
+ */
+ int max_steps;
+
+ /**
+ * Match final time requested?
+ */
+ bool match_step;
+
+ /**
+ * PETSc time step adaptor type.
+ */
+ std::string tsadapttype;
+
+ /**
+ * Minimum allowed step size for adaptive time stepping.
+ * Non-positive values indicate to use PETSc's default.
+ */
+ real_type minimum_step_size;
+
+ /**
+ * Maximum allowed step size for adaptive time stepping.
+ * Non-positive values indicate to use PETSc's default.
+ */
+ real_type maximum_step_size;
+
+ /**
+ * Absolute error tolerance for adaptive time stepping.
+ * Negative values indicate to use PETSc's default.
+ */
+ real_type absolute_tolerance;
+
+ /**
+ * Relative error tolerance for adaptive time stepping.
+ * Negative values indicate to use PETSc's default.
+ */
+ real_type relative_tolerance;
+
+ /**
+ * Ignore algebraic terms for local truncation error.
+ */
+ bool ignore_algebraic_lte;
+ };
+
+ /**
+ * Interface to the PETSc TS solver for Ordinary Differential Equations
+ * and Differential-Algebraic Equations. The TS solver is described in the
+ * [PETSc manual](https://petsc.org/release/manual/ts/).
+ *
+ * This class supports two kinds of formulations.
+ * The explicit formulation:
+ * \f[
+ * \begin{cases}
+ * \dot y = G(t,y)\, , \\
+ * y(t_0) = y_0\, , \\
+ * \end{cases}
+ * \f]
+ * and the implicit formulation:
+ * \f[
+ * \begin{cases}
+ * F(t,y,\dot y) = 0\, , \\
+ * y(t_0) = y_0\, . \\
+ * \end{cases}
+ * \f]
+ *
+ * The interface to PETSc is realized by means of std::function callbacks
+ * like in the SUNDIALS::IDA and SUNDIALS::ARKode classes.
+ *
+ * TimeStepper supports any vector and matrix type having constructors and
+ * methods:
+ *
+ * @code
+ * class VectorType : public Subscriptor
+ * ...
+ * explicit VectorType(Vec);
+ * ...
+ * Vec & petsc_vector();
+ * ...
+ * @endcode
+ *
+ * @code
+ * class MatrixType : public Subscriptor
+ * ...
+ * explicit MatrixType(Mat);
+ * ...
+ * Mat & petsc_matrix();
+ * ...
+ * @endcode
+ *
+ * In particular, the supported types are the ones that can *wrap*
+ * PETSc's native vector and matrix classes, that are able to modify
+ * them in place, and that can return PETSc native types when requested.
+ *
+ * To use explicit solvers (like for example explicit Runge-Kutta methods),
+ * the user only needs to provide the implementation of $G$ via the
+ * TimeStepper::explicit_function. For implicit solvers, users have also the
+ * alternative of providing the $F$ function via
+ * TimeStepper::implicit_function. IMEX methods are also supported by
+ * providing both callbacks.
+ *
+ * The default linearization procedure of an implicit solver instantiated with
+ * this class consists in using Jacobian-Free-Newton-Krylov; the action of
+ * tangent matrices inside a linear solver process are approximated via
+ * matrix-free finite-differencing of the nonlinear residual equations that
+ * are ODE-solver specific. For details, consult the [PETSc
+ * manual](https://petsc.org/release/manual/snes/#sec-nlmatrixfree).
+ *
+ * In alternative, users can also provide the implementations of the
+ * *Jacobians*. This can be accomplished in two ways:
+ * - a-la-PETSc style using TimeStepper::implicit_jacobian
+ * and TimeStepper::explicit_jacobian.
+ * - a-la-Deal.II style using TimeStepper::setup_jacobian and
+ * TimeStepper::solve_for_jacobian_system
+ *
+ * In case both approaches are coded, the deal.II style
+ * will be used.
+ *
+ * The second approach is more in style with the deal.II philosophy
+ * and it also allows command line customization, like for example,
+ * @code
+ * ./myApp -snes_type newtontr -ksp_type cg
+ * @endcode
+ * in case the user wants to change the default nonlinear solver to
+ * a trust region solver and iterate on the tangent system with CG,
+ * still using TimeStepper::solve_for_jacobian_system as a preconditioner.
+ *
+ * The first approach has instead the advantage that only the matrix assembly
+ * procedure has to be coded, thus allowing quicker implementations and
+ * faster turnaround for experimenting with linear solver preconditioning
+ * configurations via command line customizations, like for example,
+ * @code
+ * ./myApp -ksp_type cg -pc_type gamg
+ * @endcode
+ */
+ template <typename VectorType = PETScWrappers::VectorBase,
+ typename PMatrixType = PETScWrappers::MatrixBase,
+ typename AMatrixType = PMatrixType>
+ class TimeStepper
+ {
+ public:
+ /**
+ * Type that holds real-valued numbers.
+ *
+ * Used to represent time and norms tolerances.
+ */
+ using real_type = PetscReal;
+
+ /**
+ * Constructor.
+ */
+ TimeStepper(const TimeStepperData &data = TimeStepperData(),
+ const MPI_Comm & mpi_comm = PETSC_COMM_WORLD);
+
+ /**
+ * Destructor.
+ */
+ ~TimeStepper();
+
+ /**
+ * Return the PETSc TS object.
+ */
+ TS
+ petsc_ts();
+
+ /**
+ * Return the underlying MPI communicator.
+ */
+ MPI_Comm
+ get_mpi_communicator() const;
+
+ /**
+ * Reset the solver, it does not change the customization.
+ */
+ void
+ reinit();
+
+ /**
+ * Reset solver.
+ * Change customization according to `data`.
+ */
+ void
+ reinit(const TimeStepperData &data);
+
+ /**
+ * Set the preconditioning matrix only.
+ *
+ * When used with TimeStepper::setup_jacobian and
+ * TimeStepper::solve_for_jacobian_system, PETSc will approximate the linear
+ * system matrix-vector product using an internal matrix-free
+ * representation.
+ *
+ * When used with TimeStepper::implicit_jacobian or
+ * TimeStepper::explicit_jacobian, PETSc will use the same matrix for both
+ * preconditioning and matrix-vector products.
+ */
+ void
+ reinit_matrix(PMatrixType &P);
+
+ /**
+ * Set both the linear system matrix and the preconditioning matrix
+ * that PETSc will use.
+ */
+ void
+ reinit_matrices(AMatrixType &A, PMatrixType &P);
+
+ /**
+ * Return current time.
+ */
+ real_type
+ get_time();
+
+ /**
+ * Return current time step.
+ */
+ real_type
+ get_time_step();
+
+ /**
+ * Integrate the differential-algebraic equations starting from @p y.
+ *
+ * This function returns the final number of computed steps.
+ * Upon returning, the @p y vector contains the solution of the DAE at
+ * the end time.
+ */
+ unsigned int
+ solve(VectorType &y);
+
+ /**
+ * Integrate the differential-algebraic equations starting from @p y.
+ *
+ * This function returns the final number of computed steps.
+ * Upon returning, the @p y vector contains the solution of the DAE at
+ * the end time.
+ *
+ * Here we also pass the matrix to handle the Jacobian.
+ */
+ unsigned int
+ solve(VectorType &y, PMatrixType &P);
+
+ /**
+ * Integrate the differential-algebraic equations starting from @p y.
+ *
+ * This function returns the final number of computed steps.
+ * Upon returning, the @p y vector contains the solution of the DAE at
+ * the end time.
+ *
+ * Here we also pass the matrices to handle Jacobians.
+ */
+ unsigned int
+ solve(VectorType &y, AMatrixType &A, PMatrixType &P);
+
+ /**
+ * Callback for the computation of the implicit residual $F(t, y, \dot y)$.
+ */
+ std::function<int(const real_type t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ VectorType & res)>
+ implicit_function;
+
+ /**
+ * Callback for the computation of the implicit Jacobian
+ * $\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot
+ * y}$
+ *
+ * All implicit solvers implementations are recast to use the above
+ * linearization. The $\alpha$ parameter is time-step and solver-type
+ * specific.
+ */
+ std::function<int(const real_type t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ const real_type alpha,
+ AMatrixType & A,
+ PMatrixType & P)>
+ implicit_jacobian;
+
+ /**
+ * Callback for the computation of the explicit residual $G(t, y)$.
+ */
+ std::function<int(const real_type t, const VectorType &y, VectorType &res)>
+ explicit_function;
+
+ /**
+ * Callback for the computation of the explicit Jacobian $\dfrac{\partial
+ * G}{\partial y}$.
+ */
+ std::function<int(const real_type t,
+ const VectorType &y,
+ AMatrixType & A,
+ PMatrixType & P)>
+ explicit_jacobian;
+
+ /**
+ * Callback for monitoring the solution process.
+ *
+ * This function is called by TimeStepper at the beginning
+ * of each time step.
+ */
+ std::function<int(const real_type t,
+ const VectorType & y,
+ const unsigned int step_number)>
+ monitor;
+
+ /**
+ * Set up Jacobian callback without matrices.
+ *
+ * This callback gives full control to users to set up the linearized
+ * equations
+ * $\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot
+ * y}$.
+ *
+ * All implicit solvers implementations are recast to use the above
+ * linearization. The $\alpha$ parameter is time-step and solver-type
+ * specific.
+ *
+ * Solvers must be provided via TimeStepper::solve_for_jacobian_system.
+ */
+ std::function<int(const real_type t,
+ const VectorType &y,
+ const VectorType &ydot,
+ const real_type alpha)>
+ setup_jacobian;
+
+ /**
+ * Solution of the Jacobian system set up with TimeStepper::setup_jacobian.
+ *
+ * This is used as a preconditioner inside the Krylov process.
+ */
+ std::function<int(const VectorType &src, VectorType &dst)>
+ solve_for_jacobian_system;
+
+ /**
+ * Return an index set containing the algebraic components.
+ *
+ * Implementation of this function is optional. If your equation is also
+ * algebraic (i.e., it contains algebraic constraints, or Lagrange
+ * multipliers), you should implement this function in order to return only
+ * these components of your system.
+ */
+ std::function<IndexSet()> algebraic_components;
+
+ protected:
+ /**
+ * The PETSc object
+ */
+ TS ts;
+
+ SmartPointer<AMatrixType, TimeStepper> A;
+ SmartPointer<PMatrixType, TimeStepper> P;
+ bool need_dae_tolerances;
+ };
+
+# ifndef DOXYGEN
+ /* Inline functions */
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ inline TS
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::petsc_ts()
+ {
+ return ts;
+ }
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ inline MPI_Comm
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::get_mpi_communicator()
+ const
+ {
+ return PetscObjectComm(reinterpret_cast<PetscObject>(ts));
+ }
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ inline typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time()
+ {
+ PetscReal t;
+ PetscErrorCode ierr = TSGetTime(ts, &t);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ return t;
+ }
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ inline typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time_step()
+ {
+ PetscReal dt;
+ PetscErrorCode ierr = TSGetTimeStep(ts, &dt);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ return dt;
+ }
+# endif
+
+} // namespace PETScWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PETSC
+
+#endif
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+//---------------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#ifndef dealii_petsc_ts_templates_h
+#define dealii_petsc_ts_templates_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+# include <deal.II/base/exceptions.h>
+
+# include <deal.II/lac/petsc_precondition.h>
+# include <deal.II/lac/petsc_ts.h>
+
+# include <petscdm.h>
+# include <petscts.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Shorthand notation for PETSc error codes.
+# define AssertTS(code) \
+ do \
+ { \
+ PetscErrorCode ierr = (code); \
+ AssertThrow(ierr == 0, ExcPETScError(ierr)); \
+ } \
+ while (0)
+
+// Shorthand notation for User error codes.
+# define AssertUser(code, name) \
+ do \
+ { \
+ int ierr = (code); \
+ AssertThrow(ierr == 0, \
+ StandardExceptions::ExcFunctionNonzeroReturn(name, ierr)); \
+ } \
+ while (0)
+
+namespace PETScWrappers
+{
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::TimeStepper(
+ const TimeStepperData &data,
+ const MPI_Comm & mpi_comm)
+ {
+ AssertTS(TSCreate(mpi_comm, &ts));
+ AssertTS(TSSetApplicationContext(ts, this));
+ reinit(data);
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::~TimeStepper()
+ {
+ AssertTS(TSDestroy(&ts));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit()
+ {
+ AssertTS(TSReset(ts));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit(
+ const TimeStepperData &data)
+ {
+ reinit();
+
+ // Solver type
+ if (data.tstype.size())
+ AssertTS(TSSetType(ts, data.tstype.c_str()));
+
+ // Options prefix
+ if (data.opts_prefix.size())
+ AssertTS(TSSetOptionsPrefix(ts, data.opts_prefix.c_str()));
+
+ // Time and steps limits
+ AssertTS(TSSetTime(ts, data.initial_time));
+ if (data.final_time > data.initial_time)
+ AssertTS(TSSetMaxTime(ts, data.final_time));
+ if (data.initial_step_size > 0.0)
+ AssertTS(TSSetTimeStep(ts, data.initial_step_size));
+ if (data.max_steps >= 0)
+ AssertTS(TSSetMaxSteps(ts, data.max_steps));
+
+ // Decide how to end the integration. Either stepover the final time or
+ // match it.
+ AssertTS(TSSetExactFinalTime(ts,
+ data.match_step ? TS_EXACTFINALTIME_MATCHSTEP :
+ TS_EXACTFINALTIME_STEPOVER));
+
+ // Adaptive tolerances
+ const PetscReal atol = data.absolute_tolerance > 0.0 ?
+ data.absolute_tolerance :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+ const PetscReal rtol = data.relative_tolerance > 0.0 ?
+ data.relative_tolerance :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+ AssertTS(TSSetTolerances(ts, atol, nullptr, rtol, nullptr));
+
+ // At this point we do not know the problem size so we cannot
+ // set variable tolerances for differential and algebratic equations
+ // Store this value and use it during solve.
+ this->need_dae_tolerances = data.ignore_algebraic_lte;
+
+ // Adaptive time stepping
+ TSAdapt tsadapt;
+ AssertTS(TSGetAdapt(ts, &tsadapt));
+ AssertTS(TSAdaptSetType(tsadapt, data.tsadapttype.c_str()));
+
+ // As of 3.19, PETSc does not propagate options prefixes to the
+ // adaptors.
+ if (data.opts_prefix.size())
+ AssertTS(TSAdaptSetOptionsPrefix(tsadapt, data.opts_prefix.c_str()));
+
+ // Time step limits
+ const PetscReal hmin = data.minimum_step_size > 0.0 ?
+ data.minimum_step_size :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+ const PetscReal hmax = data.maximum_step_size > 0.0 ?
+ data.maximum_step_size :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+ AssertTS(TSAdaptSetStepLimits(tsadapt, hmin, hmax));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrix(
+ PMatrixType &P)
+ {
+ this->A = nullptr;
+ this->P = &P;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrices(
+ AMatrixType &A,
+ PMatrixType &P)
+ {
+ this->A = &A;
+ this->P = &P;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType &y)
+ {
+ auto ts_ifunction_ =
+ [](TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx)
+ -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ VectorType xdotdealii(xdot);
+ VectorType fdealii(f);
+ AssertUser(user->implicit_function(t, xdealii, xdotdealii, fdealii),
+ "implicit_function");
+ PetscFunctionReturn(0);
+ };
+
+ auto ts_ijacobian_ = [](TS ts,
+ PetscReal t,
+ Vec x,
+ Vec xdot,
+ PetscReal s,
+ Mat A,
+ Mat P,
+ void * ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ VectorType xdotdealii(xdot);
+ AMatrixType Adealii(A);
+ PMatrixType Pdealii(P);
+
+ AssertUser(
+ user->implicit_jacobian(t, xdealii, xdotdealii, s, Adealii, Pdealii),
+ "implicit_jacobian");
+
+ // Handle the Jacobian-free case
+ // This call allow to resample the linearization point
+ // of the MFFD tangent operator
+ PetscBool flg;
+ AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+ if (flg)
+ {
+ AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+ AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+ }
+ PetscFunctionReturn(0);
+ };
+
+ auto ts_ijacobian_with_setup_ = [](TS ts,
+ PetscReal t,
+ Vec x,
+ Vec xdot,
+ PetscReal s,
+ Mat A,
+ Mat P,
+ void * ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ VectorType xdotdealii(xdot);
+ AMatrixType Adealii(A);
+ PMatrixType Pdealii(P);
+
+ user->A = &Adealii;
+ user->P = &Pdealii;
+ AssertUser(user->setup_jacobian(t, xdealii, xdotdealii, s),
+ "setup_jacobian");
+
+ // Handle the Jacobian-free case
+ // This call allow to resample the linearization point
+ // of the MFFD tangent operator
+ PetscBool flg;
+ AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+ if (flg)
+ {
+ AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+ AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+ }
+ PetscFunctionReturn(0);
+ };
+
+ auto ts_rhsfunction_ =
+ [](TS ts, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ VectorType fdealii(f);
+
+ AssertUser(user->explicit_function(t, xdealii, fdealii),
+ "explicit_function");
+ PetscFunctionReturn(0);
+ };
+
+ auto ts_rhsjacobian_ =
+ [](TS ts, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ AMatrixType Adealii(A);
+ PMatrixType Pdealii(P);
+
+ AssertUser(user->explicit_jacobian(t, xdealii, Adealii, Pdealii),
+ "explicit_jacobian");
+
+ // Handle the Jacobian-free case
+ // This call allow to resample the linearization point
+ // of the MFFD tangent operator
+ PetscBool flg;
+ AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+ if (flg)
+ {
+ AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+ AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+ }
+ PetscFunctionReturn(0);
+ };
+
+ auto ts_monitor_ =
+ [](TS ts, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ (void)ts;
+ auto user = static_cast<TimeStepper *>(ctx);
+
+ VectorType xdealii(x);
+ AssertUser(user->monitor(t, xdealii, it), "monitor");
+ PetscFunctionReturn(0);
+ };
+
+ AssertThrow(explicit_function || implicit_function,
+ StandardExceptions::ExcFunctionNotProvided(
+ "explicit_function || implicit_function"));
+
+ AssertTS(TSSetSolution(ts, y.petsc_vector()));
+
+ if (explicit_function)
+ AssertTS(TSSetRHSFunction(ts, nullptr, ts_rhsfunction_, this));
+
+ if (implicit_function)
+ AssertTS(TSSetIFunction(ts, nullptr, ts_ifunction_, this));
+
+ if (setup_jacobian)
+ {
+ AssertTS(TSSetIJacobian(ts,
+ A ? A->petsc_matrix() : nullptr,
+ P ? P->petsc_matrix() : nullptr,
+ ts_ijacobian_with_setup_,
+ this));
+
+ // Tell PETSc to setup a MFFD operator for the linear system matrix
+ SNES snes;
+ AssertTS(TSGetSNES(ts, &snes));
+ if (!A)
+ AssertTS(SNESSetUseMatrixFree(snes, PETSC_TRUE, PETSC_FALSE));
+
+ // Do not waste memory by creating a dummy AIJ matrix inside PETSc.
+ if (!P)
+ {
+ DM dm;
+ AssertTS(SNESGetDM(snes, &dm));
+ AssertTS(DMSetMatType(dm, MATSHELL));
+ }
+ }
+ else
+ {
+ if (explicit_jacobian)
+ {
+ AssertTS(TSSetRHSJacobian(ts,
+ A ? A->petsc_matrix() :
+ (P ? P->petsc_matrix() : nullptr),
+ P ? P->petsc_matrix() : nullptr,
+ ts_rhsjacobian_,
+ this));
+ }
+
+ if (implicit_jacobian)
+ {
+ AssertTS(TSSetIJacobian(ts,
+ A ? A->petsc_matrix() :
+ (P ? P->petsc_matrix() : nullptr),
+ P ? P->petsc_matrix() : nullptr,
+ ts_ijacobian_,
+ this));
+ }
+
+ // The user did not set any Jacobian callback. PETSc default in this
+ // case is to use FD and thus assemble a dense operator by finite
+ // differencing the residual callbacks. Here instead we decide to
+ // use a full matrix-free approach by default. This choice can always
+ // be overriden from command line.
+ if (!explicit_jacobian && !implicit_jacobian)
+ {
+ SNES snes;
+ AssertTS(TSGetSNES(ts, &snes));
+ AssertTS(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
+ }
+ }
+
+ // In case solve_for_jacobian_system is provided, create a shell
+ // preconditioner wrapping the user call. The default internal Krylov
+ // solver only applies the preconditioner. This choice
+ // can be overriden by command line and users can use any other
+ // Krylov method if their solve is not accurate enough.
+ // Using solve_for_jacobian_system as a preconditioner allow users
+ // to provide approximate solvers and possibly iterate on a matrix-free
+ // approximation of the tangent operator.
+ PreconditionShell precond(
+ PetscObjectComm(reinterpret_cast<PetscObject>(ts)));
+ if (solve_for_jacobian_system)
+ {
+ precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
+ VectorType dst(static_cast<const Vec &>(indst));
+ const VectorType src(static_cast<const Vec &>(insrc));
+ return solve_for_jacobian_system(src, dst);
+ };
+
+ // Default Krylov solver (preconditioner only)
+ SNES snes;
+ KSP ksp;
+ AssertTS(TSGetSNES(ts, &snes));
+ AssertTS(SNESGetKSP(snes, &ksp));
+ AssertTS(KSPSetType(ksp, KSPPREONLY));
+ AssertTS(KSPSetPC(ksp, precond.get_pc()));
+ }
+
+ // Attach user monitoring routine.
+ if (monitor)
+ AssertTS(TSMonitorSet(ts, ts_monitor_, this, nullptr));
+
+ // Allow command line customization.
+ AssertTS(TSSetFromOptions(ts));
+
+ // Handle algebraic components.
+ if (algebraic_components && need_dae_tolerances)
+ {
+ PetscReal atol, rtol;
+ AssertTS(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr));
+
+ Vec av, rv;
+ AssertTS(VecDuplicate(y.petsc_vector(), &av));
+ AssertTS(VecDuplicate(y.petsc_vector(), &rv));
+
+ VectorBase avdealii(av);
+ VectorBase rvdealii(rv);
+ avdealii = atol;
+ rvdealii = rtol;
+ for (auto i : algebraic_components())
+ {
+ avdealii[i] = -1.0;
+ rvdealii[i] = -1.0;
+ }
+ avdealii.compress(VectorOperation::insert);
+ rvdealii.compress(VectorOperation::insert);
+ AssertTS(TSSetTolerances(ts, atol, av, rtol, rv));
+ AssertTS(VecDestroy(&av));
+ AssertTS(VecDestroy(&rv));
+ }
+
+ // Having set everything up, now do the actual work
+ // and let PETSc do the time stepping.
+ AssertTS(TSSolve(ts, nullptr));
+
+ // Return the number of steps taken.
+ PetscInt nt;
+ AssertTS(TSGetStepNumber(ts, &nt));
+ return nt;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType & y,
+ PMatrixType &P)
+ {
+ reinit_matrix(P);
+ return solve(y);
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType & y,
+ AMatrixType &A,
+ PMatrixType &P)
+ {
+ reinit_matrices(A, P);
+ return solve(y);
+ }
+
+} // namespace PETScWrappers
+
+# undef AssertTS
+# undef AssertUser
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PETSC
+
+#endif
if(DEAL_II_WITH_PETSC)
set(_unity_include_src
${_unity_include_src}
+ petsc_communication_pattern.cc
petsc_full_matrix.cc
petsc_matrix_base.cc
petsc_matrix_free.cc
petsc_parallel_vector.cc
petsc_precondition.cc
petsc_solver.cc
- petsc_communication_pattern.cc
+ petsc_ts.cc
petsc_sparse_matrix.cc
petsc_vector_base.cc
)
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef DOXYGEN
+
+# include <deal.II/base/config.h>
+
+# ifdef DEAL_II_WITH_PETSC
+
+# include <deal.II/lac/petsc_block_sparse_matrix.h>
+# include <deal.II/lac/petsc_ts.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template class PETScWrappers::TimeStepper<>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::BlockVector>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector,
+ PETScWrappers::MPI::SparseMatrix>;
+template class PETScWrappers::TimeStepper<
+ PETScWrappers::MPI::BlockVector,
+ PETScWrappers::MPI::BlockSparseMatrix>;
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+# endif // DEAL_II_WITH_PETSC
+#endif
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#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"
+
+/**
+ * Solves the harmonic oscillator problem
+ *
+ * u'' = -k^2 u
+ * u (t0) = sin(k * t0)
+ * u'(t0) = k cos (k* t0)
+ *
+ * using the PETScWrappers::TimeStepper class
+ * that interfaces PETSc TS ODE solver object.
+ * The ODE can be an EXPLICIT first order ode:
+ *
+ * y[0]' = y[1] \
+ * -> y' = G(y,t)
+ * y[1]' = - k^2 y[0] /
+ *
+ * or specified in IMPLICIT form:
+ *
+ * y[0]' - y[1] = 0 \
+ * -> F(y',y,t) = 0
+ * y[1]' + k^2 y[0] = 0 /
+ *
+ * The exact solution is
+ *
+ * y[0](t) = sin(k t)
+ * y[1](t) = k cos(k t)
+ *
+ * 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,
+ std::ostream & _out)
+ : time_stepper(data)
+ , out(_out)
+ , kappa(_kappa)
+ {
+ // 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) -> int {
+ res(0) = y_dot(0) - y(1);
+ res(1) = y_dot(1) + kappa * kappa * y(0);
+ res.compress(VectorOperation::insert);
+ return 0;
+ };
+
+ // 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) -> int {
+ P.set(0, 0, shift);
+ P.set(0, 1, -1);
+ P.set(1, 0, kappa * kappa);
+ P.set(1, 1, shift);
+ P.compress(VectorOperation::insert);
+ return 0;
+ };
+ }
+ 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) -> int {
+ myshift = shift;
+ return 0;
+ };
+
+ // In the solve phase we se the stored shift to solve
+ // for the implicit Jacobian system
+ time_stepper.solve_for_jacobian_system =
+ [&](const VectorType &src, VectorType &dst) -> int {
+ auto sf = 1. / (kappa * kappa + myshift * myshift);
+ dst(0) = sf * (myshift * src(0) + src(1));
+ dst(1) = sf * (-kappa * kappa * src(0) + myshift * src(1));
+ dst.compress(VectorOperation::insert);
+ return 0;
+ };
+ }
+ }
+ 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) -> int {
+ res(0) = y(1);
+ res(1) = -kappa * kappa * y(0);
+ res.compress(VectorOperation::insert);
+ return 0;
+ };
+
+ // 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) -> int {
+ P.set(0, 0, 0);
+ P.set(0, 1, 1);
+ P.set(1, 0, -kappa * kappa);
+ P.set(1, 1, 0);
+ P.compress(VectorOperation::insert);
+ return 0;
+ };
+ }
+ }
+
+ // Monitoring routine. Here we print diagnostic for the exact
+ // solution to the log file.
+ time_stepper.monitor = [&](const real_type t,
+ const VectorType & y,
+ const unsigned int step_number) -> int {
+ std::vector<real_type> exact(2);
+ exact[0] = std::sin(kappa * t);
+ exact[1] = kappa * std::cos(kappa * t);
+ out << t << " " << y(0) << " (" << exact[0] << ")"
+ << " " << y(1) << " (" << exact[1] << ")" << std::endl;
+ return 0;
+ };
+ }
+
+ 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, 2, 2);
+ y[0] = std::sin(kappa * t);
+ y[1] = kappa * std::cos(kappa * t);
+ y.compress(VectorOperation::insert);
+
+ // Integrate the ODE.
+ time_stepper.solve(y);
+ }
+
+private:
+ TimeStepper time_stepper;
+ std::ostream &out; // Used by the monitoring routine
+ real_type kappa; // Defines the oscillator
+ real_type myshift; // Used by the user solve
+};
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ std::ofstream out("output");
+
+ PETScWrappers::TimeStepperData data;
+ ParameterHandler prm;
+
+ data.add_parameters(prm);
+ out << "# Default Parameters" << std::endl;
+ prm.print_parameters(out, ParameterHandler::ShortText);
+
+ std::ifstream ifile(SOURCE_DIR "/petsc_ts_00_in.prm");
+ prm.parse_input(ifile);
+
+ out << "# Testing Parameters" << std::endl;
+ prm.print_parameters(out, ParameterHandler::ShortText);
+
+ for (int setjaci = 0; setjaci < 2; setjaci++)
+ {
+ bool setjac = setjaci ? true : false;
+
+ {
+ out << "# Test explicit interface (J " << setjac << ")" << std::endl;
+ HarmonicOscillator ode_expl(1.0, data, setjac, false, false, out);
+ ode_expl.run();
+ }
+
+ {
+ out << "# Test implicit interface (J " << setjac << ")" << std::endl;
+ HarmonicOscillator ode_impl(1.0, data, setjac, true, false, out);
+ ode_impl.run();
+ }
+ }
+
+ {
+ out << "# Test user interface" << std::endl;
+ HarmonicOscillator ode_user(1.0, data, true, true, true, out);
+ ode_user.run();
+ }
+ return 0;
+}
--- /dev/null
+# 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
+# Testing Parameters
+subsection Error control
+ set absolute error tolerance = -1
+ 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 = 6.0
+ set initial step size = 0.0625
+ set initial time = 1.0
+ set match final time = false
+ set maximum number of steps = -1
+ set options prefix =
+ set solver type = bdf
+end
+# Test explicit interface (J 0)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test implicit interface (J 0)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test explicit interface (J 1)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test implicit interface (J 1)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test user interface
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
--- /dev/null
+subsection Running parameters
+ set initial time = 1.0
+ set final time = 6.0
+ set solver type = bdf
+ set initial step size = 0.0625
+end
+subsection Error control
+ set minimum step size = 1e-7
+ set maximum step size = 1e+7
+end
--- /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.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/petsc_ts.templates.h>
+
+#include "../tests.h"
+
+/**
+ * Tests user defined Vector and Matrix types
+ * and exceptions handling for PETSCWrappers::TimeStepper.
+ */
+
+class VectorType : public Subscriptor
+{
+public:
+ explicit VectorType(Vec v)
+ : v(v)
+ {}
+
+ Vec &
+ petsc_vector()
+ {
+ return v;
+ }
+
+private:
+ Vec v;
+};
+
+class MatrixType : public Subscriptor
+{
+public:
+ explicit MatrixType(Mat A)
+ : A(A)
+ {}
+
+ Mat &
+ petsc_matrix()
+ {
+ return A;
+ }
+
+private:
+ Mat A;
+};
+
+using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ VectorType v(nullptr);
+ MatrixType A(nullptr);
+
+ TimeStepper myode;
+
+ try
+ {
+ auto ts = myode.petsc_ts();
+ auto t0 = myode.get_time();
+ auto dt = myode.get_time_step();
+ myode.solve(v, A);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << exc.what() << std::endl;
+ }
+ return 0;
+}
--- /dev/null
+
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_ts.templates.h> in function
+ unsigned int dealii::PETScWrappers::TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = VectorType; PMatrixType = MatrixType; AMatrixType = MatrixType]
+The violated condition was:
+ explicit_function || implicit_function
+Additional information:
+ Please provide an implementation for the function "explicit_function
+ || implicit_function"
+--------------------------------------------------------
+
--- /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.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+//#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_ts.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include "../tests.h"
+
+/**
+ * Solves a dummy index-1 DAE to test the solver.
+ */
+using VectorType = PETScWrappers::MPI::Vector;
+using TimeStepper = PETScWrappers::TimeStepper<VectorType>;
+using real_type = TimeStepper::real_type;
+
+
+class Dummy : public TimeStepper
+{
+public:
+ Dummy()
+ {
+ // Customize solver: use BDF and basic adaptive time stepping
+ PETScWrappers::TimeStepperData data;
+ data.tstype = "bdf";
+ data.final_time = 1.0;
+ data.tsadapttype = "basic";
+ reinit(data);
+
+ // Here we solve the two variables system:
+ // x_dot = x
+ // z = x
+ // with y = (x,z)
+ implicit_function = [&](const real_type t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ VectorType & res) -> int {
+ res(0) = y_dot(0) - y(0);
+ res(1) = y(1) - y(0);
+ res.compress(VectorOperation::insert);
+ return 0;
+ };
+
+ // Return the index set representing the
+ // algebraic components
+ algebraic_components = [&]() -> IndexSet {
+ IndexSet algebraic_is;
+
+ algebraic_is.set_size(2);
+ algebraic_is.add_index(1);
+ return algebraic_is;
+ };
+ }
+
+ unsigned int
+ solve()
+ {
+ VectorType y(MPI_COMM_SELF, 2, 2);
+ y(0) = 1.0;
+ y(1) = y(0);
+ y.compress(VectorOperation::insert);
+
+ return TimeStepper::solve(y);
+ }
+};
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ Dummy dae;
+ dae.solve();
+ return 0;
+}
projects / dealii.git / commitdiff