feat: ✨ update structure and values of input data

src/capil/calc_capillary_pressure.m

@@ -0,0 +1,23 @@
+function [pc] = calc_capillary_pressure(sw,poro,perm,sw_barrier,lambda,contact_angle,surface_tension)
+arguments
+    sw (1,1) double
+    poro double
+    perm double
+    sw_barrier (1,1) double
+    lambda (1,1) double
+    contact_angle (1,1) double
+    surface_tension (1,1) double
+end
+if sw <= sw_barrier
+    pc = Inf;
+    return;
+end
+sw = min(sw,1);
+
+    function leverett_j = calc_leverett_j(sw,sw_barrier,lambda)
+        sw_normalized = (sw-sw_barrier)/(1-sw_barrier);
+        leverett_j = sw_normalized^(-1/lambda);
+    end
+
+pc = cos(contact_angle) * surface_tension * calc_leverett_j(sw,sw_barrier,lambda) * sqrt(poro./perm);
+end

src/capil/calc_capillary_pressure_deriv.m

@@ -0,0 +1,24 @@
+function [dpc_dsw] = calc_capillary_pressure_deriv(sw,poro,perm,sw_barrier,lambda,contact_angle,surface_tension)
+arguments
+    sw (1,1) double
+    poro double
+    perm double
+    sw_barrier (1,1) double
+    lambda (1,1) double
+    contact_angle (1,1) double
+    surface_tension (1,1) double
+end
+if sw <= sw_barrier
+    dpc_dsw = Inf;
+    return;
+end
+sw = min(sw,1);
+
+    function dj_dsw = calc_leverett_j_deriv(sw,sw_barrier,lambda)
+        pow = (-1/lambda);
+        mult = (1-sw_barrier)^(-pow);
+        dj_dsw = mult * pow * (sw-sw_barrier)^(pow-1);
+    end
+
+dpc_dsw = cos(contact_angle) * surface_tension * calc_leverett_j_deriv(sw,sw_barrier,lambda) * sqrt(poro./perm);
+end

src/capil/calc_capillary_pressure_inv.m

@@ -0,0 +1,21 @@
+function sw = calc_capillary_pressure_inv(pc,poro,perm,sw_barrier,lambda,contact_angle,surface_tension)
+arguments
+    pc (1,1) double
+    poro double
+    perm double
+    sw_barrier (1,1) double
+    lambda (1,1) double
+    contact_angle (1,1) double
+    surface_tension (1,1) double
+end
+
+leverett_j = pc ./ (cos(contact_angle) * surface_tension * sqrt(poro./perm));
+
+    function sw = calc_leverett_j_inv(leverett_j,sw_barrier,lambda)
+        sw_normalized = leverett_j.^(-lambda);
+        sw = sw_normalized .* (1-sw_barrier) + sw_barrier;
+    end
+
+sw = calc_leverett_j_inv(leverett_j,sw_barrier,lambda);
+
+end

src/downscale.m

@@ -21,9 +21,10 @@
 
 %% calculate fine-scale permeability
 
-sub_permeability = rsgen3D(dr,subscale_dims,params.corr_lens,@(N)sub_permeability);
-sub_permeability = normalize(sub_permeability,log(perm_coarse(1)));
-sub_permeability = exp(sub_permeability);
+sub_permeability_log = log(sub_permeability);
+sub_permeability_log = rsgen3D(dr,subscale_dims,params.corr_lens,@(N)sub_permeability_log);
+sub_permeability_log = normalize(sub_permeability_log,log(perm_coarse));
+sub_permeability = exp(sub_permeability_log);
 
 %% calculate fine-scale sub_entry_pressures
 

src/gen_options.m

@@ -1,9 +1,8 @@
 function [options] = gen_options()
 options = struct(...
-    'subscale',       20*meter(), ...
-    'sat_num_points', 16, ...
-    'sat_min',  0.01, ...
+    'subscale',       50*meter(), ...
+    'sat_num_points', 10, ...
     'sat_tol', 1e-8, ...
-    'perm_min_threshold', 1e-8 * milli() * darcy() ...
+    'perm_min_threshold', 1e-7 * milli() * darcy() ...
     );
 end

src/gen_params.m

@@ -1,29 +1,39 @@
 function [params] = gen_params()
 params = struct(...
-    'poro',      struct(...
-        'mean',     0.23,...
-        'std_dev',  0.049,...
-        'corr_lens',[200,20,0.1].*meter() ...
-    ),...
-    'perm_corr', @(poro) 10.^(log10(poro).*4.19+5.17) * milli() * darcy(),... NOTE: produces milli Darcy
-    'capil',     struct(...
-        'lambda',         1.75, ...
-        'angle',          deg2rad(37), ...
-        'tension',        36 * milli() * newton() / meter(), ...
-        'pres_entry',     10 * kilo() * pascal(), ...
-        'sw_barrier',     0 ...
-    ),...
-    'rel_perm', struct('sw_resid', 0.01)...
-    );
+    'corr_lens', [50,50,0].*meter(), ...
+    'capil', struct(...
+        'lambda',          0.75, ...
+        'contact_angle',   deg2rad(50), ...
+        'surface_tension', 0.03098 * newton() / meter(), ...
+        'sw_barrier',      0.39 ...
+    )...
+);
 
-params.capil.pres_func = @(sw, entry_pressure) brooks_corey(sw, params.capil.sw_barrier, entry_pressure, params.capil.lambda);
+poro_gen = @(num_samples) randn(1,num_samples)*0.01 + 0.2;
+perm_gen = @(num_samples) exp(randn(1,num_samples) + log(100)) * milli * darcy;
 
-params.capil.pres_func_inv = @(pc, entry_pressure) brooks_corey_inv(pc, params.capil.sw_barrier, entry_pressure, params.capil.lambda);
+poro_perm_gen = @(num_samples) [poro_gen(num_samples); perm_gen(num_samples)];
+params.poro_perm_gen = poro_perm_gen;
 
-params.capil.pres_deriv = @(sw, entry_pressure) brooks_corey_deriv(sw, params.capil.sw_barrier, entry_pressure,params.capil.lambda);
+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);
 
-params.capil.j_entry = params.capil.pres_entry / (sqrt(params.poro.mean/params.perm_corr(params.poro.mean)) .* cos(params.capil.angle) .* params.capil.tension);
+params.capil.pres_func = @(sw,poro,perm) calc_capillary_pressure(sw,poro,perm,...
+    params.capil.sw_barrier,...
+    params.capil.lambda,...
+    params.capil.contact_angle,...
+    params.capil.surface_tension);
 
-params.rel_perm.wat_func = @(sw) max((sw - params.rel_perm.sw_resid)./(1 - params.rel_perm.sw_resid),0);
-params.rel_perm.gas_func = @(sg) min(sg./(1 - params.rel_perm.sw_resid),1);
+params.capil.pres_func_inv = @(pc,poro,perm) calc_capillary_pressure_inv(pc,poro,perm,...
+    params.capil.sw_barrier,...
+    params.capil.lambda,...
+    params.capil.contact_angle,...
+    params.capil.surface_tension);
+
+params.capil.pres_deriv = @(sw,poro,perm) calc_capillary_pressure_deriv(sw,poro,perm,...
+    params.capil.sw_barrier,...
+    params.capil.lambda,...
+    params.capil.contact_angle,...
+    params.capil.surface_tension);
 end

src/rel_perm/calc_krg.m

@@ -0,0 +1,31 @@
+function [krg] = calc_krg(sg)
+persistent data;
+
+if isempty(data)
+    data = [...
+        0.000 0.0000;
+        0.025 0.0001;
+        0.050 0.0002;
+        0.075 0.001;
+        0.100 0.002;
+        0.125 0.004;
+        0.150 0.006;
+        0.175 0.010;
+        0.200 0.015;
+        0.2250 0.020;
+        0.2500 0.025;
+        0.2750 0.034;
+        0.2800 0.035;
+        0.2950 0.041;
+        0.3000 0.042;
+        0.3500 0.065;
+        0.4000 0.10;
+        0.4500 0.13;
+        0.5000 0.17;
+        0.5500 0.23;
+        0.5800 0.26;
+        ];
+end
+
+krg = interp1(data(:,1), data(:,2), sg,"linear","extrap");
+end

src/rel_perm/calc_krw.m

@@ -0,0 +1,28 @@
+function [krw, sw_irr] = calc_krw(sw)
+persistent data sw_irr_
+
+if isempty(data)
+    data = [...
+        0.42 0.00;
+        0.45 0.01;
+        0.50 0.03;
+        0.55 0.08;
+        0.60 0.13;
+        0.65 0.20;
+        0.70 0.28;
+        0.71 0.29;
+        0.75 0.38;
+        0.80 0.50;
+        0.85 0.60;
+        0.90 0.70;
+        0.95 0.85;
+        1.00 1.00;
+        ];
+    sw_irr_ = interp1(data(:,2), data(:,1), 0);
+end
+
+krw = interp1(data(:,1), data(:,2), sw, "linear");
+sw_irr = sw_irr_;
+
+end
+