From: Stefano Zampini Date: Tue, 17 May 2022 15:45:34 +0000 (+0300) Subject: PETScWrappers: Add support for ODE solver X-Git-Tag: v9.5.0-rc1~334^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8bf193a0c314c9452688e37c725563266d15b1b8;p=dealii.git PETScWrappers: Add support for ODE solver --- diff --git a/include/deal.II/base/exceptions.h b/include/deal.II/base/exceptions.h index f59caf1852..28745c7fc9 100644 --- a/include/deal.II/base/exceptions.h +++ b/include/deal.II/base/exceptions.h @@ -860,6 +860,14 @@ namespace StandardExceptions << "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. */ diff --git a/include/deal.II/lac/petsc_ts.h b/include/deal.II/lac/petsc_ts.h new file mode 100644 index 0000000000..a5f6153bae --- /dev/null +++ b/include/deal.II/lac/petsc_ts.h @@ -0,0 +1,592 @@ +//----------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_PETSC +# include +# include +# include + +# include +# include +# include + +# include + +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 + 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 + 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 + implicit_jacobian; + + /** + * Callback for the computation of the explicit residual $G(t, y)$. + */ + std::function + explicit_function; + + /** + * Callback for the computation of the explicit Jacobian $\dfrac{\partial + * G}{\partial y}$. + */ + std::function + explicit_jacobian; + + /** + * Callback for monitoring the solution process. + * + * This function is called by TimeStepper at the beginning + * of each time step. + */ + std::function + 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 + 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 + 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 algebraic_components; + + protected: + /** + * The PETSc object + */ + TS ts; + + SmartPointer A; + SmartPointer P; + bool need_dae_tolerances; + }; + +# ifndef DOXYGEN + /* Inline functions */ + + template + inline TS + TimeStepper::petsc_ts() + { + return ts; + } + + template + inline MPI_Comm + TimeStepper::get_mpi_communicator() + const + { + return PetscObjectComm(reinterpret_cast(ts)); + } + + template + inline typename TimeStepper::real_type + TimeStepper::get_time() + { + PetscReal t; + PetscErrorCode ierr = TSGetTime(ts, &t); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + return t; + } + + template + inline typename TimeStepper::real_type + TimeStepper::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 diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h new file mode 100644 index 0000000000..93c0fb09c2 --- /dev/null +++ b/include/deal.II/lac/petsc_ts.templates.h @@ -0,0 +1,483 @@ +//----------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_PETSC +# include + +# include +# include + +# include +# include + +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 + TimeStepper::TimeStepper( + const TimeStepperData &data, + const MPI_Comm & mpi_comm) + { + AssertTS(TSCreate(mpi_comm, &ts)); + AssertTS(TSSetApplicationContext(ts, this)); + reinit(data); + } + + + + template + TimeStepper::~TimeStepper() + { + AssertTS(TSDestroy(&ts)); + } + + + + template + void + TimeStepper::reinit() + { + AssertTS(TSReset(ts)); + } + + + + template + void + TimeStepper::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(PETSC_DEFAULT); + const PetscReal rtol = data.relative_tolerance > 0.0 ? + data.relative_tolerance : + static_cast(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(PETSC_DEFAULT); + const PetscReal hmax = data.maximum_step_size > 0.0 ? + data.maximum_step_size : + static_cast(PETSC_DEFAULT); + AssertTS(TSAdaptSetStepLimits(tsadapt, hmin, hmax)); + } + + + + template + void + TimeStepper::reinit_matrix( + PMatrixType &P) + { + this->A = nullptr; + this->P = &P; + } + + + + template + void + TimeStepper::reinit_matrices( + AMatrixType &A, + PMatrixType &P) + { + this->A = &A; + this->P = &P; + } + + + + template + unsigned int + TimeStepper::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(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(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(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(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(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(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(ts))); + if (solve_for_jacobian_system) + { + precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int { + VectorType dst(static_cast(indst)); + const VectorType src(static_cast(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 + unsigned int + TimeStepper::solve(VectorType & y, + PMatrixType &P) + { + reinit_matrix(P); + return solve(y); + } + + + + template + unsigned int + TimeStepper::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 diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 1ada9bea31..8b6cdd0ded 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -98,6 +98,7 @@ set(_inst 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 @@ -107,7 +108,7 @@ if(DEAL_II_WITH_PETSC) 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 ) diff --git a/source/lac/petsc_ts.cc b/source/lac/petsc_ts.cc new file mode 100644 index 0000000000..d7018087d2 --- /dev/null +++ b/source/lac/petsc_ts.cc @@ -0,0 +1,41 @@ +// --------------------------------------------------------------------- +// +// 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 + +# ifdef DEAL_II_WITH_PETSC + +# include +# include + +DEAL_II_NAMESPACE_OPEN + + +template class PETScWrappers::TimeStepper<>; +template class PETScWrappers::TimeStepper; +template class PETScWrappers::TimeStepper; +template class PETScWrappers::TimeStepper; +template class PETScWrappers::TimeStepper< + PETScWrappers::MPI::BlockVector, + PETScWrappers::MPI::BlockSparseMatrix>; + + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_WITH_PETSC +#endif diff --git a/tests/petsc/petsc_ts_00.cc b/tests/petsc/petsc_ts_00.cc new file mode 100644 index 0000000000..d95ae8aad6 --- /dev/null +++ b/tests/petsc/petsc_ts_00.cc @@ -0,0 +1,264 @@ +//----------------------------------------------------------- +// +// 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 + +#include +#include +#include + +#include + +#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; +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 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; +} diff --git a/tests/petsc/petsc_ts_00.output b/tests/petsc/petsc_ts_00.output new file mode 100644 index 0000000000..3eba9270de --- /dev/null +++ b/tests/petsc/petsc_ts_00.output @@ -0,0 +1,446 @@ +# 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) diff --git a/tests/petsc/petsc_ts_00_in.prm b/tests/petsc/petsc_ts_00_in.prm new file mode 100644 index 0000000000..5127d37cc4 --- /dev/null +++ b/tests/petsc/petsc_ts_00_in.prm @@ -0,0 +1,10 @@ +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 diff --git a/tests/petsc/petsc_ts_01.cc b/tests/petsc/petsc_ts_01.cc new file mode 100644 index 0000000000..b2aee68d48 --- /dev/null +++ b/tests/petsc/petsc_ts_01.cc @@ -0,0 +1,88 @@ +//----------------------------------------------------------- +// +// 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 +#include + +#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; + +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; +} diff --git a/tests/petsc/petsc_ts_01.output b/tests/petsc/petsc_ts_01.output new file mode 100644 index 0000000000..761c7667f5 --- /dev/null +++ b/tests/petsc/petsc_ts_01.output @@ -0,0 +1,12 @@ + +DEAL:: +-------------------------------------------------------- +An error occurred in file in function + unsigned int dealii::PETScWrappers::TimeStepper::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" +-------------------------------------------------------- + diff --git a/tests/petsc/petsc_ts_02.cc b/tests/petsc/petsc_ts_02.cc new file mode 100644 index 0000000000..30ea46b1db --- /dev/null +++ b/tests/petsc/petsc_ts_02.cc @@ -0,0 +1,92 @@ +//----------------------------------------------------------- +// +// 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 + +#include +#include + +#include "../tests.h" + +/** + * Solves a dummy index-1 DAE to test the solver. + */ +using VectorType = PETScWrappers::MPI::Vector; +using TimeStepper = PETScWrappers::TimeStepper; +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; +} diff --git a/tests/petsc/petsc_ts_02.output b/tests/petsc/petsc_ts_02.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/petsc/petsc_ts_02.output @@ -0,0 +1 @@ +