From a7177c8019ccaba45b257427ec4f7f8319293210 Mon Sep 17 00:00:00 2001
From: Henning Scheufler <henning.scheufler@web.de>
Date: Sun, 26 Jan 2025 13:51:26 +0100
Subject: [PATCH] WIP: add implicit div calculation without boundary
 contribution

---
 .../cellCentred/operators/divOperator.hpp     | 11 ++++
 .../cellCentred/operators/gaussGreenDiv.hpp   |  7 +++
 .../cellCentred/operators/gaussGreenDiv.cpp   | 50 ++++++++++++++++++-
 3 files changed, 67 insertions(+), 1 deletion(-)

diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
index 2ed0fb16f..cf5011cfb 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
@@ -4,6 +4,7 @@
 #pragma once
 
 #include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/linearAlgebra.hpp"
 #include "NeoFOAM/core/executor/executor.hpp"
 #include "NeoFOAM/core/input.hpp"
 #include "NeoFOAM/dsl/operator.hpp"
@@ -45,6 +46,11 @@ class DivOperatorFactory :
     div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) = 0;
 
+    virtual void
+    div(la::LinearSystem<scalar, localIdx>& ls,
+        const SurfaceField<scalar>& faceFlux,
+        VolumeField<scalar>& phi) = 0;
+
     virtual void
     div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) = 0;
 
@@ -112,6 +118,11 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<scalar>>
 
     void div(Field<scalar>& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, getField()); }
 
+    void div(la::LinearSystem<scalar, localIdx>& ls)
+    {
+        divOperatorStrategy_->div(ls, faceFlux_, getField());
+    };
+
     void div(VolumeField<scalar>& divPhi)
     {
         divOperatorStrategy_->div(divPhi, faceFlux_, getField());
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
index 6d0fdd9d7..e1a67033f 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
@@ -8,6 +8,7 @@
 #include "NeoFOAM/mesh/unstructured.hpp"
 #include "NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp"
 #include "NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp"
 
 namespace NeoFOAM::finiteVolume::cellCentred
 {
@@ -28,6 +29,11 @@ class GaussGreenDiv : public DivOperatorFactory::Register<GaussGreenDiv>
     div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) override;
 
+    void
+    div(la::LinearSystem<scalar, localIdx>& ls,
+        const SurfaceField<scalar>& faceFlux,
+        VolumeField<scalar>& phi) override;
+
     void
     div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) override;
@@ -40,6 +46,7 @@ class GaussGreenDiv : public DivOperatorFactory::Register<GaussGreenDiv>
 private:
 
     SurfaceInterpolation surfaceInterpolation_;
+    const std::shared_ptr<SparsityPattern> sparsityPattern_;
 };
 
 } // namespace NeoFOAM
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
index 890bb3779..0bf38e40c 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
@@ -100,7 +100,8 @@ GaussGreenDiv::GaussGreenDiv(
     const Executor& exec, const UnstructuredMesh& mesh, const Input& inputs
 )
     : DivOperatorFactory::Register<GaussGreenDiv>(exec, mesh),
-      surfaceInterpolation_(exec, mesh, inputs) {};
+      surfaceInterpolation_(exec, mesh, inputs),
+      sparsityPattern_(SparsityPattern::readOrCreate(mesh)) {};
 
 void GaussGreenDiv::div(
     VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
@@ -109,6 +110,53 @@ void GaussGreenDiv::div(
     computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi);
 };
 
+void GaussGreenDiv::div(
+    la::LinearSystem<scalar, localIdx>& ls,
+    const SurfaceField<scalar>& faceFlux,
+    VolumeField<scalar>& phi
+)
+{
+    const UnstructuredMesh& mesh = phi.mesh();
+    const std::size_t nInternalFaces = mesh.nInternalFaces();
+    const auto exec = phi.exec();
+    const auto [sFaceFlux, owner, neighbour, diagOffs, ownOffs, neiOffs] = spans(
+        faceFlux.internalField(),
+        mesh.faceOwner(),
+        mesh.faceNeighbour(),
+        sparsityPattern_->diagOffset(),
+        sparsityPattern_->ownerOffset(),
+        sparsityPattern_->neighbourOffset()
+    );
+    const auto rowPtrs = ls.matrix().rowPtrs();
+    const auto colIdxs = ls.matrix().colIdxs();
+    const auto values = ls.matrix().values();
+
+    parallelFor(
+        exec,
+        {0, nInternalFaces},
+        KOKKOS_LAMBDA(const size_t facei) {
+            scalar flux = sFaceFlux[facei];
+            scalar weight = 0.5;
+            std::size_t own = static_cast<std::size_t>(owner[facei]);
+            std::size_t nei = static_cast<std::size_t>(neighbour[facei]);
+
+            // lower triangular part
+
+            // add neighbour contribution
+            std::size_t rowNeiStart = rowPtrs[nei];
+            values[rowNeiStart + neiOffs[facei]] += -flux * weight;
+            Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], -flux * weight);
+
+            // upper triangular part
+
+            // add owner contribution
+            std::size_t rowOwnStart = rowPtrs[own];
+            values[rowOwnStart + ownOffs[facei]] += flux * (1 - weight);
+            Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], flux * (1 - weight));
+        }
+    );
+};
+
 void GaussGreenDiv::div(
     Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
 )