#include <deal.II/base/mpi.h>
#ifdef DEAL_II_WITH_SUNDIALS
-# 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>
+# 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>
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>
-# 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
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-# include <deal.II/sundials/copy.h>
-# include <deal.II/sundials/n_vector.h>
+# 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/sundials/copy.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
-# endif
-# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-# include <deal.II/sundials/sunlinsol_newempty.h>
-# endif
-# include <iomanip>
-# include <iostream>
+
+# 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>;
-
- if (solve_jacobian_system)
- IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
- else
- AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
+ IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
# 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
+#endif
\ No newline at end of file
/**
- * Solve the Harmonic oscillator problem, using a direct solver for the
- * jacobian system.
+ * Solve the Harmonic oscillator problem.
*
* 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