From 9e2c24aacd8298284bd598f7daf881edcbf301d6 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 5 Nov 2021 11:29:44 -0400 Subject: [PATCH] Revert "Make IDA work with new sundials." This reverts commit 67c15c819d1d2f23c58c0fb7a3753dd0a917def6. --- 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 -- 7 files changed, 98 insertions(+), 611 deletions(-) delete mode 100644 tests/sundials/ida_01.with_sundials.geq.4.0.0.output delete mode 100644 tests/sundials/ida_02.cc delete mode 100644 tests/sundials/ida_02.with_sundials.geq.4.0.0.output diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 43d52e16af..7ff57dac50 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 -# 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 +# 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 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,9 +286,6 @@ 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: * @@ -313,7 +310,6 @@ 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, @@ -335,7 +331,6 @@ 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) {} /** @@ -457,9 +452,6 @@ 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(); } @@ -545,12 +537,6 @@ 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; }; /** @@ -595,25 +581,13 @@ 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. @@ -729,59 +703,10 @@ 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 @@ -903,12 +828,6 @@ 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 @@ -916,15 +835,13 @@ 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!"); @@ -933,14 +850,15 @@ 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 6788c62fca..1fe5fc9061 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -12,39 +12,38 @@ // the top level directory of deal.II. // //----------------------------------------------------------- + + #include #include #ifdef DEAL_II_WITH_SUNDIALS -# include -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# include +# include + +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif + +# include -# include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) # ifdef DEAL_II_SUNDIALS_WITH_IDAS # include # else # include # endif -# endif -# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) -# include -# endif -# include -# include + +# include +# include DEAL_II_NAMESPACE_OPEN @@ -86,7 +85,6 @@ namespace SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template int t_dae_lsetup(IDAMem IDA_mem, @@ -123,7 +121,6 @@ namespace SUNDIALS } - template int t_dae_solve(IDAMem IDA_mem, @@ -155,108 +152,10 @@ 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) @@ -275,14 +174,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 } @@ -297,14 +196,13 @@ 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(); @@ -320,7 +218,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { Assert(is_serial_vector::value, ExcInternalError( @@ -358,7 +256,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); @@ -367,7 +265,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -378,8 +276,6 @@ namespace SUNDIALS return step_number; } - - template void IDA::reset(const double current_time, @@ -395,10 +291,11 @@ 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); @@ -407,7 +304,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -419,7 +316,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(); @@ -435,7 +332,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); @@ -448,6 +345,7 @@ 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()); @@ -494,115 +392,14 @@ 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; - - if (solve_jacobian_system) - IDA_mem->ida_lsolve = t_dae_solve; - else - AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system")); + IDA_mem->ida_lsolve = t_dae_solve; # 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); @@ -660,6 +457,20 @@ 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 &, @@ -681,24 +492,26 @@ 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 899195b795..4c11264684 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 +#endif \ No newline at end of file diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index a44c679751..239424b53f 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -24,8 +24,7 @@ /** - * Solve the Harmonic oscillator problem, using a direct solver for the - * jacobian system. + * Solve the Harmonic oscillator problem. * * u'' = -k^2 u * u (0) = 0 @@ -99,22 +98,12 @@ 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 deleted file mode 100644 index 7238b14d99..0000000000 --- a/tests/sundials/ida_01.with_sundials.geq.4.0.0.output +++ /dev/null @@ -1,33 +0,0 @@ -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 deleted file mode 100644 index 2da335b998..0000000000 --- a/tests/sundials/ida_02.cc +++ /dev/null @@ -1,167 +0,0 @@ -//----------------------------------------------------------- -// -// 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 deleted file mode 100644 index 7238b14d99..0000000000 --- a/tests/sundials/ida_02.with_sundials.geq.4.0.0.output +++ /dev/null @@ -1,33 +0,0 @@ -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