From 67c15c819d1d2f23c58c0fb7a3753dd0a917def6 Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Thu, 13 May 2021 19:44:11 +0200 Subject: [PATCH] Make IDA work with new sundials. --- .../minor/20210513NicolaGiulianiLucaHeltai | 3 + include/deal.II/sundials/ida.h | 160 +++++++--- source/sundials/ida.cc | 301 ++++++++++++++---- source/sundials/n_vector.inst.in | 2 +- tests/sundials/ida_01.cc | 13 +- .../ida_01.with_sundials.geq.4.0.0.output | 33 ++ tests/sundials/ida_02.cc | 167 ++++++++++ .../ida_02.with_sundials.geq.4.0.0.output | 33 ++ 8 files changed, 614 insertions(+), 98 deletions(-) create mode 100644 doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai create mode 100644 tests/sundials/ida_01.with_sundials.geq.4.0.0.output create mode 100644 tests/sundials/ida_02.cc create mode 100644 tests/sundials/ida_02.with_sundials.geq.4.0.0.output diff --git a/doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai b/doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai new file mode 100644 index 0000000000..51815190f5 --- /dev/null +++ b/doc/news/changes/minor/20210513NicolaGiulianiLucaHeltai @@ -0,0 +1,3 @@ +New: SUNDIALS::IDA is now compatible with SUNDIALS > 4.0.0 +
+(Nicola Giuliani, Luca Heltai, 2021/05/13) diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 268bb55ca7..4a6fd13415 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -20,44 +20,44 @@ #include #ifdef DEAL_II_WITH_SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - -# include -# include -# include -# include -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif -# include -# include - -# ifdef DEAL_II_SUNDIALS_WITH_IDAS -# include -# else -# include -# endif - -# include -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) -# include -# include -# include -# endif -# include - -# include -# include -# include - -# include +# include +# include +# include +# include +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif +# include +# include + +# ifdef DEAL_II_SUNDIALS_WITH_IDAS +# include +# else +# include +# endif + +# include +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# include +# include +# include +# endif +# include + +# include + +# include +# include +# include + +# include DEAL_II_NAMESPACE_OPEN // Shorthand notation for IDA error codes. -# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code)) +# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code)) namespace SUNDIALS { @@ -286,6 +286,9 @@ namespace SUNDIALS * @param maximum_order Maximum BDF order * @param maximum_non_linear_iterations Maximum number of nonlinear * iterations + * @param ls_norm_factor Converting factor from the integrator tolerance + * to the linear solver tolerance + * iterations * * Error parameters: * @@ -310,6 +313,7 @@ namespace SUNDIALS const double minimum_step_size = 1e-6, const unsigned int maximum_order = 5, const unsigned int maximum_non_linear_iterations = 10, + const double ls_norm_factor = 0, // Error parameters const double absolute_tolerance = 1e-6, const double relative_tolerance = 1e-5, @@ -331,6 +335,7 @@ namespace SUNDIALS , reset_type(reset_type) , maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic) , maximum_non_linear_iterations(maximum_non_linear_iterations) + , ls_norm_factor(ls_norm_factor) {} /** @@ -452,6 +457,9 @@ namespace SUNDIALS }); prm.add_parameter("Maximum number of nonlinear iterations", maximum_non_linear_iterations_ic); + prm.add_parameter( + "Factor to use when converting from the integrator tolerance to the linear solver tolerance", + ls_norm_factor); prm.leave_subsection(); } @@ -537,6 +545,12 @@ namespace SUNDIALS * Maximum number of iterations for Newton method during time advancement. */ unsigned int maximum_non_linear_iterations; + + /** + * Factor to use when converting from the integrator tolerance to the + * linear solver tolerance + */ + double ls_norm_factor; }; /** @@ -581,13 +595,25 @@ namespace SUNDIALS * @param mpi_comm MPI communicator */ IDA(const AdditionalData &data = AdditionalData(), - const MPI_Comm & mpi_comm = MPI_COMM_WORLD); + const MPI_Comm mpi_comm = MPI_COMM_WORLD); /** * Destructor. */ ~IDA(); + + void + set_n_iter(const int n_iter) + { + this->n_iter = n_iter; + } + + int + get_n_iter() const + { + return this->n_iter; + } /** * Integrate differential-algebraic equations. This function returns the * final number of computed steps. @@ -703,10 +729,59 @@ namespace SUNDIALS * then last function will be attempted again * - <0: Unrecoverable error the computation will be aborted and an * assertion will be thrown. + * + * @warning Starting with SUNDIALS 4.1, SUNDIALS provides the possibility of + * specifying the tolerance for the resolution. A part from the tolerance + * only `rhs` is provided and `dst` needs to be returned. */ + DEAL_II_DEPRECATED_EARLY std::function solve_jacobian_system; + + /** + * Solve the Jacobian linear system. This function will be called by IDA + * (possibly several times) after setup_jacobian() has been called at least + * once. IDA tries to do its best to call setup_jacobian() the minimum + * amount of times. If convergence can be achieved without updating the + * Jacobian, then IDA does not call setup_jacobian() again. If, on the + * contrary, internal IDA convergence tests fail, then IDA calls again + * setup_jacobian() with updated vectors and coefficients so that successive + * calls to solve_jacobian_systems() lead to better convergence in the + * Newton process. + * + * The jacobian $J$ should be (an approximation of) the system Jacobian + * \f[ + * J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + + * \alpha \dfrac{\partial F}{\partial \dot y}. + * \f] + * + * Arguments to the function are: + * + * @param[in] rhs The system right hand side to solve for. + * @param[out] dst The solution of $J^{-1} * src$. + * @param[out] n_iter the number of iterations required to solve the + * jacobian system + * @param[in] tolerance The tolerance with which to solve the linear system + * of equations. + * + * A call to this function should store in `dst` the result of $J^{-1}$ + * applied to `src`, i.e., `J*dst = src`. It is the users responsibility + * to set up proper solvers and preconditioners inside this function. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error (IDAReinit will be called if this happens, and + * then last function will be attempted again + * - <0: Unrecoverable error the computation will be aborted and an + * assertion will be thrown. + */ + std::function + solve_with_jacobian; + /** * Process solution. This function is called by IDA at fixed time steps, * every `output_period` seconds, and it is passed a polynomial @@ -822,6 +897,12 @@ namespace SUNDIALS */ N_Vector diff_id; + /** + * Number of iteration required to solve the Jacobian system + */ + int n_iter; + + /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI @@ -829,13 +910,15 @@ namespace SUNDIALS */ MPI_Comm communicator; + + /** * Memory pool of vectors. */ GrowingVectorMemory mem; -# ifdef DEAL_II_WITH_PETSC -# ifdef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifdef PETSC_USE_COMPLEX static_assert(!std::is_same::value, "Sundials does not support complex scalar types, " "but PETSc is configured to use a complex scalar type!"); @@ -844,15 +927,14 @@ namespace SUNDIALS !std::is_same::value, "Sundials does not support complex scalar types, " "but PETSc is configured to use a complex scalar type!"); -# endif // PETSC_USE_COMPLEX -# endif // DEAL_II_WITH_PETSC +# endif // PETSC_USE_COMPLEX +# endif // DEAL_II_WITH_PETSC }; } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -# endif #endif // DEAL_II_WITH_SUNDIALS #endif diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 1fe5fc9061..6788c62fca 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -12,38 +12,39 @@ // the top level directory of deal.II. // //----------------------------------------------------------- - - #include #include #ifdef DEAL_II_WITH_SUNDIALS +# include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - -# include - -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif -# include +# include +# include +# include +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) # ifdef DEAL_II_SUNDIALS_WITH_IDAS # include # else # include # endif - -# include -# include +# endif +# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) +# include +# endif +# include +# include DEAL_II_NAMESPACE_OPEN @@ -85,6 +86,7 @@ namespace SUNDIALS +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template int t_dae_lsetup(IDAMem IDA_mem, @@ -121,6 +123,7 @@ namespace SUNDIALS } + template int t_dae_solve(IDAMem IDA_mem, @@ -152,10 +155,108 @@ namespace SUNDIALS return err; } + + +# else + template + int + t_dae_jacobian_setup(realtype tt, + realtype cj, + N_Vector yy, + N_Vector yp, + N_Vector /* residual */, + SUNMatrix /* ignored */, + void *user_data, + N_Vector /* tmp1 */, + N_Vector /* tmp2 */, + N_Vector /* tmp3 */) + { + Assert(user_data != nullptr, ExcInternalError()); + IDA &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer src_yp(mem); + solver.reinit_vector(*src_yp); + + copy(*src_yy, yy); + copy(*src_yp, yp); + + + int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj); + + + return err; + } + + + + template + int + t_dae_solve_jacobian_system(SUNLinearSolver LS, + SUNMatrix /*ignored*/, + N_Vector x, + N_Vector b, + realtype /*tol*/) + { + const IDA &solver = + *static_cast *>(LS->content); + + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); + + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); + + copy(*src_b, b); + + const int err = solver.solve_jacobian_system(*src_b, *dst_x); + + copy(x, *dst_x); + + return err; + } + + + template + int + t_dae_solve_with_jacobian(SUNLinearSolver LS, + SUNMatrix /*ignored*/, + N_Vector x, + N_Vector b, + realtype tol) + { + IDA &solver = *static_cast *>(LS->content); + + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); + + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); + + copy(*src_b, b); + int n_iter; + const int err = solver.solve_with_jacobian(*src_b, *dst_x, n_iter, tol); + solver.set_n_iter(n_iter > 0 ? n_iter : 1); + copy(x, *dst_x); + + return err; + } +# endif } // namespace + + template - IDA::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm) + IDA::IDA(const AdditionalData &data, const MPI_Comm mpi_comm) : data(data) , ida_mem(nullptr) , yy(nullptr) @@ -174,14 +275,14 @@ namespace SUNDIALS { if (ida_mem) IDAFree(&ida_mem); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const int ierr = MPI_Comm_free(&communicator); (void)ierr; AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr)); } -# endif +# endif } @@ -196,13 +297,14 @@ namespace SUNDIALS double h = data.initial_step_size; unsigned int step_number = 0; + this->n_iter = 1; int status; (void)status; // The solution is stored in // solution. Here we take only a // view of it. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = solution.locally_owned_elements(); @@ -218,7 +320,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { Assert(is_serial_vector::value, ExcInternalError( @@ -256,7 +358,7 @@ namespace SUNDIALS } // Free the vectors which are no longer used. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(yy); @@ -265,7 +367,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -276,6 +378,8 @@ namespace SUNDIALS return step_number; } + + template void IDA::reset(const double current_time, @@ -291,11 +395,10 @@ namespace SUNDIALS ida_mem = IDACreate(); - // Free the vectors which are no longer used. if (yy) { -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(yy); @@ -304,7 +407,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -316,7 +419,7 @@ namespace SUNDIALS int status; (void)status; system_size = solution.size(); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = solution.locally_owned_elements(); @@ -332,7 +435,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { yy = N_VNew_Serial(system_size); yp = N_VNew_Serial(system_size); @@ -345,7 +448,6 @@ namespace SUNDIALS status = IDAInit(ida_mem, t_dae_residual, current_time, yy, yp); AssertIDA(status); - if (get_local_tolerances) { copy(abs_tolls, get_local_tolerances()); @@ -392,14 +494,115 @@ namespace SUNDIALS AssertIDA(status); // Initialize solver +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) auto IDA_mem = static_cast(ida_mem); IDA_mem->ida_lsetup = t_dae_lsetup; - IDA_mem->ida_lsolve = t_dae_solve; + + if (solve_jacobian_system) + IDA_mem->ida_lsolve = t_dae_solve; + else + AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system")); # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) IDA_mem->ida_setupNonNull = true; # endif +# else + SUNMatrix J = nullptr; + SUNLinearSolver LS = nullptr; + + // and attach it to the SUNLinSol object. The functions that will get + // called do not actually receive the IDAMEM object, just the LS + // object, so we have to store a pointer to the current + // object in the LS object + LS = SUNLinSolNewEmpty(); + LS->content = this; + + LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type { + return SUNLINEARSOLVER_MATRIX_ITERATIVE; + }; + + LS->ops->free = [](SUNLinearSolver LS) -> int { + if (LS->content) + { + LS->content = nullptr; + } + if (LS->ops) + { + free(LS->ops); + LS->ops = nullptr; + } + free(LS); + LS = nullptr; + return 0; + }; + if (solve_with_jacobian) + { + LS->ops->solve = t_dae_solve_with_jacobian; + } + else if (solve_jacobian_system) + { + LS->ops->solve = t_dae_solve_jacobian_system; + } + else + { + AssertThrow(false, ExcFunctionNotProvided("solve_with_jacobian")); + } + // When we set an iterative solver IDA requires that resid is provided. From + // SUNDIALS docs If an iterative method computes the preconditioned initial + // residual and returns with a successful solve without performing any + // iterations (i.e., either the initial guess or the preconditioner is + // sufficiently accurate), then this optional routine may be called by the + // SUNDIALS package. This routine should return the N_Vector containing the + // preconditioned initial residual vector. + LS->ops->resid = [](SUNLinearSolver /*ignored*/) -> N_Vector { + return nullptr; + }; + // When we set an iterative solver IDA requires that last number of + // iteration is provided. Since we can't know what kind of solver the user + // has provided we set 1. This is clearly suboptimal. + LS->ops->numiters = [](SUNLinearSolver LS) -> int { + IDA &solver = *static_cast *>(LS->content); + return solver.get_n_iter(); + }; + // Even though we don't use it, IDA still wants us to set some + // kind of matrix object for the nonlinear solver. This is because + // if we don't set it, it won't call the functions that set up + // the matrix object (i.e., the argument to the 'IDASetJacFn' + // function below). + J = SUNMatNewEmpty(); + J->content = this; + + J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { + return SUNMATRIX_CUSTOM; + }; + + J->ops->destroy = [](SUNMatrix A) { + if (A->content) + { + A->content = nullptr; + } + if (A->ops) + { + free(A->ops); + A->ops = nullptr; + } + free(A); + A = nullptr; + }; + + // Now set the linear system and Jacobian objects in the solver: + status = IDASetLinearSolver(ida_mem, LS, J); + AssertIDA(status); + + status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor); + AssertIDA(status); + // Finally tell IDA about + // it as well. The manual says that this must happen *after* + // calling IDASetLinearSolver + status = IDASetJacFn(ida_mem, &t_dae_jacobian_setup); + AssertIDA(status); +# endif status = IDASetMaxOrd(ida_mem, data.maximum_order); AssertIDA(status); @@ -457,20 +660,6 @@ namespace SUNDIALS return ret; }; - setup_jacobian = [](const double, - const VectorType &, - const VectorType &, - const double) -> int { - int ret = 0; - AssertThrow(false, ExcFunctionNotProvided("setup_jacobian")); - return ret; - }; - - solve_jacobian_system = [](const VectorType &, VectorType &) -> int { - int ret = 0; - AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system")); - return ret; - }; output_step = [](const double, const VectorType &, @@ -492,26 +681,24 @@ namespace SUNDIALS template class IDA>; template class IDA>; -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI -# ifdef DEAL_II_WITH_TRILINOS +# ifdef DEAL_II_WITH_TRILINOS template class IDA; template class IDA; -# endif // DEAL_II_WITH_TRILINOS +# endif // DEAL_II_WITH_TRILINOS -# ifdef DEAL_II_WITH_PETSC -# ifndef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifndef PETSC_USE_COMPLEX template class IDA; template class IDA; -# endif // PETSC_USE_COMPLEX -# endif // DEAL_II_WITH_PETSC +# endif // PETSC_USE_COMPLEX +# endif // DEAL_II_WITH_PETSC -# endif // DEAL_II_WITH_MPI +# endif // DEAL_II_WITH_MPI } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -# endif - #endif // DEAL_II_WITH_SUNDIALS diff --git a/source/sundials/n_vector.inst.in b/source/sundials/n_vector.inst.in index 4c11264684..899195b795 100644 --- a/source/sundials/n_vector.inst.in +++ b/source/sundials/n_vector.inst.in @@ -85,4 +85,4 @@ for (V : DEAL_II_VEC_TEMPLATES) template class SUNDIALS::internal::NVectorView; template class SUNDIALS::internal::NVectorView; } -#endif \ No newline at end of file +#endif diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index 239424b53f..a44c679751 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -24,7 +24,8 @@ /** - * Solve the Harmonic oscillator problem. + * Solve the Harmonic oscillator problem, using a direct solver for the + * jacobian system. * * u'' = -k^2 u * u (0) = 0 @@ -98,12 +99,22 @@ public: return 0; }; + time_stepper.solve_jacobian_system = [&](const VectorType &src, VectorType & dst) -> int { Jinv.vmult(dst, src); return 0; }; + time_stepper.solve_with_jacobian = [&](const VectorType &src, + VectorType & dst, + int & n_iter, + const double) -> int { + Jinv.vmult(dst, src); + n_iter = 1; + return 0; + }; + time_stepper.output_step = [&](const double t, const VectorType & sol, const VectorType & sol_dot, diff --git a/tests/sundials/ida_01.with_sundials.geq.4.0.0.output b/tests/sundials/ida_01.with_sundials.geq.4.0.0.output new file mode 100644 index 0000000000..7238b14d99 --- /dev/null +++ b/tests/sundials/ida_01.with_sundials.geq.4.0.0.output @@ -0,0 +1,33 @@ +0 0 1 1 0 +0.2 0.198666 0.98005 0.980056 -0.198635 +0.4 0.389414 0.92105 0.921063 -0.3894 +0.6 0.564639 0.825325 0.825346 -0.564653 +0.8 0.717354 0.696694 0.696713 -0.717367 +1 0.841469 0.540288 0.540295 -0.841478 +1.2 0.932034 0.362343 0.362339 -0.932034 +1.4 0.985441 0.169954 0.169952 -0.985442 +1.6 0.999562 -0.0292108 -0.0292135 -0.999563 +1.8 0.973833 -0.227211 -0.227214 -0.973833 +2 0.909281 -0.416152 -0.416155 -0.909281 +2.2 0.808479 -0.588503 -0.588506 -0.808478 +2.4 0.675445 -0.737392 -0.737394 -0.675443 +2.6 0.515484 -0.856883 -0.856885 -0.515482 +2.8 0.334972 -0.942213 -0.942214 -0.334969 +3 0.141106 -0.989979 -0.98998 -0.141103 +3.2 -0.0583858 -0.998279 -0.998279 0.0583887 +3.4 -0.255549 -0.966779 -0.966779 0.255552 +3.6 -0.442525 -0.896738 -0.896737 0.442528 +3.8 -0.611858 -0.790946 -0.790945 0.61186 +4 -0.756798 -0.653622 -0.65362 0.7568 +4.2 -0.871566 -0.49024 -0.490238 0.871568 +4.4 -0.951588 -0.307314 -0.307312 0.95159 +4.6 -0.993674 -0.112137 -0.112134 0.993674 +4.8 -0.996144 0.0875108 0.0875137 0.996144 +5 -0.958901 0.28367 0.283673 0.958901 +5.2 -0.88343 0.468519 0.468522 0.883429 +5.4 -0.772739 0.63469 0.634693 0.772738 +5.6 -0.631242 0.775558 0.77556 0.63124 +5.8 -0.464579 0.885506 0.885508 0.464577 +6 -0.279395 0.960152 0.960154 0.279393 +6.2 -0.0830731 0.99652 0.996521 0.0830703 +6.3 0.016828 0.999835 0.999835 -0.0168267 diff --git a/tests/sundials/ida_02.cc b/tests/sundials/ida_02.cc new file mode 100644 index 0000000000..2da335b998 --- /dev/null +++ b/tests/sundials/ida_02.cc @@ -0,0 +1,167 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +/** + * Solve the Harmonic oscillator problem, using an iterative solver for the + * jacobian system (with the interface that is available only in version > + * 4.0.0) + * + * u'' = -k^2 u + * u (0) = 0 + * u'(0) = k + * + * write in terms of a first order ode: + * + * y[0]' - y[1] = 0 + * y[1]' + k^2 y[0] = 0 + * + * That is + * + * F(y', y, t) = y' + A y = 0 + * + * A = [ 0 , -1; k^2, 0 ] + * + * y_0 = 0, k + * y_0' = k, 0 + * + * The exact solution is + * + * y[0](t) = sin(k t) + * y[1](t) = k cos(k t) + * + * The Jacobian to assemble is the following: + * + * J = alpha I + A + */ +class HarmonicOscillator +{ +public: + HarmonicOscillator( + double _kappa, + const typename SUNDIALS::IDA>::AdditionalData &data) + : time_stepper(data) + , y(2) + , y_dot(2) + , J(2, 2) + , A(2, 2) + , Jinv(2, 2) + , kappa(_kappa) + , out("output") + { + using VectorType = Vector; + + time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(2); }; + + + time_stepper.residual = [&](const double t, + const VectorType &y, + const VectorType &y_dot, + VectorType & res) -> int { + res = y_dot; + A.vmult_add(res, y); + return 0; + }; + + time_stepper.setup_jacobian = [&](const double, + const VectorType &, + const VectorType &, + const double alpha) -> int { + A(0, 1) = -1.0; + A(1, 0) = kappa * kappa; + + J = A; + + J(0, 0) = alpha; + J(1, 1) = alpha; + + return 0; + }; + + time_stepper.solve_with_jacobian = [&](const VectorType &src, + VectorType & dst, + int & n_iter, + const double tolerance) -> int { + SolverControl solver_control(1000, tolerance); + SolverGMRES> solver(solver_control); + solver.solve(J, dst, src, PreconditionIdentity()); + n_iter = solver_control.last_step() > 0 ? solver_control.last_step() : 1; + return 0; + }; + + time_stepper.output_step = [&](const double t, + const VectorType & sol, + const VectorType & sol_dot, + const unsigned int step_number) -> int { + out << t << " " << sol[0] << " " << sol[1] << " " << sol_dot[0] << " " + << sol_dot[1] << std::endl; + return 0; + }; + } + + void + run() + { + y[1] = kappa; + y_dot[0] = kappa; + time_stepper.solve_dae(y, y_dot); + } + SUNDIALS::IDA> time_stepper; + +private: + Vector y; + Vector y_dot; + FullMatrix J; + FullMatrix A; + FullMatrix Jinv; + double kappa; + + std::ofstream out; +}; + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + SUNDIALS::IDA>::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + // std::ofstream ofile(SOURCE_DIR "/ida_01.prm"); + // prm.print_parameters(ofile, ParameterHandler::ShortText); + // ofile.close(); + + std::ifstream ifile(SOURCE_DIR "/ida_01.prm"); + prm.parse_input(ifile); + + + HarmonicOscillator ode(1.0, data); + ode.run(); + return 0; +} diff --git a/tests/sundials/ida_02.with_sundials.geq.4.0.0.output b/tests/sundials/ida_02.with_sundials.geq.4.0.0.output new file mode 100644 index 0000000000..7238b14d99 --- /dev/null +++ b/tests/sundials/ida_02.with_sundials.geq.4.0.0.output @@ -0,0 +1,33 @@ +0 0 1 1 0 +0.2 0.198666 0.98005 0.980056 -0.198635 +0.4 0.389414 0.92105 0.921063 -0.3894 +0.6 0.564639 0.825325 0.825346 -0.564653 +0.8 0.717354 0.696694 0.696713 -0.717367 +1 0.841469 0.540288 0.540295 -0.841478 +1.2 0.932034 0.362343 0.362339 -0.932034 +1.4 0.985441 0.169954 0.169952 -0.985442 +1.6 0.999562 -0.0292108 -0.0292135 -0.999563 +1.8 0.973833 -0.227211 -0.227214 -0.973833 +2 0.909281 -0.416152 -0.416155 -0.909281 +2.2 0.808479 -0.588503 -0.588506 -0.808478 +2.4 0.675445 -0.737392 -0.737394 -0.675443 +2.6 0.515484 -0.856883 -0.856885 -0.515482 +2.8 0.334972 -0.942213 -0.942214 -0.334969 +3 0.141106 -0.989979 -0.98998 -0.141103 +3.2 -0.0583858 -0.998279 -0.998279 0.0583887 +3.4 -0.255549 -0.966779 -0.966779 0.255552 +3.6 -0.442525 -0.896738 -0.896737 0.442528 +3.8 -0.611858 -0.790946 -0.790945 0.61186 +4 -0.756798 -0.653622 -0.65362 0.7568 +4.2 -0.871566 -0.49024 -0.490238 0.871568 +4.4 -0.951588 -0.307314 -0.307312 0.95159 +4.6 -0.993674 -0.112137 -0.112134 0.993674 +4.8 -0.996144 0.0875108 0.0875137 0.996144 +5 -0.958901 0.28367 0.283673 0.958901 +5.2 -0.88343 0.468519 0.468522 0.883429 +5.4 -0.772739 0.63469 0.634693 0.772738 +5.6 -0.631242 0.775558 0.77556 0.63124 +5.8 -0.464579 0.885506 0.885508 0.464577 +6 -0.279395 0.960152 0.960154 0.279393 +6.2 -0.0830731 0.99652 0.996521 0.0830703 +6.3 0.016828 0.999835 0.999835 -0.0168267 -- 2.39.5