# include <deal.II/lac/petsc_vector.h>
# endif
-# include <deal.II/sundials/copy.h>
# include <deal.II/sundials/n_vector.h>
// Make sure we #include the SUNDIALS config file...
# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
, kinsol_ctx(nullptr)
# endif
- , solution(nullptr)
- , u_scale(nullptr)
- , f_scale(nullptr)
, pending_exception(nullptr)
{
set_functions_to_trigger_an_assert();
status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
AssertKINSOL(status);
- const auto system_size = initial_guess_and_solution.size();
-
- // The solution is stored in solution. Here we take only a view of it.
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- const IndexSet &is =
- initial_guess_and_solution.locally_owned_elements();
- const auto local_system_size = is.n_elements();
-
- solution = N_VNew_Parallel(mpi_communicator,
- local_system_size,
- system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
-
- u_scale = N_VNew_Parallel(mpi_communicator,
- local_system_size,
- system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
- N_VConst_Parallel(1.e0, u_scale);
-
- f_scale = N_VNew_Parallel(mpi_communicator,
- local_system_size,
- system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
- N_VConst_Parallel(1.e0, f_scale);
- }
- else
+ // helper function to create N_Vectors compatible with different versions
+ // of SUNDIALS
+ const auto make_compatible_nvector_view = [this](auto &v) {
+ return internal::make_nvector_view(v
+# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+ ,
+ kinsol_ctx
# endif
+ );
+ };
+
+
+ VectorType ones;
+ // Prepare constant vector for scaling f or u only if we need it
+ if (!get_function_scaling || !get_solution_scaling)
{
- Assert(is_serial_vector<VectorType>::value, ExcInternalError());
- solution = N_VNew_Serial(system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
- u_scale = N_VNew_Serial(system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
- N_VConst_Serial(1.e0, u_scale);
- f_scale = N_VNew_Serial(system_size
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- kinsol_ctx
-# endif
- );
- N_VConst_Serial(1.e0, f_scale);
+ reinit_vector(ones);
+ ones = 1.0;
}
- if (get_solution_scaling)
- internal::copy(u_scale, get_solution_scaling());
-
- if (get_function_scaling)
- internal::copy(f_scale, get_function_scaling());
+ auto u_scale = make_compatible_nvector_view(
+ get_solution_scaling ? get_solution_scaling() : ones);
+ auto f_scale = make_compatible_nvector_view(
+ get_function_scaling ? get_function_scaling() : ones);
- internal::copy(solution, initial_guess_and_solution);
+ auto solution = make_compatible_nvector_view(initial_guess_and_solution);
// This must be called before KINSetMAA
status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
[](N_Vector yy, N_Vector FF, void *user_data) -> int {
KINSOL<VectorType> &solver =
*static_cast<KINSOL<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src_yy(mem);
- solver.reinit_vector(*src_yy);
-
- typename VectorMemory<VectorType>::Pointer dst_FF(mem);
- solver.reinit_vector(*dst_FF);
- internal::copy(*src_yy, yy);
+ auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+ auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
Assert(solver.iteration_function, ExcInternalError());
*src_yy,
*dst_FF);
- internal::copy(FF, *dst_FF);
-
return err;
},
solution);
[](N_Vector yy, N_Vector FF, void *user_data) -> int {
KINSOL<VectorType> &solver =
*static_cast<KINSOL<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src_yy(mem);
- solver.reinit_vector(*src_yy);
-
- typename VectorMemory<VectorType>::Pointer dst_FF(mem);
- solver.reinit_vector(*dst_FF);
- internal::copy(*src_yy, yy);
+ auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+ auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
Assert(solver.residual, ExcInternalError());
const int err = call_and_possibly_capture_exception(
solver.residual, solver.pending_exception, *src_yy, *dst_FF);
- internal::copy(FF, *dst_FF);
-
return err;
},
solution);
SUNLinearSolver LS = nullptr;
if (solve_jacobian_system ||
- solve_with_jacobian) // user assigned a function object to the solver
- // slot
+ solve_with_jacobian) // user assigned a function
+ // object to the solver slot
{
// Set the operations we care for in the sun_linear_solver object
// and attach it to the KINSOL object. The functions that will get
// signals to call. Let's first check the more modern one:
if (solver.solve_with_jacobian)
{
- // 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);
-
- internal::copy(*src_b, b);
+ auto src_b = internal::unwrap_nvector_const<VectorType>(b);
+ auto dst_x = internal::unwrap_nvector<VectorType>(x);
const int err =
call_and_possibly_capture_exception(solver.solve_with_jacobian,
*dst_x,
tol);
- internal::copy(x, *dst_x);
-
return err;
}
else
// the old signal. Check this.
Assert(solver.solve_jacobian_system, ExcInternalError());
- // Allocate temporary (deal.II-type) vectors into which to copy
- // the N_vectors
+ // Allocate temporary (deal.II-type) dummy vectors
GrowingVectorMemory<VectorType> mem;
typename VectorMemory<VectorType>::Pointer src_ycur(mem);
typename VectorMemory<VectorType>::Pointer src_fcur(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);
-
- internal::copy(*src_b, b);
+ auto src_b = internal::unwrap_nvector_const<VectorType>(b);
+ auto dst_x = internal::unwrap_nvector<VectorType>(x);
// Call the user-provided setup function with these arguments.
// Note that Sundials 4.x and later no longer provide values for
*src_b,
*dst_x);
- internal::copy(x, *dst_x);
-
return err;
}
};
const KINSOL<VectorType> &solver =
*static_cast<const KINSOL<VectorType> *>(user_data);
- // Allocate temporary (deal.II-type) vectors into which to copy the
- // N_vectors
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer ycur(mem);
- typename VectorMemory<VectorType>::Pointer fcur(mem);
- solver.reinit_vector(*ycur);
- solver.reinit_vector(*fcur);
-
- internal::copy(*ycur, u);
- internal::copy(*fcur, f);
+ auto ycur = internal::unwrap_nvector_const<VectorType>(u);
+ auto fcur = internal::unwrap_nvector<VectorType>(f);
// Call the user-provided setup function with these arguments:
return call_and_possibly_capture_exception(solver.setup_jacobian,
status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
ScopeExit upon_exit([this, &J, &LS]() mutable {
- // Free the vectors which are no longer used.
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- N_VDestroy_Parallel(solution);
- N_VDestroy_Parallel(u_scale);
- N_VDestroy_Parallel(f_scale);
- }
- else
-# endif
- {
- N_VDestroy_Serial(solution);
- N_VDestroy_Serial(u_scale);
- N_VDestroy_Serial(f_scale);
- }
-
if (J != nullptr)
SUNMatDestroy(J);
if (LS != nullptr)
}
AssertKINSOL(status);
- internal::copy(initial_guess_and_solution, solution);
-
long nniters;
status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
AssertKINSOL(status);
projects / dealii.git / commitdiff