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
parseParameters68 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
main(int argc,char * argv[])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