feat: ✨ Relative permeability via new TableFunction class

Also, `params.rel_perm` nesting was removed.

plot_result.m

@@ -47,19 +47,19 @@ function curves_plot(mask, strata_trapped, params,scale)
 
 tiles_krw = tiledlayout(figure(),3,2,TileSpacing='compact',Padding='tight');
 
-stat_plot(nexttile(tiles_krw),'Water, along X-axis','',strata_trapped.saturation,params.rel_perm.calc_krw,sub_data(strata_trapped.rel_perm_wat,mask,1));
+stat_plot(nexttile(tiles_krw),'Water, along X-axis','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,mask,1));
 yscale(scale);
-stat_plot(nexttile(tiles_krw),'Gas, along X-axis','',strata_trapped.saturation,@(sw) params.rel_perm.calc_krg(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,1));
+stat_plot(nexttile(tiles_krw),'Gas, along X-axis','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,1));
 yscale(scale);
 
-stat_plot(nexttile(tiles_krw),'Water, along Y-axis','',strata_trapped.saturation,params.rel_perm.calc_krw,sub_data(strata_trapped.rel_perm_wat,mask,2));
+stat_plot(nexttile(tiles_krw),'Water, along Y-axis','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,mask,2));
 yscale(scale);
-stat_plot(nexttile(tiles_krw),'Gas, along Y-axis','',strata_trapped.saturation,@(sw) params.rel_perm.calc_krg(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,2));
+stat_plot(nexttile(tiles_krw),'Gas, along Y-axis','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,2));
 yscale(scale);
 
-stat_plot(nexttile(tiles_krw),'Water, along Z-axis','',strata_trapped.saturation,params.rel_perm.calc_krw,sub_data(strata_trapped.rel_perm_wat,mask,3));
+stat_plot(nexttile(tiles_krw),'Water, along Z-axis','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,mask,3));
 yscale(scale);
-stat_plot(nexttile(tiles_krw),'Gas, along Z-axis','',strata_trapped.saturation,@(sw) params.rel_perm.calc_krg(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,3));
+stat_plot(nexttile(tiles_krw),'Gas, along Z-axis','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,mask,3));
 yscale(scale);
 
 xlabel(tiles_krw,'Wetting phase saturation');

src/TableFunction.m

@@ -0,0 +1,47 @@
+classdef TableFunction
+    properties
+        data (:,2) double
+    end
+
+    methods
+        function obj = TableFunction(data)
+            obj.data = data;
+        end
+
+        function out = func(obj,arg)
+            out = interp1(obj.data(:,1), obj.data(:,2), arg, "linear",'extrap');
+        end
+
+        function out = inv(obj,arg)
+            out = interp1(obj.data(:,2), obj.data(:,1), arg, "linear",'extrap');
+        end
+
+        function out = deriv(obj,arg)
+            out = compute_derivative(obj.data(:,1), obj.data(:,2), arg);
+        end
+    end
+end
+
+function df = compute_derivative(x, f, x0)
+    % x: vector of x values
+    % f: vector of f(x) values
+    % x0: point at which to compute the derivative
+
+    % Find the index of the closest point to x0
+    [~, idx] = min(abs(x - x0));
+
+    % Check if we can use central difference
+    if idx > 1 && idx < length(x)
+        % Central difference formula
+        h = x(idx+1) - x(idx-1);
+        df = (f(idx+1) - f(idx-1)) / h;
+    elseif idx == 1
+        % Forward difference formula
+        h = x(2) - x(1);
+        df = (f(2) - f(1)) / h;
+    else
+        % Backward difference formula
+        h = x(end) - x(end-1);
+        df = (f(end) - f(end-1)) / h;
+    end
+end

src/gen_params.m

@@ -15,9 +15,32 @@
 poro_perm_gen = @(num_samples) [poro_gen(num_samples); perm_gen(num_samples)];
 params.poro_perm_gen = poro_perm_gen;
 
-params.rel_perm.calc_krw = @(sw) calc_krw(sw);
-params.rel_perm.calc_krg = @(sg) calc_krg(sg);
-[~, params.rel_perm.sw_resid] = calc_krw(0);
+krw_data = [...
+     0.12	0		
+     0.15	0.05		
+     0.25	0.1		
+     0.35	0.3		
+     0.45	0.55		
+     0.65	0.85		 
+     0.8	0.97		
+     0.9	1		
+     1	    1
+     ];
+params.krw = TableFunction(krw_data);
+params.sw_resid = krw_data(1,1);
+
+krg_data = [...
+        0.0 0
+        0.1	0.01
+        0.2	0.02
+        0.35 0.04
+        0.55 0.2
+        0.65 0.4
+        0.75 0.7
+        0.85 1
+        0.88 1
+        ];
+params.krg = TableFunction(krg_data);
 
 params.capil.pres_func = @(sw,poro,perm) calc_capillary_pressure(sw,poro,perm,...
     params.capil.sw_barrier,...

src/upscale.m

@@ -49,8 +49,8 @@
 
     sw = invaded_mat_mid .* sw_mid + ~invaded_mat_mid .* 1;
 
-    kg_mat_local = params.rel_perm.calc_krg(1-sw);
-    kw_mat_local = params.rel_perm.calc_krw(sw);
+    kg_mat_local = params.krg.func(1-sw);
+    kw_mat_local = params.krw.func(sw);
 
     kg_mat_local = kg_mat_local.*permeabilities;
     kw_mat_local = kw_mat_local.*permeabilities;

strata_trapper.m

@@ -11,7 +11,7 @@
 
 perm_upscaled = zeros(grid.cells.num, 3);
 
-saturations = linspace(params.rel_perm.sw_resid,1,options.sat_num_points);
+saturations = linspace(params.sw_resid,1,options.sat_num_points);
 
 cap_pres_upscaled = zeros(grid.cells.num,length(saturations));
 krw = zeros(grid.cells.num,3,length(saturations));
@@ -49,7 +49,7 @@
     end
 end
 
-krw(:,:,saturations<=params.rel_perm.sw_resid) = 0;
+krw(:,:,saturations<=params.sw_resid) = 0;
 krg(krg<0) = 0;
 
 strata_trapped = struct(...