From cb9d6566d512e3315620914bdc1cb23a57f252a1 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 6 Oct 2022 15:50:45 +0200
Subject: [PATCH] changed: put flux assembly for tracer in separate function

---
 ebos/ecltracermodel.hh | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh
index 09dd3fbe6..3c147214e 100644
--- a/ebos/ecltracermodel.hh
+++ b/ebos/ecltracermodel.hh
@@ -290,6 +290,26 @@ protected:
         (*this->tracerMatrix_)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt;
     }
 
+    template<class TrRe>
+    void assembleTracerEquationFlux(TrRe& tr,
+                                    const ElementContext& elemCtx,
+                                    unsigned scvfIdx,
+                                    unsigned I,
+                                    unsigned J)
+    {
+        TracerEvaluation flux;
+        bool isUpF;
+        computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
+        int globalUpIdx = isUpF ? I : J;
+        for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+            tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux
+        }
+        if (isUpF) {
+            (*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
+            (*this->tracerMatrix_)[I][I][0][0] += flux.derivative(0);
+        }
+    }
+
     template <class TrRe>
     void assembleTracerEquations_(TrRe & tr)
     {
@@ -335,22 +355,11 @@ protected:
 
             size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
             for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
-                TracerEvaluation flux;
                 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
                 unsigned j = face.exteriorIndex();
                 unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
-                bool isUpF;
-                computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
-                int globalUpIdx = isUpF ? I : J;
-                for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
-                    tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux
-                }
-                if (isUpF) {
-                    (*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
-                    (*this->tracerMatrix_)[I][I][0][0] += flux.derivative(0);
-                }
+                this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J);
             }
-
         }
 
         // Wells  terms