* - reinit_vector;
* - residual;
* - setup_jacobian;
- * - solve_jacobian_system;
+ * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
*
- * Optionally, also the following functions could be rewritten. By default
+ * The function `solve_jacobian_system` should be implemented for SUNDIALS
+ * < 4.0.0. For later versions, you should use
+ * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear
+ * algorithms.
+ *
+ * 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
* that requires a not-implemented function, an Assertion will be
* thrown.
* @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();
-
+ /**
+ * Save the number of iterations of the last Jacobian solve.
+ */
void
- set_n_iter(const int n_iter)
- {
- this->n_iter = n_iter;
- }
+ set_n_iterations(const int n_iter);
+ /**
+ * Return the number of iterations of the last Jacobian solve.
+ */
int
- get_n_iter() const
- {
- return this->n_iter;
- }
+ get_n_iterations() const;
+
/**
* Integrate differential-algebraic equations. This function returns the
* final number of computed steps.
* 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().
+ * is called by IDA once, before any call to solve_jacobian_system() (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
* \f[
* If the user uses a matrix based computation of the Jacobian, than 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() can
- * assume that this function has been called at least once.
+ * Subsequent calls (possibly more than one) to solve_jacobian_system() or
+ * 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(), to
- * obtain a solution $x$ to the system $J x = b$.
+ * setup_jacobian() it is possible to call solve_jacobian_system() or
+ * solve_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
+ * system $J x = b$.
*
* This function should return:
* - 0: Success
* 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;
-
/**
- * 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
+ * Solve the Jacobian linear system up to a specified tolerance. 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_system_up_to_tolerance() lead to better convergence in the
* Newton process.
*
* The jacobian $J$ should be (an approximation of) the system Jacobian
VectorType & dst,
int & n_iter,
const double tolerance)>
- solve_with_jacobian;
+ solve_jacobian_system_up_to_tolerance;
/**
* Process solution. This function is called by IDA at fixed time steps,
/**
* IDA configuration data.
*/
- AdditionalData data;
+ const AdditionalData data;
/**
* IDA memory object.
void *ida_mem;
/**
- * Number of iteration required to solve the Jacobian system
+ * Number of iteration that were required to solve the last Jacobian system
*/
- int n_iter;
+ int n_iterations;
/**
* MPI communicator. SUNDIALS solver runs happily in
*/
MPI_Comm communicator;
-
-
/**
* Memory pool of vectors.
*/
# endif // PETSC_USE_COMPLEX
# endif // DEAL_II_WITH_PETSC
};
-
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
v->ops->nvdotprod = NVectorOperations::dot_product<VectorType>;
v->ops->nvmaxnorm = NVectorOperations::max_norm<VectorType>;
v->ops->nvwrmsnorm = NVectorOperations::weighted_rms_norm<VectorType>;
- v->ops->nvmin = NVectorOperations::min_element<VectorType>;
- v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm<VectorType>;
- v->ops->nvl1norm = NVectorOperations::l1_norm<VectorType>;
+ v->ops->nvmin = NVectorOperations::min_element<VectorType>;
+ v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm<VectorType>;
+ v->ops->nvl1norm = NVectorOperations::l1_norm<VectorType>;
v->ops->nvwrmsnormmask =
NVectorOperations::weighted_rms_norm_mask<VectorType>;
// v->ops->nvcompare = undef;
SUNMatrix /*ignored*/,
N_Vector x,
N_Vector b,
- realtype /*tol*/)
- {
- const IDA<VectorType> &solver =
- *static_cast<const IDA<VectorType> *>(LS->content);
-
- auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
- auto *dst_x = internal::unwrap_nvector<VectorType>(x);
-
- const int err = solver.solve_jacobian_system(*src_b, *dst_x);
-
- return err;
- }
-
-
- template <typename VectorType>
- int
- t_dae_solve_with_jacobian(SUNLinearSolver LS,
- SUNMatrix /*ignored*/,
- N_Vector x,
- N_Vector b,
- realtype tol)
+ realtype tol)
{
IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
auto *dst_x = internal::unwrap_nvector<VectorType>(x);
+ int err = 0;
+ if (solver.solve_jacobian_system_up_to_tolerance)
+ {
+ int n_iter = 0;
- 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);
+ 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);
+ }
+ else
+ // We have already checked this outside, so we should never get here.
+ Assert(false, ExcInternalError());
return err;
}
template <typename VectorType>
- IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm mpi_comm)
+ IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm)
: data(data)
, ida_mem(nullptr)
, communicator(is_serial_vector<VectorType>::value ?
+ 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_iter = 1;
+ this->n_iterations = 1;
int status;
(void)status;
return 0;
};
- if (solve_with_jacobian)
- {
- LS->ops->solve = t_dae_solve_with_jacobian<VectorType>;
- }
- else if (solve_jacobian_system)
- {
- LS->ops->solve = t_dae_solve_jacobian_system<VectorType>;
- }
- else
- {
- AssertThrow(false, ExcFunctionNotProvided("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
// SUNDIALS docs If an iterative method computes the preconditioned initial
// residual and returns with a successful solve without performing any
// has provided we set 1. This is clearly suboptimal.
LS->ops->numiters = [](SUNLinearSolver LS) -> int {
IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
- return solver.get_n_iter();
+ 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
GrowingVectorMemory<VectorType> mem;
typename VectorMemory<VectorType>::Pointer v(mem);
reinit_vector(*v);
- const unsigned int size = v->size();
- return complete_index_set(size);
+ return v->locally_owned_elements();
};
}
return 0;
};
-
+ // Used only in ver < 4.0.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 {
+ // Used in ver >= 4.0.0
+ time_stepper.solve_jacobian_system_up_to_tolerance =
+ [&](const VectorType &src,
+ VectorType & dst,
+ int & n_iter,
+ const double) -> int {
Jinv.vmult(dst, src);
n_iter = 1;
return 0;
return 0;
};
- time_stepper.solve_with_jacobian = [&](const VectorType &src,
- VectorType & dst,
- int & n_iter,
- const double tolerance) -> int {
+ time_stepper.solve_jacobian_system_up_to_tolerance =
+ [&](const VectorType &src,
+ VectorType & dst,
+ int & n_iter,
+ const double tolerance) -> int {
SolverControl solver_control(1000, tolerance);
SolverGMRES<Vector<double>> solver(solver_control);
solver.solve(J, dst, src, PreconditionIdentity());
projects / dealii.git / commitdiff