From bdc66c9a365b948aeb4849143c62320f1244f637 Mon Sep 17 00:00:00 2001
From: Maksim Elizarev <m.elizarev@imperial.ac.uk>
Date: Tue, 8 Oct 2024 14:04:28 +0100
Subject: [PATCH] WIP plot capillary pressure upscaling results in terms of
 Leverett J-function heterogeneity

---
 plot_result.m     | 63 +++++++++++++++++++----------------------------
 src/CapPressure.m | 14 +++++++++--
 2 files changed, 38 insertions(+), 39 deletions(-)

diff --git a/plot_result.m b/plot_result.m
index d89f504..5e218b3 100644
--- a/plot_result.m
+++ b/plot_result.m
@@ -1,51 +1,38 @@
 function plot_result(rock, mask, strata_trapped, params, kr_scale)
 arguments
-    rock 
-    mask 
-    strata_trapped 
-    params 
+    rock
+    mask
+    strata_trapped
+    params
     kr_scale = "log"
 end
 
-curves_plot(mask, strata_trapped, params, kr_scale);
-
-base_cap = strata_trapped;
-base_cap.permeability = rock.perm;
-for i = 1:length(base_cap.saturation)
-    base_cap.capillary_pressure(mask,i) = params.cap_pressure.func(base_cap.saturation(i),rock.poro(mask),rock.perm(mask,:));
-end
-
-tiles_pc = tiledlayout(figure(),'flow',TileSpacing='compact',Padding='tight');
-
-pc_base = nexttile(tiles_pc);
-y_lim_base = stat_plot(pc_base,'Base','',base_cap.saturation,[],base_cap.capillary_pressure(mask,:)./barsa());
-pc_base.YScale='log';
+leverett_j_upscaled = params.cap_pressure.inv_lj(strata_trapped.capillary_pressure(mask,:),rock.poro(mask),strata_trapped.permeability(mask,:));
 
-[y_lim_strat, pc_strat] = stat_plot(nexttile(tiles_pc),'StrataTrapped','',strata_trapped.saturation,[],strata_trapped.capillary_pressure(mask,:)./barsa());
-pc_strat.YScale='log';
+[y_lim_pc, ax_pc] =  stat_plot(axes(figure),'Leverett J-function','',strata_trapped.saturation,...
+    @(sw)params.cap_pressure.leverett_j.func(sw), leverett_j_upscaled,true);
 
-y_lim_concat = [y_lim_base;y_lim_strat];
-y_lim_both = [min(y_lim_concat(:,1)),max(y_lim_concat(:,2))];
+ax_pc.YScale='log';
+% ylim(ax_pc,y_lim_pc);
 
-ylim(pc_base,y_lim_both);
-ylim(pc_strat,y_lim_both);
+xlabel(ax_pc,'Wetting phase saturation');
+ylabel(ax_pc,'[-]');
 
-xlabel(tiles_pc,'Wetting phase saturation');
-ylabel(tiles_pc,'Capillary pressure, bar');
+curves_plot(mask, strata_trapped, params, kr_scale);
 
 end
 
 
 function curves_plot(mask, strata_trapped, params,scale)
 arguments
-    mask 
-    strata_trapped 
-    params 
+    mask
+    strata_trapped
+    params
     scale = "log"
 end
 sub_data = @(data,mask,direction) squeeze(data(mask,direction,:));
 
-tiles_krw = tiledlayout(figure(),3,2,TileSpacing='compact',Padding='tight');
+tiles_krw = tiledlayout(figure(),3,2,TileSpacing='tight',Padding='tight');
 
 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);
@@ -57,7 +44,7 @@ function curves_plot(mask, strata_trapped, params,scale)
 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,@(sw)params.krw.func(sw),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),true);
 yscale(scale);
 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);
@@ -67,7 +54,7 @@ function curves_plot(mask, strata_trapped, params,scale)
 end
 
 
-function [y_lim, ax] = stat_plot(ax, name, y_label, x_data, base_func, data,color)
+function [y_lim, ax] = stat_plot(ax, name, y_label, x_data, base_func, data,show_legend,color)
 arguments
     ax
     name char
@@ -75,6 +62,7 @@ function curves_plot(mask, strata_trapped, params,scale)
     x_data (1,:) double
     base_func
     data   (:,:) double
+    show_legend (1,1) logical = false
     color = 'blue'
 end
 
@@ -95,12 +83,14 @@ function curves_plot(mask, strata_trapped, params,scale)
 ylabel(ax,y_label);
 ax.YLimitMethod="tight";
 
-legends = {'Median','Quantiles 0.01 and 0.99',''};
-if ~isempty(base_func)
-    legends{end+1} = 'Intrinsic curve';
-end
+if show_legend
+    legends = {'Median','Quantiles 0.01 and 0.99',''};
+    if ~isempty(base_func)
+        legends{end+1} = 'Intrinsic curve';
+    end
 
-legend(ax,legends,'Location','best');
+    legend(ax,legends,'Location','best','BackgroundAlpha',0.5);
+end
 
 try
     [yu,yl,ym] = ax.Children(:).YData;
@@ -110,4 +100,3 @@ function curves_plot(mask, strata_trapped, params,scale)
     y_lim = [nan,nan];
 end
 end
-
diff --git a/src/CapPressure.m b/src/CapPressure.m
index e9773fc..831a41b 100644
--- a/src/CapPressure.m
+++ b/src/CapPressure.m
@@ -34,10 +34,10 @@
             pc = obj.mult * obj.leverett_j.func(sw) * sqrt(poro./perm);
         end
 
-        function sw = inv(obj,pc,poro,perm)
+         function lj = inv_lj(obj,pc,poro,perm)
             arguments
                 obj (1,1) CapPressure
-                pc  (1,1) double
+                pc  double
                 poro double {mustBeNonnegative}
                 perm double {mustBeNonnegative}
             end
@@ -45,7 +45,17 @@
             perm = obj.transform_perm(poro,perm);
 
             lj = pc ./ obj.mult .* sqrt(perm./poro);
+        end
+
+        function sw = inv(obj,pc,poro,perm)
+            arguments
+                obj (1,1) CapPressure
+                pc  (1,1) double
+                poro double {mustBeNonnegative}
+                perm double {mustBeNonnegative}
+            end
 
+            lj = obj.inv_lj(pc,poro,perm);
             lj = min(lj,obj.leverett_j.data(1,2));
             lj = max(lj,obj.leverett_j.data(end,2));