20template<
typename ValueType,
typename IndexType = localIdx>
23 auto [A, b, sp] = ls.
view();
24 const auto& exec = A.exec();
35template<
typename ValueType,
typename IndexType = localIdx>
46 : psi_(psi), expr_(expr), fvSchemes_(fvSchemes), fvSolution_(fvSolution),
50 expr_.build(fvSchemes_);
55 : psi_(ls.psi_), expr_(ls.expr_), fvSchemes_(ls.fvSchemes_), fvSolution_(ls.fvSolution_),
56 ls_(ls.ls_), sparsityPattern_(ls.sparsityPattern_) {};
63 if (!sparsityPattern_)
65 NF_THROW(
"fvcc:LinearSystem:sparsityPattern: sparsityPattern is null");
67 return *sparsityPattern_;
77 if (!sparsityPattern_)
79 NF_THROW(
"fvcc:LinearSystem:sparsityPattern: sparsityPattern is null");
81 return *sparsityPattern_;
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();
93 ls_ = expr_.implicitOperation();
94 expr_.implicitOperation(ls_, t, dt);
95 auto rhs = ls_.rhs().span();
100 KOKKOS_LAMBDA(
const size_t i) { rhs[i] -= expSourceSpan[i] * vol[i]; }
106 if (expr_.temporalOperators().size() == 0 && expr_.spatialOperators().size() == 0)
111 if (expr_.temporalOperators().size() > 0)
122 auto vol = psi_.mesh().cellVolumes().span();
123 auto expSource = expr_.explicitOperation(psi_.mesh().nCells());
124 auto expSourceSpan = expSource.span();
126 ls_ = expr_.implicitOperation();
127 auto rhs = ls_.rhs().span();
132 KOKKOS_LAMBDA(
const size_t i) { rhs[i] -= expSourceSpan[i] * vol[i]; }
141 if (expr_.temporalOperators().size() == 0 && expr_.spatialOperators().size() == 0)
145 if (expr_.temporalOperators().size() > 0)
155 NeoFOAM::la::ginkgo::BiCGStab<ValueType> solver(psi_.exec(), fvSolution_);
156 auto convertedLS = convertLinearSystem(ls_);
157 solver.solve(convertedLS, psi_.internalField());
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();
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;
184 const auto exec = psi_.exec();
185 const auto [owner, neighbour, surfFaceCells, ownOffs, neiOffs, internalPsi] =
spans(
189 sparsityPattern_->ownerOffset(),
190 sparsityPattern_->neighbourOffset(),
194 auto rhs = ls_.rhs().span();
195 auto [values, colIdxs, rowPtrs] = ls_.view();
200 auto resultSpan = result.
span();
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]);
209 std::size_t rowNeiStart = rowPtrs[nei];
210 std::size_t rowOwnStart = rowPtrs[own];
212 auto Upper = values[rowNeiStart + neiOffs[facei]];
213 auto Lower = values[rowOwnStart + ownOffs[facei]];
215 &resultSpan[facei], Upper * internalPsi[nei] - Lower * internalPsi[own]
229 std::shared_ptr<SparsityPattern> sparsityPattern_;
232template<
typename ValueType,
typename IndexType = localIdx>
233VolumeField<ValueType>
245 auto [result, b, x] =
247 const auto [values, colIdxs, rowPtrs] = ls.
linearSystem().view();
252 KOKKOS_LAMBDA(
const std::size_t rowi) {
253 IndexType rowStart = rowPtrs[rowi];
254 IndexType rowEnd = rowPtrs[rowi + 1];
256 for (IndexType coli = rowStart; coli < rowEnd; coli++)
258 sum += values[coli] * x[colIdxs[coli]];
260 result[rowi] =
sum - b[rowi];
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.
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
void assemble(scalar t, scalar dt)
const SparsityPattern & sparsityPattern() const
void solve(scalar t, scalar dt)
void setReference(const IndexType refCell, ValueType refValue)
const la::LinearSystem< ValueType, IndexType > & linearSystem() const
Expression(const Expression &ls)
la::LinearSystem< ValueType, IndexType > & linearSystem()
SparsityPattern & sparsityPattern()
VolumeField< ValueType > & getField()
const Executor & exec() const
Expression(dsl::Expression< ValueType > expr, VolumeField< ValueType > &psi, const Dictionary &fvSchemes, const Dictionary &fvSolution)
Field< ValueType > flux() const
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.
#define NF_ERROR_EXIT(message)
Macro for printing an error message and aborting the program.
VolumeField< ValueType > operator&(const Expression< ValueType, IndexType > ls, const VolumeField< ValueType > &psi)
la::LinearSystem< ValueType, IndexType > convert(const la::LinearSystem< scalar, IndexType > &ls)
T sum(const Field< T > &field)
void parallelFor(const Executor &exec, std::pair< size_t, size_t > range, Kernel kernel)
std::variant< SerialExecutor, CPUExecutor, GPUExecutor > Executor
auto spans(Args &... fields)