NeoN
A framework for CFD software
Loading...
Searching...
No Matches
petsc.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2025 NeoN authors
2//
3// SPDX-License-Identifier: MIT
4
5#pragma once
6
7#if NF_WITH_PETSC
8
9#include <Kokkos_Core.hpp>
10#include <petscvec_kokkos.hpp>
11#include <petscmat.h>
12#include <petscksp.h>
13
14#include "NeoN/fields/field.hpp"
20// #include "NeoN/core/database/petscSolverContextCollection.hpp"
21
23
24namespace NeoN::la::petsc
25{
26
27// template<typename ValueType>
28
29class petscSolver : public SolverFactory::template Register<petscSolver>
30{
31
32 using Base = SolverFactory::template Register<petscSolver>;
33
34private:
35
36 Dictionary solverDict_;
37 // Mat Amat_;
38 // KSP ksp_;
39
40 // Vec sol_, rhs_;
41
42 // NeoN::la::petscSolverContext::petscSolverContext<scalar> petsctx_;
43
44 /*
45 NeoN::Database& db_;
46 std::string eqnName_;
47 */
48
49public:
50
51 petscSolver(Executor exec, Dictionary solverDict) : Base(exec), solverDict_(solverDict)
52 //, Amat_(nullptr), ksp_(nullptr), sol_(nullptr),
53 // rhs_(nullptr),petsctx_(exec_)
54 //, db_(db), eqnName_(eqnName)
55 {}
56
57 //- Destructor
58 virtual ~petscSolver()
59 {
60 // MatDestroy(&Amat_);
61 // KSPDestroy(&ksp_);
62 // VecDestroy(&sol_);
63 // VecDestroy(&rhs_);
64 }
65
66 static std::string name() { return "Petsc"; }
67
68 static std::string doc() { return "TBD"; }
69
70 static std::string schema() { return "none"; }
71
72 // TODO why use a smart pointer here?
73 virtual std::unique_ptr<SolverFactory> clone() const final
74 {
75 NF_ERROR_EXIT("Not implemented");
76 return {};
77 }
78
79
80 virtual SolverStats
81 solve(const LinearSystem<scalar, localIdx>& sys, Vector<scalar>& x) const final
82 {
83
84 Mat Amat;
85 KSP ksp;
86
87 Vec sol, rhs;
88
89 NeoN::la::petscSolverContext::petscSolverContext<scalar> petsctx(exec_);
90
91 std::size_t nrows = sys.rhs().size();
92
93 petsctx.initialize(sys);
94
95 Amat = petsctx.AMat();
96 rhs = petsctx.rhs();
97 sol = petsctx.sol();
98 ksp = petsctx.ksp();
99
100 VecSetValuesCOO(rhs, sys.rhs().data(), ADD_VALUES);
101 MatSetValuesCOO(Amat, sys.matrix().values().data(), ADD_VALUES);
102
103 // MatView(Amat_, PETSC_VIEWER_STDOUT_WORLD);
104 // VecView(rhs_, PETSC_VIEWER_STDOUT_WORLD);
105
106
107 // KSPCreate(PETSC_COMM_WORLD, &ksp_);
108 KSPSetOperators(ksp, Amat, Amat);
109 // KSPSetTolerances(ksp, 1.e-9, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
110 // KSPSetFromOptions(ksp_);
111 // KSPSetUp(ksp);
112
113
114 PetscCallVoid(KSPSolve(ksp, rhs, sol));
115
116 auto numIter = 0;
117 KSPGetIterationNumber(ksp, &numIter);
118
119
120 PetscScalar* x_help = static_cast<PetscScalar*>(x.data());
121 VecGetArray(sol, &x_help);
122
123 Vector<scalar> x2(x.exec(), static_cast<scalar*>(x_help), nrows, x.exec());
124 x = x2;
125
126 // TODO residual norms are missing
127 return {numIter, 0.0, 0.0, 0.0};
128 }
129};
130
131}
132
133#endif
A template class for registering derived classes with a base class.
#define NF_ERROR_EXIT(message)
Macro for printing an error message and aborting the program.
Definition error.hpp:110
void solve(Expression< typename VectorType::ElementType > &exp, VectorType &solution, scalar t, scalar dt, const Dictionary &fvSchemes, const Dictionary &fvSolution)
Definition solver.hpp:38
float scalar
Definition scalar.hpp:16
const std::string & name(const NeoN::Document &doc)
Retrieves the name of a Document.