--- /dev/null
+Removed: The deprecated member variables
+SUNDIALS::KINSOL::solve_jacobian_system,
+and SUNDIALS::IDA::solve_jacobian_system
+have been removed.
+<br>
+(Daniel Arndt, 2023/07/02)
/**
* 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
* - 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
* 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);
* };
* 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[
* 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".
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<void(const VectorType &rhs, VectorType &dst)>
- solve_jacobian_system;
-
/**
* Solve the Jacobian linear system up to a specified tolerance. This
* function will be called by IDA (possibly several times) after
*
* 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;
*
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<void(const VectorType &ycur,
- const VectorType &fcur,
- const VectorType &rhs,
- VectorType & dst)>
- 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
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,
auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
auto *dst_x = internal::unwrap_nvector<VectorType>(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
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:
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
const KINSOL<VectorType> &solver =
*static_cast<const KINSOL<VectorType> *>(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<VectorType>(b);
- auto dst_x = internal::unwrap_nvector<VectorType>(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<VectorType>(b);
+ auto dst_x = internal::unwrap_nvector<VectorType>(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<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_ycur(mem);
- typename VectorMemory<VectorType>::Pointer src_fcur(mem);
-
- auto src_b = internal::unwrap_nvector_const<VectorType>(b);
- auto dst_x = internal::unwrap_nvector<VectorType>(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
+++ /dev/null
-//-----------------------------------------------------------
-//
-// 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 <deal.II/base/parameter_handler.h>
-
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/sundials/arkode.h>
-
-#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<double>;
-
- ParameterHandler prm;
- SUNDIALS::ARKode<VectorType>::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<VectorType> 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<double> 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<double> y(3);
- y[0] = u0;
- y[1] = v0;
- y[2] = w0;
- ode.solve_ode(y);
-}
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);
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;
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;
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;
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);
};
- 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<double> 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<double> 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<double> 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<double> J_inverse(2, 2);
+ J_inverse.invert(J);
+
+ J_inverse.vmult(dst, rhs);
+ };
VectorType v(N);
v(0) = 0.5;
};
- 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;
};
- 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;
projects / dealii.git / commitdiff