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