kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
// Robert example
- kinsol.iteration_function = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
const double dstep = 0.1;
const double y10 = 1.0;
const double y20 = 0.0;
F[0] = yd1 + y10;
F[1] = -yd1 - yd3 + y20;
F[2] = yd3 + y30;
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [&](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [&](const VectorType &u, VectorType &F) {
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;
- };
+ [&](const VectorType &rhs, VectorType &dst, double) { dst = rhs; };
kinsol.solve_jacobian_system = [&](const VectorType &,
const VectorType &,
const VectorType &rhs,
- VectorType & dst) -> int {
- dst = rhs;
- return 0;
- };
+ VectorType & dst) { dst = rhs; };
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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 {
+ kinsol.setup_jacobian = [](const VectorType &u, const VectorType &F) {
// 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 {
+ VectorType & dst) {
deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
<< ')' << std::endl;
J_inverse.invert(J);
J_inverse.vmult(dst, rhs);
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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 {
+ kinsol.setup_jacobian = [](const VectorType &u, const VectorType &F) {
// 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_with_jacobian = [](const VectorType &rhs,
- VectorType & dst,
- const double /*tolerance*/) -> int {
- deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
- << ')' << std::endl;
+ kinsol.solve_with_jacobian =
+ [](const VectorType &rhs, VectorType &dst, const double /*tolerance*/) {
+ 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;
+ // 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<double> 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<double> 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<double> J_inverse(2, 2);
- J_inverse.invert(J);
+ FullMatrix<double> J_inverse(2, 2);
+ J_inverse.invert(J);
- J_inverse.vmult(dst, rhs);
-
- return 0;
- };
+ J_inverse.vmult(dst, rhs);
+ };
VectorType v(N);
v(0) = 0.5;
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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<double> J_inverse(2, 2);
kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
- const VectorType &F) -> int {
+ const VectorType &F) {
deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ')'
<< std::endl;
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 {
+ VectorType & dst) {
deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
<< ')' << std::endl;
J_inverse.vmult(dst, rhs);
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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<double> J_inverse(2, 2);
kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
- const VectorType &F) -> int {
+ const VectorType &F) {
deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ')'
<< std::endl;
J(1, 1) = -std::cos(u[0] - u[1]) + 2;
J_inverse.invert(J);
-
- return 0;
};
kinsol.solve_with_jacobian = [&J_inverse](const VectorType &rhs,
VectorType & dst,
- const double /*tolerance*/) -> int {
+ const double /*tolerance*/) {
deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
<< ')' << std::endl;
J_inverse.vmult(dst, rhs);
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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<double> J_inverse(2, 2);
kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
- const VectorType &F) -> int {
+ const VectorType &F) {
// 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] << ')'
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 {
+ VectorType & dst) {
deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
<< ')' << std::endl;
J_inverse.vmult(dst, rhs);
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [](const VectorType &u, VectorType &F) {
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 {
+ kinsol.iteration_function = [](const VectorType &u, VectorType &F) {
// 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<double> J_inverse(2, 2);
kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
- const VectorType &F) -> int {
+ const VectorType &F) {
// 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] << ')'
J(1, 1) = -std::cos(u[0] - u[1]) + 2;
J_inverse.invert(J);
-
- return 0;
};
kinsol.solve_with_jacobian = [&J_inverse](const VectorType &rhs,
VectorType & dst,
- const double /*tolerance*/) -> int {
+ const double /*tolerance*/) {
deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
<< ')' << std::endl;
J_inverse.vmult(dst, rhs);
-
- return 0;
};
VectorType v(N);
kinsol.reinit_vector = [N](VectorType &v) { v.reinit(N); };
- kinsol.residual = [&](const VectorType &u, VectorType &F) -> int {
+ kinsol.residual = [&](const VectorType &u, VectorType &F) {
deallog << "Computing residual at " << u[0] << std::endl;
if ((u[0] < -10) || (u[0] > 20))
{
deallog << "Reporting recoverable failure." << std::endl;
- return 1;
+ throw RecoverableUserCallbackError();
}
F.reinit(u);
F[0] = std::atan(u[0]) - 0.5;
-
- return 0;
};
double J_inverse;
kinsol.setup_jacobian = [&J_inverse](const VectorType &u,
- const VectorType &F) -> int {
+ const VectorType &F) {
deallog << "Setting up Jacobian system at u=" << u[0] << std::endl;
const double J = 1. / (1 + u[0] * u[0]);
J_inverse = 1. / J;
-
- return 0;
};
- kinsol.solve_with_jacobian =
- [&](const VectorType &rhs, VectorType &dst, double) -> int {
- dst[0] = J_inverse * rhs[0];
- return 0;
- };
+ kinsol.solve_with_jacobian = [&](const VectorType &rhs,
+ VectorType & dst,
+ double) { dst[0] = J_inverse * rhs[0]; };
VectorType v(N);
v[0] = 10;
[&](const Vector<double> &evaluation_point,
Vector<double> & residual) {
compute_residual(evaluation_point, residual);
-
- return 0;
};
nonlinear_solver.setup_jacobian =
[&](const Vector<double> ¤t_u,
const Vector<double> & /*current_f*/) {
compute_and_factorize_jacobian(current_u);
-
- return 0;
};
nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
Vector<double> & dst,
const double tolerance) {
this->solve(rhs, dst, tolerance);
-
- return 0;
};
nonlinear_solver.solve(current_solution);
MinimalSurfaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run();
- return 0;
}