feat: ✨ Hydrostatic correction for MIP

Fixes #57

gen_params.m

@@ -40,6 +40,7 @@
 
 params = Params(...
     TableFunction(krw_data), TableFunction(krg_data),...
-    CapPressure(contact_angle,surface_tension,TableFunction(leverett_j_data),[1,1,0])...
+    CapPressure(contact_angle,surface_tension,TableFunction(leverett_j_data),[1,1,0]),...
+    800, 1000 ...
     );
 end

src/Options.m

@@ -2,6 +2,6 @@
     properties
         sat_num_points (1,1) uint32 {mustBePositive} = 40
         sat_tol (1,1) double {mustBeNonnegative} = 1e-4
+        hydrostatic_correction (1,1) logical = false
     end
 end
-

src/Params.m

@@ -3,17 +3,21 @@
         krw (1,1) TableFunction = TableFunction([]);
         krg (1,1) TableFunction = TableFunction([]);
         cap_pressure  (1,1) CapPressure
+        rho_gas (1,1) double
+        rho_water (1,1) double {mustBeFinite}
         sw_resid (1,1) double
     end
-    
+
     methods
-        function obj = Params(krw,krg,cap_pressure)
+        function obj = Params(krw,krg,cap_pressure, rho_gas,rho_water)
             obj.krw = krw;
             obj.krg = krg;
             obj.cap_pressure = cap_pressure;
 
             obj.sw_resid = obj.krw.data(1,1);
+
+            obj.rho_gas   = rho_gas;
+            obj.rho_water = rho_water;
         end
     end
 end
-

src/calc_percolation.m

@@ -1,6 +1,16 @@
-function [invasion,num_iter] = calc_percolation(p_boundary,p_entry)
+function [invasion,num_iter] = calc_percolation(p_boundary,p_entry, ...
+    include_gravity, Lz, rho_water, rho_gas)
+
+h_ref = Lz*0.5;
+
+Nz = size(p_entry,3);
+h(1,1,1:Nz) = linspace(Lz/Nz/2,Lz - Lz/Nz/2,Nz);
+
+hydrostatic_correction = include_gravity * isfinite(rho_gas) * std_gravity() * ...
+    (rho_water - rho_gas) * (h - h_ref);
+
 invasion = false(size(p_entry));
-invadable = p_entry < p_boundary;
+invadable = (p_entry + hydrostatic_correction) < p_boundary;
 idx_invadable = find(invadable)';
 [I,J,K] = ind2sub(size(invadable),idx_invadable);
 is_invading = true;

src/units/gram.m

@@ -0,0 +1,4 @@
+function out = gram()
+    out = 1e-3;
+end
+

src/units/sec.m

@@ -0,0 +1,4 @@
+function s = sec()
+    s = 1;
+end
+

src/units/std_gravity.m

@@ -0,0 +1,4 @@
+function g = std_gravity()
+    g = 9.80665 * meter() / sec()^2;
+end
+

src/upscale.m

@@ -37,9 +37,7 @@
         [pc_mid_tot, sw_mid, pc_mid, sub_sw, converged] = mip_iteration(...
             sw_target, dr, entry_pressures, porosities, permeabilities, pc_mid, ...
             Nz_sub, Nx_sub, Ny_sub,...
-            params.cap_pressure,...
-            options.sat_tol ...
-            );
+            params, options);
 
         if converged
             break;
@@ -84,17 +82,17 @@
 function [pc_mid_tot, sw_mid, pc_mid, sub_sw, converged] = mip_iteration(...
     sw_target, dr, entry_pressures, porosities, permeabilities, pc_mid,...
     Nz_sub, Nx_sub, Ny_sub, ...
-    cap_pressure,...
-    tol_sw)
+    params, options)
 
-invaded_mat_mid = calc_percolation(pc_mid, entry_pressures);
+invaded_mat_mid = calc_percolation(pc_mid, entry_pressures,...
+    options.hydrostatic_correction, dr(3), params.rho_water, params.rho_gas);
 
 volume = prod(dr);
 sub_volume = volume./double(Nz_sub*Nx_sub*Ny_sub);
 pore_volumes = porosities .* sub_volume;
 pore_volume = sum(pore_volumes,'all');
 
-sub_sw = invaded_mat_mid .* cap_pressure.inv(pc_mid,porosities,permeabilities) + ~invaded_mat_mid .* 1;
+sub_sw = invaded_mat_mid .* params.cap_pressure.inv(pc_mid,porosities,permeabilities) + ~invaded_mat_mid .* 1;
 sub_sw(~isfinite(sub_sw)) = 1;
 sw_mid = sum(sub_sw.*pore_volumes,'all')/pore_volume;
 
@@ -106,12 +104,12 @@
 
 sw_err = sw_target - sw_mid;
 err = abs(sw_err);
-converged = err <= tol_sw;
+converged = err <= options.sat_tol;
 if converged
     return;
 end
 
-deriv = cap_pressure.deriv(sw_mid, mean(porosities,'all'), mean(permeabilities,'all'));
+deriv = params.cap_pressure.deriv(sw_mid, mean(porosities,'all'), mean(permeabilities,'all'));
 
 dpc = sw_err*deriv;