From: Wolfgang Bangerth 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=1d98287a0cf55b2df57f92d88d6582e748ef7b5d;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 +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + +namespace Step77 +{ + using namespace dealii; + + template + class MinimalSurfaceProblem + { + public: + MinimalSurfaceProblem(); + void + run(); + + private: + void + setup_system(const bool initial_step); + void + solve(const Vector &rhs, + Vector & 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 &evaluation_point); + void + compute_residual(const Vector &evaluation_point, + Vector & residual); + + Triangulation triangulation; + + DoFHandler dof_handler; + FE_Q fe; + + AffineConstraints hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix jacobian_matrix; + std::unique_ptr jacobian_matrix_factorization; + + Vector current_solution; + + // unlike step-77, don't output timer information. + }; + + + + template + class BoundaryValues : public Function + { + public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + + template + double + BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const + { + return std::sin(2 * numbers::PI * (p[0] + p[1])); + } + + + + template + MinimalSurfaceProblem::MinimalSurfaceProblem() + : dof_handler(triangulation) + , fe(1) + {} + + + + template + void + MinimalSurfaceProblem::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 + void + MinimalSurfaceProblem::compute_and_factorize_jacobian( + const Vector &evaluation_point) + { + deallog << " Computing Jacobian matrix" << std::endl; + + { + const QGauss quadrature_formula(fe.degree + 1); + + jacobian_matrix = 0; + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + boundary_values); + Vector dummy_solution(dof_handler.n_dofs()); + Vector dummy_rhs(dof_handler.n_dofs()); + MatrixTools::apply_boundary_values(boundary_values, + jacobian_matrix, + dummy_solution, + dummy_rhs); + } + + { + jacobian_matrix_factorization = std::make_unique(); + jacobian_matrix_factorization->factorize(jacobian_matrix); + } + } + + + + template + void + MinimalSurfaceProblem::compute_residual( + const Vector &evaluation_point, + Vector & residual) + { + const QGauss quadrature_formula(fe.degree + 1); + FEValues 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 cell_residual(dofs_per_cell); + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 + void + MinimalSurfaceProblem::solve(const Vector &rhs, + Vector & solution, + const double /*tolerance*/) + { + jacobian_matrix_factorization->vmult(solution, rhs); + + hanging_node_constraints.distribute(solution); + } + + + + template + void + MinimalSurfaceProblem::refine_mesh() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + 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 solution_transfer(dof_handler); + solution_transfer.prepare_for_coarsening_and_refinement(current_solution); + + triangulation.execute_coarsening_and_refinement(); + + dof_handler.distribute_dofs(fe); + + Vector 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 + void + MinimalSurfaceProblem::set_boundary_values() + { + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + boundary_values); + for (const auto &boundary_value : boundary_values) + current_solution(boundary_value.first) = boundary_value.second; + + hanging_node_constraints.distribute(current_solution); + } + + + + template + void + MinimalSurfaceProblem::output_results( + const unsigned int refinement_cycle) + { + DataOut 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 + void + MinimalSurfaceProblem::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>::AdditionalData + additional_data; + additional_data.abs_tol = target_tolerance; + + TrilinosWrappers::NOXSolver> nonlinear_solver( + additional_data); + + nonlinear_solver.residual = + [&](const Vector &evaluation_point, + Vector & residual) { + compute_residual(evaluation_point, residual); + + return 0; + }; + + nonlinear_solver.setup_jacobian = + [&](const Vector ¤t_u) { + compute_and_factorize_jacobian(current_u); + + return 0; + }; + + nonlinear_solver.solve_with_jacobian = [&](const Vector &rhs, + Vector & 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