diff --git a/demo.m b/demo.m index 5b494b3..edf0df0 100644 --- a/demo.m +++ b/demo.m @@ -25,8 +25,7 @@ function demo(args) %% Downscaling demo - -fig_downscale = downscale_demo(params, downscale_params,visible); +downscale_demo(params, downscale_params,visible); %% Grid & rock properties @@ -34,19 +33,19 @@ function demo(args) %% Run StrataTrapper -mask = true(ceil(grid.cells.num* 0.001),1); % process only a fraction of cells +mask = true(ceil(grid.cells.num* 1e-3),1); % process only a fraction of cells sub_rock = downscale_all(grid,rock,mask,downscale_params); strata_trapped = strata_trapper(grid, sub_rock, params, ... mask=mask, options=Options(), ... enable_waitbar=args.show_progress,... - parfor_arg=args.parfor_arg ... + parfor_arg=args.parfor_arg... ); %% Visualize saturation functions -fig = plot_result(strata_trapped,visible=visible); +plot_result(strata_trapped,visible=visible); %% OGS inputs generation diff --git a/plot_result.m b/plot_result.m index 33be5a9..b55c8be 100644 --- a/plot_result.m +++ b/plot_result.m @@ -14,19 +14,21 @@ clf(fig); end -[t_all,t_kr,t_krw,t_krg,ax_pc,ax_krw_x,ax_krw_y,ax_krw_z,ax_krg_x,ax_krg_y,ax_krg_z] = nested_tiles(fig); +[t_all,t_kr,t_krw,t_krg,ax_pc,ax_krw_x,ax_krw_y,ax_krw_z,ax_krg_x,ax_krg_y,ax_krg_z] ... + = nested_tiles(fig); leverett_j_upscaled = strata_trapped.params.cap_pressure.inv_lj(... strata_trapped.capillary_pressure,... strata_trapped.porosity,... strata_trapped.permeability); -[~, ax_pc] = stat_plot(ax_pc,'Leverett J-function','',strata_trapped.saturation,... +[~, ax_pc] = stat_plot(ax_pc,strata_trapped.saturation,... @(sw)strata_trapped.params.cap_pressure.leverett_j.func(sw), leverett_j_upscaled,true); title(ax_pc,'Leverett J-function'); ylabel(ax_pc,'[-]'); ax_pc.YScale='log'; -curves_plot([ax_krw_x,ax_krw_y,ax_krw_z;ax_krg_x,ax_krg_y,ax_krg_z], strata_trapped, strata_trapped.params, args.kr_scale); +curves_plot([ax_krw_x,ax_krw_y,ax_krw_z;ax_krg_x,ax_krg_y,ax_krg_z], ... + strata_trapped, strata_trapped.params, args.kr_scale); xlabel(t_all,'Wetting phase saturation',FontSize=args.font_size); title(t_kr,'Relative permeability',FontSize=args.font_size); @@ -51,28 +53,32 @@ function curves_plot(ax_kr, strata_trapped, params,scale) end sub_data = @(data,direction) squeeze(data(:,direction,:)); -stat_plot(ax_kr(1,1),'','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,1)); +stat_plot(ax_kr(1,1),strata_trapped.saturation,... + @(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,1)); ax_kr(1,1).YScale = scale; -stat_plot(ax_kr(2,1),'','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,1)); +stat_plot(ax_kr(2,1),strata_trapped.saturation,... + @(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,1)); ax_kr(2,1).YScale = scale; -stat_plot(ax_kr(1,2),'','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,2)); +stat_plot(ax_kr(1,2),strata_trapped.saturation,... + @(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,2)); ax_kr(1,2).YScale = scale; -stat_plot(ax_kr(2,2),'','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,2)); +stat_plot(ax_kr(2,2),strata_trapped.saturation,... + @(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,2)); ax_kr(2,2).YScale = scale; -stat_plot(ax_kr(1,3),'','',strata_trapped.saturation,@(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,3)); +stat_plot(ax_kr(1,3),strata_trapped.saturation,... + @(sw)params.krw.func(sw),sub_data(strata_trapped.rel_perm_wat,3)); ax_kr(1,3).YScale = scale; -stat_plot(ax_kr(2,3),'','',strata_trapped.saturation,@(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,3)); +stat_plot(ax_kr(2,3),strata_trapped.saturation, ... + @(sw) params.krg.func(1-sw),sub_data(strata_trapped.rel_perm_gas,3)); ax_kr(2,3).YScale = scale; end -function [y_lim, ax] = stat_plot(ax, name, y_label, x_data, base_func, data,show_legend,color) +function [y_lim, ax] = stat_plot(ax, x_data, base_func, data,show_legend,color) arguments ax - name char - y_label char x_data (1,:) double base_func data (:,:) double @@ -114,7 +120,8 @@ function curves_plot(ax_kr, strata_trapped, params,scale) end end -function [t_all,t_kr,t_krw,t_krg,ax_pc,ax_krw_x,ax_krw_y,ax_krw_z,ax_krg_x,ax_krg_y,ax_krg_z] = nested_tiles(fig) +function [t_all,t_kr,t_krw,t_krg,ax_pc,ax_krw_x,ax_krw_y,ax_krw_z,ax_krg_x,ax_krg_y,ax_krg_z] ... + = nested_tiles(fig) params = {'TileSpacing','tight','Padding','tight'}; t_all = tiledlayout(fig,1,3,params{:}); diff --git a/src/downscale.m b/src/downscale.m index 2d2f882..f057faa 100644 --- a/src/downscale.m +++ b/src/downscale.m @@ -21,7 +21,8 @@ %% calculate fine-scale permeability sub_permeability_log = log(sub_perm_x); -sub_permeability_log = gen_corr_field(dr,subscale_dims,downscale_params.corr_lens,@(N)sub_permeability_log); +sub_permeability_log = gen_corr_field(dr,subscale_dims,downscale_params.corr_lens, ... + @(N)sub_permeability_log); perm_x = perm_coarse(1); sub_permeability_log = normalize(sub_permeability_log,log(perm_x)); sub_perm_x = exp(sub_permeability_log); diff --git a/src/upscale.m b/src/upscale.m index 5076977..8e78bbd 100644 --- a/src/upscale.m +++ b/src/upscale.m @@ -63,7 +63,8 @@ Kg_sub_mD = Kg_sub_mD.*(porosities~=0).*permeabilities_mD; Kw_sub_mD = Kw_sub_mD.*(porosities~=0).*permeabilities_mD; - [krg(:,index_saturation), krw(:,index_saturation)] = calc_relative_permeabilities(dr_sub, perm_upscaled_mD, Kg_sub_mD, Kw_sub_mD); + [krg(:,index_saturation), krw(:,index_saturation)] = calc_relative_permeabilities( ... + dr_sub, perm_upscaled_mD, Kg_sub_mD, Kw_sub_mD); end % NOTE: we expect only small negative values as computational errors @@ -75,7 +76,8 @@ sw_extra = [params.sw_resid,1]; [~,unique_idx] = unique(sw_upscaled); -pc_upscaled_extra = exp(interp1(sw_upscaled(unique_idx), log(pc_upscaled(unique_idx)), sw_extra, "linear","extrap")); +pc_upscaled_extra = exp(interp1( ... + sw_upscaled(unique_idx), log(pc_upscaled(unique_idx)), sw_extra, "linear","extrap")); pc_upscaled = [pc_upscaled,pc_upscaled_extra]; sw_upscaled(end+1) = sw_extra(1); @@ -121,7 +123,8 @@ pore_volumes = porosities .* sub_volume; pore_volume = sum(pore_volumes,'all'); -sub_sw = invaded_mat_mid .* params.cap_pressure.inv(pc_mid,porosities,permeabilities) + ~invaded_mat_mid .* 1; +sub_sw = invaded_mat_mid .* params.cap_pressure.inv(pc_mid,porosities,permeabilities) ... + + ~invaded_mat_mid .* 1; sub_sw(~isfinite(sub_sw)) = 1; sw_mid = sum(sub_sw.*pore_volumes,'all')/pore_volume;