DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
petsc.h
1// Copyright (C) 2004-2018 Johan Hoffman, Johan Jansson, Anders Logg and
2// Garth N. Wells
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#ifdef HAS_PETSC
11
12#include "Vector.h"
13#include "utils.h"
14#include <boost/lexical_cast.hpp>
15#include <functional>
16#include <optional>
17#include <petscksp.h>
18#include <petscmat.h>
19#include <petscoptions.h>
20#include <petscvec.h>
21#include <span>
22#include <string>
23#include <vector>
24
25namespace dolfinx::common
26{
27class IndexMap;
28} // namespace dolfinx::common
29
30namespace dolfinx::la
31{
32class SparsityPattern;
33
35namespace petsc
36{
38void error(int error_code, std::string filename, std::string petsc_function);
39
47std::vector<Vec>
48create_vectors(MPI_Comm comm,
49 const std::vector<std::span<const PetscScalar>>& x);
50
56Vec create_vector(const common::IndexMap& map, int bs);
57
66Vec create_vector(MPI_Comm comm, std::array<std::int64_t, 2> range,
67 std::span<const std::int64_t> ghosts, int bs);
68
78Vec create_vector_wrap(const common::IndexMap& map, int bs,
79 std::span<const PetscScalar> x);
80
84template <class V>
85Vec create_vector_wrap(const la::Vector<V>& x)
86{
87 assert(x.index_map());
88 return create_vector_wrap(*x.index_map(), x.bs(), x.array());
89}
90
103std::vector<IS> create_index_sets(
104 const std::vector<
105 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
106
108std::vector<std::vector<PetscScalar>> get_local_vectors(
109 const Vec x,
110 const std::vector<
111 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
112
114void scatter_local_vectors(
115 Vec x, const std::vector<std::span<const PetscScalar>>& x_b,
116 const std::vector<
117 std::pair<std::reference_wrapper<const common::IndexMap>, int>>& maps);
118
121Mat create_matrix(MPI_Comm comm, const SparsityPattern& sp,
122 std::optional<std::string> type = std::nullopt);
123
129MatNullSpace create_nullspace(MPI_Comm comm, std::span<const Vec> basis);
130
136namespace options
137{
139void set(std::string option);
140
142template <typename T>
143void set(std::string option, const T value)
144{
145 if (option[0] != '-')
146 option = '-' + option;
147
148 PetscErrorCode ierr;
149 ierr = PetscOptionsSetValue(nullptr, option.c_str(),
150 boost::lexical_cast<std::string>(value).c_str());
151 if (ierr != 0)
152 petsc::error(ierr, __FILE__, "PetscOptionsSetValue");
153}
154
156void clear(std::string option);
157
159void clear();
160} // namespace options
161
167class Vector
168{
169public:
174 Vector(const common::IndexMap& map, int bs);
175
176 // Delete copy constructor to avoid accidental copying of 'heavy' data
177 Vector(const Vector& x) = delete;
178
180 Vector(Vector&& x);
181
193 Vector(Vec x, bool inc_ref_count);
194
196 virtual ~Vector();
197
198 // Assignment operator (disabled)
199 Vector& operator=(const Vector& x) = delete;
200
202 Vector& operator=(Vector&& x);
203
206 Vector copy() const;
207
209 std::int64_t size() const;
210
212 std::int32_t local_size() const;
213
215 std::array<std::int64_t, 2> local_range() const;
216
218 MPI_Comm comm() const;
219
221 void set_options_prefix(std::string options_prefix);
222
225 std::string get_options_prefix() const;
226
228 void set_from_options();
229
231 Vec vec() const;
232
233private:
234 // PETSc Vec pointer
235 Vec _x;
236};
237
240class Operator
241{
242public:
244 Operator(Mat A, bool inc_ref_count);
245
246 // Copy constructor (deleted)
247 Operator(const Operator& A) = delete;
248
250 Operator(Operator&& A);
251
253 virtual ~Operator();
254
256 Operator& operator=(const Operator& A) = delete;
257
259 Operator& operator=(Operator&& A);
260
263 std::array<std::int64_t, 2> size() const;
264
270 Vec create_vector(std::size_t dim) const;
271
273 Mat mat() const;
274
275protected:
276 // PETSc Mat pointer
277 Mat _matA;
278};
279
285class Matrix : public Operator
286{
287public:
292 static auto set_fn(Mat A, InsertMode mode)
293 {
294 return [A, mode, cache = std::vector<PetscInt>()](
295 std::span<const std::int32_t> rows,
296 std::span<const std::int32_t> cols,
297 std::span<const PetscScalar> vals) mutable -> int
298 {
299 PetscErrorCode ierr;
300#ifdef PETSC_USE_64BIT_INDICES
301 cache.resize(rows.size() + cols.size());
302 std::ranges::copy(rows, cache.begin());
303 std::ranges::copy(cols, std::next(cache.begin(), rows.size()));
304 const PetscInt* _rows = cache.data();
305 const PetscInt* _cols = cache.data() + rows.size();
306 ierr = MatSetValuesLocal(A, rows.size(), _rows, cols.size(), _cols,
307 vals.data(), mode);
308#else
309 ierr = MatSetValuesLocal(A, rows.size(), rows.data(), cols.size(),
310 cols.data(), vals.data(), mode);
311#endif
312
313#ifndef NDEBUG
314 if (ierr != 0)
315 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
316#endif
317 return ierr;
318 };
319 }
320
326 static auto set_block_fn(Mat A, InsertMode mode)
327 {
328 return [A, mode, cache = std::vector<PetscInt>()](
329 std::span<const std::int32_t> rows,
330 std::span<const std::int32_t> cols,
331 std::span<const PetscScalar> vals) mutable -> int
332 {
333 PetscErrorCode ierr;
334#ifdef PETSC_USE_64BIT_INDICES
335 cache.resize(rows.size() + cols.size());
336 std::ranges::copy(rows, cache.begin());
337 std::ranges::copy(cols, std::next(cache.begin(), rows.size()));
338 const PetscInt* _rows = cache.data();
339 const PetscInt* _cols = cache.data() + rows.size();
340 ierr = MatSetValuesBlockedLocal(A, rows.size(), _rows, cols.size(), _cols,
341 vals.data(), mode);
342#else
343 ierr = MatSetValuesBlockedLocal(A, rows.size(), rows.data(), cols.size(),
344 cols.data(), vals.data(), mode);
345#endif
346
347#ifndef NDEBUG
348 if (ierr != 0)
349 petsc::error(ierr, __FILE__, "MatSetValuesBlockedLocal");
350#endif
351 return ierr;
352 };
353 }
354
363 static auto set_block_expand_fn(Mat A, int bs0, int bs1, InsertMode mode)
364 {
365 return [A, bs0, bs1, mode, cache0 = std::vector<PetscInt>(),
366 cache1 = std::vector<PetscInt>()](
367 std::span<const std::int32_t> rows,
368 std::span<const std::int32_t> cols,
369 std::span<const PetscScalar> vals) mutable -> int
370 {
371 PetscErrorCode ierr;
372 cache0.resize(bs0 * rows.size());
373 cache1.resize(bs1 * cols.size());
374 for (std::size_t i = 0; i < rows.size(); ++i)
375 for (int k = 0; k < bs0; ++k)
376 cache0[bs0 * i + k] = bs0 * rows[i] + k;
377
378 for (std::size_t i = 0; i < cols.size(); ++i)
379 for (int k = 0; k < bs1; ++k)
380 cache1[bs1 * i + k] = bs1 * cols[i] + k;
381
382 ierr = MatSetValuesLocal(A, cache0.size(), cache0.data(), cache1.size(),
383 cache1.data(), vals.data(), mode);
384#ifndef NDEBUG
385 if (ierr != 0)
386 petsc::error(ierr, __FILE__, "MatSetValuesLocal");
387#endif
388 return ierr;
389 };
390 }
391
393 Matrix(MPI_Comm comm, const SparsityPattern& sp,
394 std::optional<std::string> type = std::nullopt);
395
400 Matrix(Mat A, bool inc_ref_count);
401
402 // Copy constructor (deleted)
403 Matrix(const Matrix& A) = delete;
404
406 Matrix(Matrix&& A) = default;
407
409 ~Matrix() = default;
410
412 Matrix& operator=(const Matrix& A) = delete;
413
415 Matrix& operator=(Matrix&& A) = default;
416
420 enum class AssemblyType : std::int32_t
421 {
422 FINAL,
423 FLUSH
424 };
425
431 void apply(AssemblyType type);
432
434 double norm(Norm norm_type) const;
435
436 //--- Special PETSc Functions ---
437
440 void set_options_prefix(std::string options_prefix);
441
444 std::string get_options_prefix() const;
445
447 void set_from_options();
448};
449
452class KrylovSolver
453{
454public:
457 explicit KrylovSolver(MPI_Comm comm);
458
462 KrylovSolver(KSP ksp, bool inc_ref_count);
463
464 // Copy constructor (deleted)
465 KrylovSolver(const KrylovSolver& solver) = delete;
466
468 KrylovSolver(KrylovSolver&& solver);
469
471 ~KrylovSolver();
472
473 // Assignment operator (deleted)
474 KrylovSolver& operator=(const KrylovSolver&) = delete;
475
477 KrylovSolver& operator=(KrylovSolver&& solver);
478
480 void set_operator(const Mat A);
481
483 void set_operators(const Mat A, const Mat P);
484
487 int solve(Vec x, const Vec b, bool transpose = false) const;
488
491 void set_options_prefix(std::string options_prefix);
492
495 std::string get_options_prefix() const;
496
498 void set_from_options() const;
499
501 KSP ksp() const;
502
503private:
504 // PETSc solver pointer
505 KSP _ksp;
506};
507} // namespace petsc
508} // namespace dolfinx::la
509
510#endif
Definition IndexMap.h:94
Definition SparsityPattern.h:26
int size(MPI_Comm comm)
Definition MPI.cpp:72
constexpr std::array< std::int64_t, 2 > local_range(int rank, std::int64_t N, int size)
Return local range for the calling process, partitioning the global [0, N - 1] range across all ranks...
Definition MPI.h:90
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
Linear algebra interface.
Definition sparsitybuild.h:15
auto norm(const V &x, Norm type=Norm::l2)
Definition Vector.h:268