From 3b2c4c2123e14fd319555bc388d06182a89a46a6 Mon Sep 17 00:00:00 2001
From: ftomei <ftomei@arpae.it>
Date: Mon, 13 Jan 2025 17:47:39 +0100
Subject: [PATCH] update 3D

---
 soilFluxes3D/header/solver.h |  2 ++
 soilFluxes3D/solver.cpp      | 36 ++++++++++++++++++++++++++++++++++++
 soilFluxes3D/water.cpp       |  1 -
 3 files changed, 38 insertions(+), 1 deletion(-)

diff --git a/soilFluxes3D/header/solver.h b/soilFluxes3D/header/solver.h
index 52195724a..9c0129ea3 100644
--- a/soilFluxes3D/header/solver.h
+++ b/soilFluxes3D/header/solver.h
@@ -13,5 +13,7 @@
 
     bool GaussSeidelRelaxation (int myApproximation, double myResidualTolerance, int myProcess);
 
+    bool JacobiRelaxation (int approximation, double toleranceThreshold, int process);
+
 #endif  // SOLVER_H
 
diff --git a/soilFluxes3D/solver.cpp b/soilFluxes3D/solver.cpp
index dbcb97363..cc17334f9 100644
--- a/soilFluxes3D/solver.cpp
+++ b/soilFluxes3D/solver.cpp
@@ -233,3 +233,39 @@ bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int pr
 
     return true;
 }
+
+
+bool JacobiRelaxation (int approximation, double toleranceThreshold, int process)
+{
+    double currentNorm = 1.0;
+    double bestNorm = currentNorm;
+
+    int maxIterationsNr = getMaxIterationsNr(approximation);
+    int iteration = 0;
+
+    while ((currentNorm > toleranceThreshold) && (iteration < maxIterationsNr))
+    {
+        if (process == PROCESS_HEAT)
+        {
+            currentNorm = GaussSeidelIterationHeat();
+        }
+        if (process == PROCESS_WATER)
+        {
+            // TODO Jacobi method
+
+            if (currentNorm < bestNorm)
+            {
+                bestNorm = currentNorm;
+            }
+            else if (currentNorm > (bestNorm * 10.0))
+            {
+                // non-convergent system
+                return false;
+            }
+        }
+
+        iteration++;
+    }
+
+    return true;
+}
diff --git a/soilFluxes3D/water.cpp b/soilFluxes3D/water.cpp
index 13bb92d22..709c590b3 100644
--- a/soilFluxes3D/water.cpp
+++ b/soilFluxes3D/water.cpp
@@ -201,7 +201,6 @@ double redistribution(long i, TlinkedNode *link, int linkType)
 }
 
 
-
 bool computeFlux(long i, int matrixIndex, TlinkedNode *link, double deltaT, unsigned long myApprox, int linkType)
 {
     if ((*link).index == NOLINK) return false;