NeoFOAM
WIP Prototype of a modern OpenFOAM core
Loading...
Searching...
No Matches
gaussGreenLaplacian.hpp
Go to the documentation of this file.
1// SPDX-License-Identifier: MIT
2// SPDX-FileCopyrightText: 2024 NeoFOAM authors
3
4#pragma once
5
13
15{
16
17template<typename ValueType>
19 const FaceNormalGradient<ValueType>& faceNormalGradient,
20 const SurfaceField<scalar>& gamma,
22 Field<ValueType>& lapPhi,
23 const dsl::Coeff operatorScaling
24)
25{
26 const UnstructuredMesh& mesh = phi.mesh();
27 const auto exec = phi.exec();
28
29 SurfaceField<ValueType> faceNormalGrad = faceNormalGradient.faceNormalGrad(phi);
30
31 const auto [owner, neighbour, surfFaceCells] =
32 spans(mesh.faceOwner(), mesh.faceNeighbour(), mesh.boundaryMesh().faceCells());
33
34 auto [refGradient, value, valueFraction, refValue] = spans(
35 phi.boundaryField().refGrad(),
36 phi.boundaryField().value(),
39 );
40
41 const auto [result, faceArea, fnGrad, vol] =
42 spans(lapPhi, mesh.magFaceAreas(), faceNormalGrad.internalField(), mesh.cellVolumes());
43
44
45 size_t nInternalFaces = mesh.nInternalFaces();
46
47 // check if the executor is GPU
48 if (std::holds_alternative<SerialExecutor>(exec))
49 {
50 for (size_t i = 0; i < nInternalFaces; i++)
51 {
52 ValueType flux = faceArea[i] * fnGrad[i];
53 result[static_cast<size_t>(owner[i])] += flux;
54 result[static_cast<size_t>(neighbour[i])] -= flux;
55 }
56
57 for (size_t i = nInternalFaces; i < fnGrad.size(); i++)
58 {
59 auto own = static_cast<size_t>(surfFaceCells[i - nInternalFaces]);
60 ValueType valueOwn = faceArea[i] * fnGrad[i];
61 result[own] += valueOwn;
62 }
63
64 for (size_t celli = 0; celli < mesh.nCells(); celli++)
65 {
66 result[celli] *= operatorScaling[celli] / vol[celli];
67 }
68 }
69 else
70 {
72 exec,
73 {0, nInternalFaces},
74 KOKKOS_LAMBDA(const size_t i) {
75 ValueType flux = faceArea[i] * fnGrad[i];
76 Kokkos::atomic_add(&result[static_cast<size_t>(owner[i])], flux);
77 Kokkos::atomic_sub(&result[static_cast<size_t>(neighbour[i])], flux);
78 }
79 );
80
82 exec,
83 {nInternalFaces, fnGrad.size()},
84 KOKKOS_LAMBDA(const size_t i) {
85 auto own = static_cast<size_t>(surfFaceCells[i - nInternalFaces]);
86 ValueType valueOwn = faceArea[i] * fnGrad[i];
87 Kokkos::atomic_add(&result[own], valueOwn);
88 }
89 );
90
92 exec,
93 {0, mesh.nCells()},
94 KOKKOS_LAMBDA(const size_t celli) {
95 result[celli] *= operatorScaling[celli] / vol[celli];
96 }
97 );
98 }
99}
100
101template<typename ValueType>
103 public LaplacianOperatorFactory<ValueType>::template Register<GaussGreenLaplacian<ValueType>>
104{
105 using Base =
107
108public:
109
110 static std::string name() { return "Gauss"; }
111
112 static std::string doc() { return "Gauss-Green Divergence"; }
113
114 static std::string schema() { return "none"; }
115
116 GaussGreenLaplacian(const Executor& exec, const UnstructuredMesh& mesh, const Input& inputs)
117 : Base(exec, mesh), surfaceInterpolation_(exec, mesh, inputs),
118 faceNormalGradient_(exec, mesh, inputs),
119 sparsityPattern_(SparsityPattern::readOrCreate(mesh)) {};
120
122 {
123 la::LinearSystem<scalar, localIdx> ls(sparsityPattern_->linearSystem());
124 auto [A, b, sp] = ls.view();
125 const auto& exec = A.exec();
126
127 Field<ValueType> values(exec, A.nNonZeros(), zero<ValueType>());
128 Field<localIdx> mColIdxs(exec, A.colIdxs().data(), A.nNonZeros());
129 Field<localIdx> mRowPtrs(exec, A.rowPtrs().data(), A.rowPtrs().size());
130
131 la::CSRMatrix<ValueType, localIdx> matrix(values, mColIdxs, mRowPtrs);
132 Field<ValueType> rhs(exec, b.size(), zero<ValueType>());
133
134 return {matrix, rhs, ls.sparsityPattern()};
135 };
136
137 virtual void laplacian(
139 const SurfaceField<scalar>& gamma,
141 const dsl::Coeff operatorScaling
142 ) override
143 {
144 computeLaplacian<ValueType>(
145 faceNormalGradient_, gamma, phi, lapPhi.internalField(), operatorScaling
146 );
147 };
148
149 virtual void laplacian(
150 Field<ValueType>& lapPhi,
151 const SurfaceField<scalar>& gamma,
153 const dsl::Coeff operatorScaling
154 ) override
155 {
156 computeLaplacian<ValueType>(faceNormalGradient_, gamma, phi, lapPhi, operatorScaling);
157 };
158
159 virtual void laplacian(
161 const SurfaceField<scalar>& gamma,
163 const dsl::Coeff operatorScaling
164 ) override
165 {
166 const UnstructuredMesh& mesh = phi.mesh();
167 const std::size_t nInternalFaces = mesh.nInternalFaces();
168 const auto exec = phi.exec();
169 const auto [owner, neighbour, surfFaceCells, diagOffs, ownOffs, neiOffs] = spans(
170 mesh.faceOwner(),
171 mesh.faceNeighbour(),
172 mesh.boundaryMesh().faceCells(),
173 sparsityPattern_->diagOffset(),
174 sparsityPattern_->ownerOffset(),
175 sparsityPattern_->neighbourOffset()
176 );
177
178 const auto [sGamma, deltaCoeffs, magFaceArea] = spans(
179 gamma.internalField(),
180 faceNormalGradient_.deltaCoeffs().internalField(),
181 mesh.magFaceAreas()
182 );
183
184 auto [refGradient, value, valueFraction, refValue] = spans(
185 phi.boundaryField().refGrad(),
186 phi.boundaryField().value(),
188 phi.boundaryField().refValue()
189 );
190
191 const auto rowPtrs = ls.matrix().rowPtrs();
192 const auto colIdxs = ls.matrix().colIdxs();
193 auto values = ls.matrix().values();
194 auto rhs = ls.rhs().span();
195
197 exec,
198 {0, nInternalFaces},
199 KOKKOS_LAMBDA(const size_t facei) {
200 scalar flux = deltaCoeffs[facei] * sGamma[facei] * magFaceArea[facei];
201
202 std::size_t own = static_cast<std::size_t>(owner[facei]);
203 std::size_t nei = static_cast<std::size_t>(neighbour[facei]);
204
205 // add neighbour contribution upper
206 std::size_t rowNeiStart = rowPtrs[nei];
207 std::size_t rowOwnStart = rowPtrs[own];
208
209 // scalar valueNei = (1 - weight) * flux;
210 values[rowNeiStart + neiOffs[facei]] += flux * one<ValueType>();
211 Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], flux * one<ValueType>());
212
213 // upper triangular part
214
215 // add owner contribution lower
216 values[rowOwnStart + ownOffs[facei]] += flux * one<ValueType>();
217 Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], flux * one<ValueType>());
218 }
219 );
220
222 exec,
223 {nInternalFaces, sGamma.size()},
224 KOKKOS_LAMBDA(const size_t facei) {
225 std::size_t bcfacei = facei - nInternalFaces;
226 scalar flux = sGamma[facei] * magFaceArea[facei];
227
228 std::size_t own = static_cast<std::size_t>(surfFaceCells[bcfacei]);
229 std::size_t rowOwnStart = rowPtrs[own];
230
231 values[rowOwnStart + diagOffs[own]] -=
232 flux * valueFraction[bcfacei] * deltaCoeffs[facei] * one<ValueType>();
233 rhs[own] -=
234 (flux
235 * (valueFraction[bcfacei] * deltaCoeffs[facei] * refValue[bcfacei]
236 + (1.0 - valueFraction[bcfacei]) * refGradient[bcfacei]));
237 }
238 );
239
241 exec,
242 {0, rhs.size()},
243 KOKKOS_LAMBDA(const size_t celli) {
244 rhs[celli] *= operatorScaling[celli];
245 for (size_t i = rowPtrs[celli]; i < rowPtrs[celli + 1]; i++)
246 {
247 values[i] *= operatorScaling[celli];
248 }
249 }
250 );
251 };
252
253 std::unique_ptr<LaplacianOperatorFactory<ValueType>> clone() const override
254 {
255 return std::make_unique<GaussGreenLaplacian<ValueType>>(*this);
256 };
257
258private:
259
260 SurfaceInterpolation<scalar> surfaceInterpolation_;
261 FaceNormalGradient<ValueType> faceNormalGradient_;
262 const std::shared_ptr<SparsityPattern> sparsityPattern_;
263};
264
265} // namespace NeoFOAM
const NeoFOAM::Field< T > & refGrad() const
Get the view storing the Neumann boundary values.
const NeoFOAM::Field< T > & value() const
Get the view storing the computed values from the boundary condition.
const NeoFOAM::Field< T > & refValue() const
Get the view storing the Dirichlet boundary values.
const NeoFOAM::Field< scalar > & valueFraction() const
Get the view storing the fraction of the boundary value.
const labelField & faceCells() const
Get the field of face cells.
A class to contain the data and executors for a field and define some basic operations.
Definition field.hpp:49
Represents an unstructured mesh in NeoFOAM.
const BoundaryMesh & boundaryMesh() const
Get the boundary mesh.
const scalarField & cellVolumes() const
Get the field of cell volumes in the mesh.
size_t nInternalFaces() const
Get the number of internal faces in the mesh.
const scalarField & magFaceAreas() const
Get the field of magnitudes of face areas.
size_t nCells() const
Get the number of cells 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.
A class that represents a coefficient for the NeoFOAM dsl.
Definition coeff.hpp:22
void faceNormalGrad(const VolumeField< ValueType > &volField, SurfaceField< ValueType > &surfaceField) const
la::LinearSystem< ValueType, localIdx > createEmptyLinearSystem() const override
virtual void laplacian(VolumeField< ValueType > &lapPhi, const SurfaceField< scalar > &gamma, VolumeField< ValueType > &phi, const dsl::Coeff operatorScaling) override
std::unique_ptr< LaplacianOperatorFactory< ValueType > > clone() const override
virtual void laplacian(la::LinearSystem< ValueType, localIdx > &ls, const SurfaceField< scalar > &gamma, VolumeField< ValueType > &phi, const dsl::Coeff operatorScaling) override
GaussGreenLaplacian(const Executor &exec, const UnstructuredMesh &mesh, const Input &inputs)
virtual void laplacian(Field< ValueType > &lapPhi, const SurfaceField< scalar > &gamma, VolumeField< ValueType > &phi, const dsl::Coeff operatorScaling) override
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 surface field in a finite volume method.
Represents a volume field in a finite volume method.
A class representing a linear system of equations.
LinearSystemView< ValueType, IndexType > view()
Field< ValueType > & rhs()
std::string sparsityPattern() const
CSRMatrix< ValueType, IndexType > & matrix()
A template class for registering derived classes with a base class.
void computeLaplacian(const FaceNormalGradient< ValueType > &faceNormalGradient, const SurfaceField< scalar > &gamma, VolumeField< ValueType > &phi, Field< ValueType > &lapPhi, const dsl::Coeff operatorScaling)
float scalar
Definition scalar.hpp:14
std::variant< Dictionary, TokenList > Input
Definition input.hpp:13
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)