NeoFOAM
WIP Prototype of a modern OpenFOAM core
Loading...
Searching...
No Matches
expression.hpp
Go to the documentation of this file.
1// SPDX-License-Identifier: MIT
2// SPDX-FileCopyrightText: 2025 NeoFOAM authors
3
4#pragma once
5
11
12
14
15namespace dsl = NeoFOAM::dsl;
16
18{
19// TODO extend sparsity pattern to return the correct type
20template<typename ValueType, typename IndexType = localIdx>
22{
23 auto [A, b, sp] = ls.view();
24 const auto& exec = A.exec();
25
26 Field<ValueType> values(exec, A.nNonZeros(), zero<ValueType>());
27 Field<localIdx> mColIdxs(exec, A.colIdxs().data(), A.nNonZeros());
28 Field<localIdx> mRowPtrs(exec, A.rowPtrs().data(), A.rowPtrs().size());
29
30 la::CSRMatrix<ValueType, localIdx> matrix(values, mColIdxs, mRowPtrs);
31 Field<ValueType> rhs(exec, b.size(), zero<ValueType>());
32 return {matrix, rhs, ls.sparsityPattern()};
33}
34
35template<typename ValueType, typename IndexType = localIdx>
37{
38public:
39
43 [[maybe_unused]] const Dictionary& fvSchemes,
44 [[maybe_unused]] const Dictionary& fvSolution
45 )
46 : psi_(psi), expr_(expr), fvSchemes_(fvSchemes), fvSolution_(fvSolution),
47 ls_(convert<ValueType>(SparsityPattern::readOrCreate(psi.mesh())->linearSystem())),
48 sparsityPattern_(SparsityPattern::readOrCreate(psi.mesh()))
49 {
50 expr_.build(fvSchemes_);
51 assemble();
52 };
53
55 : psi_(ls.psi_), expr_(ls.expr_), fvSchemes_(ls.fvSchemes_), fvSolution_(ls.fvSolution_),
56 ls_(ls.ls_), sparsityPattern_(ls.sparsityPattern_) {};
57
58 ~Expression() = default;
59
62 {
63 if (!sparsityPattern_)
64 {
65 NF_THROW("fvcc:LinearSystem:sparsityPattern: sparsityPattern is null");
66 }
67 return *sparsityPattern_;
68 }
69
70 VolumeField<ValueType>& getField() { return this->psi_; }
71
72 const VolumeField<ValueType>& getField() const { return this->psi_; }
73
74 [[nodiscard]] const la::LinearSystem<ValueType, IndexType>& linearSystem() const { return ls_; }
75 [[nodiscard]] const SparsityPattern& sparsityPattern() const
76 {
77 if (!sparsityPattern_)
78 {
79 NF_THROW("fvcc:LinearSystem:sparsityPattern: sparsityPattern is null");
80 }
81 return *sparsityPattern_;
82 }
83
84 const Executor& exec() const { return ls_.exec(); }
85
87 {
88 auto vol = psi_.mesh().cellVolumes().span();
89 auto expSource = expr_.explicitOperation(psi_.mesh().nCells());
90 expr_.explicitOperation(expSource, t, dt);
91 auto expSourceSpan = expSource.span();
92
93 ls_ = expr_.implicitOperation();
94 expr_.implicitOperation(ls_, t, dt);
95 auto rhs = ls_.rhs().span();
96 // we subtract the explicit source term from the rhs
98 exec(),
99 {0, rhs.size()},
100 KOKKOS_LAMBDA(const size_t i) { rhs[i] -= expSourceSpan[i] * vol[i]; }
101 );
102 }
103
104 void assemble()
105 {
106 if (expr_.temporalOperators().size() == 0 && expr_.spatialOperators().size() == 0)
107 {
108 NF_ERROR_EXIT("No temporal or implicit terms to solve.");
109 }
110
111 if (expr_.temporalOperators().size() > 0)
112 {
113 // integrate equations in time
114 // NeoFOAM::timeIntegration::TimeIntegration<VolumeField<ValueType>> timeIntegrator(
115 // fvSchemes_.subDict("ddtSchemes"), fvSolution_
116 // );
117 // timeIntegrator.solve(expr_, psi_, t, dt);
118 }
119 else
120 {
121 // solve sparse matrix system
122 auto vol = psi_.mesh().cellVolumes().span();
123 auto expSource = expr_.explicitOperation(psi_.mesh().nCells());
124 auto expSourceSpan = expSource.span();
125
126 ls_ = expr_.implicitOperation();
127 auto rhs = ls_.rhs().span();
128 // we subtract the explicit source term from the rhs
130 exec(),
131 {0, rhs.size()},
132 KOKKOS_LAMBDA(const size_t i) { rhs[i] -= expSourceSpan[i] * vol[i]; }
133 );
134 }
135 }
136
137 void solve(scalar t, scalar dt)
138 {
139 // dsl::solve(expr_, psi_, t, dt, fvSchemes_, fvSolution_);
140 // FIXME:
141 if (expr_.temporalOperators().size() == 0 && expr_.spatialOperators().size() == 0)
142 {
143 NF_ERROR_EXIT("No temporal or implicit terms to solve.");
144 }
145 if (expr_.temporalOperators().size() > 0)
146 {
147 // // integrate equations in time
148 // NeoFOAM::timeIntegration::TimeIntegration<VolumeField<ValueType>> timeIntegrator(
149 // fvSchemes_.subDict("ddtSchemes"), fvSolution_
150 // );
151 // timeIntegrator.solve(expr_, psi_, t, dt);
152 }
153 else
154 {
155 NeoFOAM::la::ginkgo::BiCGStab<ValueType> solver(psi_.exec(), fvSolution_);
156 auto convertedLS = convertLinearSystem(ls_);
157 solver.solve(convertedLS, psi_.internalField());
158 }
159 }
160
161 void setReference(const IndexType refCell, ValueType refValue)
162 {
163 // TODO currently assumes that matrix is already assembled
164 const auto diagOffset = sparsityPattern_->diagOffset().span();
165 const auto rowPtrs = ls_.matrix().rowPtrs();
166 auto rhs = ls_.rhs().span();
167 auto values = ls_.matrix().values();
169 ls_.exec(),
170 {refCell, refCell + 1},
171 KOKKOS_LAMBDA(const std::size_t refCelli) {
172 auto diagIdx = rowPtrs[refCelli] + diagOffset[refCelli];
173 auto diagValue = values[diagIdx];
174 rhs[refCelli] += diagValue * refValue;
175 values[diagIdx] += diagValue;
176 }
177 );
178 }
179
181 {
182 const UnstructuredMesh& mesh = psi_.mesh();
183 const std::size_t nInternalFaces = mesh.nInternalFaces();
184 const auto exec = psi_.exec();
185 const auto [owner, neighbour, surfFaceCells, ownOffs, neiOffs, internalPsi] = spans(
186 mesh.faceOwner(),
187 mesh.faceNeighbour(),
188 mesh.boundaryMesh().faceCells(),
189 sparsityPattern_->ownerOffset(),
190 sparsityPattern_->neighbourOffset(),
191 psi_.internalField()
192 );
193
194 auto rhs = ls_.rhs().span();
195 auto [values, colIdxs, rowPtrs] = ls_.view();
196
197 Field<ValueType> result(exec, neighbour.size(), 0.0);
198
199
200 auto resultSpan = result.span();
201
203 exec,
204 {0, nInternalFaces},
205 KOKKOS_LAMBDA(const size_t facei) {
206 std::size_t own = static_cast<std::size_t>(owner[facei]);
207 std::size_t nei = static_cast<std::size_t>(neighbour[facei]);
208
209 std::size_t rowNeiStart = rowPtrs[nei];
210 std::size_t rowOwnStart = rowPtrs[own];
211
212 auto Upper = values[rowNeiStart + neiOffs[facei]];
213 auto Lower = values[rowOwnStart + ownOffs[facei]];
214 Kokkos::atomic_add(
215 &resultSpan[facei], Upper * internalPsi[nei] - Lower * internalPsi[own]
216 );
217 }
218 );
219 return result;
220 }
221
222private:
223
226 const Dictionary& fvSchemes_;
227 const Dictionary& fvSolution_;
229 std::shared_ptr<SparsityPattern> sparsityPattern_;
230};
231
232template<typename ValueType, typename IndexType = localIdx>
233VolumeField<ValueType>
235{
236 VolumeField<ValueType> resultField(
237 psi.exec(),
238 "ls_" + psi.name,
239 psi.mesh(),
240 psi.internalField(),
241 psi.boundaryField(),
243 );
244
245 auto [result, b, x] =
246 spans(resultField.internalField(), ls.linearSystem().rhs(), psi.internalField());
247 const auto [values, colIdxs, rowPtrs] = ls.linearSystem().view();
248
250 resultField.exec(),
251 {0, result.size()},
252 KOKKOS_LAMBDA(const std::size_t rowi) {
253 IndexType rowStart = rowPtrs[rowi];
254 IndexType rowEnd = rowPtrs[rowi + 1];
255 ValueType sum = 0.0;
256 for (IndexType coli = rowStart; coli < rowEnd; coli++)
257 {
258 sum += values[coli] * x[colIdxs[coli]];
259 }
260 result[rowi] = sum - b[rowi];
261 }
262 );
263
264 return resultField;
265}
266
267}
const labelField & faceCells() const
Get the field of face cells.
A class representing a dictionary that stores key-value pairs.
A class to contain the data and executors for a field and define some basic operations.
Definition field.hpp:49
std::span< ValueType > span() &&=delete
Represents an unstructured mesh in NeoFOAM.
const BoundaryMesh & boundaryMesh() const
Get the boundary mesh.
size_t nInternalFaces() const
Get the number of internal faces in the mesh.
const labelField & faceOwner() const
Get the field of face owner cells.
const labelField & faceNeighbour() const
Get the field of face neighbour cells.
const VolumeField< ValueType > & getField() const
const SparsityPattern & sparsityPattern() const
void setReference(const IndexType refCell, ValueType refValue)
const la::LinearSystem< ValueType, IndexType > & linearSystem() const
la::LinearSystem< ValueType, IndexType > & linearSystem()
Expression(dsl::Expression< ValueType > expr, VolumeField< ValueType > &psi, const Dictionary &fvSchemes, const Dictionary &fvSolution)
const UnstructuredMesh & mesh() const
Returns a const reference to the unstructured mesh object.
const Executor & exec() const
Returns a const reference to the executor object.
const Field< ValueType > & internalField() const
Returns a const reference to the internal field.
const BoundaryFields< ValueType > & boundaryField() const
Returns a const reference to the boundary field.
Represents a volume field in a finite volume method.
std::vector< VolumeBoundary< ValueType > > boundaryConditions() const
A class representing a linear system of equations.
LinearSystemView< ValueType, IndexType > view()
std::string sparsityPattern() const
#define NF_THROW(message)
Macro for throwing a NeoFOAMException with the specified error message.
Definition error.hpp:127
#define NF_ERROR_EXIT(message)
Macro for printing an error message and aborting the program.
Definition error.hpp:108
VolumeField< ValueType > operator&(const Expression< ValueType, IndexType > ls, const VolumeField< ValueType > &psi)
la::LinearSystem< ValueType, IndexType > convert(const la::LinearSystem< scalar, IndexType > &ls)
float scalar
Definition scalar.hpp:14
T sum(const Field< T > &field)
Definition sum.hpp:58
void parallelFor(const Executor &exec, std::pair< size_t, size_t > range, Kernel kernel)
std::variant< SerialExecutor, CPUExecutor, GPUExecutor > Executor
Definition executor.hpp:16
auto spans(Args &... fields)