*
* 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.
+ * 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.
*
* 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 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).
+ * 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).
*
- * If the use of a Newton or Picard method is desired, then the user should
- * also supply
- * - solve_jacobian_system or solve_with_jacobian;
+ * If the use of a Newton method is desired, then the user should also supply
+ * - solve_jacobian_system;
* and optionally
* - setup_jacobian;
*
- * Fixed point iteration does not require the solution of any linear system.
+ * 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.
*
* 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_callback(N_Vector yy, N_Vector FF, void *user_data)
+ residual_or_iteration_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
- 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)
+ else 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 &)
+ KINSOL<VectorType>::KINSOL(const AdditionalData &data,
+ const MPI_Comm & 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 = 0;
+ int status =
+ KINInit(kinsol_mem, residual_or_iteration_callback<VectorType>, solution);
(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);
- // 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>;
+ auto KIN_mem = static_cast<KINMem>(kinsol_mem);
+ 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);
- if (J != nullptr)
- SUNMatDestroy(J);
- if (LS != nullptr)
- SUNLinSolFree(LS);
+ SUNMatDestroy(J);
+ 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 using fixed point iteration, and Anderson
-// acceleration
+// 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.
/**
- * 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.
+ * Solve the non linear problem
+ *
+ * F(u) = 0 , where f_i(u) = u_i^2 - i^2, 0 <= i < N
*/
int
main(int argc, char **argv)
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_fixed_point.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Size of the problem
- unsigned int N = 3;
+ unsigned int N = 10;
SUNDIALS::KINSOL<VectorType> kinsol(data);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- // 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];
+ 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;
+ };
- 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[0] = 1;
+ v = 1.0;
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 in the form accepted by Picard iteration.
+// 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.
//
+// Compared to the _01 test, this is simply a more complicated function:
// We solve the nonlinear problem
//
-// F(u) = L u + N(u) = 0
+// F(u) = 0
//
-// where L is a constant matrix, and N(u) is non linear.
+// with a 2-dimensional vector u and where
//
-// We set L = id and
-//
-// N_i(u) = .1*u_i^2 - i - 1
+// F(u) = [ cos(u1 + u2) - 1 ] -> u1=-u2
+// [ sin(u1 - u2) ] -> u1=u2
//
+// 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_picard.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.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 = u;
-
- F[0] += .1 * u[0] * u[0] - 1;
- F[1] += .1 * u[1] * u[1] - 2;
+ 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]);
return 0;
};
- kinsol.solve_with_jacobian =
- [&](const VectorType &rhs, VectorType &dst, double) -> int {
- dst = rhs;
- return 0;
- };
- kinsol.solve_jacobian_system = [&](const VectorType &,
- const VectorType &,
- const VectorType &rhs,
- VectorType & dst) -> int {
- dst = rhs;
+ 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];
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_newton.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Size of the problem
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm");
prm.parse_input(ifile);
// Update the Jacobian in each iteration:
ParameterHandler prm;
data.add_parameters(prm);
- std::ifstream ifile(SOURCE_DIR "/kinsol_linesearch.prm");
+ std::ifstream ifile(SOURCE_DIR "/kinsol_01.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