From: Wolfgang Bangerth Date: Fri, 5 Mar 2021 01:36:59 +0000 (+0100) Subject: Add more KINSOL tests. X-Git-Tag: v9.3.0-rc1~295^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=57ec5b495950ff26e8e66b2ab41a3aba887c2040;p=dealii.git Add more KINSOL tests. --- diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index 6738c0cc01..28a8a3ea17 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -442,7 +442,7 @@ namespace SUNDIALS // called do not actually receive the KINSOL object, just the LS // object, so we have to store a pointer to the current // object in the LS object - LS = SUNLinSolNewEmpty(); + LS = SUNLinSolNewEmpty(); LS->content = this; LS->ops->gettype = @@ -472,8 +472,8 @@ namespace SUNDIALS // if we don't set it, it won't call the functions that set up // the matrix object (i.e., the argument to the 'KINSetJacFn' // function below). - J = SUNMatNewEmpty(); - J->content = this; + J = SUNMatNewEmpty(); + J->content = this; J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { return SUNMATRIX_CUSTOM; diff --git a/tests/sundials/kinsol_02.cc b/tests/sundials/kinsol_02.cc new file mode 100644 index 0000000000..8ef0bf95e2 --- /dev/null +++ b/tests/sundials/kinsol_02.cc @@ -0,0 +1,91 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#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. +// +// Compared to the _01 test, this is simply a more complicated function: +// We solve the nonlinear problem +// +// F(u) = 0 +// +// with a 2-dimensional vector u and where +// +// 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) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + using VectorType = Vector; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(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(0) = std::cos(u[0] + u[1]) - 1; + F(1) = std::sin(u[0] - u[1]); + 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]; + 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.output b/tests/sundials/kinsol_02.with_sundials=true.output new file mode 100644 index 0000000000..ee3d8fab69 --- /dev/null +++ b/tests/sundials/kinsol_02.with_sundials=true.output @@ -0,0 +1,3 @@ + +9.761e-04 9.761e-04 +DEAL::Converged in 27 iterations. diff --git a/tests/sundials/kinsol_03.cc b/tests/sundials/kinsol_03.cc new file mode 100644 index 0000000000..da8cdc8f1b --- /dev/null +++ b/tests/sundials/kinsol_03.cc @@ -0,0 +1,137 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#include "../tests.h" + +// Solve a nonlinear system. +// +// Similar to the _02 test, but we're now actually providing a solver with +// the Jacobian matrix. For the current case, +// +// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ] +// [ sin(u1 - u2) + 2*u2 ] +// +// the Jacobian is the 2x2 matrix +// +// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)] +// [ cos(u1 - u2) -cos(u1 - u2) + 2] +// +// The addition of the +2u_i to the function F does not move the solution +// (it is still u=0) but it makes sure that the Jacobian at the solution +// remains non-singular. +// +// This test case has a flaw in that starting at SUNDIALS 4.x, the solve +// function no longer receives the current 'u' vector. This means that +// one can't compute a proper Jacobian in that function because we don't +// know what to linearize around. The _04 test fixes this by computing +// the Jacobian in the setup function. +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + using VectorType = Vector; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(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 { + deallog << "Evaluating the solution at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + + F(0) = std::cos(u[0] + u[1]) - 1 + 2 * u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1]; + 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 + 2 * u[0] - u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1] - u[1]; + return 0; + }; + + + kinsol.setup_jacobian = [](const VectorType &u, const VectorType &F) -> int { + // We don't do any kind of set-up in this program, but we can at least + // say that we're here + deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + return 0; + }; + + + kinsol.solve_jacobian_system = [](const VectorType &, + const VectorType &, + const VectorType &rhs, + VectorType & dst) -> int { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ")" << std::endl; + + // This isn't right for SUNDIALS >4.0: We don't actually get a valid + // 'u' vector, and so do the linearization of the problem around + // the zero vector. This *happens* to converge, but it isn't the + // right approach. Check the _04 test for a better approach. + VectorType u(2); + u[0] = u[1] = 0; + + FullMatrix J(2, 2); + J(0, 0) = -std::sin(u[0] + u[1]) + 2; + J(0, 1) = -std::sin(u[0] + u[1]); + J(1, 0) = std::cos(u[0] - u[1]); + J(1, 1) = -std::cos(u[0] - u[1]) + 2; + + FullMatrix J_inverse(2, 2); + J_inverse.invert(J); + + J_inverse.vmult(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; +} diff --git a/tests/sundials/kinsol_03.with_sundials=true.output b/tests/sundials/kinsol_03.with_sundials=true.output new file mode 100644 index 0000000000..e8a1815705 --- /dev/null +++ b/tests/sundials/kinsol_03.with_sundials=true.output @@ -0,0 +1,20 @@ + +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Setting up Jacobian system at u=(0.500000,1.23400) +DEAL::Solving Jacobian system with rhs=(0.162480,-1.79816) +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Evaluating the solution at u=(0.581240,-0.645395) +DEAL::Solving Jacobian system with rhs=(-1.16042,0.349431) +DEAL::Evaluating the solution at u=(0.581240,-0.645395) +DEAL::Evaluating the solution at u=(0.00102861,0.284248) +DEAL::Solving Jacobian system with rhs=(0.0383589,-0.289047) +DEAL::Evaluating the solution at u=(0.00102861,0.284248) +DEAL::Evaluating the solution at u=(0.0202080,-0.0239792) +DEAL::Solving Jacobian system with rhs=(-0.0404090,0.00378553) +DEAL::Evaluating the solution at u=(0.0202081,-0.0239792) +DEAL::Evaluating the solution at u=(3.55540e-06,1.08225e-05) +DEAL::Solving Jacobian system with rhs=(-7.11070e-06,-1.43779e-05) +DEAL::Evaluating the solution at u=(3.56127e-06,1.08404e-05) +DEAL::Evaluating the solution at u=(5.16813e-11,-5.16814e-11) +5.168e-11 -5.168e-11 +DEAL::Converged in 5 iterations. diff --git a/tests/sundials/kinsol_04.cc b/tests/sundials/kinsol_04.cc new file mode 100644 index 0000000000..b2b97a0586 --- /dev/null +++ b/tests/sundials/kinsol_04.cc @@ -0,0 +1,142 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#include "../tests.h" + +// Solve a nonlinear system. +// +// Similar to the _02 test, but we're now actually providing a solver with +// the Jacobian matrix. For the current case, +// +// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ] +// [ sin(u1 - u2) + 2*u2 ] +// +// the Jacobian is the 2x2 matrix +// +// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)] +// [ cos(u1 - u2) -cos(u1 - u2) + 2] +// +// The addition of the +2u_i to the function F does not move the solution +// (it is still u=0) but it makes sure that the Jacobian at the solution +// remains non-singular +// +// The _03 test case has a flaw in that starting at SUNDIALS 4.x, the solve +// function no longer receives the current 'u' vector. This means that +// one can't compute a proper Jacobian in that function because we don't +// know what to linearize around. This test fixes this by computing +// the Jacobian in the setup function. +// +// It turns out that doing this leads to a *larger* number of +// iterations than the kludge in _03, but that is because _03 +// linearizes around the zero vector which just so also happens to be +// solution of the problem. In other words, it accidentally has the +// perfect Jacobian matrix to use, and consequently converges rapidly, +// whereas we here have to deal with a poor Jacobian for the first few +// iterations. Furthermore, because we don't update the Jacobian very +// frequently, we are stuck with the poor Jacobian for numerous +// iterations. The _05 test therefore forces updates in every +// iteration and, unsurprisingly, converges much quicker. + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + using VectorType = Vector; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(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 { + deallog << "Evaluating the solution at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + + F(0) = std::cos(u[0] + u[1]) - 1 + 2 * u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1]; + 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 + 2 * u[0] - u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1] - u[1]; + return 0; + }; + + FullMatrix J_inverse(2, 2); + + kinsol.setup_jacobian = [&J_inverse](const VectorType &u, + const VectorType &F) -> int { + deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + + FullMatrix J(2, 2); + J(0, 0) = -std::sin(u[0] + u[1]) + 2; + J(0, 1) = -std::sin(u[0] + u[1]); + J(1, 0) = std::cos(u[0] - u[1]); + J(1, 1) = -std::cos(u[0] - u[1]) + 2; + + J_inverse.invert(J); + + return 0; + }; + + + kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u, + const VectorType &, + const VectorType &rhs, + VectorType & dst) -> int { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ")" << std::endl; + + J_inverse.vmult(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; +} diff --git a/tests/sundials/kinsol_04.with_sundials=true.output b/tests/sundials/kinsol_04.with_sundials=true.output new file mode 100644 index 0000000000..6d54d87cb9 --- /dev/null +++ b/tests/sundials/kinsol_04.with_sundials=true.output @@ -0,0 +1,39 @@ + +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Setting up Jacobian system at u=(0.500000,1.23400) +DEAL::Solving Jacobian system with rhs=(0.162480,-1.79816) +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Evaluating the solution at u=(-0.282294,0.265967) +DEAL::Solving Jacobian system with rhs=(0.564722,-0.0107297) +DEAL::Evaluating the solution at u=(-0.282294,0.265967) +DEAL::Evaluating the solution at u=(0.0662880,0.0516109) +DEAL::Solving Jacobian system with rhs=(-0.125634,-0.117898) +DEAL::Evaluating the solution at u=(0.0662880,0.0516109) +DEAL::Evaluating the solution at u=(-0.0704025,0.0385647) +DEAL::Solving Jacobian system with rhs=(0.141312,0.0316222) +DEAL::Evaluating the solution at u=(-0.0704025,0.0385647) +DEAL::Evaluating the solution at u=(0.0336920,0.00224815) +DEAL::Solving Jacobian system with rhs=(-0.0667383,-0.0359350) +DEAL::Evaluating the solution at u=(0.0336921,0.00224815) +DEAL::Evaluating the solution at u=(-0.0257948,0.00879611) +DEAL::Solving Jacobian system with rhs=(0.0517342,0.0169918) +DEAL::Evaluating the solution at u=(-0.0257949,0.00879612) +DEAL::Evaluating the solution at u=(0.0149765,-0.00176527) +DEAL::Solving Jacobian system with rhs=(-0.0298657,-0.0132104) +DEAL::Evaluating the solution at u=(0.0149765,-0.00176527) +DEAL::Evaluating the solution at u=(-0.0102328,0.00261442) +DEAL::Solving Jacobian system with rhs=(0.0204945,0.00761799) +DEAL::Evaluating the solution at u=(-0.0102328,0.00261442) +DEAL::Evaluating the solution at u=(0.00635479,-0.00112180) +DEAL::Solving Jacobian system with rhs=(-0.0126959,-0.00523293) +DEAL::Evaluating the solution at u=(0.00635481,-0.00112180) +DEAL::Evaluating the solution at u=(-0.00417342,0.000933297) +DEAL::Solving Jacobian system with rhs=(0.00835208,0.00324010) +DEAL::Evaluating the solution at u=(-0.00417343,0.000933300) +DEAL::Evaluating the solution at u=(0.00265311,-0.000520866) +DEAL::Setting up Jacobian system at u=(0.00265311,-0.000520866) +DEAL::Solving Jacobian system with rhs=(-0.00530395,-0.00213224) +DEAL::Evaluating the solution at u=(0.00265313,-0.000520869) +DEAL::Evaluating the solution at u=(-1.13663e-06,1.12596e-06) +-1.137e-06 1.126e-06 +DEAL::Converged in 11 iterations. diff --git a/tests/sundials/kinsol_05.cc b/tests/sundials/kinsol_05.cc new file mode 100644 index 0000000000..898a95ca7a --- /dev/null +++ b/tests/sundials/kinsol_05.cc @@ -0,0 +1,142 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#include "../tests.h" + +// Solve a nonlinear system. +// +// Similar to the _02 test, but we're now actually providing a solver with +// the Jacobian matrix. For the current case, +// +// F(u) = [ cos(u1 + u2) - 1 + 2*u1 ] +// [ sin(u1 - u2) + 2*u2 ] +// +// the Jacobian is the 2x2 matrix +// +// J(u) = [ -sin(u1 + u2) + 2 -sin(u1 + u2)] +// [ cos(u1 - u2) -cos(u1 - u2) + 2] +// +// The addition of the +2u_i to the function F does not move the solution +// (it is still u=0) but it makes sure that the Jacobian at the solution +// remains non-singular +// +// The _03 test case has a flaw in that starting at SUNDIALS 4.x, the solve +// function no longer receives the current 'u' vector. This means that +// one can't compute a proper Jacobian in that function because we don't +// know what to linearize around. This test fixes this by computing +// the Jacobian in the setup function. +// +// The _04 test fixes this by setting up the Jacobian in the setup +// function, which has the current iterate. But it converges slowly +// because it doesn't update the Jacobian very often. This test +// finally updates it in every iteration. + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + using VectorType = Vector; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm"); + prm.parse_input(ifile); + + // Update the Jacobian in each iteration: + data.maximum_setup_calls = 1; + + + // 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 { + deallog << "Evaluating the solution at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + + F(0) = std::cos(u[0] + u[1]) - 1 + 2 * u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1]; + 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 + 2 * u[0] - u[0]; + F(1) = std::sin(u[0] - u[1]) + 2 * u[1] - u[1]; + return 0; + }; + + FullMatrix J_inverse(2, 2); + + kinsol.setup_jacobian = [&J_inverse](const VectorType &u, + const VectorType &F) -> int { + // We don't do any kind of set-up in this program, but we can at least + // say that we're here + deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ")" + << std::endl; + + FullMatrix J(2, 2); + J(0, 0) = -std::sin(u[0] + u[1]) + 2; + J(0, 1) = -std::sin(u[0] + u[1]); + J(1, 0) = std::cos(u[0] - u[1]); + J(1, 1) = -std::cos(u[0] - u[1]) + 2; + + J_inverse.invert(J); + + return 0; + }; + + + kinsol.solve_jacobian_system = [&J_inverse](const VectorType &u, + const VectorType &, + const VectorType &rhs, + VectorType & dst) -> int { + deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1] + << ")" << std::endl; + + J_inverse.vmult(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; +} diff --git a/tests/sundials/kinsol_05.with_sundials=true.output b/tests/sundials/kinsol_05.with_sundials=true.output new file mode 100644 index 0000000000..431b21399d --- /dev/null +++ b/tests/sundials/kinsol_05.with_sundials=true.output @@ -0,0 +1,20 @@ + +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Setting up Jacobian system at u=(0.500000,1.23400) +DEAL::Solving Jacobian system with rhs=(0.162480,-1.79816) +DEAL::Evaluating the solution at u=(0.500000,1.23400) +DEAL::Evaluating the solution at u=(-0.282294,0.265967) +DEAL::Setting up Jacobian system at u=(-0.282294,0.265967) +DEAL::Solving Jacobian system with rhs=(0.564722,-0.0107297) +DEAL::Evaluating the solution at u=(-0.282294,0.265967) +DEAL::Evaluating the solution at u=(-0.000445202,0.0468183) +DEAL::Setting up Jacobian system at u=(-0.000445202,0.0468183) +DEAL::Solving Jacobian system with rhs=(0.00196544,-0.0463907) +DEAL::Evaluating the solution at u=(-0.000445202,0.0468183) +DEAL::Evaluating the solution at u=(-0.000536539,0.000570488) +DEAL::Setting up Jacobian system at u=(-0.000536539,0.000570488) +DEAL::Solving Jacobian system with rhs=(0.00107308,-3.39491e-05) +DEAL::Evaluating the solution at u=(-0.000536554,0.000570504) +DEAL::Evaluating the solution at u=(-2.88123e-10,7.40347e-10) +-2.881e-10 7.403e-10 +DEAL::Converged in 4 iterations.