From: Daniel Arndt Date: Fri, 5 Nov 2021 15:24:41 +0000 (-0400) Subject: Revert "Fix Kinsol interface." X-Git-Tag: v9.3.2~1^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=713b795316b34d4377051b7ea5aba26d05ba136d;p=dealii.git Revert "Fix Kinsol interface." This reverts commit 1032025f3a076f97fc42e43fbd927278b6aec327. --- diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 9bf6c55180..939215b05f 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -66,7 +66,10 @@ namespace SUNDIALS * * 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: @@ -133,6 +136,22 @@ namespace SUNDIALS * 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: @@ -161,17 +180,20 @@ namespace SUNDIALS * 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 @@ -689,6 +711,11 @@ namespace SUNDIALS */ void *kinsol_mem; + /** + * MPI communicator. SUNDIALS solver runs happily in parallel. + */ + MPI_Comm communicator; + /** * Memory pool of vectors. */ diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index bdbac29035..cb12781b0e 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -133,7 +133,7 @@ namespace SUNDIALS { template 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 &solver = *static_cast *>(user_data); @@ -144,26 +144,7 @@ namespace SUNDIALS int err = 0; if (solver.residual) err = solver.residual(*src_yy, *dst_FF); - else - Assert(false, ExcInternalError()); - - return err; - } - - - - template - int - iteration_callback(N_Vector yy, N_Vector FF, void *user_data) - { - KINSOL &solver = - *static_cast *>(user_data); - - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *dst_FF = internal::unwrap_nvector(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()); @@ -300,15 +281,20 @@ namespace SUNDIALS return err; } } + # endif } // namespace template - KINSOL::KINSOL(const AdditionalData &data, const MPI_Comm &) + KINSOL::KINSOL(const AdditionalData &data, + const MPI_Comm & mpi_comm) : data(data) , kinsol_mem(nullptr) + , communicator(is_serial_vector::value ? + MPI_COMM_SELF : + Utilities::MPI::duplicate_communicator(mpi_comm)) { set_functions_to_trigger_an_assert(); } @@ -320,6 +306,15 @@ namespace SUNDIALS { if (kinsol_mem) KINFree(&kinsol_mem); + +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + const int ierr = MPI_Comm_free(&communicator); + (void)ierr; + AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr)); + } +# endif } @@ -328,6 +323,8 @@ namespace SUNDIALS unsigned int KINSOL::solve(VectorType &initial_guess_and_solution) { + unsigned int system_size = initial_guess_and_solution.size(); + NVectorView u_scale, f_scale; VectorType u_scale_temp, f_scale_temp; @@ -350,20 +347,6 @@ namespace SUNDIALS 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) @@ -371,26 +354,17 @@ namespace SUNDIALS kinsol_mem = KINCreate(); - int status = 0; + int status = + KINInit(kinsol_mem, residual_or_iteration_callback, solution); (void)status; + AssertKINSOL(status); status = KINSetUserData(kinsol_mem, static_cast(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, solution); - else - status = KINInit(kinsol_mem, residual_callback, solution); - AssertKINSOL(status); - status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance); AssertKINSOL(status); @@ -409,6 +383,9 @@ namespace SUNDIALS 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); @@ -421,11 +398,8 @@ namespace SUNDIALS { /* interface up to and including 4.0 */ # if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0) - auto KIN_mem = static_cast(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; + auto KIN_mem = static_cast(kinsol_mem); + KIN_mem->kin_lsolve = solve_with_jacobian_callback; if (setup_jacobian) // user assigned a function object to the Jacobian // set-up slot KIN_mem->kin_lsetup = setup_jacobian_callback; @@ -512,6 +486,21 @@ namespace SUNDIALS } # 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); @@ -521,10 +510,8 @@ namespace SUNDIALS 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(nniters); @@ -537,8 +524,6 @@ namespace SUNDIALS reinit_vector = [](VectorType &) { AssertThrow(false, ExcFunctionNotProvided("reinit_vector")); }; - - setup_jacobian = [](const VectorType &, const VectorType &) { return 0; }; } template class KINSOL>; diff --git a/tests/sundials/kinsol_01.cc b/tests/sundials/kinsol_01.cc index 028491fcea..15e9005bd7 100644 --- a/tests/sundials/kinsol_01.cc +++ b/tests/sundials/kinsol_01.cc @@ -22,22 +22,15 @@ #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) @@ -53,35 +46,31 @@ 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 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; diff --git a/tests/sundials/kinsol_linesearch.prm b/tests/sundials/kinsol_01.prm similarity index 100% rename from tests/sundials/kinsol_linesearch.prm rename to tests/sundials/kinsol_01.prm diff --git a/tests/sundials/kinsol_01.with_sundials=true.expect=run.output b/tests/sundials/kinsol_01.with_sundials=true.expect=run.output new file mode 100644 index 0000000000..2a13a0ef6e --- /dev/null +++ b/tests/sundials/kinsol_01.with_sundials=true.expect=run.output @@ -0,0 +1,3 @@ + +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. diff --git a/tests/sundials/kinsol_01.with_sundials=true.output b/tests/sundials/kinsol_01.with_sundials=true.output deleted file mode 100644 index 61ccde6ace..0000000000 --- a/tests/sundials/kinsol_01.with_sundials=true.output +++ /dev/null @@ -1,3 +0,0 @@ - -9.968e-01 2.953e-03 2.616e-04 -DEAL::Converged in 8 iterations. diff --git a/tests/sundials/kinsol_02.cc b/tests/sundials/kinsol_02.cc index dc56917859..5a01ab35c6 100644 --- a/tests/sundials/kinsol_02.cc +++ b/tests/sundials/kinsol_02.cc @@ -22,18 +22,23 @@ #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) { @@ -44,54 +49,43 @@ main(int argc, char **argv) using VectorType = Vector; - // Size of the problem - unsigned int N = 2; - - FullMatrix L(N, N); - L(0, 0) = 1; - L(1, 1) = 1; - L(0, 1) = 1; - - FullMatrix Linv(N, N); - Linv.invert(L); - SUNDIALS::KINSOL::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 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; } diff --git a/tests/sundials/kinsol_02.with_sundials=true.expect=run.output b/tests/sundials/kinsol_02.with_sundials=true.expect=run.output new file mode 100644 index 0000000000..ee3d8fab69 --- /dev/null +++ b/tests/sundials/kinsol_02.with_sundials=true.expect=run.output @@ -0,0 +1,3 @@ + +9.761e-04 9.761e-04 +DEAL::Converged in 27 iterations. diff --git a/tests/sundials/kinsol_02.with_sundials=true.output b/tests/sundials/kinsol_02.with_sundials=true.output deleted file mode 100644 index 38ba409016..0000000000 --- a/tests/sundials/kinsol_02.with_sundials=true.output +++ /dev/null @@ -1,3 +0,0 @@ - -9.161e-01 1.708e+00 -DEAL::Converged in 7 iterations. diff --git a/tests/sundials/kinsol_03.cc b/tests/sundials/kinsol_03.cc index 1c2c5ebee6..d5368565cd 100644 --- a/tests/sundials/kinsol_03.cc +++ b/tests/sundials/kinsol_03.cc @@ -58,7 +58,7 @@ main(int argc, char **argv) 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 diff --git a/tests/sundials/kinsol_03_new_interface.cc b/tests/sundials/kinsol_03_new_interface.cc index a0d79e3993..f473627aa8 100644 --- a/tests/sundials/kinsol_03_new_interface.cc +++ b/tests/sundials/kinsol_03_new_interface.cc @@ -58,7 +58,7 @@ main(int argc, char **argv) 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 diff --git a/tests/sundials/kinsol_04.cc b/tests/sundials/kinsol_04.cc index b0f079685d..f7c8baaf41 100644 --- a/tests/sundials/kinsol_04.cc +++ b/tests/sundials/kinsol_04.cc @@ -70,7 +70,7 @@ main(int argc, char **argv) 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 diff --git a/tests/sundials/kinsol_04_new_interface.cc b/tests/sundials/kinsol_04_new_interface.cc index da932551ed..73ab99866a 100644 --- a/tests/sundials/kinsol_04_new_interface.cc +++ b/tests/sundials/kinsol_04_new_interface.cc @@ -70,7 +70,7 @@ main(int argc, char **argv) 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 diff --git a/tests/sundials/kinsol_05.cc b/tests/sundials/kinsol_05.cc index a5dfbed054..5c5c3464f0 100644 --- a/tests/sundials/kinsol_05.cc +++ b/tests/sundials/kinsol_05.cc @@ -64,7 +64,7 @@ main(int argc, char **argv) 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: diff --git a/tests/sundials/kinsol_05_new_interface.cc b/tests/sundials/kinsol_05_new_interface.cc index 4d2e05ecf5..83648603e1 100644 --- a/tests/sundials/kinsol_05_new_interface.cc +++ b/tests/sundials/kinsol_05_new_interface.cc @@ -64,7 +64,7 @@ main(int argc, char **argv) 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: diff --git a/tests/sundials/kinsol_fixed_point.prm b/tests/sundials/kinsol_fixed_point.prm deleted file mode 100644 index a64b3cac28..0000000000 --- a/tests/sundials/kinsol_fixed_point.prm +++ /dev/null @@ -1,16 +0,0 @@ -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 diff --git a/tests/sundials/kinsol_newton.prm b/tests/sundials/kinsol_newton.prm deleted file mode 100644 index 3e5ce272a8..0000000000 --- a/tests/sundials/kinsol_newton.prm +++ /dev/null @@ -1,16 +0,0 @@ -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 diff --git a/tests/sundials/kinsol_picard.prm b/tests/sundials/kinsol_picard.prm deleted file mode 100644 index f684cc3e2b..0000000000 --- a/tests/sundials/kinsol_picard.prm +++ /dev/null @@ -1,16 +0,0 @@ -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