1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2023 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // Authors: Peter Munch, Martin Kronbichler 15 // 16 // --------------------------------------------------------------------- 17 18 // deal.II includes 19 #include <deal.II/base/conditional_ostream.h> 20 #include <deal.II/base/mpi.h> 21 #include <deal.II/base/parameter_handler.h> 22 #include <deal.II/base/quadrature_lib.h> 23 24 #include <deal.II/distributed/shared_tria.h> 25 #include <deal.II/distributed/tria.h> 26 27 #include <deal.II/dofs/dof_handler.h> 28 #include <deal.II/dofs/dof_renumbering.h> 29 30 #include <deal.II/fe/fe_nothing.h> 31 #include <deal.II/fe/fe_q.h> 32 #include <deal.II/fe/fe_system.h> 33 #include <deal.II/fe/fe_tools.h> 34 #include <deal.II/fe/fe_values.h> 35 #include <deal.II/fe/mapping_q1.h> 36 37 #include <deal.II/grid/grid_generator.h> 38 39 #include <deal.II/lac/affine_constraints.h> 40 #include <deal.II/lac/precondition.h> 41 #include <deal.II/lac/solver_cg.h> 42 43 // boost 44 #include <boost/algorithm/string.hpp> 45 46 #include <sstream> 47 48 // include operators 49 #include "bps-ceed.h" 50 #include "bps-cpu.h" 51 52 // Test cases 53 //TESTARGS(name="BP1") --resource {ceed_resource} --bp BP1 --fe_degree 2 --print_timings 0 54 //TESTARGS(name="BP4") --resource {ceed_resource} --bp BP4 --fe_degree 3 --print_timings 0 55 56 /** 57 * Relevant parameters. 58 */ 59 struct Parameters 60 { 61 BPType bp = BPType::BP5; 62 unsigned int n_global_refinements = 1; 63 unsigned int fe_degree = 2; 64 bool print_timings = true; 65 std::string libCEED_resource = "/cpu/self"; 66 67 bool 68 parse(int argc, char *argv[]) 69 { 70 if (argc == 1 && (std::string(argv[0]) == "--help")) 71 { 72 std::cout << "Usage: ./bp [OPTION]..." << std::endl; 73 std::cout << std::endl; 74 std::cout << "--bp name of benchmark (BP1-BP6)" << std::endl; 75 std::cout << "--n_refinements number of refinements (0-)" << std::endl; 76 std::cout << "--fe_degree polynomial degree (1-)" << std::endl; 77 std::cout << "--print_timings name of benchmark (0, 1)" << std::endl; 78 std::cout << "--resource name of resource (e.g., /cpu/self/avx/blocked)" << std::endl; 79 80 return true; 81 } 82 83 AssertThrow(argc % 2 == 0, ExcInternalError()); 84 85 while (argc > 0) 86 { 87 std::string label(argv[0]); 88 89 if ("--bp" == label) 90 { 91 std::string bp_string(argv[1]); 92 93 if (bp_string == "BP1") 94 bp = BPType::BP1; 95 else if (bp_string == "BP2") 96 bp = BPType::BP2; 97 else if (bp_string == "BP3") 98 bp = BPType::BP3; 99 else if (bp_string == "BP4") 100 bp = BPType::BP4; 101 else if (bp_string == "BP5") 102 bp = BPType::BP5; 103 else if (bp_string == "BP6") 104 bp = BPType::BP6; 105 else 106 AssertThrow(false, ExcInternalError()); 107 } 108 else if ("--n_refinements" == label) 109 { 110 n_global_refinements = std::atoi(argv[1]); 111 } 112 else if ("--fe_degree" == label) 113 { 114 fe_degree = std::atoi(argv[1]); 115 } 116 else if ("--print_timings" == label) 117 { 118 print_timings = std::atoi(argv[1]); 119 } 120 else if ("--resource" == label) 121 { 122 libCEED_resource = std::string(argv[1]); 123 } 124 else 125 { 126 AssertThrow(false, ExcNotImplemented()); 127 } 128 129 130 argc -= 2; 131 argv += 2; 132 } 133 134 return false; 135 } 136 }; 137 138 139 140 int 141 main(int argc, char *argv[]) 142 { 143 Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); 144 145 Parameters params; 146 if (params.parse(argc - 1, argv + 1)) 147 return 0; 148 149 ConditionalOStream pout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0); 150 151 // configuration 152 const BPType bp = params.bp; 153 154 using Number = double; 155 using VectorType = LinearAlgebra::distributed::Vector<Number>; 156 const unsigned int dim = 2; 157 const unsigned int fe_degree = params.fe_degree; 158 const unsigned int n_q_points = (bp <= BPType::BP4) ? (fe_degree + 2) : (fe_degree + 1); 159 const unsigned int n_refinements = params.n_global_refinements; 160 const unsigned int n_components = 161 (bp == BPType::BP1 || bp == BPType::BP3 || bp == BPType::BP5) ? 1 : dim; 162 163 // create mapping, quadrature, fe, mesh, ... 164 MappingQ1<dim> mapping; 165 QGauss<dim> quadrature(n_q_points); 166 FESystem<dim> fe(FE_Q<dim>(fe_degree), n_components); 167 168 #ifdef DEAL_II_WITH_P4EST 169 parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD); 170 #else 171 Triangulation<dim> tria; 172 #endif 173 174 GridGenerator::hyper_cube(tria); 175 tria.refine_global(n_refinements); 176 177 DoFHandler<dim> dof_handler(tria); 178 dof_handler.distribute_dofs(fe); 179 180 DoFRenumbering::support_point_wise(dof_handler); 181 182 AffineConstraints<Number> constraints; 183 184 if (!(bp == BPType::BP1 || bp == BPType::BP2)) 185 { 186 // for stiffness matrix 187 DoFTools::make_zero_boundary_constraints(dof_handler, constraints); 188 constraints.close(); 189 } 190 191 const auto test = [&](const std::string &label, const auto &op) { 192 (void)label; 193 194 // initialize vector 195 VectorType u, v; 196 op.initialize_dof_vector(u); 197 op.initialize_dof_vector(v); 198 u = 1.0; 199 200 constraints.set_zero(u); 201 202 // perform matrix-vector product 203 op.vmult(v, u); 204 205 // create solver 206 ReductionControl reduction_control(100, 1e-20, 1e-6); 207 208 // create preconditioner 209 DiagonalMatrix<VectorType> diagonal_matrix; 210 op.compute_inverse_diagonal(diagonal_matrix.get_vector()); 211 212 std::chrono::time_point<std::chrono::system_clock> now; 213 214 bool not_converged = false; 215 216 try 217 { 218 // solve problem 219 SolverCG<VectorType> solver(reduction_control); 220 now = std::chrono::system_clock::now(); 221 solver.solve(op, v, u, diagonal_matrix); 222 } 223 catch (const SolverControl::NoConvergence &) 224 { 225 pout << "Error: solver failed to converge with" << std::endl; 226 not_converged = true; 227 } 228 229 230 const auto time = 231 std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now() - now) 232 .count() / 233 1e9; 234 235 236 if (params.print_timings || not_converged) 237 { 238 pout << label << ": " << reduction_control.last_step() << " " << v.l2_norm() << " " 239 << (params.print_timings ? time : 0.0) << std::endl; 240 } 241 }; 242 243 // create and test the libCEED operator 244 OperatorCeed<dim, Number> op_ceed( 245 mapping, dof_handler, constraints, quadrature, bp, params.libCEED_resource); 246 test("ceed", op_ceed); 247 248 // create and test a native deal.II operator 249 OperatorDealii<dim, Number> op_dealii(mapping, dof_handler, constraints, quadrature, bp); 250 test("dealii", op_dealii); 251 } 252