* - reinit_vector;
* - residual;
* - setup_jacobian;
- * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance;
+ * - solve_jacobian_system;
*
- * 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
+ * Optionally, also the following functions could be rewritten. 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_iterations(const int n_iter);
+ set_n_iter(const int n_iter)
+ {
+ this->n_iter = n_iter;
+ }
- /**
- * Return the number of iterations of the last Jacobian solve.
- */
int
- get_n_iterations() const;
-
+ get_n_iter() const
+ {
+ return this->n_iter;
+ }
/**
* 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() (for
- * SUNDIALS < 4.0.0) or solve_jacobian_system_up_to_tolerance() (for
- * SUNDIALS >= 4.0.0).
+ * is called by IDA once, before any call to solve_jacobian_system().
*
* 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() or
- * solve_jacobian_system_up_to_tolerance() can assume that this function has
- * been called at least once.
+ * Subsequent calls (possibly more than one) to solve_jacobian_system() 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_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the
- * system $J x = b$.
+ * setup_jacobian() it is possible to call solve_jacobian_system(), 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 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
+ * 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
VectorType & dst,
int & n_iter,
const double tolerance)>
- solve_jacobian_system_up_to_tolerance;
+ solve_with_jacobian;
/**
* Process solution. This function is called by IDA at fixed time steps,
/**
* IDA configuration data.
*/
- const AdditionalData data;
+ AdditionalData data;
/**
* IDA memory object.
void *ida_mem;
/**
- * Number of iteration that were required to solve the last Jacobian system
+ * Number of iteration required to solve the Jacobian system
*/
- int n_iterations;
+ int n_iter;
/**
* 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)
+ 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)
{
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;
- 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());
+ 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);
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_iterations = 1;
+ this->n_iter = 1;
int status;
(void)status;
return 0;
};
- 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>;
-
+ 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"));
+ }
// 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_iterations();
+ 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
GrowingVectorMemory<VectorType> mem;
typename VectorMemory<VectorType>::Pointer v(mem);
reinit_vector(*v);
- return v->locally_owned_elements();
+ const unsigned int size = v->size();
+ return complete_index_set(size);
};
}
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;
};
- // 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 {
+ 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;
return 0;
};
- time_stepper.solve_jacobian_system_up_to_tolerance =
- [&](const VectorType &src,
- VectorType & dst,
- int & n_iter,
- const double tolerance) -> int {
+ time_stepper.solve_with_jacobian = [&](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