* - reinit_vector;
* - residual;
* - setup_jacobian;
- * - solve_jacobian_system/solve_with_jacobian;
+ * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
*
* The function `solve_jacobian_system` should be implemented for SUNDIALS
* < 4.0.0. For later versions, you should use
- * `solve_with_jacobian` to leverage better non-linear
+ * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear
* algorithms.
*
* Optionally, also the following functions could be provided. By default
*/
~IDA();
+ /**
+ * Save the number of iterations of the last Jacobian solve.
+ */
+ void
+ set_n_iterations(const int n_iter);
+
+ /**
+ * Return the number of iterations of the last Jacobian solve. This
+ * piece of information corresponds to the what the
+ * solve_jacobian_system_up_to_tolerance() function returns through
+ * its third argument.
+ */
+ int
+ get_n_iterations() const;
+
/**
* Integrate differential-algebraic equations. This function returns the
* final number of computed steps.
* 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() (for
- * SUNDIALS < 4.0.0) or solve_with_jacobian() (for
+ * SUNDIALS < 4.0.0) or solve_jacobian_system_up_to_tolerance() (for
* SUNDIALS >= 4.0.0).
*
* The Jacobian $J$ should be a (possibly inexact) computation of
* 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
+ * solve_jacobian_system_up_to_tolerance() 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
+ * solve_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
* system $J x = b$.
*
* This function should return:
* specifying the tolerance for the resolution. A part from the tolerance
* only `rhs` is provided and `dst` needs to be returned.
*/
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
DEAL_II_DEPRECATED_EARLY
+# endif
std::function<int(const VectorType &rhs, VectorType &dst)>
solve_jacobian_system;
* 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_with_jacobian() lead to better convergence in the
+ * solve_jacobian_system_up_to_tolerance() lead to better convergence in the
* Newton process.
*
* The Jacobian $J$ should be (an approximation of) the system Jacobian
*
* @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. This is an output argument through which the
+ * function can communicate how many iterations it took to solve the
+ * linear system, and that can then be queried from the outside using the
+ * get_n_iterations() function.
* @param[in] tolerance The tolerance with which to solve the linear system
* of equations.
*
* - <0: Unrecoverable error the computation will be aborted and an
* assertion will be thrown.
*/
- std::function<
- int(const VectorType &rhs, VectorType &dst, const double tolerance)>
- solve_with_jacobian;
+ std::function<int(const VectorType &rhs,
+ VectorType & dst,
+ int & n_iter,
+ const double tolerance)>
+ solve_jacobian_system_up_to_tolerance;
/**
* Process solution. This function is called by IDA at fixed time steps,
*/
void *ida_mem;
+ /**
+ * Number of iteration that were required to solve the last
+ * Jacobian system. This variable is set by the wrapper that calls
+ * `solve_jacobian_system_up_to_tolerance`.
+ */
+ int n_iterations;
+
/**
* MPI communicator. SUNDIALS solver runs happily in
* parallel. Note that if the library is compiled without MPI
auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
auto *dst_x = internal::unwrap_nvector<VectorType>(x);
int err = 0;
- if (solver.solve_with_jacobian)
- err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
+ if (solver.solve_jacobian_system_up_to_tolerance)
+ {
+ int n_iter = 0;
+
+ err = solver.solve_jacobian_system_up_to_tolerance(*src_b,
+ *dst_x,
+ n_iter,
+ tol);
+ solver.set_n_iterations(n_iter > 0 ? n_iter : 1);
+ }
else if (solver.solve_jacobian_system)
- err = solver.solve_jacobian_system(*src_b, *dst_x);
+ {
+ err = solver.solve_jacobian_system(*src_b, *dst_x);
+ }
else
// We have already checked this outside, so we should never get here.
Assert(false, ExcInternalError());
+ template <typename VectorType>
+ void
+ IDA<VectorType>::set_n_iterations(const int n_iter)
+ {
+ n_iterations = n_iter;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ IDA<VectorType>::get_n_iterations() const
+ {
+ return n_iterations;
+ }
+
+
+
template <typename VectorType>
unsigned int
IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
double h = data.initial_step_size;
unsigned int step_number = 0;
+ this->n_iterations = 1;
int status;
(void)status;
return 0;
};
- AssertThrow(solve_jacobian_system || solve_with_jacobian,
- ExcFunctionNotProvided(
- "solve_jacobian_system or solve_with_jacobian"));
+ AssertThrow(
+ solve_jacobian_system || solve_jacobian_system_up_to_tolerance,
+ ExcFunctionNotProvided(
+ "solve_jacobian_system or solve_jacobian_system_up_to_tolerance"));
LS->ops->solve = t_dae_solve_jacobian_system<VectorType>;
// When we set an iterative solver IDA requires that resid is provided. From
// 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 /*ignored*/) -> int { return 1; };
+ LS->ops->numiters = [](SUNLinearSolver LS) -> int {
+ IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
+ return solver.get_n_iterations();
+ };
// 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
projects / dealii.git / commitdiff