refactor: ♻️ make upscale MEX-compatible

Author: djmaxus

src/CapPressure.m

@@ -26,7 +26,7 @@
                 obj (1,1) CapPressure
                 sw  (1,1) double
                 poro double {mustBeNonnegative}
-                perm double {mustBeNonnegative}
+                perm (:,:,:,:) double {mustBeNonnegative}
             end
 
             perm = obj.transform_perm(poro,perm);
@@ -43,7 +43,7 @@
                 obj (1,1) CapPressure
                 pc  (1,1) double
                 poro double {mustBeNonnegative}
-                perm double {mustBeNonnegative}
+                perm (:,:,:,:) double {mustBeNonnegative}
             end
 
             lj = obj.inv_lj(pc,poro,perm);
@@ -58,7 +58,7 @@
                 obj (1,1) CapPressure
                 pc  double
                 poro double {mustBeNonnegative}
-                perm double {mustBeNonnegative}
+                perm (:,:,:,:) double {mustBeNonnegative}
             end
 
             perm = obj.transform_perm(poro,perm);
@@ -74,7 +74,7 @@
                 obj (1,1) CapPressure
                 sw  (1,1) double
                 poro double {mustBeNonnegative}
-                perm double {mustBePositive}
+                perm (:,:,:,:) double {mustBePositive}
             end
 
             perm = obj.transform_perm(poro,perm);
@@ -85,17 +85,24 @@
             dpc_dsw = obj.mult * obj.leverett_j.deriv(sw) * sqrt(poresize_mult);
         end
 
-        function perm = transform_perm(obj,poro,perm)
+        function perm_transformed = transform_perm(obj,poro,perm)
+        arguments
+            obj (1,1) CapPressure
+            poro double
+            perm (:,:,:,:) double 
+        end
             dims_poro = size(poro);
             dims_perm = size(perm);
             if (length(dims_poro) == 3) && (length(dims_perm) == 4) && (dims_perm(4) == 3)
-                perm = sum(perm .* reshape(obj.perm_weights, 1, 1, 1, 3), 4);
+                perm_transformed = sum(perm .* reshape(obj.perm_weights, 1, 1, 1, 3), 4);
                 return;
             end
 
             if (length(dims_poro) == 2) && (length(dims_perm) == 2) && (dims_perm(2) == 3)
-                perm = sum(perm .* reshape(obj.perm_weights, 1, 3), 2);
+                perm_transformed = sum(perm .* reshape(obj.perm_weights, 1, 3), 2);
+                return;
             end
+            perm_transformed = perm;
         end
     end
 end

src/Options.m

@@ -3,7 +3,7 @@
         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;
+        m_save_mip_step (1,1) {mustBeA(m_save_mip_step,'logical')} = false;
         m_perm_threshold_mD (1,3) double {mustBeNonnegative} = [0,0,0];
         m_poro_threshold    (1,1) double {mustBeNonnegative} = 0;
     end

src/calc_percolation.m

@@ -4,6 +4,7 @@
 h_ref = Lz*0.5;
 
 Nz = size(p_entry,3);
+h = zeros(1,1,Nz);
 h(1,1,1:Nz) = linspace(Lz/Nz/2,Lz - Lz/Nz/2,Nz);
 
 hydrostatic_correction = 0;
@@ -28,10 +29,10 @@
 function [invasion,has_changed] = calc_percolation_iter(invasion,Idx)
 [Nx,Ny,Nz] = size(invasion);
 has_changed = false;
-for idx = Idx
-    i = idx(1);
-    j = idx(2);
-    k = idx(3);
+for num_col = 1:size(Idx,2)
+    i = Idx(1,num_col);
+    j = Idx(2,num_col);
+    k = Idx(3,num_col);
 
     is_boundary = i==1 || i==Nx || j==1 || j==Ny || k==1 || k==Nz;
 

src/upscale.m

@@ -1,11 +1,19 @@
 function [perm_upscaled, pc_upscaled, krw, krg, mip] = upscale(...
     dr, saturations, params, options, porosities, permeabilities)
+arguments
+    dr (1,3) double
+    saturations (1, :) double
+    params (1,1) Params
+    options (1,1) Options
+    porosities (:,:,:) double
+    permeabilities (:,:,:,3) double
+end
 
 if max(porosities,[],'all') <= 0
     error('inactive cell');
 end
 
-mip(1:length(saturations)) = struct('sw',nan,'sub_sw',[]);
+mip = repmat(construct_mip_cell_data(),1,numel(saturations));
 
 % apply thresholds before proceeding
 small_poro = porosities < options.m_poro_threshold;
@@ -49,7 +57,10 @@
 
 entry_pressures = params.cap_pressure.func(1,porosities,permeabilities);
 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));
+is_pc_max_finite = isfinite(pc_max(:));
+pc_max_finite = pc_max(:);
+pc_max_finite = pc_max_finite(is_pc_max_finite);
+pc_points = linspace(max(pc_max_finite),min(entry_pressures(:)),length(saturations));
 
 for index_saturation = 1:length(saturations)
 
@@ -61,6 +72,8 @@
     calc_endpoint = index_saturation == 1 || index_saturation == length(saturations);
     max_iterations = calc_endpoint*1000 + ~calc_endpoint*100;
     err_prev = Inf;
+    pc_mid_tot = 0;
+    sub_sw = zeros(0,0,0); coder.varsize('sub_sw');
     for iteration_num=1:max_iterations
         [pc_mid_tot, sw_mid, pc_mid, sub_sw, converged, err] = mip_iteration(...
             sw_target, dr, entry_pressures, porosities, permeabilities, pc_mid, ...
@@ -109,12 +122,12 @@
 pc_upscaled = [pc_upscaled,pc_upscaled_extra];
 
 sw_upscaled(end+1) = sw_extra(1);
-krw(:,end+1) = 0;
-krg(:,end+1) = max(params.krg.data(:,2));
+krw = extend_cols_with(krw,0);
+krg = extend_cols_with(krg,max(params.krg.data(:,2)));
 
 sw_upscaled(end+1) = sw_extra(2);
-krw(:,end+1) = max(params.krw.data(:,2));
-krg(:,end+1) = 0;
+krw = extend_cols_with(krw,max(params.krw.data(:,2)));
+krg = extend_cols_with(krg,0);
 
 [sw_upscaled,unique_idx] = unique(sw_upscaled);
 pc_upscaled = pc_upscaled(unique_idx);
@@ -205,3 +218,20 @@
 krw = K_phase_upscaled(4:6)';
 
 end
+
+function kr = extend_cols_with(kr,value)
+krx = kr(1,:);
+kry = kr(2,:);
+krz = kr(3,:);
+
+krx(end+1) = value;
+kry(end+1) = value;
+krz(end+1) = value;
+kr = [krx;kry;krz];
+end
+
+function mip_cell_data = construct_mip_cell_data()
+sub_sw = zeros(0,0,0);
+coder.varsize('sub_sw');
+mip_cell_data = struct('sw',nan,'sub_sw',sub_sw);
+end

src/upscale_permeability.m

@@ -39,6 +39,7 @@
 
 idx_conductive = find(cond_sum>0);
 
+idx_adj = zeros(numel(idx_conductive),3,2);
 idx_adj(:,1,1) = idx_conductive - 1;
 idx_adj(:,1,2) = idx_conductive + 1;
 idx_adj(:,2,1) = idx_conductive - x;
@@ -53,11 +54,11 @@
 A_diag(:,2) = 1:num_eqs;
 
 num_cells_ext = x*y*z;
-
+map_to_cond = zeros(num_cells_ext,1);
 map_to_cond(1:num_cells_ext,1)=0;
 map_to_cond(idx_conductive) = 1:length(idx_conductive);
 
-A_ndiag{3,2} = [];
+A_ndiag = cell(3,2);
 
 for dir=dirs
     cond_dir = cond(:,dir);
@@ -75,7 +76,7 @@
 end
 
 A_ndiag = cell2mat(reshape(A_ndiag,[],1));
-
+lin_to_sub = zeros(num_cells_ext,3);
 [lin_to_sub(:,1),lin_to_sub(:,2),lin_to_sub(:,3)] = ind2sub([x,y,z],1:num_cells_ext);
 is_inner = false(num_cells_ext,3);