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/dofs/dof_tools.h> 20 21 #include <deal.II/fe/mapping.h> 22 23 #include <deal.II/lac/la_parallel_vector.h> 24 25 #include <deal.II/matrix_free/fe_evaluation.h> 26 #include <deal.II/matrix_free/matrix_free.h> 27 #include <deal.II/matrix_free/shape_info.h> 28 #include <deal.II/matrix_free/tools.h> 29 30 // libCEED includes 31 #include <ceed.h> 32 #include <ceed/backend.h> 33 34 // QFunction source 35 #include "bps-qfunctions.h" 36 37 using namespace dealii; 38 39 /** 40 * BP types. For more details, see https://ceed.exascaleproject.org/bps/. 41 */ 42 enum class BPType : unsigned int 43 { 44 BP1, 45 BP2, 46 BP3, 47 BP4, 48 BP5, 49 BP6 50 }; 51 52 53 54 /** 55 * Struct storing relevant information regarding each BP. 56 */ 57 struct BPInfo 58 { 59 BPInfo(const BPType type, const int dim, const int fe_degree) 60 : type(type) 61 , dim(dim) 62 , fe_degree(fe_degree) 63 { 64 if (type == BPType::BP1) 65 type_string = "BP1"; 66 else if (type == BPType::BP2) 67 type_string = "BP2"; 68 else if (type == BPType::BP3) 69 type_string = "BP3"; 70 else if (type == BPType::BP4) 71 type_string = "BP4"; 72 else if (type == BPType::BP5) 73 type_string = "BP5"; 74 else if (type == BPType::BP6) 75 type_string = "BP6"; 76 77 this->n_q_points_1d = (type <= BPType::BP4) ? (fe_degree + 2) : (fe_degree + 1); 78 79 this->n_components = 80 (type == BPType::BP1 || type == BPType::BP3 || type == BPType::BP5) ? 1 : dim; 81 } 82 83 84 BPType type; 85 std::string type_string; 86 unsigned int dim; 87 unsigned int fe_degree; 88 unsigned int n_q_points_1d; 89 unsigned int n_components; 90 }; 91 92 93 94 /** 95 * Base class of operators. 96 */ 97 template <typename Number> 98 class OperatorBase 99 { 100 public: 101 /** 102 * deal.II vector type 103 */ 104 using VectorType = LinearAlgebra::distributed::Vector<Number>; 105 106 /** 107 * Initialize vector. 108 */ 109 virtual void 110 reinit() = 0; 111 112 /** 113 * Perform matrix-vector product 114 */ 115 virtual void 116 vmult(VectorType &dst, const VectorType &src) const = 0; 117 118 /** 119 * Initialize vector. 120 */ 121 virtual void 122 initialize_dof_vector(VectorType &vec) const = 0; 123 124 /** 125 * Compute inverse of diagonal. 126 */ 127 virtual void 128 compute_inverse_diagonal(VectorType &diagonal) const = 0; 129 }; 130 131 132 /** 133 * Operator implementation using libCEED. 134 */ 135 template <int dim, typename Number> 136 class OperatorCeed : public OperatorBase<Number> 137 { 138 public: 139 using VectorType = typename OperatorBase<Number>::VectorType; 140 141 /** 142 * Constructor. 143 */ 144 OperatorCeed(const Mapping<dim> &mapping, 145 const DoFHandler<dim> &dof_handler, 146 const AffineConstraints<Number> &constraints, 147 const Quadrature<dim> &quadrature, 148 const BPType &bp, 149 const std::string &resource) 150 : mapping(mapping) 151 , dof_handler(dof_handler) 152 , constraints(constraints) 153 , quadrature(quadrature) 154 , bp(bp) 155 , resource(resource) 156 { 157 reinit(); 158 } 159 160 /** 161 * Destructor. 162 */ 163 ~OperatorCeed() 164 { 165 CeedVectorDestroy(&src_ceed); 166 CeedVectorDestroy(&dst_ceed); 167 CeedOperatorDestroy(&op_apply); 168 CeedDestroy(&ceed); 169 } 170 171 /** 172 * Initialized internal data structures, particularly, libCEED. 173 */ 174 void 175 reinit() override 176 { 177 CeedVector q_data; 178 CeedBasis sol_basis; 179 CeedElemRestriction sol_restriction; 180 CeedElemRestriction q_data_restriction; 181 BuildContext build_ctx_data; 182 CeedQFunctionContext build_ctx; 183 CeedQFunction qf_apply; 184 185 const auto &tria = dof_handler.get_triangulation(); 186 const auto &fe = dof_handler.get_fe(); 187 188 const auto n_components = fe.n_components(); 189 190 if (bp == BPType::BP1 || bp == BPType::BP3 || bp == BPType::BP5) 191 { 192 AssertThrow(n_components == 1, ExcInternalError()); 193 } 194 else 195 { 196 AssertThrow(n_components == dim, ExcInternalError()); 197 } 198 199 // 1) create CEED instance -> "MatrixFree" 200 const char *ceed_spec = resource.c_str(); 201 CeedInit(ceed_spec, &ceed); 202 203 // 2) create shape functions -> "ShapeInfo" 204 const unsigned int fe_degree = fe.tensor_degree(); 205 const unsigned int n_q_points = quadrature.get_tensor_basis()[0].size(); 206 { 207 const dealii::internal::MatrixFreeFunctions::ShapeInfo<double> shape_info(quadrature, fe, 0); 208 const auto &shape_data = shape_info.get_shape_data(); 209 std::vector<CeedScalar> q_ref_1d; 210 for (const auto q : shape_data.quadrature.get_points()) 211 q_ref_1d.push_back(q(0)); 212 213 // transpose bases for compatibility with restriction 214 std::vector<CeedScalar> interp_1d(shape_data.shape_values.size()); 215 std::vector<CeedScalar> grad_1d(shape_data.shape_gradients.size()); 216 for (unsigned int i = 0; i < n_q_points; ++i) 217 for (unsigned int j = 0; j < fe_degree + 1; ++j) 218 { 219 interp_1d[j + i * (fe_degree + 1)] = shape_data.shape_values[j * n_q_points + i]; 220 grad_1d[j + i * (fe_degree + 1)] = shape_data.shape_gradients[j * n_q_points + i]; 221 } 222 223 CeedBasisCreateTensorH1(ceed, 224 dim, 225 n_components, 226 fe_degree + 1, 227 n_q_points, 228 interp_1d.data(), 229 grad_1d.data(), 230 q_ref_1d.data(), 231 quadrature.get_tensor_basis()[0].get_weights().data(), 232 &sol_basis); 233 } 234 235 // 3) create restriction matrix -> DoFInfo 236 unsigned int n_local_active_cells = 0; 237 238 for (const auto &cell : dof_handler.active_cell_iterators()) 239 if (cell->is_locally_owned()) 240 n_local_active_cells++; 241 242 partitioner = 243 std::make_shared<Utilities::MPI::Partitioner>(dof_handler.locally_owned_dofs(), 244 DoFTools::extract_locally_active_dofs( 245 dof_handler), 246 dof_handler.get_communicator()); 247 248 std::vector<CeedInt> indices; 249 indices.reserve(n_local_active_cells * fe.n_dofs_per_cell() / n_components); 250 251 const auto dof_mapping = FETools::lexicographic_to_hierarchic_numbering<dim>(fe_degree); 252 253 std::vector<types::global_dof_index> local_indices(fe.n_dofs_per_cell()); 254 255 for (const auto &cell : dof_handler.active_cell_iterators()) 256 if (cell->is_locally_owned()) 257 { 258 cell->get_dof_indices(local_indices); 259 260 for (const auto i : dof_mapping) 261 indices.emplace_back( 262 partitioner->global_to_local(local_indices[fe.component_to_system_index(0, i)])); 263 } 264 265 CeedElemRestrictionCreate(ceed, 266 n_local_active_cells, 267 fe.n_dofs_per_cell() / n_components, 268 n_components, 269 1, 270 this->extended_local_size(), 271 CEED_MEM_HOST, 272 CEED_COPY_VALUES, 273 indices.data(), 274 &sol_restriction); 275 276 // 4) create mapping -> MappingInfo 277 const unsigned int n_components_metric = (bp <= BPType::BP2) ? 1 : (dim * (dim + 1) / 2); 278 279 this->weights = compute_metric_data(ceed, mapping, tria, quadrature, bp); 280 281 strides = {{1, 282 static_cast<int>(quadrature.size()), 283 static_cast<int>(quadrature.size() * n_components_metric)}}; 284 CeedVectorCreate(ceed, weights.size(), &q_data); 285 CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_USE_POINTER, weights.data()); 286 CeedElemRestrictionCreateStrided(ceed, 287 n_local_active_cells, 288 quadrature.size(), 289 n_components_metric, 290 weights.size(), 291 strides.data(), 292 &q_data_restriction); 293 294 build_ctx_data.dim = dim; 295 build_ctx_data.space_dim = dim; 296 297 CeedQFunctionContextCreate(ceed, &build_ctx); 298 CeedQFunctionContextSetData( 299 build_ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(build_ctx_data), &build_ctx_data); 300 301 // 5) create q operation 302 if (bp == BPType::BP1) 303 CeedQFunctionCreateInterior(ceed, 1, f_apply_mass, f_apply_mass_loc, &qf_apply); 304 else if (bp == BPType::BP2) 305 CeedQFunctionCreateInterior(ceed, 1, f_apply_mass_vec, f_apply_mass_vec_loc, &qf_apply); 306 else if (bp == BPType::BP3 || bp == BPType::BP5) 307 CeedQFunctionCreateInterior(ceed, 1, f_apply_poisson, f_apply_poisson_loc, &qf_apply); 308 else if (bp == BPType::BP4 || bp == BPType::BP6) 309 CeedQFunctionCreateInterior(ceed, 1, f_apply_poisson_vec, f_apply_poisson_vec_loc, &qf_apply); 310 else 311 AssertThrow(false, ExcInternalError()); 312 313 if (bp <= BPType::BP2) 314 CeedQFunctionAddInput(qf_apply, "u", n_components, CEED_EVAL_INTERP); 315 else 316 CeedQFunctionAddInput(qf_apply, "u", dim * n_components, CEED_EVAL_GRAD); 317 318 CeedQFunctionAddInput(qf_apply, "qdata", n_components_metric, CEED_EVAL_NONE); 319 320 if (bp <= BPType::BP2) 321 CeedQFunctionAddOutput(qf_apply, "v", n_components, CEED_EVAL_INTERP); 322 else 323 CeedQFunctionAddOutput(qf_apply, "v", dim * n_components, CEED_EVAL_GRAD); 324 325 CeedQFunctionSetContext(qf_apply, build_ctx); 326 327 // 6) put everything together 328 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply); 329 330 CeedOperatorSetField(op_apply, "u", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE); 331 CeedOperatorSetField(op_apply, "qdata", q_data_restriction, CEED_BASIS_NONE, q_data); 332 CeedOperatorSetField(op_apply, "v", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE); 333 334 // 7) libCEED vectors 335 CeedElemRestrictionCreateVector(sol_restriction, &src_ceed, NULL); 336 CeedElemRestrictionCreateVector(sol_restriction, &dst_ceed, NULL); 337 338 // 8) cleanup 339 CeedVectorDestroy(&q_data); 340 CeedElemRestrictionDestroy(&q_data_restriction); 341 CeedElemRestrictionDestroy(&sol_restriction); 342 CeedBasisDestroy(&sol_basis); 343 CeedQFunctionContextDestroy(&build_ctx); 344 CeedQFunctionDestroy(&qf_apply); 345 } 346 347 /** 348 * Perform matrix-vector product. 349 */ 350 void 351 vmult(VectorType &dst, const VectorType &src) const override 352 { 353 // communicate: update ghost values 354 src.update_ghost_values(); 355 356 // pass memory buffers to libCEED 357 VectorTypeCeed x(src_ceed); 358 VectorTypeCeed y(dst_ceed); 359 x.import_array(src, CEED_MEM_HOST); 360 y.import_array(dst, CEED_MEM_HOST); 361 362 // apply operator 363 CeedOperatorApply(op_apply, x(), y(), CEED_REQUEST_IMMEDIATE); 364 365 // pull arrays back to deal.II 366 x.take_array(); 367 y.take_array(); 368 369 // communicate: compress 370 src.zero_out_ghost_values(); 371 dst.compress(VectorOperation::add); 372 373 // apply constraints: we assume homogeneous DBC 374 constraints.set_zero(dst); 375 } 376 377 /** 378 * Initialized vector. 379 */ 380 void 381 initialize_dof_vector(VectorType &vec) const override 382 { 383 vec.reinit(partitioner); 384 } 385 386 /** 387 * Compute inverse of diagonal. 388 */ 389 void 390 compute_inverse_diagonal(VectorType &diagonal) const override 391 { 392 this->initialize_dof_vector(diagonal); 393 394 // pass memory buffer to libCEED 395 VectorTypeCeed y(dst_ceed); 396 y.import_array(diagonal, CEED_MEM_HOST); 397 398 CeedOperatorLinearAssembleDiagonal(op_apply, y(), CEED_REQUEST_IMMEDIATE); 399 400 // pull array back to deal.II 401 y.take_array(); 402 403 diagonal.compress(VectorOperation::add); 404 405 for (auto &i : diagonal) 406 i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; 407 } 408 409 private: 410 /** 411 * Wrapper around a deal.II vector to create a libCEED vector view. 412 */ 413 class VectorTypeCeed 414 { 415 public: 416 /** 417 * Constructor. 418 */ 419 VectorTypeCeed(const CeedVector &vec_orig) 420 { 421 vec_ceed = NULL; 422 CeedVectorReferenceCopy(vec_orig, &vec_ceed); 423 } 424 425 /** 426 * Return libCEED vector view. 427 */ 428 CeedVector & 429 operator()() 430 { 431 return vec_ceed; 432 } 433 434 /** 435 * Set deal.II memory in libCEED vector. 436 */ 437 void 438 import_array(const VectorType &vec, const CeedMemType space) 439 { 440 mem_space = space; 441 CeedVectorSetArray(vec_ceed, mem_space, CEED_USE_POINTER, vec.get_values()); 442 } 443 444 /** 445 * Sync memory from device to host. 446 */ 447 void 448 sync_array() 449 { 450 CeedVectorSyncArray(vec_ceed, mem_space); 451 } 452 453 /** 454 * Take previously set deal.II array from libCEED vector 455 */ 456 void 457 take_array() 458 { 459 CeedScalar *ptr; 460 CeedVectorTakeArray(vec_ceed, mem_space, &ptr); 461 } 462 463 /** 464 * Destructor: destroy vector view. 465 */ 466 ~VectorTypeCeed() 467 { 468 bool has_array; 469 CeedVectorHasBorrowedArrayOfType(vec_ceed, mem_space, &has_array); 470 if (has_array) 471 { 472 CeedScalar *ptr; 473 CeedVectorTakeArray(vec_ceed, mem_space, &ptr); 474 } 475 CeedVectorDestroy(&vec_ceed); 476 } 477 478 private: 479 /** 480 * libCEED vector view. 481 */ 482 CeedMemType mem_space; 483 CeedVector vec_ceed; 484 }; 485 486 /** 487 * Number of locally active DoFs. 488 */ 489 unsigned int 490 extended_local_size() const 491 { 492 return partitioner->locally_owned_size() + partitioner->n_ghost_indices(); 493 } 494 495 /** 496 * Compute metric data: Jacobian, ... 497 */ 498 static std::vector<double> 499 compute_metric_data(const Ceed &ceed, 500 const Mapping<dim> &mapping, 501 const Triangulation<dim> &tria, 502 const Quadrature<dim> &quadrature, 503 const BPType bp) 504 { 505 std::vector<double> weights; 506 507 if (false) 508 { 509 FE_Nothing<dim> dummy_fe; 510 FEValues<dim> fe_values(mapping, dummy_fe, quadrature, update_JxW_values); 511 512 for (const auto &cell : tria.active_cell_iterators()) 513 if (cell->is_locally_owned()) 514 { 515 fe_values.reinit(cell); 516 517 for (const auto q : fe_values.quadrature_point_indices()) 518 weights.emplace_back(fe_values.JxW(q)); 519 } 520 521 return weights; 522 } 523 524 CeedBasis geo_basis; 525 CeedVector q_data; 526 CeedElemRestriction q_data_restriction; 527 CeedVector node_coords; 528 CeedElemRestriction geo_restriction; 529 CeedQFunctionContext build_ctx; 530 CeedQFunction qf_build; 531 CeedOperator op_build; 532 533 const unsigned int n_q_points = quadrature.get_tensor_basis()[0].size(); 534 535 const unsigned int n_components = (bp <= BPType::BP2) ? 1 : (dim * (dim + 1) / 2); 536 537 const auto mapping_q = dynamic_cast<const MappingQ<dim> *>(&mapping); 538 539 AssertThrow(mapping_q, ExcMessage("Wrong mapping!")); 540 541 const unsigned int fe_degree = mapping_q->get_degree(); 542 543 FE_Q<dim> geo_fe(fe_degree); 544 545 { 546 const dealii::internal::MatrixFreeFunctions::ShapeInfo<double> shape_info(quadrature, 547 geo_fe, 548 0); 549 const auto &shape_data = shape_info.get_shape_data(); 550 std::vector<CeedScalar> q_ref_1d; 551 for (const auto q : shape_data.quadrature.get_points()) 552 q_ref_1d.push_back(q(0)); 553 554 // transpose bases for compatibility with restriction 555 std::vector<CeedScalar> interp_1d(shape_data.shape_values.size()); 556 std::vector<CeedScalar> grad_1d(shape_data.shape_gradients.size()); 557 for (unsigned int i = 0; i < n_q_points; ++i) 558 for (unsigned int j = 0; j < fe_degree + 1; ++j) 559 { 560 interp_1d[j + i * (fe_degree + 1)] = shape_data.shape_values[j * n_q_points + i]; 561 grad_1d[j + i * (fe_degree + 1)] = shape_data.shape_gradients[j * n_q_points + i]; 562 } 563 564 CeedBasisCreateTensorH1(ceed, 565 dim, 566 dim, 567 fe_degree + 1, 568 n_q_points, 569 interp_1d.data(), 570 grad_1d.data(), 571 q_ref_1d.data(), 572 quadrature.get_tensor_basis()[0].get_weights().data(), 573 &geo_basis); 574 } 575 576 unsigned int n_local_active_cells = 0; 577 578 for (const auto &cell : tria.active_cell_iterators()) 579 if (cell->is_locally_owned()) 580 n_local_active_cells++; 581 582 std::vector<double> geo_support_points; 583 std::vector<CeedInt> geo_indices; 584 585 DoFHandler<dim> geo_dof_handler(tria); 586 geo_dof_handler.distribute_dofs(geo_fe); 587 588 const auto geo_partitioner = 589 std::make_shared<Utilities::MPI::Partitioner>(geo_dof_handler.locally_owned_dofs(), 590 DoFTools::extract_locally_active_dofs( 591 geo_dof_handler), 592 geo_dof_handler.get_communicator()); 593 594 geo_indices.reserve(n_local_active_cells * geo_fe.n_dofs_per_cell()); 595 596 const auto dof_mapping = FETools::lexicographic_to_hierarchic_numbering<dim>(fe_degree); 597 598 FEValues<dim> fe_values(mapping, 599 geo_fe, 600 geo_fe.get_unit_support_points(), 601 update_quadrature_points); 602 603 std::vector<types::global_dof_index> local_indices(geo_fe.n_dofs_per_cell()); 604 605 const unsigned int n_points = 606 geo_partitioner->locally_owned_size() + geo_partitioner->n_ghost_indices(); 607 608 geo_support_points.resize(dim * n_points); 609 610 for (const auto &cell : geo_dof_handler.active_cell_iterators()) 611 if (cell->is_locally_owned()) 612 { 613 fe_values.reinit(cell); 614 cell->get_dof_indices(local_indices); 615 616 for (const auto i : dof_mapping) 617 { 618 const auto index = geo_partitioner->global_to_local(local_indices[i]); 619 geo_indices.emplace_back(index * dim); 620 621 const auto point = fe_values.quadrature_point(i); 622 623 for (unsigned int d = 0; d < dim; ++d) 624 geo_support_points[index * dim + d] = point[d]; 625 } 626 } 627 628 weights.resize(n_local_active_cells * quadrature.size() * n_components); 629 630 CeedInt strides[3] = {1, 631 static_cast<int>(quadrature.size()), 632 static_cast<int>(quadrature.size() * n_components)}; 633 634 CeedVectorCreate(ceed, weights.size(), &q_data); 635 CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_USE_POINTER, weights.data()); 636 CeedElemRestrictionCreateStrided(ceed, 637 n_local_active_cells, 638 quadrature.size(), 639 n_components, 640 weights.size(), 641 strides, 642 &q_data_restriction); 643 644 CeedVectorCreate(ceed, geo_support_points.size(), &node_coords); 645 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, geo_support_points.data()); 646 647 CeedElemRestrictionCreate(ceed, 648 n_local_active_cells, 649 geo_fe.n_dofs_per_cell(), 650 dim, 651 1, 652 geo_support_points.size(), 653 CEED_MEM_HOST, 654 CEED_COPY_VALUES, 655 geo_indices.data(), 656 &geo_restriction); 657 658 BuildContext build_ctx_data; 659 build_ctx_data.dim = dim; 660 build_ctx_data.space_dim = dim; 661 662 CeedQFunctionContextCreate(ceed, &build_ctx); 663 CeedQFunctionContextSetData( 664 build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(build_ctx_data), &build_ctx_data); 665 666 // 5) create q operation 667 if (bp <= BPType::BP2) 668 CeedQFunctionCreateInterior(ceed, 1, f_build_mass, f_build_mass_loc, &qf_build); 669 else 670 CeedQFunctionCreateInterior(ceed, 1, f_build_poisson, f_build_poisson_loc, &qf_build); 671 672 CeedQFunctionAddInput(qf_build, "geo", dim * dim, CEED_EVAL_GRAD); 673 CeedQFunctionAddInput(qf_build, "weights", 1, CEED_EVAL_WEIGHT); 674 CeedQFunctionAddOutput(qf_build, "qdata", n_components, CEED_EVAL_NONE); 675 CeedQFunctionSetContext(qf_build, build_ctx); 676 677 // 6) put everything together 678 CeedOperatorCreate(ceed, qf_build, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_build); 679 CeedOperatorSetField(op_build, "geo", geo_restriction, geo_basis, CEED_VECTOR_ACTIVE); 680 CeedOperatorSetField( 681 op_build, "weights", CEED_ELEMRESTRICTION_NONE, geo_basis, CEED_VECTOR_NONE); 682 CeedOperatorSetField( 683 op_build, "qdata", q_data_restriction, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 684 685 CeedOperatorApply(op_build, node_coords, q_data, CEED_REQUEST_IMMEDIATE); 686 687 CeedVectorDestroy(&node_coords); 688 CeedVectorSyncArray(q_data, CEED_MEM_HOST); 689 CeedVectorDestroy(&q_data); 690 CeedElemRestrictionDestroy(&geo_restriction); 691 CeedElemRestrictionDestroy(&q_data_restriction); 692 CeedBasisDestroy(&geo_basis); 693 CeedQFunctionContextDestroy(&build_ctx); 694 CeedQFunctionDestroy(&qf_build); 695 CeedOperatorDestroy(&op_build); 696 697 return weights; 698 } 699 700 /** 701 * Mapping object passed to the constructor. 702 */ 703 const Mapping<dim> &mapping; 704 705 /** 706 * DoFHandler object passed to the constructor. 707 */ 708 const DoFHandler<dim> &dof_handler; 709 710 /** 711 * Constraints object passed to the constructor. 712 */ 713 const AffineConstraints<Number> &constraints; 714 715 /** 716 * Quadrature rule object passed to the constructor. 717 */ 718 const Quadrature<dim> &quadrature; 719 720 /** 721 * Selected BP. 722 */ 723 const BPType bp; 724 725 /** 726 * Resource name. 727 */ 728 const std::string resource; 729 730 /** 731 * Partitioner for distributed vectors. 732 */ 733 std::shared_ptr<Utilities::MPI::Partitioner> partitioner; 734 735 /** 736 * libCEED data structures. 737 */ 738 Ceed ceed; 739 std::vector<double> weights; 740 std::array<CeedInt, 3> strides; 741 CeedVector src_ceed; 742 CeedVector dst_ceed; 743 CeedOperator op_apply; 744 745 /** 746 * Temporal (tempral) vectors. 747 * 748 * @note Only needed for multiple components. 749 */ 750 mutable VectorType src_tmp; 751 mutable VectorType dst_tmp; 752 }; 753 754 755 756 template <int dim, typename Number> 757 class OperatorDealii : public OperatorBase<Number> 758 { 759 public: 760 using VectorType = typename OperatorBase<Number>::VectorType; 761 762 /** 763 * Constructor. 764 */ 765 OperatorDealii(const Mapping<dim> &mapping, 766 const DoFHandler<dim> &dof_handler, 767 const AffineConstraints<Number> &constraints, 768 const Quadrature<dim> &quadrature, 769 const BPType &bp) 770 : mapping(mapping) 771 , dof_handler(dof_handler) 772 , constraints(constraints) 773 , quadrature(quadrature) 774 , bp(bp) 775 { 776 reinit(); 777 } 778 779 /** 780 * Destructor. 781 */ 782 ~OperatorDealii() = default; 783 784 /** 785 * Initialized internal data structures, particularly, MatrixFree. 786 */ 787 void 788 reinit() override 789 { 790 // configure MatrixFree 791 typename MatrixFree<dim, Number>::AdditionalData additional_data; 792 additional_data.tasks_parallel_scheme = 793 MatrixFree<dim, Number>::AdditionalData::TasksParallelScheme::none; 794 795 // create MatrixFree 796 matrix_free.reinit(mapping, dof_handler, constraints, quadrature, additional_data); 797 } 798 799 /** 800 * Matrix-vector product. 801 */ 802 void 803 vmult(VectorType &dst, const VectorType &src) const override 804 { 805 if (dof_handler.get_fe().n_components() == 1) 806 { 807 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<1>, this, dst, src, true); 808 } 809 else 810 { 811 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 812 813 matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<dim>, this, dst, src, true); 814 } 815 } 816 817 /** 818 * Initialize vector. 819 */ 820 void 821 initialize_dof_vector(VectorType &vec) const override 822 { 823 matrix_free.initialize_dof_vector(vec); 824 } 825 826 /** 827 * Compute inverse of diagonal. 828 */ 829 void 830 compute_inverse_diagonal(VectorType &diagonal) const override 831 { 832 this->initialize_dof_vector(diagonal); 833 834 if (dof_handler.get_fe().n_components() == 1) 835 { 836 MatrixFreeTools::compute_diagonal(matrix_free, 837 diagonal, 838 &OperatorDealii::do_cell_integral_local<1>, 839 this); 840 } 841 else 842 { 843 AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError()); 844 845 MatrixFreeTools::compute_diagonal(matrix_free, 846 diagonal, 847 &OperatorDealii::do_cell_integral_local<dim>, 848 this); 849 } 850 851 for (auto &i : diagonal) 852 i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0; 853 } 854 855 private: 856 /** 857 * Cell integral without vector access. 858 */ 859 template <int n_components> 860 void 861 do_cell_integral_local(FEEvaluation<dim, -1, 0, n_components, Number> &phi) const 862 { 863 if (bp <= BPType::BP2) // mass matrix 864 { 865 phi.evaluate(EvaluationFlags::values); 866 for (const auto q : phi.quadrature_point_indices()) 867 phi.submit_value(phi.get_value(q), q); 868 phi.integrate(EvaluationFlags::values); 869 } 870 else // Poisson operator 871 { 872 phi.evaluate(EvaluationFlags::gradients); 873 for (const auto q : phi.quadrature_point_indices()) 874 phi.submit_gradient(phi.get_gradient(q), q); 875 phi.integrate(EvaluationFlags::gradients); 876 } 877 } 878 879 /** 880 * Cell integral on a range of cells. 881 */ 882 template <int n_components> 883 void 884 do_cell_integral_range(const MatrixFree<dim, Number> &matrix_free, 885 VectorType &dst, 886 const VectorType &src, 887 const std::pair<unsigned int, unsigned int> &range) const 888 { 889 FEEvaluation<dim, -1, 0, n_components, Number> phi(matrix_free, range); 890 891 for (unsigned cell = range.first; cell < range.second; ++cell) 892 { 893 phi.reinit(cell); 894 phi.read_dof_values(src); // read source vector 895 do_cell_integral_local(phi); // cell integral 896 phi.distribute_local_to_global(dst); // write to destination vector 897 } 898 } 899 900 /** 901 * Mapping object passed to the constructor. 902 */ 903 const Mapping<dim> &mapping; 904 905 /** 906 * DoFHandler object passed to the constructor. 907 */ 908 const DoFHandler<dim> &dof_handler; 909 910 /** 911 * Constraints object passed to the constructor. 912 */ 913 const AffineConstraints<Number> &constraints; 914 915 /** 916 * Quadrature rule object passed to the constructor. 917 */ 918 const Quadrature<dim> &quadrature; 919 920 /** 921 * Selected BP. 922 */ 923 const BPType bp; 924 925 /** 926 * MatrixFree object. 927 */ 928 MatrixFree<dim, Number> matrix_free; 929 }; 930