From: Daniel Arndt Date: Sun, 2 Jul 2023 19:30:12 +0000 (-0400) Subject: Remove deprecated solve_jacobian_system X-Git-Tag: relicensing~795^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15596%2Fhead;p=dealii.git Remove deprecated solve_jacobian_system --- diff --git a/doc/news/changes/incompatibilities/20230702DanielArndt-2 b/doc/news/changes/incompatibilities/20230702DanielArndt-2 new file mode 100644 index 0000000000..d459dec102 --- /dev/null +++ b/doc/news/changes/incompatibilities/20230702DanielArndt-2 @@ -0,0 +1,6 @@ +Removed: The deprecated member variables +SUNDIALS::KINSOL::solve_jacobian_system, +and SUNDIALS::IDA::solve_jacobian_system +have been removed. +
+(Daniel Arndt, 2023/07/02) diff --git a/include/deal.II/numerics/nonlinear.h b/include/deal.II/numerics/nonlinear.h index cf4dc0c441..90b009991a 100644 --- a/include/deal.II/numerics/nonlinear.h +++ b/include/deal.II/numerics/nonlinear.h @@ -302,7 +302,7 @@ public: /** * A function object that users may supply and that is intended to * prepare the linear solver for subsequent calls to - * solve_jacobian_system(). + * solve_with_jacobian). * * The job of setup_jacobian() is to prepare the linear solver for * subsequent calls to solve_with_jacobian(), in the solution of linear diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index b038a458ad..d9ae292715 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -69,10 +69,7 @@ namespace SUNDIALS * - reinit_vector; * - residual; * - setup_jacobian; - * - solve_jacobian_system/solve_with_jacobian; - * - * The function solve_jacobian_system() is deprecated. You should use - * solve_with_jacobian() to leverage better non-linear algorithms. + * - solve_with_jacobian; * * Optionally, also the following functions could be provided. By default * they do nothing, or are not required. If you call the constructor in a way @@ -217,8 +214,8 @@ namespace SUNDIALS * Jinv.invert(J); * }; * - * time_stepper.solve_jacobian_system = [&](const VectorType &src, - * VectorType &dst) + * time_stepper.solve_with_jacobian_system = [&](const VectorType &src, + * VectorType &dst, double) * { * Jinv.vmult(dst,src); * }; @@ -665,8 +662,7 @@ namespace SUNDIALS * Compute Jacobian. This function is called by IDA any time a Jacobian * update is required. The user should compute the Jacobian (or update all * the variables that allow the application of the Jacobian). This function - * is called by IDA once, before any call to solve_jacobian_system() or - * solve_with_jacobian(). + * is called by IDA once, before any call to solve_with_jacobian(). * * The Jacobian $J$ should be a (possibly inexact) computation of * \f[ @@ -677,15 +673,13 @@ namespace SUNDIALS * If the user uses a matrix based computation of the Jacobian, then this * is the right place where an assembly routine should be called to * assemble both a matrix and a preconditioner for the Jacobian system. - * Subsequent calls (possibly more than one) to solve_jacobian_system() or - * solve_with_jacobian() can assume that this function has - * been called at least once. + * Subsequent calls (possibly more than one) to solve_with_jacobian() can + * assume that this function has been called at least once. * * Notice that no assumption is made by this interface on what the user * should do in this function. IDA only assumes that after a call to - * setup_jacobian() it is possible to call solve_jacobian_system() or - * solve_with_jacobian() to obtain a solution $x$ to the - * system $J x = b$. + * setup_jacobian() it is possible to call solve_with_jacobian() to obtain a + * solution $x$ to the system $J x = b$. * * @note This variable represents a * @ref GlossUserProvidedCallBack "user provided callback". @@ -700,41 +694,6 @@ namespace SUNDIALS const double alpha)> setup_jacobian; - /** - * 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] - * - * 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. - * - * @note This variable represents a - * @ref GlossUserProvidedCallBack "user provided callback". - * See there for a description of how to deal with errors and other - * requirements and conventions. In particular, IDA can deal - * with "recoverable" errors in some circumstances, so callbacks - * can throw exceptions of type RecoverableUserCallbackError. - * - * @deprecated Use solve_with_jacobian() instead which also uses a numerical - * tolerance. - */ - DEAL_II_DEPRECATED - std::function - solve_jacobian_system; - /** * Solve the Jacobian linear system up to a specified tolerance. This * function will be called by IDA (possibly several times) after diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index afb50ce1d0..35721ce801 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -168,7 +168,7 @@ namespace SUNDIALS * * If the use of a Newton or Picard method is desired, then the user should * also supply - * - solve_jacobian_system or solve_with_jacobian; + * - solve_with_jacobian; * and optionally * - setup_jacobian; * @@ -505,78 +505,6 @@ namespace SUNDIALS const VectorType ¤t_f)> setup_jacobian; - /** - * @deprecated Versions of SUNDIALS after 4.0 no longer provide all - * of the information necessary for this callback (see below). Use the - * `solve_with_jacobian` callback described below. - * - * A function object that users may supply and that is intended to solve - * a linear system with the Jacobian matrix. This function will be called by - * KINSOL (possibly several times) after setup_jacobian() has been called at - * least once. KINSOL tries to do its best to call setup_jacobian() the - * minimum number of times. If convergence can be achieved without updating - * the Jacobian, then KINSOL does not call setup_jacobian() again. If, on - * the contrary, internal KINSOL convergence tests fail, then KINSOL calls - * setup_jacobian() again with updated vectors and coefficients so that - * successive calls to solve_jacobian_system() lead to better convergence - * in the Newton process. - * - * If you do not specify a `solve_jacobian_system` or `solve_with_jacobian` - * function, then only a fixed point iteration strategy can be used. Notice - * that this may not converge, or may converge very slowly. - * - * A call to this function should store in `dst` the result of $J^{-1}$ - * applied to `rhs`, i.e., $J \cdot dst = rhs$. It is the user's - * responsibility to set up proper solvers and preconditioners inside this - * function (or in the `setup_jacobian` callback above). - * - * - * Arguments to the function are: - * - * @param[in] ycur The current $y$ vector for the current KINSOL - * internal step. In the documentation above, this $y$ vector is generally - * denoted by $u$. - * @param[in] fcur The current value of the implicit right-hand side at - * `ycur`, $f_I (t_n, ypred)$. - * @param[in] rhs The system right hand side to solve for - * @param[out] dst The solution of $J^{-1} * src$ - * - * This function should return: - * - 0: Success - * - >0: Recoverable error (KINSOL will try to change its internal - * parameters and attempt a new solution step) - * - <0: Unrecoverable error the computation will be aborted and an - * assertion will be thrown. - * - * @warning Starting with SUNDIALS 4.1, SUNDIALS no longer provides the - * `ycur` and `fcur` variables -- only `rhs` is provided and `dst` - * needs to be returned. The first two arguments will therefore be - * empty vectors if you use a SUNDIALS version newer than 4.1. - * In practice, that means that one - * can no longer compute a Jacobian matrix for the current iterate - * within this function. Rather, this has to happen inside the - * `setup_jacobian` function above that receives this information. - * If it is important that the Jacobian corresponds to the *current* - * iterate (rather than a re-used Jacobian matrix that had been - * computed in a previous iteration and that therefore corresponds - * to a *previous* iterate), then you will also have to set the - * AdditionalData::maximum_newton_step variable to one, indicating - * that the Jacobian should be re-computed in every iteration. - * - * @note This variable represents a - * @ref GlossUserProvidedCallBack "user provided callback". - * See there for a description of how to deal with errors and other - * requirements and conventions. In particular, KINSOL can deal - * with "recoverable" errors in some circumstances, so callbacks - * can throw exceptions of type RecoverableUserCallbackError. - */ - DEAL_II_DEPRECATED - std::function - solve_jacobian_system; - /** * A function object that users may supply and that is intended to solve * a linear system with the Jacobian matrix. This function will be called by diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 98b04474ca..d2bf1a7474 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -337,9 +337,8 @@ namespace SUNDIALS return 0; }; - AssertThrow(solve_jacobian_system || solve_with_jacobian, - ExcFunctionNotProvided( - "solve_jacobian_system or solve_with_jacobian")); + AssertThrow(solve_with_jacobian, + ExcFunctionNotProvided("solve_with_jacobian")); LS->ops->solve = [](SUNLinearSolver LS, SUNMatrix /*ignored*/, N_Vector x, @@ -349,25 +348,12 @@ namespace SUNDIALS auto *src_b = internal::unwrap_nvector_const(b); auto *dst_x = internal::unwrap_nvector(x); - if (solver.solve_with_jacobian) - return Utilities::call_and_possibly_capture_exception( - solver.solve_with_jacobian, - solver.pending_exception, - *src_b, - *dst_x, - tol); - else if (solver.solve_jacobian_system) - return Utilities::call_and_possibly_capture_exception( - solver.solve_jacobian_system, - solver.pending_exception, - *src_b, - *dst_x); - else - { - // We have already checked this outside, so we should never get here. - Assert(false, ExcInternalError()); - return -1; - } + return Utilities::call_and_possibly_capture_exception( + solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); }; // When we set an iterative solver IDA requires that resid is provided. From diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index e819db7c1a..fcee8cc703 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -195,9 +195,8 @@ namespace SUNDIALS else { Assert(residual, ExcFunctionNotProvided("residual")); - Assert(solve_jacobian_system || solve_with_jacobian, - ExcFunctionNotProvided( - "solve_jacobian_system || solve_with_jacobian")); + Assert(solve_with_jacobian, + ExcFunctionNotProvided("solve_with_jacobian")); } // Create a new solver object: @@ -326,9 +325,8 @@ namespace SUNDIALS SUNMatrix J = nullptr; SUNLinearSolver LS = nullptr; - if (solve_jacobian_system || - solve_with_jacobian) // user assigned a function - // object to the solver slot + // user assigned a function object to the solver slot + if (solve_with_jacobian) { // Set the operations we care for in the sun_linear_solver object // and attach it to the KINSOL object. The functions that will get @@ -373,52 +371,19 @@ namespace SUNDIALS const KINSOL &solver = *static_cast *>(LS->content); - // This is where we have to make a decision about which of the two - // signals to call. Let's first check the more modern one: - if (solver.solve_with_jacobian) - { - auto src_b = internal::unwrap_nvector_const(b); - auto dst_x = internal::unwrap_nvector(x); + Assert(solver.solve_with_jacobian, ExcInternalError()); - const int err = Utilities::call_and_possibly_capture_exception( - solver.solve_with_jacobian, - solver.pending_exception, - *src_b, - *dst_x, - tol); + auto src_b = internal::unwrap_nvector_const(b); + auto dst_x = internal::unwrap_nvector(x); - return err; - } - else - { - // User has not provided the modern callback, so the fact that - // we are here means that they must have given us something for - // the old signal. Check this. - Assert(solver.solve_jacobian_system, ExcInternalError()); - - // Allocate temporary (deal.II-type) dummy vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_ycur(mem); - typename VectorMemory::Pointer src_fcur(mem); - - auto src_b = internal::unwrap_nvector_const(b); - auto dst_x = internal::unwrap_nvector(x); - - // Call the user-provided setup function with these arguments. - // Note that Sundials 4.x and later no longer provide values for - // src_ycur and src_fcur, and so we simply pass dummy vector in. - // These vectors will have zero lengths because we don't reinit - // them above. - const int err = Utilities::call_and_possibly_capture_exception( - solver.solve_jacobian_system, - solver.pending_exception, - *src_ycur, - *src_fcur, - *src_b, - *dst_x); - - return err; - } + const int err = Utilities::call_and_possibly_capture_exception( + solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); + + return err; }; // Even though we don't use it, KINSOL still wants us to set some diff --git a/tests/sundials/arkode_05.cc b/tests/sundials/arkode_05.cc deleted file mode 100644 index 047b3be604..0000000000 --- a/tests/sundials/arkode_05.cc +++ /dev/null @@ -1,123 +0,0 @@ -//----------------------------------------------------------- -// -// Copyright (C) 2017 - 2023 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE.md at -// the top level directory of deal.II. -// -//----------------------------------------------------------- - -#include - -#include -#include - -#include - -#include "../tests.h" - - -// Test implicit-explicit time stepper. Both setup and solve_jacobian_system. -// Brusselator benchmark - -/** - * This test problem is called "brusselator", and is a typical benchmark for - * ODE solvers. This problem has 3 dependent variables u, v and w, that depend - * on the independent variable t via the IVP system - * - * du/dt = a - (w + 1)u + v u^2 - * dv/dt = w u - v u^2 - * dw/dt = (b - w)/eps - w u - * - * We integrate over the interval 0 <= t <= 10, with the initial conditions - * - * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8, - * - * and parameters - * - * a = 1.2, b = 2.5, and eps = 10^−5 - * - * The implicit part only contains the stiff part of the problem (the part with - * eps in right hand side of the third equation). - */ -int -main() -{ - initlog(); - - using VectorType = Vector; - - ParameterHandler prm; - SUNDIALS::ARKode::AdditionalData data; - data.add_parameters(prm); - - if (false) - { - std::ofstream ofile(SOURCE_DIR "/arkode_05_in.prm"); - prm.print_parameters(ofile, ParameterHandler::ShortText); - ofile.close(); - } - - std::ifstream ifile(SOURCE_DIR "/arkode_05_in.prm"); - prm.parse_input(ifile); - - SUNDIALS::ARKode ode(data); - - // Parameters - double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5; - // Explicit jacobian. - FullMatrix J(3, 3); - - ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { - ydot[0] = 0; - ydot[1] = 0; - ydot[2] = -y[2] / eps; - }; - - - ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { - ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; - ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; - ydot[2] = b / eps - y[2] * y[0]; - }; - - - ode.setup_jacobian = [&](const int, - const double, - const double gamma, - const VectorType &, - const VectorType &, - bool &j_is_current) { - J = 0; - J(0, 0) = 1; - J(1, 1) = 1; - J(2, 2) = 1 + gamma / eps; - J.gauss_jordan(); - j_is_current = true; - }; - - ode.solve_jacobian_system = [&](const double t, - const double gamma, - const VectorType &, - const VectorType &, - const VectorType &src, - VectorType & dst) { J.vmult(dst, src); }; - - ode.output_step = - [&](const double t, const VectorType &sol, const unsigned int step_number) { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] - << std::endl; - }; - - Vector y(3); - y[0] = u0; - y[1] = v0; - y[2] = w0; - ode.solve_ode(y); -} diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index d795143d0a..928fddb70e 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -96,11 +96,6 @@ public: Jinv.invert(J); }; - // Used only in ver < 4.0.0 - time_stepper.solve_jacobian_system = - [&](const VectorType &src, VectorType &dst) { Jinv.vmult(dst, src); }; - - // Used in ver >= 4.0.0 time_stepper.solve_with_jacobian = [&](const VectorType &src, VectorType &dst, const double) { Jinv.vmult(dst, src); diff --git a/tests/sundials/ida_03.cc b/tests/sundials/ida_03.cc index 9a5e301038..a2b576808d 100644 --- a/tests/sundials/ida_03.cc +++ b/tests/sundials/ida_03.cc @@ -74,11 +74,6 @@ public: J = kappa + alpha; }; - // Used only in ver < 4.0.0 - time_stepper.solve_jacobian_system = - [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; }; - - // Used in ver >= 4.0.0 time_stepper.solve_with_jacobian = [&](const VectorType &src, VectorType &dst, const double) { dst[0] = src[0] / J; diff --git a/tests/sundials/ida_04.cc b/tests/sundials/ida_04.cc index 8b95e5b5f0..0f00b4e75f 100644 --- a/tests/sundials/ida_04.cc +++ b/tests/sundials/ida_04.cc @@ -89,11 +89,6 @@ public: J = kappa + alpha; }; - // Used only in ver < 4.0.0 - time_stepper.solve_jacobian_system = - [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; }; - - // Used in ver >= 4.0.0 time_stepper.solve_with_jacobian = [&](const VectorType &src, VectorType &dst, const double) { dst[0] = src[0] / J; diff --git a/tests/sundials/ida_05.cc b/tests/sundials/ida_05.cc index 4dfb63da80..7b4d6c38d1 100644 --- a/tests/sundials/ida_05.cc +++ b/tests/sundials/ida_05.cc @@ -90,11 +90,6 @@ public: J = kappa + alpha; }; - // Used only in ver < 4.0.0 - time_stepper.solve_jacobian_system = - [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; }; - - // Used in ver >= 4.0.0 time_stepper.solve_with_jacobian = [&](const VectorType &src, VectorType &dst, const double) { dst[0] = src[0] / J; diff --git a/tests/sundials/kinsol_02.cc b/tests/sundials/kinsol_02.cc index 2b6fe97de7..3392c988f3 100644 --- a/tests/sundials/kinsol_02.cc +++ b/tests/sundials/kinsol_02.cc @@ -73,11 +73,6 @@ main() kinsol.solve_with_jacobian = [&](const VectorType &rhs, VectorType &dst, double) { dst = rhs; }; - kinsol.solve_jacobian_system = [&](const VectorType &, - const VectorType &, - const VectorType &rhs, - VectorType & dst) { dst = rhs; }; - VectorType v(N); auto niter = kinsol.solve(v); diff --git a/tests/sundials/kinsol_03.cc b/tests/sundials/kinsol_03.cc index f9f1564852..782c21a729 100644 --- a/tests/sundials/kinsol_03.cc +++ b/tests/sundials/kinsol_03.cc @@ -93,31 +93,29 @@ main() }; - kinsol.solve_jacobian_system = [](const VectorType &, - const VectorType &, - const VectorType &rhs, - VectorType & dst) { - deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] - << ')' << std::endl; - - // This isn't right for SUNDIALS >4.0: We don't actually get a valid - // 'u' vector, and so do the linearization of the problem around - // the zero vector. This *happens* to converge, but it isn't the - // right approach. Check the _04 test for a better approach. - VectorType u(2); - u[0] = u[1] = 0; - - FullMatrix J(2, 2); - J(0, 0) = -std::sin(u[0] + u[1]) + 2; - J(0, 1) = -std::sin(u[0] + u[1]); - J(1, 0) = std::cos(u[0] - u[1]); - J(1, 1) = -std::cos(u[0] - u[1]) + 2; - - FullMatrix J_inverse(2, 2); - J_inverse.invert(J); - - J_inverse.vmult(dst, rhs); - }; + kinsol.solve_with_jacobian = + [](const VectorType &rhs, VectorType &dst, double) { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ')' << std::endl; + + // This isn't right for SUNDIALS >4.0: We don't actually get a valid + // 'u' vector, and so do the linearization of the problem around + // the zero vector. This *happens* to converge, but it isn't the + // right approach. Check the _04 test for a better approach. + VectorType u(2); + u[0] = u[1] = 0; + + FullMatrix J(2, 2); + J(0, 0) = -std::sin(u[0] + u[1]) + 2; + J(0, 1) = -std::sin(u[0] + u[1]); + J(1, 0) = std::cos(u[0] - u[1]); + J(1, 1) = -std::cos(u[0] - u[1]) + 2; + + FullMatrix J_inverse(2, 2); + J_inverse.invert(J); + + J_inverse.vmult(dst, rhs); + }; VectorType v(N); v(0) = 0.5; diff --git a/tests/sundials/kinsol_04.cc b/tests/sundials/kinsol_04.cc index ed6bc6720c..aef6665761 100644 --- a/tests/sundials/kinsol_04.cc +++ b/tests/sundials/kinsol_04.cc @@ -113,15 +113,13 @@ main() }; - kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u, - const VectorType &, - const VectorType &rhs, - VectorType & dst) { - deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] - << ')' << std::endl; - - J_inverse.vmult(dst, rhs); - }; + kinsol.solve_with_jacobian = + [&J_inverse](const VectorType &rhs, VectorType &dst, double) { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ')' << std::endl; + + J_inverse.vmult(dst, rhs); + }; VectorType v(N); v(0) = 0.5; diff --git a/tests/sundials/kinsol_05.cc b/tests/sundials/kinsol_05.cc index 23b8ce1397..ebb00cafc9 100644 --- a/tests/sundials/kinsol_05.cc +++ b/tests/sundials/kinsol_05.cc @@ -113,15 +113,13 @@ main() }; - kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u, - const VectorType &, - const VectorType &rhs, - VectorType & dst) { - deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] - << ')' << std::endl; - - J_inverse.vmult(dst, rhs); - }; + kinsol.solve_with_jacobian = + [&J_inverse](const VectorType &rhs, VectorType &dst, double) { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ')' << std::endl; + + J_inverse.vmult(dst, rhs); + }; VectorType v(N); v(0) = 0.5;