17 #include <deal.II/base/config.h> 19 #include <deal.II/sundials/kinsol.h> 21 #ifdef DEAL_II_WITH_SUNDIALS 23 # include <deal.II/base/utilities.h> 25 # include <deal.II/lac/block_vector.h> 26 # ifdef DEAL_II_WITH_TRILINOS 27 # include <deal.II/lac/trilinos_parallel_block_vector.h> 28 # include <deal.II/lac/trilinos_vector.h> 30 # ifdef DEAL_II_WITH_PETSC 31 # include <deal.II/lac/petsc_block_vector.h> 32 # include <deal.II/lac/petsc_vector.h> 34 # include <deal.II/base/utilities.h> 36 # include <deal.II/sundials/copy.h> 38 # include <sundials/sundials_config.h> 39 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) 40 # include <kinsol/kinsol_direct.h> 41 # include <sunlinsol/sunlinsol_dense.h> 42 # include <sunmatrix/sunmatrix_dense.h> 44 # include <kinsol/kinsol_dense.h> 50 DEAL_II_NAMESPACE_OPEN
58 template <
typename VectorType>
60 t_kinsol_function(N_Vector yy, N_Vector FF,
void *user_data)
62 KINSOL<VectorType> &solver =
63 *
static_cast<KINSOL<VectorType> *
>(user_data);
67 solver.reinit_vector(*src_yy);
70 solver.reinit_vector(*dst_FF);
76 err = solver.residual(*src_yy, *dst_FF);
77 else if (solver.iteration_function)
78 err = solver.iteration_function(*src_yy, *dst_FF);
89 template <
typename VectorType>
91 t_kinsol_setup_jacobian(KINMem kinsol_mem)
93 KINSOL<VectorType> &solver =
94 *
static_cast<KINSOL<VectorType> *
>(kinsol_mem->kin_user_data);
98 solver.reinit_vector(*src_ycur);
101 solver.reinit_vector(*src_fcur);
103 copy(*src_ycur, kinsol_mem->kin_uu);
104 copy(*src_fcur, kinsol_mem->kin_fval);
106 int err = solver.setup_jacobian(*src_ycur, *src_fcur);
112 template <
typename VectorType>
114 t_kinsol_solve_jacobian(KINMem kinsol_mem,
120 KINSOL<VectorType> &solver =
121 *
static_cast<KINSOL<VectorType> *
>(kinsol_mem->kin_user_data);
125 solver.reinit_vector(*src_ycur);
128 solver.reinit_vector(*src_fcur);
130 copy(*src_ycur, kinsol_mem->kin_uu);
131 copy(*src_fcur, kinsol_mem->kin_fval);
134 solver.reinit_vector(*src);
137 solver.reinit_vector(*dst);
141 int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst);
144 *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale);
145 N_VProd(b, kinsol_mem->kin_fscale, b);
146 N_VProd(b, kinsol_mem->kin_fscale, b);
147 *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b);
153 template <
typename VectorType>
155 const MPI_Comm mpi_comm)
157 , kinsol_mem(nullptr)
163 Utilities::MPI::duplicate_communicator(mpi_comm))
170 template <
typename VectorType>
174 KINFree(&kinsol_mem);
175 # ifdef DEAL_II_WITH_MPI 187 template <
typename VectorType>
191 unsigned int system_size = initial_guess_and_solution.size();
196 # ifdef DEAL_II_WITH_MPI 199 const IndexSet is = initial_guess_and_solution.locally_owned_elements();
200 const unsigned int local_system_size = is.
n_elements();
203 N_VNew_Parallel(
communicator, local_system_size, system_size);
206 N_VConst_Parallel(1.e0,
u_scale);
209 N_VConst_Parallel(1.e0,
f_scale);
216 "Trying to use a serial code with a parallel vector."));
217 solution = N_VNew_Serial(system_size);
218 u_scale = N_VNew_Serial(system_size);
219 N_VConst_Serial(1.e0,
u_scale);
220 f_scale = N_VNew_Serial(system_size);
221 N_VConst_Serial(1.e0,
f_scale);
230 copy(
solution, initial_guess_and_solution);
233 KINFree(&kinsol_mem);
235 kinsol_mem = KINCreate();
237 int status = KINInit(kinsol_mem, t_kinsol_function<VectorType>,
solution);
239 AssertKINSOL(status);
241 status = KINSetUserData(kinsol_mem, (
void *)
this);
242 AssertKINSOL(status);
245 AssertKINSOL(status);
248 AssertKINSOL(status);
251 AssertKINSOL(status);
254 AssertKINSOL(status);
257 AssertKINSOL(status);
260 AssertKINSOL(status);
263 AssertKINSOL(status);
266 AssertKINSOL(status);
269 AssertKINSOL(status);
271 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) 272 SUNMatrix J =
nullptr;
273 SUNLinearSolver LS =
nullptr;
278 KINMem KIN_mem = (KINMem)kinsol_mem;
279 KIN_mem->kin_lsolve = t_kinsol_solve_jacobian<VectorType>;
282 KIN_mem->kin_lsetup = t_kinsol_setup_jacobian<VectorType>;
283 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) 284 KIN_mem->kin_setupNonNull =
true;
290 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) 291 J = SUNDenseMatrix(system_size, system_size);
292 LS = SUNDenseLinearSolver(
u_scale, J);
293 status = KINDlsSetLinearSolver(kinsol_mem, LS, J);
295 status = KINDense(kinsol_mem, system_size);
297 AssertKINSOL(status);
310 AssertKINSOL(status);
312 copy(initial_guess_and_solution,
solution);
315 # ifdef DEAL_II_WITH_MPI 331 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
332 AssertKINSOL(status);
334 # if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) 338 KINFree(&kinsol_mem);
340 return (
unsigned int)nniters;
343 template <
typename VectorType>
355 # ifdef DEAL_II_WITH_MPI 357 # ifdef DEAL_II_WITH_TRILINOS 362 # ifdef DEAL_II_WITH_PETSC 363 # ifndef PETSC_USE_COMPLEX 373 DEAL_II_NAMESPACE_CLOSE
#define AssertNothrow(cond, exc)
unsigned int maximum_non_linear_iterations
std::function< int(const VectorType &src, VectorType &dst)> residual
unsigned int maximum_setup_calls
size_type n_elements() const
double function_tolerance
std::function< int(const VectorType &src, VectorType &dst)> iteration_function
void set_functions_to_trigger_an_assert()
#define AssertThrow(cond, exc)
SolutionStrategy strategy
std::function< VectorType &()> get_function_scaling
unsigned int anderson_subspace_size
double maximum_newton_step
#define Assert(cond, exc)
std::function< VectorType &()> get_solution_scaling
static::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
unsigned int maximum_beta_failures
std::function< int(const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system
std::function< void(VectorType &)> reinit_vector
KINSOL(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
std::function< int(const VectorType ¤t_u, const VectorType ¤t_f)> setup_jacobian
unsigned int solve(VectorType &initial_guess_and_solution)
static::ExceptionBase & ExcInternalError()