feat(options): ✨ save sub-saturations from MIP step

Adds support for saving MIP step sub-saturations

Introduces an option to save sub-saturations during the MIP step.
Updates the `Options` class to include a `save_mip_step` method and a corresponding flag.
Modifies related functions to handle and store sub-saturation data when the option is enabled.
Enhances debugging and analysis capabilities by retaining intermediate MIP results.
Enables a posteriori downscaling capabilities.

Author: djmaxus

.vscode/settings.json

@@ -8,6 +8,7 @@
         "plot",
         "matlab",
         "export",
-        "matlab-repo-init"
+        "matlab-repo-init",
+        "options"
     ]
 }

demo.m

@@ -37,8 +37,10 @@ function demo(args)
 
 sub_rock = downscale_all(grid,rock,mask,downscale_params);
 
+options = Options().save_mip_step(true);
+
 strata_trapped = strata_trapper(grid, sub_rock, params, ...
-    mask=mask, options=Options(), ...
+    mask=mask, options=options, ...
     enable_waitbar=args.show_progress,...
     parfor_arg=args.parfor_arg...
     );

src/Options.m

@@ -3,5 +3,12 @@
         sat_num_points (1,1) uint32 {mustBePositive} = 40
         sat_tol (1,1) double {mustBeNonnegative} = 1e-4
         hydrostatic_correction (1,1) logical = false
+        m_save_mip_step (1,1) {mustBeOfClass(m_save_mip_step,'logical')} = false;
+    end
+
+    methods % builders
+        function self = save_mip_step(self,flag)
+            self.m_save_mip_step = flag;
+        end
     end
 end

src/upscale.m

@@ -1,4 +1,4 @@
-function [perm_upscaled, pc_upscaled, krw, krg] = upscale(...
+function [perm_upscaled, pc_upscaled, krw, krg, mip] = upscale(...
     dr, saturations, params, options, porosities, permeabilities)
 
 if max(porosities,[],'all') <= 0
@@ -28,6 +28,8 @@
 pc_max = params.cap_pressure.func(params.sw_resid,porosities,permeabilities);
 pc_points = linspace(max(pc_max(isfinite(pc_max))),min(entry_pressures(:)),length(saturations));
 
+mip(1:length(saturations)) = struct('sw',nan,'sub_sw',[]);
+
 for index_saturation = 1:length(saturations)
 
     sw_target = saturations(index_saturation);
@@ -65,6 +67,11 @@
 
     [krg(:,index_saturation), krw(:,index_saturation)] = calc_relative_permeabilities( ...
         dr_sub, perm_upscaled_mD, Kg_sub_mD, Kw_sub_mD);
+
+    if options.m_save_mip_step
+        mip(index_saturation).sw = sw_mid;
+        mip(index_saturation).sub_sw = sub_sw;
+    end
 end
 
 % NOTE: we expect only small negative values as computational errors
@@ -128,6 +135,7 @@
 sub_sw(~isfinite(sub_sw)) = 1;
 sw_mid = sum(sub_sw.*pore_volumes,'all')/pore_volume;
 
+% FIXME: Pc should converge as well as Sw
 pc_mid_tot = sum((1-sub_sw).*pore_volumes.*pc_mid,"all")/(pore_volume*(1-sw_mid));
 
 if sw_mid >=1

strata_trapper.m

@@ -35,16 +35,26 @@
 num_sub = zeros(subset_len,1);
 elapsed = zeros(subset_len,1);
 
+if options.m_save_mip_step
+    mip(1:subset_len,1:numel(saturations)) = struct('sw',nan,'sub_sw',[]);
+else
+    mip = struct([]);
+end
+
 % for cell_index = 1:subset_len
 parfor (cell_index = 1:subset_len, args.parfor_arg)
     sub_porosity = sub_rock(cell_index).poro;
     sub_permeability = sub_rock(cell_index).perm;
 
     timer_start = tic;
 
-    [perm_upscaled_cell, pc_upscaled, krw_cell, krg_cell] = upscale(...
+    [perm_upscaled_cell, pc_upscaled, krw_cell, krg_cell, mip_cell] = upscale(...
         DR(cell_index,:), saturations, params, options, sub_porosity, sub_permeability);
 
+    if options.m_save_mip_step
+        mip(cell_index,:) = mip_cell;
+    end
+
     for i = 1:3
         krg_cell(i,:) = monotonize(saturations, krg_cell(i,:), -1);
     end
@@ -78,7 +88,8 @@
     'porosity',poro_upscaled, ...
     'params', params, ...
     'options', args.options, ...
-    'grid', grid ...
+    'grid', grid, ...
+    'mip', mip ...
     );
 
 perf.num_coarse = subset_len;