From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 4 May 2023 19:54:22 +0000 (-0600)
Subject: Add a version of step-77 solved with NOX.
X-Git-Tag: v9.5.0-rc1~264^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15180%2Fhead;p=dealii.git

Add a version of step-77 solved with NOX.
---

diff --git a/tests/trilinos/step-77-with-nox.cc b/tests/trilinos/step-77-with-nox.cc
new file mode 100644
index 0000000000..cfc82d5436
--- /dev/null
+++ b/tests/trilinos/step-77-with-nox.cc
@@ -0,0 +1,470 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2021 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+
+// step-77 for the test suite - verifies KINSOL.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/trilinos/nox.h>
+
+#include "../tests.h"
+
+namespace Step77
+{
+  using namespace dealii;
+
+  template <int dim>
+  class MinimalSurfaceProblem
+  {
+  public:
+    MinimalSurfaceProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system(const bool initial_step);
+    void
+    solve(const Vector<double> &rhs,
+          Vector<double> &      solution,
+          const double          tolerance);
+    void
+    refine_mesh();
+    void
+    output_results(const unsigned int refinement_cycle);
+    void
+    set_boundary_values();
+    void
+    compute_and_factorize_jacobian(const Vector<double> &evaluation_point);
+    void
+    compute_residual(const Vector<double> &evaluation_point,
+                     Vector<double> &      residual);
+
+    Triangulation<dim> triangulation;
+
+    DoFHandler<dim> dof_handler;
+    FE_Q<dim>       fe;
+
+    AffineConstraints<double> hanging_node_constraints;
+
+    SparsityPattern                      sparsity_pattern;
+    SparseMatrix<double>                 jacobian_matrix;
+    std::unique_ptr<SparseDirectUMFPACK> jacobian_matrix_factorization;
+
+    Vector<double> current_solution;
+
+    // unlike step-77, don't output timer information.
+  };
+
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+  };
+
+
+  template <int dim>
+  double
+  BoundaryValues<dim>::value(const Point<dim> &p,
+                             const unsigned int /*component*/) const
+  {
+    return std::sin(2 * numbers::PI * (p[0] + p[1]));
+  }
+
+
+
+  template <int dim>
+  MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
+    : dof_handler(triangulation)
+    , fe(1)
+  {}
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
+  {
+    if (initial_step)
+      {
+        dof_handler.distribute_dofs(fe);
+        current_solution.reinit(dof_handler.n_dofs());
+
+        hanging_node_constraints.clear();
+        DoFTools::make_hanging_node_constraints(dof_handler,
+                                                hanging_node_constraints);
+        hanging_node_constraints.close();
+      }
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+    hanging_node_constraints.condense(dsp);
+
+    sparsity_pattern.copy_from(dsp);
+    jacobian_matrix.reinit(sparsity_pattern);
+    jacobian_matrix_factorization.reset();
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
+    const Vector<double> &evaluation_point)
+  {
+    deallog << "  Computing Jacobian matrix" << std::endl;
+
+    {
+      const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+      jacobian_matrix = 0;
+
+      FEValues<dim> fe_values(fe,
+                              quadrature_formula,
+                              update_gradients | update_quadrature_points |
+                                update_JxW_values);
+
+      const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+      const unsigned int n_q_points    = quadrature_formula.size();
+
+      FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+      std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+      std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          cell_matrix = 0;
+
+          fe_values.reinit(cell);
+
+          fe_values.get_function_gradients(evaluation_point,
+                                           evaluation_point_gradients);
+
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              const double coeff =
+                1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                      evaluation_point_gradients[q]);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) +=
+                      (((fe_values.shape_grad(i, q)    // ((\nabla \phi_i
+                         * coeff                       //   * a_n
+                         * fe_values.shape_grad(j, q)) //   * \nabla \phi_j)
+                        -                              //  -
+                        (fe_values.shape_grad(i, q)    //  (\nabla \phi_i
+                         * coeff * coeff * coeff       //   * a_n^3
+                         *
+                         (fe_values.shape_grad(j, q)       //   * (\nabla \phi_j
+                          * evaluation_point_gradients[q]) //      * \nabla u_n)
+                         * evaluation_point_gradients[q])) //   * \nabla u_n)))
+                       * fe_values.JxW(q));                // * dx
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          hanging_node_constraints.distribute_local_to_global(cell_matrix,
+                                                              local_dof_indices,
+                                                              jacobian_matrix);
+        }
+
+      std::map<types::global_dof_index, double> boundary_values;
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               0,
+                                               Functions::ZeroFunction<dim>(),
+                                               boundary_values);
+      Vector<double> dummy_solution(dof_handler.n_dofs());
+      Vector<double> dummy_rhs(dof_handler.n_dofs());
+      MatrixTools::apply_boundary_values(boundary_values,
+                                         jacobian_matrix,
+                                         dummy_solution,
+                                         dummy_rhs);
+    }
+
+    {
+      jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
+      jacobian_matrix_factorization->factorize(jacobian_matrix);
+    }
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_residual(
+    const Vector<double> &evaluation_point,
+    Vector<double> &      residual)
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+    FEValues<dim>     fe_values(fe,
+                            quadrature_formula,
+                            update_gradients | update_quadrature_points |
+                              update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    Vector<double>              cell_residual(dofs_per_cell);
+    std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        cell_residual = 0;
+        fe_values.reinit(cell);
+
+        fe_values.get_function_gradients(evaluation_point,
+                                         evaluation_point_gradients);
+
+
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const double coeff =
+              1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                    evaluation_point_gradients[q]);
+
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              cell_residual(i) +=
+                (fe_values.shape_grad(i, q)      // \nabla \phi_i
+                 * coeff                         // * a_n
+                 * evaluation_point_gradients[q] // * \nabla u_n
+                 * fe_values.JxW(q));            // * dx
+          }
+
+        cell->get_dof_indices(local_dof_indices);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          residual(local_dof_indices[i]) += cell_residual(i);
+      }
+
+    hanging_node_constraints.condense(residual);
+
+    for (const types::global_dof_index i :
+         DoFTools::extract_boundary_dofs(dof_handler))
+      residual(i) = 0;
+
+    for (const types::global_dof_index i :
+         DoFTools::extract_hanging_node_dofs(dof_handler))
+      residual(i) = 0;
+
+    deallog << "  Computed residual vector with norm " << residual.l2_norm()
+            << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
+                                    Vector<double> &      solution,
+                                    const double /*tolerance*/)
+  {
+    jacobian_matrix_factorization->vmult(solution, rhs);
+
+    hanging_node_constraints.distribute(solution);
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::refine_mesh()
+  {
+    Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+
+    KellyErrorEstimator<dim>::estimate(
+      dof_handler,
+      QGauss<dim - 1>(fe.degree + 1),
+      std::map<types::boundary_id, const Function<dim> *>(),
+      current_solution,
+      estimated_error_per_cell);
+
+    GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3,
+                                                    0.03);
+
+    triangulation.prepare_coarsening_and_refinement();
+
+    SolutionTransfer<dim> solution_transfer(dof_handler);
+    solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
+
+    triangulation.execute_coarsening_and_refinement();
+
+    dof_handler.distribute_dofs(fe);
+
+    Vector<double> tmp(dof_handler.n_dofs());
+    solution_transfer.interpolate(current_solution, tmp);
+    current_solution = std::move(tmp);
+
+    hanging_node_constraints.clear();
+
+    DoFTools::make_hanging_node_constraints(dof_handler,
+                                            hanging_node_constraints);
+    hanging_node_constraints.close();
+
+    hanging_node_constraints.distribute(current_solution);
+
+    set_boundary_values();
+
+    setup_system(/*initial_step=*/false);
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::set_boundary_values()
+  {
+    std::map<types::global_dof_index, double> boundary_values;
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             BoundaryValues<dim>(),
+                                             boundary_values);
+    for (const auto &boundary_value : boundary_values)
+      current_solution(boundary_value.first) = boundary_value.second;
+
+    hanging_node_constraints.distribute(current_solution);
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::output_results(
+    const unsigned int refinement_cycle)
+  {
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(current_solution, "solution");
+    data_out.build_patches();
+
+    deallog << "Current refinement cycle = " << refinement_cycle << std::endl;
+    data_out.write_vtk(deallog.get_file_stream());
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::run()
+  {
+    GridGenerator::hyper_ball(triangulation);
+    triangulation.refine_global(2);
+
+    setup_system(/*initial_step=*/true);
+    set_boundary_values();
+
+    for (unsigned int refinement_cycle = 0; refinement_cycle < 3;
+         ++refinement_cycle)
+      {
+        deallog << "Mesh refinement step " << refinement_cycle << std::endl;
+
+        if (refinement_cycle != 0)
+          refine_mesh();
+
+        const double target_tolerance = 1e-10 * std::pow(0.1, refinement_cycle);
+        deallog << "  Target_tolerance: " << target_tolerance << std::endl
+                << std::endl;
+
+        {
+          typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData
+            additional_data;
+          additional_data.abs_tol = target_tolerance;
+
+          TrilinosWrappers::NOXSolver<Vector<double>> nonlinear_solver(
+            additional_data);
+
+          nonlinear_solver.residual =
+            [&](const Vector<double> &evaluation_point,
+                Vector<double> &      residual) {
+              compute_residual(evaluation_point, residual);
+
+              return 0;
+            };
+
+          nonlinear_solver.setup_jacobian =
+            [&](const Vector<double> &current_u) {
+              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);
+        }
+
+        //        output_results(refinement_cycle);
+      }
+  }
+} // namespace Step77
+
+
+int
+main()
+{
+  initlog();
+  using namespace Step77;
+
+  MinimalSurfaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run();
+  return 0;
+}
diff --git a/tests/trilinos/step-77-with-nox.with_trilinos_with_nox=true.output b/tests/trilinos/step-77-with-nox.with_trilinos_with_nox=true.output
new file mode 100644
index 0000000000..0625ba83c8
--- /dev/null
+++ b/tests/trilinos/step-77-with-nox.with_trilinos_with_nox=true.output
@@ -0,0 +1,37 @@
+
+DEAL::Mesh refinement step 0
+DEAL::  Target_tolerance: 1.00000e-10
+DEAL::
+DEAL::  Computed residual vector with norm 0.867975
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.212073
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.0189603
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.000314854
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 1.14048e-07
+DEAL::Mesh refinement step 1
+DEAL::  Target_tolerance: 1.00000e-11
+DEAL::
+DEAL::  Computed residual vector with norm 0.246479
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.0359168
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.00208943
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 8.26297e-06
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 1.43448e-10
+DEAL::Mesh refinement step 2
+DEAL::  Target_tolerance: 1.00000e-12
+DEAL::
+DEAL::  Computed residual vector with norm 0.504420
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.171685
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.0241205
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 0.000550324
+DEAL::  Computing Jacobian matrix
+DEAL::  Computed residual vector with norm 5.73528e-07