From 16b361c85b4f263d1c7459839a17bf5789e9634c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@uni-potsdam.de>
Date: Thu, 6 Feb 2025 16:18:04 +0100
Subject: [PATCH] fix(velocities): prevent redundant velocity calculations

---
 include/tug/Advection/Advection.hpp  | 10 +++++++---
 include/tug/Advection/Velocities.hpp |  6 ++++++
 test/testVelocities.cpp              | 12 +++---------
 3 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp
index ff448c5..6fd52b7 100644
--- a/include/tug/Advection/Advection.hpp
+++ b/include/tug/Advection/Advection.hpp
@@ -119,13 +119,17 @@ public:
   void run() {
     this->setDomain(velocities.domainX(), velocities.domainY());
 
-    velocities.run();
+    if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
+      if (!velocities.isSteady()) {
+        velocities.run();
+      }
+    }
 
     for (int i = 0; i < this->getIterations(); i++) {
       if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) {
         velocities.run();
       }
-      adv();
+      startAdvection();
     }
   }
 
@@ -133,7 +137,7 @@ private:
   /**
    * @brief function calculating material transport for one timestep
    */
-  void adv() {
+  void startAdvection() {
     const std::size_t rows = this->rows();
     const std::size_t cols = this->cols();
     const T volume = this->deltaCol() * this->deltaRow();
diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp
index 8cf2fc3..b98a195 100644
--- a/include/tug/Advection/Velocities.hpp
+++ b/include/tug/Advection/Velocities.hpp
@@ -43,6 +43,7 @@ private:
   T timestep{-1};
   T epsilon{1E-5};
   int numThreads{omp_get_num_procs()};
+  bool steady{false};
 
   RowMajMat<T> velocitiesX;
   RowMajMat<T> velocitiesY;
@@ -118,6 +119,8 @@ public:
     this->permKY = alphaY;
   }
 
+  bool isSteady() const { return steady; }
+
   /**
    * @brief Set the timestep per iteration
    *
@@ -219,6 +222,9 @@ public:
         oldConcentrations = this->getConcentrationMatrix();
         (void)calculate_hydraulic_flow(input);
       } while (!checkConvergance(oldConcentrations));
+
+      steady = true;
+
     } else {
       if (timestep == -1) {
         throw_invalid_argument("Timestep is not set");
diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp
index 0c8af66..40f16d2 100644
--- a/test/testVelocities.cpp
+++ b/test/testVelocities.cpp
@@ -20,9 +20,7 @@ VELOCITIES_TEST(SteadyStateCenter) {
   tug::RowMajMat<double> hydHeads =
       tug::RowMajMat<double>::Constant(rows, cols, 1);
 
-  tug::RowMajMat<double> permKX =
-      tug::RowMajMat<double>::Constant(rows, cols, K);
-  tug::RowMajMat<double> permKY =
+  tug::RowMajMat<double> permK =
       tug::RowMajMat<double>::Constant(rows, cols, K);
 
   tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
@@ -30,14 +28,10 @@ VELOCITIES_TEST(SteadyStateCenter) {
       velo(hydHeads);
 
   velo.setDomain(100, 100);
-  velo.setPermKX(permKX);
-  velo.setPermKY(permKY);
+  velo.setPermKX(permK);
+  velo.setPermKY(permK);
 
   tug::Boundary<double> &bcH = velo.getBoundaryConditions();
-  bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1);
-  bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
-  bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1);
-  bcH.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1);
 
   bcH.setInnerBoundary(center_row, center_col, 10);