DOLFINx
DOLFINx C++ interface
Loading...
Searching...
No Matches
petsc.h
1// Copyright (C) 2018-2021 Garth N. Wells
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "Form.h"
10#include "assembler.h"
11#include "utils.h"
12#include <concepts>
13#include <dolfinx/la/petsc.h>
14#include <map>
15#include <memory>
16#include <petscmat.h>
17#include <petscvec.h>
18#include <span>
19#include <utility>
20#include <vector>
21
22namespace dolfinx::common
23{
24class IndexMap;
25}
26
27namespace dolfinx::fem
28{
29template <dolfinx::scalar T, std::floating_point U>
30class DirichletBC;
31
33namespace petsc
34{
41template <std::floating_point T>
43 std::string type = std::string())
44{
46 pattern.finalize();
47 return la::petsc::create_matrix(a.mesh()->comm(), pattern, type);
48}
49
58template <std::floating_point T>
60 const std::vector<std::vector<const Form<PetscScalar, T>*>>& a,
61 std::string type = std::string())
62{
63 // Extract and check row/column ranges
64 std::array<std::vector<std::shared_ptr<const FunctionSpace<T>>>, 2> V
66 std::array<std::vector<int>, 2> bs_dofs;
67 for (std::size_t i = 0; i < 2; ++i)
68 {
69 for (auto& _V : V[i])
70 bs_dofs[i].push_back(_V->dofmap()->bs());
71 }
72
73 // Build sparsity pattern for each block
74 std::shared_ptr<const mesh::Mesh<T>> mesh;
75 std::vector<std::vector<std::unique_ptr<la::SparsityPattern>>> patterns(
76 V[0].size());
77 for (std::size_t row = 0; row < V[0].size(); ++row)
78 {
79 for (std::size_t col = 0; col < V[1].size(); ++col)
80 {
81 if (const Form<PetscScalar, T>* form = a[row][col]; form)
82 {
83 patterns[row].push_back(std::make_unique<la::SparsityPattern>(
85 if (!mesh)
86 mesh = form->mesh();
87 }
88 else
89 patterns[row].push_back(nullptr);
90 }
91 }
92
93 if (!mesh)
94 throw std::runtime_error("Could not find a Mesh.");
95
96 // Compute offsets for the fields
97 std::array<std::vector<std::pair<
98 std::reference_wrapper<const common::IndexMap>, int>>,
99 2>
100 maps;
101 for (std::size_t d = 0; d < 2; ++d)
102 {
103 for (auto& space : V[d])
104 {
105 maps[d].emplace_back(*space->dofmap()->index_map,
106 space->dofmap()->index_map_bs());
107 }
108 }
109
110 // Create merged sparsity pattern
111 std::vector<std::vector<const la::SparsityPattern*>> p(V[0].size());
112 for (std::size_t row = 0; row < V[0].size(); ++row)
113 for (std::size_t col = 0; col < V[1].size(); ++col)
114 p[row].push_back(patterns[row][col].get());
115
116 la::SparsityPattern pattern(mesh->comm(), p, maps, bs_dofs);
117 pattern.finalize();
118
119 // FIXME: Add option to pass customised local-to-global map to PETSc
120 // Mat constructor
121
122 // Initialise matrix
123 Mat A = la::petsc::create_matrix(mesh->comm(), pattern, type);
124
125 // Create row and column local-to-global maps (field0, field1, field2,
126 // etc), i.e. ghosts of field0 appear before owned indices of field1
127 std::array<std::vector<PetscInt>, 2> _maps;
128 for (int d = 0; d < 2; ++d)
129 {
130 // TODO: Index map concatenation has already been computed inside
131 // the SparsityPattern constructor, but we also need it here to
132 // build the PETSc local-to-global map. Compute outside and pass
133 // into SparsityPattern constructor.
134 // TODO: avoid concatenating the same maps twice in case that V[0]
135 // == V[1].
136
137 const std::vector<
138 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& map
139 = maps[d];
140 std::vector<PetscInt>& _map = _maps[d];
141
142 // Concatenate the block index map in the row and column directions
143 const auto [rank_offset, local_offset, ghosts, _]
145 for (std::size_t f = 0; f < map.size(); ++f)
146 {
147 auto offset = local_offset[f];
148 const common::IndexMap& imap = map[f].first.get();
149 int bs = map[f].second;
150 for (std::int32_t i = 0; i < bs * imap.size_local(); ++i)
151 _map.push_back(i + rank_offset + offset);
152 for (std::int32_t i = 0; i < bs * imap.num_ghosts(); ++i)
153 _map.push_back(ghosts[f][i]);
154 }
155 }
156
157 // Create PETSc local-to-global map/index sets and attach to matrix
158 ISLocalToGlobalMapping petsc_local_to_global0;
159 ISLocalToGlobalMappingCreate(MPI_COMM_SELF, 1, _maps[0].size(),
160 _maps[0].data(), PETSC_COPY_VALUES,
161 &petsc_local_to_global0);
162 if (V[0] == V[1])
163 {
164 MatSetLocalToGlobalMapping(A, petsc_local_to_global0,
165 petsc_local_to_global0);
166 ISLocalToGlobalMappingDestroy(&petsc_local_to_global0);
167 }
168 else
169 {
170
171 ISLocalToGlobalMapping petsc_local_to_global1;
172 ISLocalToGlobalMappingCreate(MPI_COMM_SELF, 1, _maps[1].size(),
173 _maps[1].data(), PETSC_COPY_VALUES,
174 &petsc_local_to_global1);
175 MatSetLocalToGlobalMapping(A, petsc_local_to_global0,
176 petsc_local_to_global1);
177 ISLocalToGlobalMappingDestroy(&petsc_local_to_global0);
178 ISLocalToGlobalMappingDestroy(&petsc_local_to_global1);
179 }
180
181 return A;
182}
183
187template <std::floating_point T>
189 const std::vector<std::vector<const Form<PetscScalar, T>*>>& a,
190 const std::vector<std::vector<std::string>>& types)
191{
192 // Extract and check row/column ranges
194
195 std::vector<std::vector<std::string>> _types(
196 a.size(), std::vector<std::string>(a.front().size()));
197 if (!types.empty())
198 _types = types;
199
200 // Loop over each form and create matrix
201 int rows = a.size();
202 int cols = a.front().size();
203 std::vector<Mat> mats(rows * cols, nullptr);
204 std::shared_ptr<const mesh::Mesh<T>> mesh;
205 for (int i = 0; i < rows; ++i)
206 {
207 for (int j = 0; j < cols; ++j)
208 {
209 if (const Form<PetscScalar, T>* form = a[i][j]; form)
210 {
211 mats[i * cols + j] = create_matrix(*form, _types[i][j]);
212 mesh = form->mesh();
213 }
214 }
215 }
216
217 if (!mesh)
218 throw std::runtime_error("Could not find a Mesh.");
219
220 // Initialise block (MatNest) matrix
221 Mat A;
222 MatCreate(mesh->comm(), &A);
223 MatSetType(A, MATNEST);
224 MatNestSetSubMats(A, rows, nullptr, cols, nullptr, mats.data());
225 MatSetUp(A);
226
227 // De-reference Mat objects
228 for (std::size_t i = 0; i < mats.size(); ++i)
229 {
230 if (mats[i])
231 MatDestroy(&mats[i]);
232 }
233
234 return A;
235}
236
241 const std::vector<
242 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
243
246 const std::vector<
247 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
248
249// -- Vectors ----------------------------------------------------------------
250
263template <std::floating_point T>
265 Vec b, const Form<PetscScalar, T>& L,
266 std::span<const PetscScalar> constants,
267 const std::map<std::pair<IntegralType, int>,
268 std::pair<std::span<const PetscScalar>, int>>& coeffs)
269{
270 Vec b_local;
271 VecGhostGetLocalForm(b, &b_local);
272 PetscInt n = 0;
273 VecGetSize(b_local, &n);
274 PetscScalar* array = nullptr;
275 VecGetArray(b_local, &array);
276 std::span<PetscScalar> _b(array, n);
277 fem::assemble_vector(_b, L, constants, coeffs);
278 VecRestoreArray(b_local, &array);
279 VecGhostRestoreLocalForm(b, &b_local);
280}
281
291template <std::floating_point T>
293{
294 Vec b_local;
295 VecGhostGetLocalForm(b, &b_local);
296 PetscInt n = 0;
297 VecGetSize(b_local, &n);
298 PetscScalar* array = nullptr;
299 VecGetArray(b_local, &array);
300 std::span<PetscScalar> _b(array, n);
302 VecRestoreArray(b_local, &array);
303 VecGhostRestoreLocalForm(b, &b_local);
304}
305
306// FIXME: clarify how x0 is used
307// FIXME: if bcs entries are set
308
309// FIXME: need to pass an array of Vec for x0?
310// FIXME: clarify zeroing of vector
311
327template <std::floating_point T>
329 Vec b, const std::vector<std::shared_ptr<const Form<PetscScalar, T>>>& a,
330 const std::vector<std::span<const PetscScalar>>& constants,
331 const std::vector<std::map<std::pair<IntegralType, int>,
332 std::pair<std::span<const PetscScalar>, int>>>&
333 coeffs,
334 const std::vector<
335 std::vector<std::shared_ptr<const DirichletBC<PetscScalar, T>>>>& bcs1,
336 const std::vector<Vec>& x0, PetscScalar scale)
337{
338 Vec b_local;
339 VecGhostGetLocalForm(b, &b_local);
340 PetscInt n = 0;
341 VecGetSize(b_local, &n);
342 PetscScalar* array = nullptr;
343 VecGetArray(b_local, &array);
344 std::span<PetscScalar> _b(array, n);
345
346 if (x0.empty())
347 fem::apply_lifting(_b, a, constants, coeffs, bcs1, {}, scale);
348 else
349 {
350 std::vector<std::span<const PetscScalar>> x0_ref;
351 std::vector<Vec> x0_local(a.size());
352 std::vector<const PetscScalar*> x0_array(a.size());
353 for (std::size_t i = 0; i < a.size(); ++i)
354 {
355 assert(x0[i]);
356 VecGhostGetLocalForm(x0[i], &x0_local[i]);
357 PetscInt n = 0;
358 VecGetSize(x0_local[i], &n);
359 VecGetArrayRead(x0_local[i], &x0_array[i]);
360 x0_ref.emplace_back(x0_array[i], n);
361 }
362
363 std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
364 fem::apply_lifting(_b, a, constants, coeffs, bcs1, x0_tmp, scale);
365
366 for (std::size_t i = 0; i < x0_local.size(); ++i)
367 {
368 VecRestoreArrayRead(x0_local[i], &x0_array[i]);
369 VecGhostRestoreLocalForm(x0[i], &x0_local[i]);
370 }
371 }
372
373 VecRestoreArray(b_local, &array);
374 VecGhostRestoreLocalForm(b, &b_local);
375}
376
377// FIXME: clarify how x0 is used
378// FIXME: if bcs entries are set
379
380// FIXME: need to pass an array of Vec for x0?
381// FIXME: clarify zeroing of vector
382
395template <std::floating_point T>
397 Vec b,
398 const std::vector<std::shared_ptr<const Form<PetscScalar, double>>>& a,
399 const std::vector<
400 std::vector<std::shared_ptr<const DirichletBC<PetscScalar, double>>>>&
401 bcs1,
402 const std::vector<Vec>& x0, PetscScalar scale)
403{
404 Vec b_local;
405 VecGhostGetLocalForm(b, &b_local);
406 PetscInt n = 0;
407 VecGetSize(b_local, &n);
408 PetscScalar* array = nullptr;
409 VecGetArray(b_local, &array);
410 std::span<PetscScalar> _b(array, n);
411
412 if (x0.empty())
413 fem::apply_lifting<PetscScalar>(_b, a, bcs1, {}, scale);
414 else
415 {
416 std::vector<std::span<const PetscScalar>> x0_ref;
417 std::vector<Vec> x0_local(a.size());
418 std::vector<const PetscScalar*> x0_array(a.size());
419 for (std::size_t i = 0; i < a.size(); ++i)
420 {
421 assert(x0[i]);
422 VecGhostGetLocalForm(x0[i], &x0_local[i]);
423 PetscInt n = 0;
424 VecGetSize(x0_local[i], &n);
425 VecGetArrayRead(x0_local[i], &x0_array[i]);
426 x0_ref.emplace_back(x0_array[i], n);
427 }
428
429 std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
430 fem::apply_lifting<PetscScalar>(_b, a, bcs1, x0_tmp, scale);
431
432 for (std::size_t i = 0; i < x0_local.size(); ++i)
433 {
434 VecRestoreArrayRead(x0_local[i], &x0_array[i]);
435 VecGhostRestoreLocalForm(x0[i], &x0_local[i]);
436 }
437 }
438
439 VecRestoreArray(b_local, &array);
440 VecGhostRestoreLocalForm(b, &b_local);
441}
442
443// -- Setting bcs ------------------------------------------------------------
444
445// FIXME: Move these function elsewhere?
446
447// FIXME: clarify x0
448// FIXME: clarify what happens with ghosts
449
453template <std::floating_point T>
455 Vec b,
456 const std::vector<std::shared_ptr<const DirichletBC<PetscScalar, T>>>& bcs,
457 const Vec x0, PetscScalar scale = 1)
458{
459 PetscInt n = 0;
460 VecGetLocalSize(b, &n);
461 PetscScalar* array = nullptr;
462 VecGetArray(b, &array);
463 std::span<PetscScalar> _b(array, n);
464 if (x0)
465 {
466 Vec x0_local;
467 VecGhostGetLocalForm(x0, &x0_local);
468 PetscInt n = 0;
469 VecGetSize(x0_local, &n);
470 const PetscScalar* array = nullptr;
471 VecGetArrayRead(x0_local, &array);
472 std::span<const PetscScalar> _x0(array, n);
473 fem::set_bc(_b, bcs, _x0, scale);
474 VecRestoreArrayRead(x0_local, &array);
475 VecGhostRestoreLocalForm(x0, &x0_local);
476 }
477 else
478 fem::set_bc(_b, bcs, scale);
479
480 VecRestoreArray(b, &array);
481}
482
483} // namespace petsc
484} // namespace dolfinx::fem
This class represents the distribution index arrays across processes. An index array is a contiguous ...
Definition IndexMap.h:94
std::int32_t num_ghosts() const noexcept
Number of ghost indices on this process.
Definition IndexMap.cpp:866
std::int32_t size_local() const noexcept
Number of indices owned by this process.
Definition IndexMap.cpp:868
Object for setting (strong) Dirichlet boundary conditions.
Definition DirichletBC.h:259
A representation of finite element variational forms.
Definition Form.h:120
Sparsity pattern data structure that can be used to initialize sparse matrices. After assembly,...
Definition SparsityPattern.h:26
void finalize()
Finalize sparsity pattern and communicate off-process entries.
Definition SparsityPattern.cpp:226
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
std::tuple< std::int64_t, std::vector< std::int32_t >, std::vector< std::vector< std::int64_t > >, std::vector< std::vector< int > > > stack_index_maps(const std::vector< std::pair< std::reference_wrapper< const IndexMap >, int > > &maps)
Compute layout data and ghost indices for a stacked (concatenated) index map, i.e....
Definition IndexMap.cpp:583
void set_bc(Vec b, const std::vector< std::shared_ptr< const DirichletBC< PetscScalar, T > > > &bcs, const Vec x0, PetscScalar scale=1)
Set bc values in owned (local) part of the PETSc vector, multiplied by 'scale'. The vectors b and x0 ...
Definition petsc.h:454
void apply_lifting(Vec b, const std::vector< std::shared_ptr< const Form< PetscScalar, T > > > &a, const std::vector< std::span< const PetscScalar > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const PetscScalar >, int > > > &coeffs, const std::vector< std::vector< std::shared_ptr< const DirichletBC< PetscScalar, T > > > > &bcs1, const std::vector< Vec > &x0, PetscScalar scale)
Modify RHS vector to account for Dirichlet boundary conditions.
Definition petsc.h:328
Mat create_matrix(const Form< PetscScalar, T > &a, std::string type=std::string())
Create a matrix.
Definition petsc.h:42
void assemble_vector(Vec b, const Form< PetscScalar, T > &L, std::span< const PetscScalar > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const PetscScalar >, int > > &coeffs)
Assemble linear form into an already allocated PETSc vector.
Definition petsc.h:264
Mat create_matrix_block(const std::vector< std::vector< const Form< PetscScalar, T > * > > &a, std::string type=std::string())
Initialise a monolithic matrix for an array of bilinear forms.
Definition petsc.h:59
Mat create_matrix_nest(const std::vector< std::vector< const Form< PetscScalar, T > * > > &a, const std::vector< std::vector< std::string > > &types)
Create nested (MatNest) matrix.
Definition petsc.h:188
Vec create_vector_block(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Initialise monolithic vector. Vector is not zeroed.
Definition petsc.cpp:15
Vec create_vector_nest(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Create nested (VecNest) vector. Vector is not zeroed.
Definition petsc.cpp:60
Finite element method functionality.
Definition assemble_matrix_impl.h:26
void assemble_vector(std::span< T > b, const Form< T, U > &L, std::span< const T > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble linear form into a vector.
Definition assembler.h:109
void apply_lifting(std::span< T > b, const std::vector< std::shared_ptr< const Form< T, U > > > &a, const std::vector< std::span< const T > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > > &coeffs, const std::vector< std::vector< std::shared_ptr< const DirichletBC< T, U > > > > &bcs1, const std::vector< std::span< const T > > &x0, T scale)
Modify b such that:
Definition assembler.h:150
std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< U > >, 2 > > > extract_function_spaces(const std::vector< std::vector< const Form< T, U > * > > &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition utils.h:125
void set_bc(std::span< T > b, const std::vector< std::shared_ptr< const DirichletBC< T, U > > > &bcs, std::span< const T > x0, T scale=1)
Set bc values in owned (local) part of the vector, multiplied by 'scale'. The vectors b and x0 must h...
Definition assembler.h:438
la::SparsityPattern create_sparsity_pattern(const Form< T, U > &a)
Create a sparsity pattern for a given form.
Definition utils.h:150
std::array< std::vector< std::shared_ptr< const FunctionSpace< T > > >, 2 > common_function_spaces(const std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< T > >, 2 > > > &V)
Extract FunctionSpaces for (0) rows blocks and (1) columns blocks from a rectangular array of (test,...
Definition FunctionSpace.h:374
Mat create_matrix(MPI_Comm comm, const SparsityPattern &sp, std::string type=std::string())
Create a PETSc Mat. Caller is responsible for destroying the returned object.
Definition petsc.cpp:232