From 37079e412df069f078a4ddb79658f9dd06cafc79 Mon Sep 17 00:00:00 2001
From: Maziar Veyskarami <maziar.veyskarami@iws.uni-stuttgart.de>
Date: Tue, 16 Apr 2024 15:09:45 +0000
Subject: [PATCH] Cleanup/pnm 2pnc

---
 CHANGELOG.md                                  |  2 +-
 test/porenetwork/2pnc/CMakeLists.txt          |  6 +-
 .../2pnc/grids/{1dGrid.dgf => 1d_grid.dgf}    | 18 ++---
 test/porenetwork/2pnc/main.cc                 |  8 --
 test/porenetwork/2pnc/params.input            |  7 +-
 test/porenetwork/2pnc/params_ni.input         |  8 +-
 test/porenetwork/2pnc/problem.hh              | 78 +++++++------------
 test/references/test_pnm_2pnc-reference.vtp   | 60 +++++++-------
 .../references/test_pnm_2pnc_ni-reference.vtp | 62 +++++++--------
 9 files changed, 106 insertions(+), 143 deletions(-)
 rename test/porenetwork/2pnc/grids/{1dGrid.dgf => 1d_grid.dgf} (64%)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f2a57c7103..5b915e08e4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,7 +14,7 @@ The term is added as a source term in the problem using the new helper function
 The function uses new spatial parameter interface implemented in the new `BrinkmanSpatialParams` class (`dumux/freeflow/spatialparams.hh`). The helper function can deal with isotropic and anisotropic permeabilites.
 - __Facet-Coupling__: Fixed the handling of duplicate degrees of freedom in the box facet-coupling model in the corner case that an internal fracture turns into a boundary fracture (see [merge request 3748](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/merge_requests/3748) for images).
 - __Periodic Boundaries__: Fixed an issue for vector-valued unknowns. Other schemes that provide a periodic map at boundaries now also support periodicity.
-
+- __Pore network__: Added a model and a test case for two-phase compositional fluid flow.
 ### Immediate interface changes not allowing/requiring a deprecation period:
 - __RichardsNewtonSolver__: It is now possible to select the MPICommunicator used by the RichardsNewtonSolver (e.g., either real or dummy communicator)
 - __CompositionalFluidState__: setRelativeHumidity was removed. Use the other setters. This setter was removed because it was very specific with a lot of specific prerequisites not fitting the general concept of the class. It was also outdated and not used in any example or test and didn't fit the index convention used in the fluid systems anymore.
diff --git a/test/porenetwork/2pnc/CMakeLists.txt b/test/porenetwork/2pnc/CMakeLists.txt
index 04791d2e1a..9b8735774e 100644
--- a/test/porenetwork/2pnc/CMakeLists.txt
+++ b/test/porenetwork/2pnc/CMakeLists.txt
@@ -12,7 +12,7 @@ dumux_add_test(NAME test_pnm_2pnc
                CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
                CMD_ARGS      --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2pnc-reference.vtp
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc-00010.vtp
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc-00013.vtp
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc")
 
 dumux_add_test(NAME test_pnm_2pnc_ni
@@ -23,5 +23,5 @@ dumux_add_test(NAME test_pnm_2pnc_ni
                CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
                CMD_ARGS      --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_2pnc_ni-reference.vtp
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni-00004.vtp
-                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni params_ni.input -Problem.Name test_pnm_2pnc_ni")
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni-00014.vtp
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_2pnc_ni params_ni.input")
diff --git a/test/porenetwork/2pnc/grids/1dGrid.dgf b/test/porenetwork/2pnc/grids/1d_grid.dgf
similarity index 64%
rename from test/porenetwork/2pnc/grids/1dGrid.dgf
rename to test/porenetwork/2pnc/grids/1d_grid.dgf
index e3aa06359f..17c29c0198 100644
--- a/test/porenetwork/2pnc/grids/1dGrid.dgf
+++ b/test/porenetwork/2pnc/grids/1d_grid.dgf
@@ -5,16 +5,16 @@ DGF
 % Element parameters: ThroatInscribedRadius ThroatLength ThroatLabel
 Vertex % Coordinates and volumes of the pore bodies
 parameters 2
-0 0 0 0.0002263 2.0
-1 0 0 0.0002263 -1.0
-2 0 0 0.0002263 -1.0
-3 0 0 0.0002263 -1.0
-4 0 0 0.0002263 3.0
+0 0 0 2e-3 2.0
+1 0 0 2e-3 -1.0
+2 0 0 2e-3 -1.0
+3 0 0 2e-3 -1.0
+4 0 0 2e-3 3.0
 #
 SIMPLEX % Connections of the pore bodies (pore throats)
 parameters 3
-0 1 3.3304e-05 6.6609e-05 2
-1 2 3.3304e-05 6.6609e-05 -1
-2 3 3.3304e-05 6.6609e-05 -1
-3 4 3.3304e-05 6.6609e-05 3
+0 1 1e-04 3e-05 2
+1 2 1e-04 3e-05 -1
+2 3 1e-04 3e-05 -1
+3 4 1e-04 3e-05 3
 #
diff --git a/test/porenetwork/2pnc/main.cc b/test/porenetwork/2pnc/main.cc
index b12fb47681..3d70f36923 100644
--- a/test/porenetwork/2pnc/main.cc
+++ b/test/porenetwork/2pnc/main.cc
@@ -11,14 +11,6 @@
  */
 #include <config.h>
 
-#include <ctime>
-#include <iostream>
-
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/timer.hh>
-#include <dune/grid/io/file/dgfparser/dgfexception.hh>
-#include <dune/grid/io/file/vtk.hh>
-
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/common/initialize.hh>
 #include <dumux/common/properties.hh>
diff --git a/test/porenetwork/2pnc/params.input b/test/porenetwork/2pnc/params.input
index 42f13ccf19..4ad23fe3c1 100644
--- a/test/porenetwork/2pnc/params.input
+++ b/test/porenetwork/2pnc/params.input
@@ -1,9 +1,9 @@
 [TimeLoop]
 DtInitial = 1e-5 # [s]
-TEnd = 0.01 # [s]
+TEnd = 1e2 # [s]
 
 [Grid]
-File = ./grids/1dGrid.dgf
+File = ./grids/1d_grid.dgf
 PoreGeometry = Cube
 ThroatCrossSectionShape = Square
 
@@ -11,9 +11,6 @@ ThroatCrossSectionShape = Square
 Name = test_pnm_2pnc
 VtpOutputFrequency = 10 # Write every n-th time step. 0 only writes a file if an invasion / snap-off occurred. -1 writes every step
 EnableGravity = false
-CapillaryPressure = 5000
-UseFixedPressureAndSaturationBoundary = false
-Source = 1e-5
 
 [Vtk]
 AddVelocity = 1
diff --git a/test/porenetwork/2pnc/params_ni.input b/test/porenetwork/2pnc/params_ni.input
index 90c30695d6..83b222700c 100644
--- a/test/porenetwork/2pnc/params_ni.input
+++ b/test/porenetwork/2pnc/params_ni.input
@@ -1,18 +1,16 @@
 [TimeLoop]
 DtInitial = 1e-5 # [s]
-TEnd = 0.01 # [s]
+TEnd = 1e2 # [s]
 
 [Grid]
-File = ./grids/1dGrid.dgf
+File = ./grids/1d_grid.dgf
 PoreGeometry = Cube
 ThroatCrossSectionShape = Square
 
 [Problem]
-Name = test_pnm2pnc_ni
+Name = test_pnm_2pnc_ni
 VtpOutputFrequency = 10 # Write every n-th time step. 0 only writes a file if an invasion / snap-off occurred
-CapillaryPressure = 5000
 EnableGravity = false
-UseFixedPressureAndSaturationBoundary = true
 InletTemperature = 288.15
 OutletTemperature = 283.15
 Source = 1e-5
diff --git a/test/porenetwork/2pnc/problem.hh b/test/porenetwork/2pnc/problem.hh
index bf71e5a8ab..5b34d54332 100644
--- a/test/porenetwork/2pnc/problem.hh
+++ b/test/porenetwork/2pnc/problem.hh
@@ -45,10 +45,7 @@ public:
     : ParentType(gridGeometry, spatialParams)
     {
         vtpOutputFrequency_ = getParam<int>("Problem.VtpOutputFrequency");
-        useFixedPressureAndSaturationBoundary_ = getParam<bool>("Problem.UseFixedPressureAndSaturationBoundary", false);
-        pc_ = getParam<Scalar>("Problem.CapillaryPressure");
-        source_ = getParam<Scalar>("Problem.Source");
-        inletPressure_ = getParam<Scalar>("Problem.InletPressure", 1e5);
+        inletPressure_ = getParam<Scalar>("Problem.InletPressure", 1.1e5);
         outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1e5);
 #if !ISOTHERMAL
         inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 288.15);
@@ -84,14 +81,12 @@ public:
     {
         BoundaryTypes bcTypes;
 
-        // If a global phase pressure difference (pn,inlet - pw,outlet) with fixed saturations is specified, use a Dirichlet BC here
-        if (useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+        // Use Dirichlet BCs for both inlet and outlet
+        if (isInletPore_(scv) || isOutletPore_(scv))
             bcTypes.setAllDirichlet();
-        else if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
-            bcTypes.setAllNeumann();
-        else if (isOutletPore_(scv))
-            bcTypes.setAllDirichlet();
-
+#if !ISOTHERMAL
+        bcTypes.setDirichlet(Indices::temperatureIdx);
+#endif
         return bcTypes;
     }
 
@@ -101,30 +96,26 @@ public:
                                const SubControlVolume& scv) const
     {
         PrimaryVariables values(0.0);
-        values[Indices::pressureIdx] = 1e5;
-        values[Indices::switchIdx] = 0.0;
 
-        // If a global phase pressure difference (pn,inlet - pw,outlet) is specified and the saturation shall also be fixed, apply:
-        // pw,inlet = pw,outlet = 1e5; pn,outlet = pw,outlet + pc(S=0) = pw,outlet; pn,inlet = pw,inlet + pc_
-        if (useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
+        if (isInletPore_(scv))
         {
             values.setState(Indices::bothPhases);
             values[Indices::pressureIdx] = inletPressure_;
-            values[Indices::switchIdx] = 1.0 - this->spatialParams().fluidMatrixInteraction(element, scv, int()/*dummyElemsol*/).sw(pc_);
-#if !ISOTHERMAL
-            values[Indices::temperatureIdx] = inletTemperature_;
-#endif
+            values[Indices::switchIdx] = 1.0;
         }
-        else if (isOutletPore_(scv))
+        else
         {
-            values.setState(Indices::firstPhaseOnly);
+            values.setState(Indices::bothPhases);
             values[Indices::pressureIdx] = outletPressure_;
             values[Indices::switchIdx] = 0.0;
+        }
+
 #if !ISOTHERMAL
+        if (isInletPore_(scv))
+            values[Indices::temperatureIdx] = inletTemperature_;
+        else
             values[Indices::temperatureIdx] = outletTemperature_;
 #endif
-        }
-
         return values;
     }
 
@@ -142,19 +133,6 @@ public:
                             const SubControlVolume& scv) const
     {
         PrimaryVariables values(0.0);
-
-        // for isothermal case, we fix injection rate of non-wetting phase at inlet
-        // for non-isothermal case, we fix injection of air enthalpy at inlet
-        if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv))
-        {
-            values[Indices::conti0EqIdx + 1] = source_/scv.volume();
-#if !ISOTHERMAL
-            const auto pressure = elemVolVars[scv].pressure(1);
-            const auto airEnthalpy = Components::Air<Scalar>::gasEnthalpy(inletTemperature_, pressure);
-            values[Indices::temperatureIdx] = airEnthalpy * source_ * Components::Air<Scalar>::molarMass()/scv.volume();
-#endif
-        }
-
         return values;
     }
     // \}
@@ -163,24 +141,25 @@ public:
     PrimaryVariables initial(const Vertex& vertex) const
     {
         PrimaryVariables values(0.0);
+
+        values.setState(Indices::bothPhases);
         values[Indices::pressureIdx] = outletPressure_;
+        values[Indices::switchIdx] = 0.0;
+
+#if !ISOTHERMAL
+        values[Indices::temperatureIdx] = outletTemperature_;
+#endif
 
-        // get global index of pore
         const auto dofIdxGlobal = this->gridGeometry().vertexMapper().index(vertex);
         if (isInletPore_(dofIdxGlobal))
         {
-            values.setState(Indices::firstPhaseOnly);
-            values[Indices::switchIdx] = 0.0;
-        }
-        else
-        {
-            values.setState(Indices::firstPhaseOnly);
-            values[Indices::switchIdx] = 0.0;
-        }
-
+            values.setState(Indices::bothPhases);
+            values[Indices::pressureIdx] = inletPressure_;
+            values[Indices::switchIdx] = 1.0;
 #if !ISOTHERMAL
-        values[Indices::temperatureIdx] = outletTemperature_;
+            values[Indices::temperatureIdx] = inletTemperature_;
 #endif
+        }
 
         return values;
     }
@@ -209,9 +188,6 @@ private:
     }
 
     int vtpOutputFrequency_;
-    bool useFixedPressureAndSaturationBoundary_;
-    Scalar pc_;
-    Scalar source_;
     Scalar inletPressure_;
     Scalar outletPressure_;
 #if !ISOTHERMAL
diff --git a/test/references/test_pnm_2pnc-reference.vtp b/test/references/test_pnm_2pnc-reference.vtp
index 5ad168b997..2470cf92d9 100644
--- a/test/references/test_pnm_2pnc-reference.vtp
+++ b/test/references/test_pnm_2pnc-reference.vtp
@@ -4,70 +4,70 @@
     <Piece NumberOfLines="4" NumberOfPoints="5">
       <PointData Scalars="S_liq">
         <DataArray type="Float32" Name="S_liq" NumberOfComponents="1" format="ascii">
-          0.0374111 0.0474591 0.0343665 0.0391619 1
+          0 0 0 0.0020251 1
         </DataArray>
         <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
-          99081.6 99126.6 97908.1 97768.2 100000
+          110000 107090 104104 101247 100000
         </DataArray>
         <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
-          999.712 999.712 999.711 999.711 999.701
+          999.718 973.791 960.601 999.713 999.712
         </DataArray>
         <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
-          765.754 765.754 765.753 765.753 765.754
+          765.759 765.758 765.756 765.755 765.754
         </DataArray>
         <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
-          0.962589 0.952541 0.965634 0.960838 0
+          1 1 1 0.997975 0
         </DataArray>
         <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
-          101923 101441 100971 100498 100000
+          112159 109250 106263 103192 100000
         </DataArray>
         <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
-          1.24808 1.24214 1.23636 1.23055 0.00940657
+          1.37399 1.33835 1.30168 1.26368 1.22442
         </DataArray>
         <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
-          56885.7 56887.5 56889.3 56891.1 103741
+          56849.8 56851.5 56857.4 56880.9 56893.1
         </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
-          2841.63 2314.36 3062.68 2730.15 0
+          2159.24 2159.24 2159.24 1944.36 0
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_liq" NumberOfComponents="1" format="ascii">
-          0.999981 0.999982 0.999982 0.999982 1
+          0.99998 0.974047 0.960855 0.999981 0.999982
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_liq" NumberOfComponents="1" format="ascii">
-          0.99997 0.99997 0.99997 0.999971 1
+          0.999967 0.974034 0.960844 0.99997 0.999971
         </DataArray>
         <DataArray type="Float32" Name="x^Air_liq" NumberOfComponents="1" format="ascii">
-          1.8543e-05 1.84542e-05 1.83676e-05 1.82806e-05 0
+          2.0428e-05 1.98981e-05 1.9351e-05 1.87766e-05 1.81889e-05
         </DataArray>
         <DataArray type="Float32" Name="X^Air_liq" NumberOfComponents="1" format="ascii">
-          2.98082e-05 2.96654e-05 2.95262e-05 2.93864e-05 0
+          3.28382e-05 3.19864e-05 3.1107e-05 3.01837e-05 2.92389e-05
         </DataArray>
         <DataArray type="Float32" Name="rhoMolar_liq" NumberOfComponents="1" format="ascii">
-          55492.1 55492.1 55492.1 55492.1 55492.2
+          55492.4 55492.3 55492.3 55492.2 55492.2
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.0120499 0.0121072 0.0121635 0.0122207 0.0122818
+          0.0109501 0.0109502 0.0111055 0.0119017 0.0122816
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.00753016 0.00756613 0.00760153 0.00763743 0.0122818
+          0.00684007 0.00684012 0.00693756 0.00743717 0.00767567
         </DataArray>
         <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
-          0.98795 0.987893 0.987836 0.987779 0
+          0.98905 0.98905 0.988894 0.988098 0.987718
         </DataArray>
         <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
-          0.99247 0.992434 0.992399 0.992363 0
+          0.99316 0.99316 0.993062 0.992563 0.992324
         </DataArray>
         <DataArray type="Float32" Name="rhoMolar_gas" NumberOfComponents="1" format="ascii">
-          43.2939 43.089 42.8893 42.6886 0.522147
+          47.6418 46.4059 45.1371 43.8327 42.477
         </DataArray>
         <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
-          3 3 3 3 1
+          3 2 2 3 3
         </DataArray>
         <DataArray type="Float32" Name="poreInscribedRadius" NumberOfComponents="1" format="ascii">
-          0.0002263 0.0002263 0.0002263 0.0002263 0.0002263
+          0.002 0.002 0.002 0.002 0.002
         </DataArray>
         <DataArray type="Float32" Name="coordinationNumber" NumberOfComponents="1" format="ascii">
           1 2 2 2 1
@@ -78,10 +78,10 @@
       </PointData>
       <CellData Scalars="process rank" Vectors="velocity_liq (m/s)">
         <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
-          -0.00296711 -0 -0 0.0692308 0 0 0.00794862 0 0 -0.159574 -0 -0
+          0.738437 0 0 0.758079 0 0 0.72489 0 0 0.390444 0 0
         </DataArray>
         <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
-          60.2882 0 0 59.377 0 0 59.6628 0 0 61.9152 0 0
+          7704.64 0 0 7909.83 0 0 8133.44 0 0 8430.99 0 0
         </DataArray>
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
           0 0 0 0
@@ -90,25 +90,25 @@
           2 -1 -1 3
         </DataArray>
         <DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
-          3.3304e-05 3.3304e-05 3.3304e-05 3.3304e-05
+          0.0001 0.0001 0.0001 0.0001
         </DataArray>
         <DataArray type="Float32" Name="throatLength" NumberOfComponents="1" format="ascii">
-          6.6609e-05 6.6609e-05 6.6609e-05 6.6609e-05
+          3e-05 3e-05 3e-05 3e-05
         </DataArray>
         <DataArray type="Float32" Name="pcEntry" NumberOfComponents="1" format="ascii">
-          4106.16 4106.16 4106.16 4106.16
+          1367.51 1367.51 1367.51 1367.51
         </DataArray>
         <DataArray type="Float32" Name="transmissibilityW" NumberOfComponents="1" format="ascii">
-          4.81607e-17 3.56908e-17 3.56908e-17 5.65214e-17
+          3.20753e-16 3.20753e-16 3.20753e-16 4.87829e-16
         </DataArray>
         <DataArray type="Float32" Name="transmissibilityN" NumberOfComponents="1" format="ascii">
-          8.5211e-15 8.78119e-15 8.78119e-15 8.36682e-15
+          1.81814e-12 1.81814e-12 1.81814e-12 1.80211e-12
         </DataArray>
         <DataArray type="Float32" Name="volumeFluxW" NumberOfComponents="1" format="ascii">
-          1.65794e-12 3.33017e-11 3.82348e-12 9.65956e-11
+          7.14632e-10 7.3364e-10 7.01522e-10 4.6599e-10
         </DataArray>
         <DataArray type="Float32" Name="volumeFluxN" NumberOfComponents="1" format="ascii">
-          2.33789e-07 2.34872e-07 2.36002e-07 2.37215e-07
+          0.000300729 0.000308738 0.000317467 0.000327177
         </DataArray>
       </CellData>
       <Points>
diff --git a/test/references/test_pnm_2pnc_ni-reference.vtp b/test/references/test_pnm_2pnc_ni-reference.vtp
index be60712cea..23c6593dac 100644
--- a/test/references/test_pnm_2pnc_ni-reference.vtp
+++ b/test/references/test_pnm_2pnc_ni-reference.vtp
@@ -4,73 +4,73 @@
     <Piece NumberOfLines="4" NumberOfPoints="5">
       <PointData Scalars="S_liq">
         <DataArray type="Float32" Name="S_liq" NumberOfComponents="1" format="ascii">
-          0.0200785 0.0261951 0.948722 1 1
+          0 0 0 0.00293691 1
         </DataArray>
         <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
-          100000 101052 104283 102142 100000
+          110000 107089 104100 101346 100000
         </DataArray>
         <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
-          999.111 999.711 999.715 999.705 999.701
+          999.117 974.141 986.478 999.227 999.712
         </DataArray>
         <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
-          879.067 766.229 765.823 765.781 765.754
+          879.071 878.683 869.94 860.79 765.754
         </DataArray>
         <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
-          0.979922 0.973805 0.0512775 0 0
+          1 1 1 0.997063 0
         </DataArray>
         <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
-          105000 104963 104925 102142 100000
+          112159 109248 106259 103193 100000
         </DataArray>
         <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
-          1.26143 1.28536 1.28499 0.379949 0.00940657
+          1.34797 1.31306 1.27869 1.24339 1.22442
         </DataArray>
         <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
-          56207.5 56871.4 56874.1 57633.3 103741
+          56175.8 56179.9 56246.7 56317 56893.1
         </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
-          5000 3911.23 641.727 0 -0
+          2159.24 2159.24 2159.24 1847.61 0
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_liq" NumberOfComponents="1" format="ascii">
-          0.999983 0.999981 0.999981 0.999994 1
+          0.999982 0.974984 0.987279 0.999983 0.999982
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_liq" NumberOfComponents="1" format="ascii">
-          0.999973 0.999969 0.999969 0.999991 1
+          0.999971 0.974973 0.987268 0.999973 0.999971
         </DataArray>
         <DataArray type="Float32" Name="x^Air_liq" NumberOfComponents="1" format="ascii">
-          1.70916e-05 1.90931e-05 1.90944e-05 5.54692e-06 0
+          1.82763e-05 1.78078e-05 1.74492e-05 1.70791e-05 1.81889e-05
         </DataArray>
         <DataArray type="Float32" Name="X^Air_liq" NumberOfComponents="1" format="ascii">
-          2.74751e-05 3.06925e-05 3.06946e-05 8.91683e-06 0
+          2.93794e-05 2.86264e-05 2.80498e-05 2.74549e-05 2.92389e-05
         </DataArray>
         <DataArray type="Float32" Name="rhoMolar_liq" NumberOfComponents="1" format="ascii">
-          55458.8 55492.1 55492.3 55492.2 55492.2
+          55459.1 55459.1 55462.2 55465.3 55492.2
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.0162449 0.0117178 0.0117075 0.0120252 0.0122818
+          0.015208 0.0152067 0.0154517 0.0157083 0.0122816
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.0101679 0.00732173 0.00731526 0.00759294 0.0122818
+          0.00951512 0.0095143 0.00966851 0.00983002 0.00767567
         </DataArray>
         <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
-          0.983755 0.988282 0.988293 0.294909 0
+          0.984792 0.984793 0.984548 0.984292 0.987718
         </DataArray>
         <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
-          0.989832 0.992678 0.992685 0.299341 0
+          0.990485 0.990486 0.990331 0.99017 0.992324
         </DataArray>
         <DataArray type="Float32" Name="rhoMolar_gas" NumberOfComponents="1" format="ascii">
-          43.8272 44.5816 44.5685 13.3171 0.522147
+          46.8154 45.6028 44.4133 43.1913 42.477
         </DataArray>
         <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
-          3 3 3 1 1
+          3 2 2 3 3
         </DataArray>
         <DataArray type="Float32" Name="poreInscribedRadius" NumberOfComponents="1" format="ascii">
-          0.0002263 0.0002263 0.0002263 0.0002263 0.0002263
+          0.002 0.002 0.002 0.002 0.002
         </DataArray>
         <DataArray type="Float32" Name="T" NumberOfComponents="1" format="ascii">
-          288.15 283.172 283.153 283.151 283.15
+          288.15 288.133 287.757 287.361 283.15
         </DataArray>
         <DataArray type="Float32" Name="coordinationNumber" NumberOfComponents="1" format="ascii">
           1 2 2 2 1
@@ -81,10 +81,10 @@
       </PointData>
       <CellData Scalars="process rank" Vectors="velocity_liq (m/s)">
         <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
-          -0.0224304 -0 -0 -0.112597 -0 -0 3.83914 0 0 3.83908 0 0
+          0.848257 0 0 0.870363 0 0 0.794098 0 0 0.524429 0 0
         </DataArray>
         <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
-          4.82473 0 0 4.88488 0 0 0 0 0 0 0 0
+          7618.25 0 0 7820.81 0 0 8032.27 0 0 8337.29 0 0
         </DataArray>
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
           0 0 0 0
@@ -93,25 +93,25 @@
           2 -1 -1 3
         </DataArray>
         <DataArray type="Float32" Name="throatInscribedRadius" NumberOfComponents="1" format="ascii">
-          3.3304e-05 3.3304e-05 3.3304e-05 3.3304e-05
+          0.0001 0.0001 0.0001 0.0001
         </DataArray>
         <DataArray type="Float32" Name="throatLength" NumberOfComponents="1" format="ascii">
-          6.6609e-05 6.6609e-05 6.6609e-05 6.6609e-05
+          3e-05 3e-05 3e-05 3e-05
         </DataArray>
         <DataArray type="Float32" Name="pcEntry" NumberOfComponents="1" format="ascii">
-          4106.16 4106.16 4106.16 4106.16
+          1367.51 1367.51 1367.51 1367.51
         </DataArray>
         <DataArray type="Float32" Name="transmissibilityW" NumberOfComponents="1" format="ascii">
-          5.02435e-18 1.34185e-17 1.03853e-14 1.03853e-14
+          3.20753e-16 3.20753e-16 3.20753e-16 5.98319e-16
         </DataArray>
         <DataArray type="Float32" Name="transmissibilityN" NumberOfComponents="1" format="ascii">
-          9.81542e-15 9.41615e-15 0 0
+          1.81814e-12 1.81814e-12 1.81814e-12 1.79303e-12
         </DataArray>
         <DataArray type="Float32" Name="volumeFluxW" NumberOfComponents="1" format="ascii">
-          4.04824e-12 3.32099e-11 1.70328e-08 1.70325e-08
+          8.20912e-10 8.42305e-10 7.68498e-10 6.93165e-10
         </DataArray>
         <DataArray type="Float32" Name="volumeFluxN" NumberOfComponents="1" format="ascii">
-          2.05347e-08 2.02316e-08 0 0
+          0.000297357 0.000305264 0.000313517 0.000322472
         </DataArray>
       </CellData>
       <Points>
-- 
GitLab