]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some code to where it is actually used. 15179/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 4 May 2023 18:46:54 +0000 (12:46 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 4 May 2023 23:38:06 +0000 (17:38 -0600)
While there, just use '{...}' instead of 'if (true) {...}'.

tests/trilinos/nox_solver_01.cc

index ba1ff2acb14b23a18ff0d29f9daabe8b0c2ca78f..ab26231f83964616d12f0b78076f69fc590d7168 100644 (file)
 #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<Epetra_CrsMatrix *>(&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<VectorType> 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<VectorType> 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<Epetra_CrsMatrix *>(&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<NoxInterface> interface = Teuchos::rcp(new NoxInterface());
 
-      Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = interface;
-      Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = interface;
+    // setup linear system and group
 
-      Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> 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<NOX::Epetra::Group> 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<NoxInterface> 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<NOX::Epetra::Interface::Required> iReq = interface;
+    Teuchos::RCP<NOX::Epetra::Interface::Jacobian> 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<NOX::Epetra::LinearSystemAztecOO> linSys =
+      Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(
+        printParams, lsParams, iReq, iJac, A, InitialGuess));
 
-      Teuchos::RCP<NOX::Solver::Generic> solver =
-        NOX::Solver::buildSolver(group, check, non_linear_parameters);
-      auto status = solver->solve();
+    NOX::Epetra::Vector              noxInitGuess(InitialGuess, NOX::DeepCopy);
+    Teuchos::RCP<NOX::Epetra::Group> 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<NOX::Solver::Generic> 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;
+  }
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.