From ae48085a58ad70346aa4b4db0e29e48ecb68ca0f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 4 May 2023 12:46:54 -0600 Subject: [PATCH] Move some code to where it is actually used. While there, just use '{...}' instead of 'if (true) {...}'. --- tests/trilinos/nox_solver_01.cc | 253 ++++++++++++++++---------------- 1 file changed, 125 insertions(+), 128 deletions(-) diff --git a/tests/trilinos/nox_solver_01.cc b/tests/trilinos/nox_solver_01.cc index ba1ff2acb1..ab26231f83 100644 --- a/tests/trilinos/nox_solver_01.cc +++ b/tests/trilinos/nox_solver_01.cc @@ -48,37 +48,6 @@ #include "../tests.h" -class NoxInterface : public NOX::Epetra::Interface::Required, - public NOX::Epetra::Interface::Jacobian -{ -public: - bool - computeF(const Epetra_Vector & x, - Epetra_Vector & f, - NOX::Epetra::Interface::Required::FillType F) override - { - (void)F; - - f[0] = x[0] * x[0]; - - return true; - } - - bool - computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac) override - { - auto jac = dynamic_cast(&Jac); - - AssertThrow(jac, ExcNotImplemented()); - - jac->PutScalar(2.0 * x[0]); - - return true; - } -}; - - - int main(int argc, char **argv) { @@ -114,125 +83,153 @@ main(int argc, char **argv) /* * First check: Use the deal.II-based wrappers for the test. */ - if (true) - { - // set up solver - TrilinosWrappers::NOXSolver solver(additional_data, - non_linear_parameters); - - // ... helper functions - double J = 0.0; - - solver.residual = [](const auto &src, auto &dst) { - // compute residual - dst[0] = src[0] * src[0]; - return 0; - }; - - solver.setup_jacobian = [&](const auto &src) { - // compute Jacobian - J = 2.0 * src[0]; - return 0; - }; - - solver.apply_jacobian = [&](const auto &src, auto &dst) { - // solve with Jacobian - dst[0] = src[0] * J; - return 0; - }; - - solver.solve_with_jacobian = [&](const auto &src, auto &dst, const auto) { - // solve with Jacobian - dst[0] = src[0] / J; - return 0; - }; - - // initial guess - VectorType solution(1); - solution[0] = 2.0; - - // solve with the given initial guess - solver.solve(solution); - - deallog << "The solution is: " << solution[0] << std::endl; - } + { + // set up solver + TrilinosWrappers::NOXSolver solver(additional_data, + non_linear_parameters); + + // ... helper functions + double J = 0.0; + + solver.residual = [](const auto &src, auto &dst) { + // compute residual + dst[0] = src[0] * src[0]; + return 0; + }; + + solver.setup_jacobian = [&](const auto &src) { + // compute Jacobian + J = 2.0 * src[0]; + return 0; + }; + + solver.apply_jacobian = [&](const auto &src, auto &dst) { + // solve with Jacobian + dst[0] = src[0] * J; + return 0; + }; + + solver.solve_with_jacobian = [&](const auto &src, auto &dst, const auto) { + // solve with Jacobian + dst[0] = src[0] / J; + return 0; + }; + + // initial guess + VectorType solution(1); + solution[0] = 2.0; + + // solve with the given initial guess + solver.solve(solution); + + deallog << "The solution is: " << solution[0] << std::endl; + } /* * Second check: Run the same test through the native NOX interfaces. */ - if (true) + { + class NoxInterface : public NOX::Epetra::Interface::Required, + public NOX::Epetra::Interface::Jacobian { - // convert data structures to Epetra structures - IndexSet is(1); - is.add_index(0); + public: + bool + computeF(const Epetra_Vector & x, + Epetra_Vector & f, + NOX::Epetra::Interface::Required::FillType F) override + { + (void)F; + + f[0] = x[0] * x[0]; + + return true; + } + + bool + computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac) override + { + auto jac = dynamic_cast(&Jac); - TrilinosWrappers::MPI::Vector solution(is, MPI_COMM_WORLD); - solution[0] = 2.0; + AssertThrow(jac, ExcNotImplemented()); - TrilinosWrappers::SparsityPattern dsp(is, MPI_COMM_WORLD); - dsp.add(0, 0); - dsp.compress(); + jac->PutScalar(2.0 * x[0]); - TrilinosWrappers::SparseMatrix system_matrix; - system_matrix.reinit(dsp); + return true; + } + }; - // setup linear system and group + // convert data structures to Epetra structures + IndexSet is(1); + is.add_index(0); - Epetra_Vector InitialGuess(Copy, solution.trilinos_vector(), 0); + TrilinosWrappers::MPI::Vector solution(is, MPI_COMM_WORLD); + solution[0] = 2.0; - auto A = - Teuchos::rcp(new Epetra_CrsMatrix(system_matrix.trilinos_matrix())); + TrilinosWrappers::SparsityPattern dsp(is, MPI_COMM_WORLD); + dsp.add(0, 0); + dsp.compress(); - Teuchos::ParameterList &printParams = - non_linear_parameters->sublist("Printing"); - Teuchos::ParameterList &lsParams = - non_linear_parameters->sublist("Newton").sublist("Linear Solver"); + TrilinosWrappers::SparseMatrix system_matrix; + system_matrix.reinit(dsp); - Teuchos::RCP interface = Teuchos::rcp(new NoxInterface()); - Teuchos::RCP iReq = interface; - Teuchos::RCP iJac = interface; + // setup linear system and group - Teuchos::RCP linSys = - Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO( - printParams, lsParams, iReq, iJac, A, InitialGuess)); + Epetra_Vector InitialGuess(Copy, solution.trilinos_vector(), 0); - NOX::Epetra::Vector noxInitGuess(InitialGuess, NOX::DeepCopy); - Teuchos::RCP group = Teuchos::rcp( - new NOX::Epetra::Group(printParams, iReq, noxInitGuess, linSys)); + auto A = + Teuchos::rcp(new Epetra_CrsMatrix(system_matrix.trilinos_matrix())); - // setup solver control - auto check = - Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR)); + Teuchos::ParameterList &printParams = + non_linear_parameters->sublist("Printing"); + Teuchos::ParameterList &lsParams = + non_linear_parameters->sublist("Newton").sublist("Linear Solver"); - if (additional_data.abs_tol > 0.0) - { - const auto additional_data_norm_f_abs = - Teuchos::rcp(new NOX::StatusTest::NormF(additional_data.abs_tol)); - check->addStatusTest(additional_data_norm_f_abs); - } + Teuchos::RCP interface = Teuchos::rcp(new NoxInterface()); - if (additional_data.rel_tol > 0.0) - { - const auto additional_data_norm_f_rel = Teuchos::rcp( - new NOX::StatusTest::RelativeNormF(additional_data.rel_tol)); - check->addStatusTest(additional_data_norm_f_rel); - } + Teuchos::RCP iReq = interface; + Teuchos::RCP iJac = interface; - if (additional_data.max_iter > 0) - { - const auto additional_data_max_iterations = Teuchos::rcp( - new NOX::StatusTest::MaxIters(additional_data.max_iter)); - check->addStatusTest(additional_data_max_iterations); - } + Teuchos::RCP linSys = + Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO( + printParams, lsParams, iReq, iJac, A, InitialGuess)); - Teuchos::RCP solver = - NOX::Solver::buildSolver(group, check, non_linear_parameters); - auto status = solver->solve(); + NOX::Epetra::Vector noxInitGuess(InitialGuess, NOX::DeepCopy); + Teuchos::RCP group = Teuchos::rcp( + new NOX::Epetra::Group(printParams, iReq, noxInitGuess, linSys)); - AssertThrow(status == NOX::StatusTest::Converged, ExcInternalError()); + // setup solver control + auto check = + Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR)); - deallog << "The solution is: " << group->getX().norm() << std::endl; - } + if (additional_data.abs_tol > 0.0) + { + const auto additional_data_norm_f_abs = + Teuchos::rcp(new NOX::StatusTest::NormF(additional_data.abs_tol)); + check->addStatusTest(additional_data_norm_f_abs); + } + + if (additional_data.rel_tol > 0.0) + { + const auto additional_data_norm_f_rel = Teuchos::rcp( + new NOX::StatusTest::RelativeNormF(additional_data.rel_tol)); + check->addStatusTest(additional_data_norm_f_rel); + } + + if (additional_data.max_iter > 0) + { + const auto additional_data_max_iterations = + Teuchos::rcp(new NOX::StatusTest::MaxIters(additional_data.max_iter)); + check->addStatusTest(additional_data_max_iterations); + } + + Teuchos::RCP solver = + NOX::Solver::buildSolver(group, check, non_linear_parameters); + auto status = solver->solve(); + + AssertThrow(status == NOX::StatusTest::Converged, ExcInternalError()); + + deallog << "The solution is: " << group->getX().norm() << std::endl; + } } -- 2.39.5