--- /dev/null
+New: SUNDIALS::IDA is now compatible with SUNDIALS > 4.0.0
+<br>
+(Nicola Giuliani, Luca Heltai, 2021/05/13)
#include <deal.II/base/mpi.h>
#ifdef DEAL_II_WITH_SUNDIALS
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-
-# include <deal.II/base/conditional_ostream.h>
-# include <deal.II/base/exceptions.h>
-# include <deal.II/base/logstream.h>
-# include <deal.II/base/parameter_handler.h>
-# ifdef DEAL_II_WITH_PETSC
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# endif
-# include <deal.II/lac/vector.h>
-# include <deal.II/lac/vector_memory.h>
-
-# ifdef DEAL_II_SUNDIALS_WITH_IDAS
-# include <idas/idas.h>
-# else
-# include <ida/ida.h>
-# endif
-
-# include <sundials/sundials_config.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
-# include <ida/ida_spbcgs.h>
-# include <ida/ida_spgmr.h>
-# include <ida/ida_sptfqmr.h>
-# endif
-# include <boost/signals2.hpp>
-
-# include <nvector/nvector_serial.h>
-# include <sundials/sundials_math.h>
-# include <sundials/sundials_types.h>
-
-# include <memory>
+# include <deal.II/base/conditional_ostream.h>
+# include <deal.II/base/exceptions.h>
+# include <deal.II/base/logstream.h>
+# include <deal.II/base/parameter_handler.h>
+# ifdef DEAL_II_WITH_PETSC
+# include <deal.II/lac/petsc_block_vector.h>
+# include <deal.II/lac/petsc_vector.h>
+# endif
+# include <deal.II/lac/vector.h>
+# include <deal.II/lac/vector_memory.h>
+
+# ifdef DEAL_II_SUNDIALS_WITH_IDAS
+# include <idas/idas.h>
+# else
+# include <ida/ida.h>
+# endif
+
+# include <sundials/sundials_config.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
+# include <ida/ida_spbcgs.h>
+# include <ida/ida_spgmr.h>
+# include <ida/ida_sptfqmr.h>
+# endif
+# include <deal.II/sundials/sunlinsol_wrapper.h>
+
+# include <boost/signals2.hpp>
+
+# include <nvector/nvector_serial.h>
+# include <sundials/sundials_math.h>
+# include <sundials/sundials_types.h>
+
+# include <memory>
DEAL_II_NAMESPACE_OPEN
// Shorthand notation for IDA error codes.
-# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
+# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code))
namespace SUNDIALS
{
* @param maximum_order Maximum BDF order
* @param maximum_non_linear_iterations Maximum number of nonlinear
* iterations
+ * @param ls_norm_factor Converting factor from the integrator tolerance
+ * to the linear solver tolerance
+ * iterations
*
* Error parameters:
*
const double minimum_step_size = 1e-6,
const unsigned int maximum_order = 5,
const unsigned int maximum_non_linear_iterations = 10,
+ const double ls_norm_factor = 0,
// Error parameters
const double absolute_tolerance = 1e-6,
const double relative_tolerance = 1e-5,
, reset_type(reset_type)
, maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic)
, maximum_non_linear_iterations(maximum_non_linear_iterations)
+ , ls_norm_factor(ls_norm_factor)
{}
/**
});
prm.add_parameter("Maximum number of nonlinear iterations",
maximum_non_linear_iterations_ic);
+ prm.add_parameter(
+ "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
+ ls_norm_factor);
prm.leave_subsection();
}
* Maximum number of iterations for Newton method during time advancement.
*/
unsigned int maximum_non_linear_iterations;
+
+ /**
+ * Factor to use when converting from the integrator tolerance to the
+ * linear solver tolerance
+ */
+ double ls_norm_factor;
};
/**
* @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();
+
+ void
+ set_n_iter(const int n_iter)
+ {
+ this->n_iter = n_iter;
+ }
+
+ int
+ get_n_iter() const
+ {
+ return this->n_iter;
+ }
/**
* Integrate differential-algebraic equations. This function returns the
* final number of computed steps.
* then last function will be attempted again
* - <0: Unrecoverable error the computation will be aborted and an
* assertion will be thrown.
+ *
+ * @warning Starting with SUNDIALS 4.1, SUNDIALS provides the possibility of
+ * specifying the tolerance for the resolution. A part from the tolerance
+ * only `rhs` is provided and `dst` needs to be returned.
*/
+ DEAL_II_DEPRECATED_EARLY
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
+ * 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]
+ *
+ * Arguments to the function are:
+ *
+ * @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
+ * @param[in] tolerance The tolerance with which to solve the linear system
+ * of equations.
+ *
+ * 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.
+ *
+ * This function should return:
+ * - 0: Success
+ * - >0: Recoverable error (IDAReinit will be called if this happens, and
+ * then last function will be attempted again
+ * - <0: Unrecoverable error the computation will be aborted and an
+ * assertion will be thrown.
+ */
+ std::function<int(const VectorType &rhs,
+ VectorType & dst,
+ int & n_iter,
+ const double tolerance)>
+ solve_with_jacobian;
+
/**
* Process solution. This function is called by IDA at fixed time steps,
* every `output_period` seconds, and it is passed a polynomial
*/
N_Vector diff_id;
+ /**
+ * Number of iteration required to solve the Jacobian system
+ */
+ int n_iter;
+
+
/**
* MPI communicator. SUNDIALS solver runs happily in
* parallel. Note that if the library is compiled without MPI
*/
MPI_Comm communicator;
+
+
/**
* Memory pool of vectors.
*/
GrowingVectorMemory<VectorType> mem;
-# ifdef DEAL_II_WITH_PETSC
-# ifdef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifdef PETSC_USE_COMPLEX
static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
"Sundials does not support complex scalar types, "
"but PETSc is configured to use a complex scalar type!");
!std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
"Sundials does not support complex scalar types, "
"but PETSc is configured to use a complex scalar type!");
-# endif // PETSC_USE_COMPLEX
-# endif // DEAL_II_WITH_PETSC
+# endif // PETSC_USE_COMPLEX
+# endif // DEAL_II_WITH_PETSC
};
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
-# endif
#endif // DEAL_II_WITH_SUNDIALS
#endif
// the top level directory of deal.II.
//
//-----------------------------------------------------------
-
-
#include <deal.II/base/config.h>
#include <deal.II/sundials/ida.h>
#ifdef DEAL_II_WITH_SUNDIALS
+# include <deal.II/base/utilities.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-
-# include <deal.II/base/utilities.h>
-
-# include <deal.II/lac/block_vector.h>
-# ifdef DEAL_II_WITH_TRILINOS
-# include <deal.II/lac/trilinos_parallel_block_vector.h>
-# include <deal.II/lac/trilinos_vector.h>
-# endif
-# ifdef DEAL_II_WITH_PETSC
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# endif
+# include <deal.II/lac/block_vector.h>
+# ifdef DEAL_II_WITH_TRILINOS
+# include <deal.II/lac/trilinos_parallel_block_vector.h>
+# include <deal.II/lac/trilinos_vector.h>
+# endif
+# ifdef DEAL_II_WITH_PETSC
+# include <deal.II/lac/petsc_block_vector.h>
+# include <deal.II/lac/petsc_vector.h>
+# endif
-# include <deal.II/sundials/copy.h>
+# include <deal.II/sundials/copy.h>
+# include <deal.II/sundials/n_vector.h>
+# include <sundials/sundials_config.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
# ifdef DEAL_II_SUNDIALS_WITH_IDAS
# include <idas/idas_impl.h>
# else
# include <ida/ida_impl.h>
# endif
-
-# include <iomanip>
-# include <iostream>
+# endif
+# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
+# include <deal.II/sundials/sunlinsol_newempty.h>
+# endif
+# include <iomanip>
+# include <iostream>
DEAL_II_NAMESPACE_OPEN
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
template <typename VectorType>
int
t_dae_lsetup(IDAMem IDA_mem,
}
+
template <typename VectorType>
int
t_dae_solve(IDAMem IDA_mem,
return err;
}
+
+
+# else
+ template <typename VectorType>
+ int
+ t_dae_jacobian_setup(realtype tt,
+ realtype cj,
+ N_Vector yy,
+ N_Vector yp,
+ N_Vector /* residual */,
+ SUNMatrix /* ignored */,
+ void *user_data,
+ N_Vector /* tmp1 */,
+ N_Vector /* tmp2 */,
+ N_Vector /* tmp3 */)
+ {
+ Assert(user_data != nullptr, ExcInternalError());
+ IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
+ GrowingVectorMemory<VectorType> mem;
+
+ typename VectorMemory<VectorType>::Pointer src_yy(mem);
+ solver.reinit_vector(*src_yy);
+
+ typename VectorMemory<VectorType>::Pointer src_yp(mem);
+ solver.reinit_vector(*src_yp);
+
+ copy(*src_yy, yy);
+ copy(*src_yp, yp);
+
+
+ int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj);
+
+
+ return err;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ t_dae_solve_jacobian_system(SUNLinearSolver LS,
+ SUNMatrix /*ignored*/,
+ N_Vector x,
+ N_Vector b,
+ realtype /*tol*/)
+ {
+ const IDA<VectorType> &solver =
+ *static_cast<const IDA<VectorType> *>(LS->content);
+
+ // Allocate temporary (deal.II-type) vectors into which to copy the
+ // N_vectors
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_b(mem);
+ typename VectorMemory<VectorType>::Pointer dst_x(mem);
+
+ solver.reinit_vector(*src_b);
+ solver.reinit_vector(*dst_x);
+
+ copy(*src_b, b);
+
+ const int err = solver.solve_jacobian_system(*src_b, *dst_x);
+
+ copy(x, *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);
+
+ // Allocate temporary (deal.II-type) vectors into which to copy the
+ // N_vectors
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_b(mem);
+ typename VectorMemory<VectorType>::Pointer dst_x(mem);
+
+ solver.reinit_vector(*src_b);
+ solver.reinit_vector(*dst_x);
+
+ copy(*src_b, b);
+ 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);
+ copy(x, *dst_x);
+
+ return err;
+ }
+# endif
} // namespace
+
+
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)
, yy(nullptr)
{
if (ida_mem)
IDAFree(&ida_mem);
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const int ierr = MPI_Comm_free(&communicator);
(void)ierr;
AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
}
-# endif
+# endif
}
double h = data.initial_step_size;
unsigned int step_number = 0;
+ this->n_iter = 1;
int status;
(void)status;
// The solution is stored in
// solution. Here we take only a
// view of it.
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const IndexSet is = solution.locally_owned_elements();
N_VNew_Parallel(communicator, local_system_size, system_size);
}
else
-# endif
+# endif
{
Assert(is_serial_vector<VectorType>::value,
ExcInternalError(
}
// Free the vectors which are no longer used.
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
N_VDestroy_Parallel(yy);
N_VDestroy_Parallel(diff_id);
}
else
-# endif
+# endif
{
N_VDestroy_Serial(yy);
N_VDestroy_Serial(yp);
return step_number;
}
+
+
template <typename VectorType>
void
IDA<VectorType>::reset(const double current_time,
ida_mem = IDACreate();
-
// Free the vectors which are no longer used.
if (yy)
{
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
N_VDestroy_Parallel(yy);
N_VDestroy_Parallel(diff_id);
}
else
-# endif
+# endif
{
N_VDestroy_Serial(yy);
N_VDestroy_Serial(yp);
int status;
(void)status;
system_size = solution.size();
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const IndexSet is = solution.locally_owned_elements();
N_VNew_Parallel(communicator, local_system_size, system_size);
}
else
-# endif
+# endif
{
yy = N_VNew_Serial(system_size);
yp = N_VNew_Serial(system_size);
status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
AssertIDA(status);
-
if (get_local_tolerances)
{
copy(abs_tolls, get_local_tolerances());
AssertIDA(status);
// Initialize solver
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
auto IDA_mem = static_cast<IDAMem>(ida_mem);
IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
- IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
+
+ if (solve_jacobian_system)
+ IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
+ else
+ AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
IDA_mem->ida_setupNonNull = true;
# endif
+# else
+ SUNMatrix J = nullptr;
+ SUNLinearSolver LS = nullptr;
+
+ // and attach it to the SUNLinSol object. The functions that will get
+ // called do not actually receive the IDAMEM object, just the LS
+ // object, so we have to store a pointer to the current
+ // object in the LS object
+ LS = SUNLinSolNewEmpty();
+ LS->content = this;
+
+ LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
+ return SUNLINEARSOLVER_MATRIX_ITERATIVE;
+ };
+
+ LS->ops->free = [](SUNLinearSolver LS) -> int {
+ if (LS->content)
+ {
+ LS->content = nullptr;
+ }
+ if (LS->ops)
+ {
+ free(LS->ops);
+ LS->ops = nullptr;
+ }
+ free(LS);
+ LS = nullptr;
+ 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"));
+ }
+ // 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
+ // iterations (i.e., either the initial guess or the preconditioner is
+ // sufficiently accurate), then this optional routine may be called by the
+ // SUNDIALS package. This routine should return the N_Vector containing the
+ // preconditioned initial residual vector.
+ LS->ops->resid = [](SUNLinearSolver /*ignored*/) -> N_Vector {
+ return nullptr;
+ };
+ // 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 LS) -> int {
+ IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
+ 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
+ // if we don't set it, it won't call the functions that set up
+ // the matrix object (i.e., the argument to the 'IDASetJacFn'
+ // function below).
+ J = SUNMatNewEmpty();
+ J->content = this;
+
+ J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
+ return SUNMATRIX_CUSTOM;
+ };
+
+ J->ops->destroy = [](SUNMatrix A) {
+ if (A->content)
+ {
+ A->content = nullptr;
+ }
+ if (A->ops)
+ {
+ free(A->ops);
+ A->ops = nullptr;
+ }
+ free(A);
+ A = nullptr;
+ };
+
+ // Now set the linear system and Jacobian objects in the solver:
+ status = IDASetLinearSolver(ida_mem, LS, J);
+ AssertIDA(status);
+
+ status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
+ AssertIDA(status);
+ // Finally tell IDA about
+ // it as well. The manual says that this must happen *after*
+ // calling IDASetLinearSolver
+ status = IDASetJacFn(ida_mem, &t_dae_jacobian_setup<VectorType>);
+ AssertIDA(status);
+# endif
status = IDASetMaxOrd(ida_mem, data.maximum_order);
AssertIDA(status);
return ret;
};
- setup_jacobian = [](const double,
- const VectorType &,
- const VectorType &,
- const double) -> int {
- int ret = 0;
- AssertThrow(false, ExcFunctionNotProvided("setup_jacobian"));
- return ret;
- };
-
- solve_jacobian_system = [](const VectorType &, VectorType &) -> int {
- int ret = 0;
- AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
- return ret;
- };
output_step = [](const double,
const VectorType &,
template class IDA<Vector<double>>;
template class IDA<BlockVector<double>>;
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
-# ifdef DEAL_II_WITH_TRILINOS
+# ifdef DEAL_II_WITH_TRILINOS
template class IDA<TrilinosWrappers::MPI::Vector>;
template class IDA<TrilinosWrappers::MPI::BlockVector>;
-# endif // DEAL_II_WITH_TRILINOS
+# endif // DEAL_II_WITH_TRILINOS
-# ifdef DEAL_II_WITH_PETSC
-# ifndef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifndef PETSC_USE_COMPLEX
template class IDA<PETScWrappers::MPI::Vector>;
template class IDA<PETScWrappers::MPI::BlockVector>;
-# endif // PETSC_USE_COMPLEX
-# endif // DEAL_II_WITH_PETSC
+# endif // PETSC_USE_COMPLEX
+# endif // DEAL_II_WITH_PETSC
-# endif // DEAL_II_WITH_MPI
+# endif // DEAL_II_WITH_MPI
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
-# endif
-
#endif // DEAL_II_WITH_SUNDIALS
template class SUNDIALS::internal::NVectorView<PETScWrappers::MPI::V>;
template class SUNDIALS::internal::NVectorView<const PETScWrappers::MPI::V>;
}
-#endif
\ No newline at end of file
+#endif
/**
- * Solve the Harmonic oscillator problem.
+ * Solve the Harmonic oscillator problem, using a direct solver for the
+ * jacobian system.
*
* u'' = -k^2 u
* u (0) = 0
return 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 {
+ Jinv.vmult(dst, src);
+ n_iter = 1;
+ return 0;
+ };
+
time_stepper.output_step = [&](const double t,
const VectorType & sol,
const VectorType & sol_dot,
--- /dev/null
+0 0 1 1 0
+0.2 0.198666 0.98005 0.980056 -0.198635
+0.4 0.389414 0.92105 0.921063 -0.3894
+0.6 0.564639 0.825325 0.825346 -0.564653
+0.8 0.717354 0.696694 0.696713 -0.717367
+1 0.841469 0.540288 0.540295 -0.841478
+1.2 0.932034 0.362343 0.362339 -0.932034
+1.4 0.985441 0.169954 0.169952 -0.985442
+1.6 0.999562 -0.0292108 -0.0292135 -0.999563
+1.8 0.973833 -0.227211 -0.227214 -0.973833
+2 0.909281 -0.416152 -0.416155 -0.909281
+2.2 0.808479 -0.588503 -0.588506 -0.808478
+2.4 0.675445 -0.737392 -0.737394 -0.675443
+2.6 0.515484 -0.856883 -0.856885 -0.515482
+2.8 0.334972 -0.942213 -0.942214 -0.334969
+3 0.141106 -0.989979 -0.98998 -0.141103
+3.2 -0.0583858 -0.998279 -0.998279 0.0583887
+3.4 -0.255549 -0.966779 -0.966779 0.255552
+3.6 -0.442525 -0.896738 -0.896737 0.442528
+3.8 -0.611858 -0.790946 -0.790945 0.61186
+4 -0.756798 -0.653622 -0.65362 0.7568
+4.2 -0.871566 -0.49024 -0.490238 0.871568
+4.4 -0.951588 -0.307314 -0.307312 0.95159
+4.6 -0.993674 -0.112137 -0.112134 0.993674
+4.8 -0.996144 0.0875108 0.0875137 0.996144
+5 -0.958901 0.28367 0.283673 0.958901
+5.2 -0.88343 0.468519 0.468522 0.883429
+5.4 -0.772739 0.63469 0.634693 0.772738
+5.6 -0.631242 0.775558 0.77556 0.63124
+5.8 -0.464579 0.885506 0.885508 0.464577
+6 -0.279395 0.960152 0.960154 0.279393
+6.2 -0.0830731 0.99652 0.996521 0.0830703
+6.3 0.016828 0.999835 0.999835 -0.0168267
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2017 - 2021 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/precondition.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/ida.h>
+
+#include "../tests.h"
+
+
+/**
+ * Solve the Harmonic oscillator problem, using an iterative solver for the
+ * jacobian system (with the interface that is available only in version >
+ * 4.0.0)
+ *
+ * u'' = -k^2 u
+ * u (0) = 0
+ * u'(0) = k
+ *
+ * write in terms of a first order ode:
+ *
+ * y[0]' - y[1] = 0
+ * y[1]' + k^2 y[0] = 0
+ *
+ * That is
+ *
+ * F(y', y, t) = y' + A y = 0
+ *
+ * A = [ 0 , -1; k^2, 0 ]
+ *
+ * y_0 = 0, k
+ * y_0' = k, 0
+ *
+ * The exact solution is
+ *
+ * y[0](t) = sin(k t)
+ * y[1](t) = k cos(k t)
+ *
+ * The Jacobian to assemble is the following:
+ *
+ * J = alpha I + A
+ */
+class HarmonicOscillator
+{
+public:
+ HarmonicOscillator(
+ double _kappa,
+ const typename SUNDIALS::IDA<Vector<double>>::AdditionalData &data)
+ : time_stepper(data)
+ , y(2)
+ , y_dot(2)
+ , J(2, 2)
+ , A(2, 2)
+ , Jinv(2, 2)
+ , kappa(_kappa)
+ , out("output")
+ {
+ using VectorType = Vector<double>;
+
+ time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(2); };
+
+
+ time_stepper.residual = [&](const double t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ VectorType & res) -> int {
+ res = y_dot;
+ A.vmult_add(res, y);
+ return 0;
+ };
+
+ time_stepper.setup_jacobian = [&](const double,
+ const VectorType &,
+ const VectorType &,
+ const double alpha) -> int {
+ A(0, 1) = -1.0;
+ A(1, 0) = kappa * kappa;
+
+ J = A;
+
+ J(0, 0) = alpha;
+ J(1, 1) = alpha;
+
+ return 0;
+ };
+
+ 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());
+ n_iter = solver_control.last_step() > 0 ? solver_control.last_step() : 1;
+ return 0;
+ };
+
+ time_stepper.output_step = [&](const double t,
+ const VectorType & sol,
+ const VectorType & sol_dot,
+ const unsigned int step_number) -> int {
+ out << t << " " << sol[0] << " " << sol[1] << " " << sol_dot[0] << " "
+ << sol_dot[1] << std::endl;
+ return 0;
+ };
+ }
+
+ void
+ run()
+ {
+ y[1] = kappa;
+ y_dot[0] = kappa;
+ time_stepper.solve_dae(y, y_dot);
+ }
+ SUNDIALS::IDA<Vector<double>> time_stepper;
+
+private:
+ Vector<double> y;
+ Vector<double> y_dot;
+ FullMatrix<double> J;
+ FullMatrix<double> A;
+ FullMatrix<double> Jinv;
+ double kappa;
+
+ std::ofstream out;
+};
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, numbers::invalid_unsigned_int);
+
+ SUNDIALS::IDA<Vector<double>>::AdditionalData data;
+ ParameterHandler prm;
+ data.add_parameters(prm);
+
+ // std::ofstream ofile(SOURCE_DIR "/ida_01.prm");
+ // prm.print_parameters(ofile, ParameterHandler::ShortText);
+ // ofile.close();
+
+ std::ifstream ifile(SOURCE_DIR "/ida_01.prm");
+ prm.parse_input(ifile);
+
+
+ HarmonicOscillator ode(1.0, data);
+ ode.run();
+ return 0;
+}
--- /dev/null
+0 0 1 1 0
+0.2 0.198666 0.98005 0.980056 -0.198635
+0.4 0.389414 0.92105 0.921063 -0.3894
+0.6 0.564639 0.825325 0.825346 -0.564653
+0.8 0.717354 0.696694 0.696713 -0.717367
+1 0.841469 0.540288 0.540295 -0.841478
+1.2 0.932034 0.362343 0.362339 -0.932034
+1.4 0.985441 0.169954 0.169952 -0.985442
+1.6 0.999562 -0.0292108 -0.0292135 -0.999563
+1.8 0.973833 -0.227211 -0.227214 -0.973833
+2 0.909281 -0.416152 -0.416155 -0.909281
+2.2 0.808479 -0.588503 -0.588506 -0.808478
+2.4 0.675445 -0.737392 -0.737394 -0.675443
+2.6 0.515484 -0.856883 -0.856885 -0.515482
+2.8 0.334972 -0.942213 -0.942214 -0.334969
+3 0.141106 -0.989979 -0.98998 -0.141103
+3.2 -0.0583858 -0.998279 -0.998279 0.0583887
+3.4 -0.255549 -0.966779 -0.966779 0.255552
+3.6 -0.442525 -0.896738 -0.896737 0.442528
+3.8 -0.611858 -0.790946 -0.790945 0.61186
+4 -0.756798 -0.653622 -0.65362 0.7568
+4.2 -0.871566 -0.49024 -0.490238 0.871568
+4.4 -0.951588 -0.307314 -0.307312 0.95159
+4.6 -0.993674 -0.112137 -0.112134 0.993674
+4.8 -0.996144 0.0875108 0.0875137 0.996144
+5 -0.958901 0.28367 0.283673 0.958901
+5.2 -0.88343 0.468519 0.468522 0.883429
+5.4 -0.772739 0.63469 0.634693 0.772738
+5.6 -0.631242 0.775558 0.77556 0.63124
+5.8 -0.464579 0.885506 0.885508 0.464577
+6 -0.279395 0.960152 0.960154 0.279393
+6.2 -0.0830731 0.99652 0.996521 0.0830703
+6.3 0.016828 0.999835 0.999835 -0.0168267
projects / dealii.git / commitdiff