From dcf3c42a336e10038e37949b3301316fa33ad1ad Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Wed, 5 Jan 2022 09:55:16 +0100
Subject: [PATCH] add perm multiplier aquifer

---
 .../SingleNumericalAquifer.cpp                | 16 +++++++++++++-
 tests/parser/AquiferTests.cpp                 | 21 +++++++++++++++++++
 2 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.cpp b/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.cpp
index 811eb77e5..93221081d 100644
--- a/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.cpp
+++ b/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.cpp
@@ -105,29 +105,43 @@ namespace Opm {
             const auto& cell_dims = grid.getCellDims(gc2);
             double face_area = 0;
             std::string perm_string;
+            std::string mult_string;
             double d = 0.;
             if (con.face_dir == FaceDir::XPlus || con.face_dir == FaceDir::XMinus) {
                 face_area = cell_dims[1] * cell_dims[2];
                 perm_string = "PERMX";
                 d = cell_dims[0];
+                mult_string = "MULTX";
+                if (con.face_dir == FaceDir::XMinus)
+                    mult_string = mult_string + "-";
             }
             if (con.face_dir == FaceDir::YMinus || con.face_dir == FaceDir::YPlus) {
                 face_area = cell_dims[0] * cell_dims[2];
                 perm_string = "PERMY";
                 d = cell_dims[1];
+                mult_string = "MULTY";
+                if (con.face_dir == FaceDir::YMinus)
+                    mult_string = mult_string + "-";
             }
 
             if (con.face_dir == FaceDir::ZMinus || con.face_dir == FaceDir::ZPlus) {
                 face_area = cell_dims[0] * cell_dims[1];
                 perm_string = "PERMZ";
                 d = cell_dims[2];
+                mult_string = "MULTZ";
+                if (con.face_dir == FaceDir::ZMinus)
+                    mult_string = mult_string + "-";
             }
 
             const double trans_cell = (con.trans_option == 0) ?
                                       cell1.transmissiblity() : (2 * cell1.permeability * face_area / cell1.length);
 
             const double cell_perm = (fp.get_double(perm_string))[grid.activeIndex(gc2)];
-            const double trans_con = 2 * cell_perm * face_area * ntg[grid.activeIndex(con.global_index)] / d;
+            double cell_multxyz = 1.0;
+            if (fp.has_double(mult_string))
+                cell_multxyz = (fp.get_double(mult_string))[grid.activeIndex(gc2)];
+
+            const double trans_con = 2 * cell_multxyz * cell_perm * face_area * ntg[grid.activeIndex(con.global_index)] / d;
 
             const double tran = trans_con * trans_cell / (trans_con + trans_cell) * con.trans_multipler;
             if (gc1 < gc2) {
diff --git a/tests/parser/AquiferTests.cpp b/tests/parser/AquiferTests.cpp
index 130017fe9..d3c747874 100644
--- a/tests/parser/AquiferTests.cpp
+++ b/tests/parser/AquiferTests.cpp
@@ -684,6 +684,15 @@ PERMY
     360*1000./
 PERMZ
     360*10./
+
+BOX
+1	8 15 15 3 3 /
+
+MULTY
+ 1e9 1e-9 1.0 2.0 3.0 4.0 5.0 6.0/
+
+ENDBOX
+
 -- setting the three cells for numerical aquifer to be inactive
 ACTNUM
 0 1 0 0 356*1 /
@@ -790,6 +799,18 @@ BOOST_AUTO_TEST_CASE(NumericalAquiferTest)
     BOOST_CHECK(aquifer.numCells() == 3);
     BOOST_CHECK(aquifer.numConnections() == 8 );
 
+    const auto& nncs = aquifer.aquiferConnectionNNCs(grid, ecl_state.fieldProps());
+    BOOST_CHECK(nncs.size() == 8 );
+    // get the half transmissibilites by using small/large multipler m
+    // mab/(ma+b) -> b for ma >> b and -> ma for ma << b
+    const double taq = nncs[0].trans;
+    const double tcell = nncs[1].trans*1.0e9;
+    // now check the multiplier for the rest of the connection i > 2 where m = 1->6
+    for (int i = 2; i < 8; ++i) {
+      const double mult = (i - 1);
+      const double t = mult*tcell*taq / (taq + mult*tcell);
+      BOOST_CHECK_CLOSE(t, nncs[i].trans, 1.0e-6);
+    }
     BOOST_CHECK(grid.getNumActive() == 360);
     // the three aquifer cells are active
     BOOST_CHECK(grid.cellActive(0, 0, 0));