*
* KINSOL's Newton solver employs the inexact Newton method. As this solver
* is intended mainly for large systems, the user is required to provide
- * their own solver function. If a solver function is not provided, the
- * internal dense solver of KINSOL is used. Be warned that this solver
- * computes the Jacobian approximately, and may be efficient only for small
- * systems.
+ * their own solver function.
*
* At the highest level, KINSOL implements the following iteration
* scheme:
* algorithm is that the full Newton step tends to be taken close to the
* solution.
*
- * As a user option, KINSOL permits the application of inequality
- * constraints, $u_i > 0$ and $u_i < 0$, as well as $u_i \geq 0$ and $u_i
- * \leq 0$, where $u_i$ is the $i$-th component of $u$. Any such constraint,
- * or no constraint, may be imposed on each component by providing the
- * optional functions
- * - get_lower_than_zero_constrained_entries()
- * - get_greater_than_zero_constrained_entries()
- * - get_lower_equal_than_zero_constrained_entries()
- * - get_greater_or_equal_than_zero_constrained_entries()
- *
- * KINSOL will reduce step lengths in order to ensure that no constraint is
- * violated. Specifically, if a new Newton iterate will violate a constraint,
- * the maximum step length along the Newton direction that will satisfy all
- * constraints is found, and $\delta_n$ is scaled to take a step of that
- * length.
- *
* The basic fixed-point iteration scheme implemented in KINSOL is given by:
* - Set $u_0 =$ an initial guess
* - For $n = 0, 1, 2, \dots$ until convergence do:
* or
* - iteration_function;
*
- * Specifying residual() allows the user to use Newton strategies (i.e.,
- * $F(u)=0$ will be solved), while specifying iteration_function(), fixed
- * point iteration or Picard iteration will be used (i.e., $G(u)=u$ will be
- * solved).
+ * Specifying residual() allows the user to use Newton and Picard strategies
+ * (i.e., $F(u)=0$ will be solved), while specifying iteration_function(), a
+ * fixed point iteration will be used (i.e., $G(u)=u$ will be solved).
*
- * If the use of a Newton method is desired, then the user should also supply
- * - solve_jacobian_system;
+ * If the use of a Newton or Picard method is desired, then the user should
+ * also supply
+ * - solve_jacobian_system or solve_with_jacobian;
* and optionally
* - setup_jacobian;
*
- * If the solve_jacobian_system() function is not supplied, then KINSOL will
- * use its internal dense solver for Newton methods, with approximate
- * Jacobian. This may be very expensive for large systems. Fixed point
- * iteration does not require the solution of any linear system.
+ * Fixed point iteration does not require the solution of any linear system.
*
* Also the following functions could be rewritten, to provide additional
* scaling factors for both the solution and the residual evaluation during
*/
void *kinsol_mem;
- /**
- * MPI communicator. SUNDIALS solver runs happily in parallel.
- */
- MPI_Comm communicator;
-
/**
* Memory pool of vectors.
*/
{
template <typename VectorType>
int
- residual_or_iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
+ residual_callback(N_Vector yy, N_Vector FF, void *user_data)
{
KINSOL<VectorType> &solver =
*static_cast<KINSOL<VectorType> *>(user_data);
int err = 0;
if (solver.residual)
err = solver.residual(*src_yy, *dst_FF);
- else if (solver.iteration_function)
+ else
+ Assert(false, ExcInternalError());
+
+ return err;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
+ {
+ KINSOL<VectorType> &solver =
+ *static_cast<KINSOL<VectorType> *>(user_data);
+
+ auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+ auto *dst_FF = internal::unwrap_nvector<VectorType>(FF);
+
+ int err = 0;
+ if (solver.iteration_function)
err = solver.iteration_function(*src_yy, *dst_FF);
else
Assert(false, ExcInternalError());
return err;
}
}
-
# endif
} // namespace
template <typename VectorType>
- KINSOL<VectorType>::KINSOL(const AdditionalData &data,
- const MPI_Comm & mpi_comm)
+ KINSOL<VectorType>::KINSOL(const AdditionalData &data, const MPI_Comm &)
: data(data)
, kinsol_mem(nullptr)
- , communicator(is_serial_vector<VectorType>::value ?
- MPI_COMM_SELF :
- Utilities::MPI::duplicate_communicator(mpi_comm))
{
set_functions_to_trigger_an_assert();
}
{
if (kinsol_mem)
KINFree(&kinsol_mem);
-
-# 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
}
unsigned int
KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
{
- unsigned int system_size = initial_guess_and_solution.size();
-
NVectorView<VectorType> u_scale, f_scale;
VectorType u_scale_temp, f_scale_temp;
f_scale = internal::make_nvector_view(f_scale_temp);
}
+ // Make sure we have what we need
+ if (data.strategy == AdditionalData::fixed_point)
+ {
+ Assert(iteration_function,
+ ExcFunctionNotProvided("iteration_function"));
+ }
+ else
+ {
+ Assert(residual, ExcFunctionNotProvided("residual"));
+ Assert(solve_jacobian_system || solve_with_jacobian,
+ ExcFunctionNotProvided(
+ "solve_jacobian_system || solve_with_jacobian"));
+ }
+
auto solution = internal::make_nvector_view(initial_guess_and_solution);
if (kinsol_mem)
kinsol_mem = KINCreate();
- int status =
- KINInit(kinsol_mem, residual_or_iteration_callback<VectorType>, solution);
+ int status = 0;
(void)status;
- AssertKINSOL(status);
status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
AssertKINSOL(status);
+ // This must be called before KINSetMAA
status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
AssertKINSOL(status);
+ // From the manual: this must be called BEFORE KINInit
+ status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
+ AssertKINSOL(status);
+
+ if (data.strategy == AdditionalData::fixed_point)
+ status = KINInit(kinsol_mem, iteration_callback<VectorType>, solution);
+ else
+ status = KINInit(kinsol_mem, residual_callback<VectorType>, solution);
+ AssertKINSOL(status);
+
status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
AssertKINSOL(status);
status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
AssertKINSOL(status);
- status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
- AssertKINSOL(status);
-
status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
AssertKINSOL(status);
{
/* interface up to and including 4.0 */
# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
- auto KIN_mem = static_cast<KINMem>(kinsol_mem);
- KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
+ auto KIN_mem = static_cast<KINMem>(kinsol_mem);
+ // Old version only works with solve_jacobian_system
+ Assert(solve_jacobian_system,
+ ExcFunctionNotProvided("solve_jacobian_system"))
+ KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
if (setup_jacobian) // user assigned a function object to the Jacobian
// set-up slot
KIN_mem->kin_lsetup = setup_jacobian_callback<VectorType>;
}
# endif
}
- else
- {
- J = SUNDenseMatrix(system_size, system_size);
- LS = SUNDenseLinearSolver(u_scale, J);
- status = KINDlsSetLinearSolver(kinsol_mem, LS, J);
- AssertKINSOL(status);
- }
-
- if (data.strategy == AdditionalData::newton ||
- data.strategy == AdditionalData::linesearch)
- Assert(residual, ExcFunctionNotProvided("residual"));
-
- if (data.strategy == AdditionalData::fixed_point ||
- data.strategy == AdditionalData::picard)
- Assert(iteration_function, ExcFunctionNotProvided("iteration_function"));
// call to KINSol
status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
AssertKINSOL(status);
- SUNMatDestroy(J);
- SUNLinSolFree(LS);
+ if (J != nullptr)
+ SUNMatDestroy(J);
+ if (LS != nullptr)
+ SUNLinSolFree(LS);
KINFree(&kinsol_mem);
return static_cast<unsigned int>(nniters);
reinit_vector = [](VectorType &) {
AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
};
+
+ setup_jacobian = [](const VectorType &, const VectorType &) { return 0; };
}
template class KINSOL<Vector<double>>;
#include "../tests.h"
-// Solve a nonlinear system but provide only residual function. KINSOL
-// then uses its internal solvers which are based on a
-// finite-difference approximation to the Jacobian and a direct
-// solver.
+// Solve a nonlinear system using fixed point iteration, and Anderson
+// acceleration
/**
- * Solve the non linear problem
- *
- * F(u) = 0 , where f_i(u) = u_i^2 - i^2, 0 <= i < N
+ * The following is a simple example problem, with the coding
+ * needed for its solution by the accelerated fixed point solver in
+ * KINSOL.
+ * The problem is from chemical kinetics, and consists of solving
+ * the first time step in a Backward Euler solution for the
+ * following three rate equations:
+ * dy1/dt = -.04*y1 + 1.e4*y2*y3
+ * dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e2*(y2)^2
+ * dy3/dt = 3.e2*(y2)^2
+ * on the interval from t = 0.0 to t = 0.1, with initial
+ * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
+ * Run statistics (optional outputs) are printed at the end.
*/
int
main(int argc, char **argv)
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_fixed_point.prm");
prm.parse_input(ifile);
// Size of the problem
- unsigned int N = 10;
+ unsigned int N = 3;
SUNDIALS::KINSOL<VectorType> kinsol(data);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
- for (unsigned int i = 0; i < u.size(); ++i)
- F[i] = u[i] * u[i] - (i + 1) * (i + 1);
- return 0;
- };
+ // Robert example
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+ const double dstep = 0.1;
+ const double y10 = 1.0;
+ const double y20 = 0.0;
+ const double y30 = 0.0;
+ const double yd1 = dstep * (-0.04 * u[0] + 1.0e4 * u[1] * u[2]);
+ const double yd3 = dstep * 3.0e2 * u[1] * u[1];
+
+ F[0] = yd1 + y10;
+ F[1] = -yd1 - yd3 + y20;
+ F[2] = yd3 + y30;
- kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
- for (unsigned int i = 0; i < u.size(); ++i)
- F[i] = u[i] * u[i] - i * i - u[i];
return 0;
};
VectorType v(N);
- v = 1.0;
+ v[0] = 1;
auto niter = kinsol.solve(v);
v.print(deallog.get_file_stream());
deallog << "Converged in " << niter << " iterations." << std::endl;
+++ /dev/null
-
-1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000 10.0000
-DEAL::Converged in 23 iterations.
--- /dev/null
+
+9.968e-01 2.953e-03 2.616e-04
+DEAL::Converged in 8 iterations.
#include "../tests.h"
-// Solve a nonlinear system but provide only residual function. KINSOL
-// then uses its internal solvers which are based on a
-// finite-difference approximation to the Jacobian and a direct
-// solver.
+// Solve a nonlinear system in the form accepted by Picard iteration.
//
-// Compared to the _01 test, this is simply a more complicated function:
// We solve the nonlinear problem
//
-// F(u) = 0
+// F(u) = L u + N(u) = 0
//
-// with a 2-dimensional vector u and where
+// where L is a constant matrix, and N(u) is non linear.
//
-// F(u) = [ cos(u1 + u2) - 1 ] -> u1=-u2
-// [ sin(u1 - u2) ] -> u1=u2
+// We set L = id and
+//
+// N_i(u) = .1*u_i^2 - i - 1
//
-// In other words, we need to find the solution u1=u2=0.
-
int
main(int argc, char **argv)
{
using VectorType = Vector<double>;
+ // Size of the problem
+ unsigned int N = 2;
+
+ FullMatrix<double> L(N, N);
+ L(0, 0) = 1;
+ L(1, 1) = 1;
+ L(0, 1) = 1;
+
+ FullMatrix<double> Linv(N, N);
+ Linv.invert(L);
+
SUNDIALS::KINSOL<VectorType>::AdditionalData data;
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_picard.prm");
prm.parse_input(ifile);
- // Size of the problem
- unsigned int N = 2;
-
SUNDIALS::KINSOL<VectorType> kinsol(data);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
- F(0) = std::cos(u[0] + u[1]) - 1;
- F(1) = std::sin(u[0] - u[1]);
+ kinsol.residual = [&](const VectorType &u, VectorType &F) -> int {
+ F = u;
+
+ F[0] += .1 * u[0] * u[0] - 1;
+ F[1] += .1 * u[1] * u[1] - 2;
return 0;
};
+ kinsol.solve_with_jacobian =
+ [&](const VectorType &rhs, VectorType &dst, double) -> int {
+ dst = rhs;
+ return 0;
+ };
- kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
- // We want a Newton-type scheme, not a fixed point iteration. So we
- // shouldn't get into this function.
- std::abort();
-
- // But if anyone wanted to see how it would look like:
- F(0) = std::cos(u[0] + u[1]) - 1 - u[0];
- F(1) = std::sin(u[0] - u[1]) - u[1];
+ kinsol.solve_jacobian_system = [&](const VectorType &,
+ const VectorType &,
+ const VectorType &rhs,
+ VectorType & dst) -> int {
+ dst = rhs;
return 0;
};
VectorType v(N);
- v(0) = 0.5;
- v(1) = 1.234;
auto niter = kinsol.solve(v);
+
v.print(deallog.get_file_stream());
deallog << "Converged in " << niter << " iterations." << std::endl;
}
+++ /dev/null
-
-9.761e-04 9.761e-04
-DEAL::Converged in 27 iterations.
--- /dev/null
+
+9.161e-01 1.708e+00
+DEAL::Converged in 7 iterations.
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_newton.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
prm.parse_input(ifile);
// Update the Jacobian in each iteration:
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
prm.parse_input(ifile);
// Update the Jacobian in each iteration:
--- /dev/null
+set Function norm stopping tolerance = 1e-10
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance = 1e-10
+set Solution strategy = fixed_point
+subsection Fixed point and Picard parameters
+ set Anderson acceleration subspace size = 2
+end
+subsection Linesearch parameters
+ set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+ set Maximum allowable scaled length of the Newton step = 0.000000
+ set Maximum iterations without matrix setup = 0
+ set No initial matrix setup = false
+ set Relative error for different quotient computation = 0.000000
+end
--- /dev/null
+set Function norm stopping tolerance = 0.000000
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance = 0.000000
+set Solution strategy = newton
+subsection Fixed point and Picard parameters
+ set Anderson acceleration subspace size = 5
+end
+subsection Linesearch parameters
+ set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+ set Maximum allowable scaled length of the Newton step = 0.000000
+ set Maximum iterations without matrix setup = 0
+ set No initial matrix setup = false
+ set Relative error for different quotient computation = 0.000000
+end
--- /dev/null
+set Function norm stopping tolerance = 1e-10
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance = 1e-10
+set Solution strategy = picard
+subsection Fixed point and Picard parameters
+ set Anderson acceleration subspace size = 2
+end
+subsection Linesearch parameters
+ set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+ set Maximum allowable scaled length of the Newton step = 0.000000
+ set Maximum iterations without matrix setup = 0
+ set No initial matrix setup = false
+ set Relative error for different quotient computation = 0.000000
+end
projects / dealii.git / commitdiff