From 5993a2a56543d6f740c28179db960b7382b2f749 Mon Sep 17 00:00:00 2001 From: JonathanMatthewsFCAS Date: Wed, 10 Oct 2018 21:32:59 -0700 Subject: [PATCH] Added test for project_boundary_values_div_conforming usage on a distributed triangulation. --- tests/mpi/project_bv_div_conf.cc | 283 ++++++++++++++++++ ...v_div_conf.with_p4est=true.mpirun=2.output | 15 + ...v_div_conf.with_p4est=true.mpirun=4.output | 23 ++ 3 files changed, 321 insertions(+) create mode 100644 tests/mpi/project_bv_div_conf.cc create mode 100644 tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=2.output create mode 100644 tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=4.output diff --git a/tests/mpi/project_bv_div_conf.cc b/tests/mpi/project_bv_div_conf.cc new file mode 100644 index 0000000000..c27328143b --- /dev/null +++ b/tests/mpi/project_bv_div_conf.cc @@ -0,0 +1,283 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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. +// +// --------------------------------------------------------------------- + + +// Check that project_boundary_values_div_conforming works for a +// distributed triangulation + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +namespace ResFlow +{ + using namespace dealii; + + + template + class FluxBoundaryValues : public Function + { + public: + FluxBoundaryValues() + : Function(dim) + {} + + virtual void + vector_value(const Point &p, Vector &value) const override; + + virtual ~FluxBoundaryValues(){}; + }; + + template + void + FluxBoundaryValues::vector_value(const Point &p, + Vector & bdry_flux) const + { + Assert(bdry_flux.size() == dim, + ExcDimensionMismatch(bdry_flux.size(), dim)); + + const double alpha = 0.3; + const double beta = 1; + + bdry_flux(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2; + bdry_flux(1) = alpha * p[0] * p[1]; + } + + template + class ResFlowProblem + { + public: + ResFlowProblem(const unsigned int degree); + ~ResFlowProblem(); + + void + run(); + + private: + void + make_grid(); + void + setup_system(); + + const unsigned int degree; + + MPI_Comm mpi_communicator; + + FESystem fe; + parallel::distributed::Triangulation triangulation; + DoFHandler dof_handler; + + std::vector owned_partitioning; + std::vector relevant_partitioning; + + AffineConstraints constraints; + + ConditionalOStream pcout; + }; + + template + ResFlowProblem::ResFlowProblem(const unsigned int degree) + : degree(degree) + , mpi_communicator(MPI_COMM_WORLD) + , fe(FE_RaviartThomas(degree), 1, FE_DGQ(degree), 1) + , triangulation(mpi_communicator) + , dof_handler(triangulation) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + {} + + template + ResFlowProblem::~ResFlowProblem() + { + dof_handler.clear(); + } + + template + void + ResFlowProblem::make_grid() + { + unsigned int nx = 1; + unsigned int ny = nx, nz = nx; + bool colorize = true; + + std::vector elem_dim; + if (dim == 2) + elem_dim = {nx, nz}; + else + elem_dim = {nx, ny, nz}; + + Point p1, p2; + + if (dim == 2) + { + p1 = Point(-1., -1.); + p2 = Point(1., 1.); + } + else + { + p1 = Point(-1., -1., -1.); + p2 = Point(1., 1., 1.); + } + + GridGenerator::subdivided_hyper_rectangle( + triangulation, elem_dim, p1, p2, colorize); + triangulation.refine_global(2); + } // end make_grid + + template + void + ResFlowProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + std::vector resflow_sub_blocks(2, 0); + resflow_sub_blocks[1] = 1; + // Renumber to yield block structure + DoFRenumbering::block_wise(dof_handler); + + std::vector dofs_per_block(2); + DoFTools::count_dofs_per_block(dof_handler, + dofs_per_block, + resflow_sub_blocks); + + const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() << " (" + << n_u << '+' << n_p << ')' << std::endl; + + // We split up the IndexSet for locally owned and locally relevant DoFs + // into two IndexSets based on how we want to create the block matrices + // and vectors. + owned_partitioning.resize(2); + owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_u); + owned_partitioning[1] = + dof_handler.locally_owned_dofs().get_view(n_u, n_u + n_p); + + pcout << "Number of Owned dofs: " << dof_handler.n_locally_owned_dofs() + << " (" << owned_partitioning[0].n_elements() << '+' + << owned_partitioning[1].n_elements() << ')' << std::endl + << std::endl; + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + relevant_partitioning.resize(2); + relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_u); + relevant_partitioning[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p); + + pcout << "Number of Relevant dofs: " << locally_relevant_dofs.n_elements() + << " (" << relevant_partitioning[0].n_elements() << '+' + << relevant_partitioning[1].n_elements() << ')' << std::endl + << std::endl; + + { + constraints.reinit(locally_relevant_dofs); + + FEValuesExtractors::Vector velocities(0); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + VectorTools::project_boundary_values_div_conforming( + dof_handler, 0, FluxBoundaryValues(), 0, constraints); + + deallog << "Constraints" << std::endl; + constraints.print(deallog.get_file_stream()); + constraints.close(); + } + } // end setup_system + + + template + void + ResFlowProblem::run() + { + pcout << "Running with " + << Utilities::MPI::n_mpi_processes(mpi_communicator) + << " MPI rank(s)..." << std::endl; + + make_grid(); + pcout << "Grid made." << std::endl; + setup_system(); + pcout << "System set up." << std::endl; + } + +} // end namespace ResFlow + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + const unsigned int dim = 2; + + try + { + using namespace dealii; + using namespace ResFlow; + + Assert(dim == 2 || dim == 3, ExcNotImplemented()); + const unsigned int press_order = 1; + ResFlowProblem resflow(press_order); + resflow.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + deallog << "OK" << std::endl; + return 0; +} diff --git a/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=2.output b/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..e0037702a2 --- /dev/null +++ b/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=2.output @@ -0,0 +1,15 @@ + +DEAL:0::Constraints + 0 = 0.468750 + 1 = -0.0162380 + 22 = 0.431250 + 23 = -0.00541266 +DEAL:0::OK + +DEAL:1::Constraints + 76 = 0.431250 + 77 = 0.00541266 + 94 = 0.468750 + 95 = 0.0162380 +DEAL:1::OK + diff --git a/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=4.output b/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..31a50baeae --- /dev/null +++ b/tests/mpi/project_bv_div_conf.with_p4est=true.mpirun=4.output @@ -0,0 +1,23 @@ + +DEAL:0::Constraints + 0 = 0.468750 + 1 = -0.0162380 + 22 = 0.431250 + 23 = -0.00541266 +DEAL:0::OK + +DEAL:1::Constraints +DEAL:1::OK + + +DEAL:2::Constraints + 76 = 0.431250 + 77 = 0.00541266 + 94 = 0.468750 + 95 = 0.0162380 +DEAL:2::OK + + +DEAL:3::Constraints +DEAL:3::OK + -- 2.39.5