448 lines
14 KiB
Matlab
448 lines
14 KiB
Matlab
function [r, Times, OutputRs, Wellpara, trun] = run_case(config, progress_callback)
|
|
%RUN_CASE Execute a simulation from unified config.
|
|
%
|
|
% This is the first migration step from hard-coded main1.m scripts toward a
|
|
% unified GUI entry point. The current implementation is intended to
|
|
% reproduce case-template runs while the remaining flow-model internals are
|
|
% still being normalized.
|
|
|
|
arguments
|
|
config (1,1) struct
|
|
progress_callback = []
|
|
end
|
|
|
|
project_root = fileparts(fileparts(fileparts(mfilename('fullpath'))));
|
|
source_case_folder = resolve_source_case_folder(project_root, config.meta);
|
|
if isempty(source_case_folder)
|
|
error('run_case:MissingSourceCaseFolder', ...
|
|
'config.meta.source_case_folder is required.');
|
|
end
|
|
|
|
if ~isfolder(source_case_folder)
|
|
error('run_case:SourceCaseFolderNotFound', ...
|
|
'Source case folder not found: %s', source_case_folder);
|
|
end
|
|
|
|
addpath(genpath(project_root));
|
|
|
|
old_dir = pwd;
|
|
cleanup_obj = onCleanup(@() cd(old_dir));
|
|
cd(source_case_folder);
|
|
|
|
overall_tic = tic;
|
|
|
|
config.grid.nx = numel(config.grid.dx);
|
|
config.grid.ny = numel(config.grid.dy);
|
|
config.grid.nz = numel(config.grid.dz);
|
|
|
|
r = GridProp_pre( ...
|
|
config.grid.dx, ...
|
|
config.grid.dy, ...
|
|
config.grid.dz, ...
|
|
config.grid.nx, ...
|
|
config.grid.ny, ...
|
|
config.grid.nz, ...
|
|
config.grid.NTG);
|
|
|
|
[f, fellip, frac_information, nf] = build_fracture_inputs(config.fracture);
|
|
|
|
modelflag = config.model.modelflag;
|
|
grid_model = config.model.grid_model;
|
|
flow_model = config.model.flow_model;
|
|
|
|
if grid_model == 1
|
|
r = grid_discretization_SP_model(modelflag, r, f, frac_information, fellip, nf);
|
|
elseif grid_model == 2
|
|
r = grid_discretization_DP_model(modelflag, r, f, frac_information, fellip, nf);
|
|
else
|
|
error('run_case:UnsupportedGridModel', 'Unsupported grid model: %d', grid_model);
|
|
end
|
|
|
|
os = OperatorRS(r.N, r.nex, r.nc);
|
|
|
|
[fluid_model, state0] = build_flow_inputs(flow_model, r, config);
|
|
|
|
well1 = handle_well1(config.wells.well1, r);
|
|
[well2, welloc] = build_fracture_wells(config.wells, r);
|
|
Wellc = handle_well1_well2(well1, well2, welloc, r);
|
|
|
|
[Times, OutputRs, Wellpara, trun] = solver_NR( ...
|
|
r, ...
|
|
flow_model, ...
|
|
fluid_model, ...
|
|
os, ...
|
|
state0, ...
|
|
config.solver.yitap, ...
|
|
config.solver.yitas, ...
|
|
config.solver.omega, ...
|
|
config.solver.Nmax, ...
|
|
config.solver.epsave, ...
|
|
config.solver.epsmax, ...
|
|
config.schedule.dtmin, ...
|
|
config.schedule.dtmax, ...
|
|
Wellc, ...
|
|
config.schedule.time, ...
|
|
config.schedule.well_schedules, ...
|
|
progress_callback);
|
|
|
|
trun.run_time = toc(overall_tic);
|
|
end
|
|
|
|
function source_case_folder = resolve_source_case_folder(project_root, meta_config)
|
|
source_case_folder = meta_config.source_case_folder;
|
|
|
|
if isfolder(source_case_folder)
|
|
return;
|
|
end
|
|
|
|
case_id = lower(string(meta_config.case_id));
|
|
switch case_id
|
|
case "case01"
|
|
source_case_folder = find_case_folder_by_index(project_root, 1);
|
|
case "case02"
|
|
source_case_folder = find_case_folder_by_index(project_root, 2);
|
|
case "case03"
|
|
source_case_folder = find_case_folder_by_index(project_root, 3);
|
|
case "case04"
|
|
source_case_folder = find_case_folder_by_index(project_root, 4);
|
|
case "case05"
|
|
source_case_folder = find_case_folder_by_index(project_root, 5);
|
|
case "case06"
|
|
source_case_folder = find_case_folder_by_index(project_root, 6);
|
|
case "case07"
|
|
source_case_folder = find_case_folder_by_index(project_root, 7);
|
|
otherwise
|
|
if isempty(source_case_folder) || startsWith(string(source_case_folder), "CASE")
|
|
source_case_folder = '';
|
|
end
|
|
end
|
|
end
|
|
|
|
function source_case_folder = find_case_folder_by_index(project_root, case_index)
|
|
folder_candidates = dir(fullfile(project_root, sprintf('*%d-*', case_index)));
|
|
folder_candidates = folder_candidates([folder_candidates.isdir]);
|
|
|
|
if isempty(folder_candidates)
|
|
source_case_folder = '';
|
|
else
|
|
source_case_folder = fullfile(project_root, folder_candidates(1).name);
|
|
end
|
|
end
|
|
|
|
function [f, fellip, frac_information, nf] = build_fracture_inputs(fracture_config)
|
|
input_style = fracture_config.input_style;
|
|
|
|
switch input_style
|
|
case 1
|
|
[f, fellip] = sort_fracture(fracture_config.input_content);
|
|
m = size(f, 1);
|
|
n = size(fellip, 1);
|
|
nf = (m + n) / 5;
|
|
frac_information = fractureInformation_input_engineering_vector( ...
|
|
f, fellip, fracture_config.flowBarrierFlags);
|
|
case 2
|
|
f = fracture_config.f;
|
|
fellip = fracture_config.fellip;
|
|
m = size(f, 1);
|
|
n = size(fellip, 1);
|
|
nf = (m + n) / 5;
|
|
frac_information = fractureInformation_input_engineering_vector(f, fellip);
|
|
case 3
|
|
fellip = fracture_config.fellip;
|
|
[f, frac_information] = input_fracture_2D( ...
|
|
fracture_config.fractureLines, ...
|
|
fracture_config.fractureHeights, ...
|
|
fracture_config.flowBarrierFlags);
|
|
m = size(f, 1);
|
|
n = size(fellip, 1);
|
|
nf = (m + n) / 5;
|
|
otherwise
|
|
error('run_case:UnsupportedFractureInputStyle', ...
|
|
'Unsupported fracture input style: %d', input_style);
|
|
end
|
|
end
|
|
|
|
function [fluid_model, state0] = build_flow_inputs(flow_model, r, config)
|
|
switch flow_model
|
|
case 1
|
|
if has_complete_gas_water_config(config.flow.gas_water)
|
|
[fluid_model, state0] = build_gas_water_flow_from_config( ...
|
|
config.flow.gas_water, config.initial, r);
|
|
else
|
|
[fluid_model, state0] = gas_water_flow(r);
|
|
end
|
|
case 2
|
|
if has_complete_oil_water_config(config.flow.oil_water)
|
|
[fluid_model, state0] = build_oil_water_flow_from_config( ...
|
|
config.flow.oil_water, config.initial, r);
|
|
else
|
|
[fluid_model, state0] = oil_water_flow(r);
|
|
end
|
|
case 3
|
|
if has_complete_multicomponent_config(config.flow.multi_component)
|
|
[fluid_model, state0] = build_multi_component_flow_from_config( ...
|
|
config.flow.multi_component, config.initial, r);
|
|
else
|
|
[fluid_model, state0] = multi_component_flow(r);
|
|
end
|
|
otherwise
|
|
error('run_case:UnsupportedFlowModel', ...
|
|
'Unsupported flow model: %d', flow_model);
|
|
end
|
|
end
|
|
|
|
function tf = has_complete_gas_water_config(gw)
|
|
required_fields = { ...
|
|
'density_g_sc', 'density_w_sc', 'gas_model', 'prw', 'Bwi', 'cw', ...
|
|
'vwi', 'cvw', 'ifpcgl', 'matrix_relperm_table', 'fracture_relperm_table', ...
|
|
'p_grad_threshold'};
|
|
|
|
tf = true;
|
|
for i = 1:numel(required_fields)
|
|
field_name = required_fields{i};
|
|
if ~isfield(gw, field_name) || isempty(gw.(field_name))
|
|
tf = false;
|
|
return;
|
|
end
|
|
end
|
|
|
|
if ~isfield(gw, 'gas_prop') || isempty(gw.gas_prop)
|
|
tf = false;
|
|
return;
|
|
end
|
|
|
|
if gw.gas_model == 1
|
|
tf = all(isfield(gw, {'prg', 'Bgi', 'cg', 'vgi', 'cvg'})) && ...
|
|
~isempty(gw.prg) && ~isempty(gw.Bgi) && ~isempty(gw.cg) && ...
|
|
~isempty(gw.vgi) && ~isempty(gw.cvg);
|
|
else
|
|
tf = all(isfield(gw, {'Ppr', 'BG', 'MUG'})) && ...
|
|
~isempty(gw.Ppr) && ~isempty(gw.BG) && ~isempty(gw.MUG);
|
|
end
|
|
end
|
|
|
|
function tf = has_complete_oil_water_config(ow)
|
|
required_fields = { ...
|
|
'density_o_sc', 'density_w_sc', 'oil_model', 'prw', 'Bwi', 'cw', ...
|
|
'vwi', 'cvw', 'ifpcow', 'matrix_relperm_table', 'fracture_relperm_table', ...
|
|
'beta_non_darcy_flow', 'p_grad_threshold'};
|
|
|
|
tf = true;
|
|
for i = 1:numel(required_fields)
|
|
field_name = required_fields{i};
|
|
if ~isfield(ow, field_name) || isempty(ow.(field_name))
|
|
tf = false;
|
|
return;
|
|
end
|
|
end
|
|
|
|
if ow.oil_model == 1
|
|
tf = all(isfield(ow, {'pro', 'Boi', 'co', 'voi', 'cvo'})) && ...
|
|
~isempty(ow.pro) && ~isempty(ow.Boi) && ~isempty(ow.co) && ...
|
|
~isempty(ow.voi) && ~isempty(ow.cvo);
|
|
else
|
|
tf = all(isfield(ow, {'Ppr', 'BO', 'MUO'})) && ...
|
|
~isempty(ow.Ppr) && ~isempty(ow.BO) && ~isempty(ow.MUO);
|
|
end
|
|
end
|
|
|
|
function [f, state0] = build_gas_water_flow_from_config(gw, initial_config, r)
|
|
if gw.gas_model == 1
|
|
[Ppr, BG, MUG] = cal_gas_prop(gw.prg, gw.Bgi, gw.cg, gw.vgi, gw.cvg);
|
|
else
|
|
Ppr = gw.Ppr;
|
|
BG = gw.BG;
|
|
MUG = gw.MUG;
|
|
end
|
|
|
|
matrix_table = gw.matrix_relperm_table;
|
|
fracture_table = gw.fracture_relperm_table;
|
|
|
|
f = fluidPVT_gas_water_flow( ...
|
|
Ppr, BG, MUG, gw.Bwi, gw.prw, gw.cw, gw.vwi, gw.cvw, ...
|
|
matrix_table(:, 1), matrix_table(:, 3), matrix_table(:, 2), matrix_table(:, 4), ...
|
|
fracture_table(:, 1), fracture_table(:, 3), fracture_table(:, 2), fracture_table(:, 4), ...
|
|
gw.density_g_sc, gw.ifpcgl, r.rpt);
|
|
|
|
gas_prop = gw.gas_prop;
|
|
c1 = 1;
|
|
c2 = 8 / 9;
|
|
gas_prop.Kn_modified_factor = 1 + 8 * c1 * gas_prop.Kn + 16 * c2 * gas_prop.Kn^2;
|
|
|
|
f.Dwsi = gw.density_w_sc;
|
|
f.Dgsi = gw.density_g_sc;
|
|
f.gas_prop = gas_prop;
|
|
f.p_grad_threshold = gw.p_grad_threshold;
|
|
|
|
total_cells = r.nmc + r.nfc;
|
|
pressure0 = expand_initial_value(initial_config.pressure, total_cells);
|
|
sw0 = expand_initial_value(initial_config.sw, total_cells);
|
|
state0 = initialRS_gas_water_flow(pressure0, sw0);
|
|
end
|
|
|
|
function [f, state0] = build_oil_water_flow_from_config(ow, initial_config, r)
|
|
if ow.oil_model == 1
|
|
[Ppr, BO, MUO] = cal_oil_prop(ow.pro, ow.Boi, ow.co, ow.voi, ow.cvo);
|
|
else
|
|
Ppr = ow.Ppr;
|
|
BO = ow.BO;
|
|
MUO = ow.MUO;
|
|
end
|
|
|
|
matrix_table = ow.matrix_relperm_table;
|
|
fracture_table = ow.fracture_relperm_table;
|
|
|
|
f = fluidPVT_oil_water_flow( ...
|
|
Ppr, BO, MUO, ow.Bwi, ow.prw, ow.cw, ow.vwi, ow.cvw, ...
|
|
matrix_table(:, 1), matrix_table(:, 3), matrix_table(:, 2), matrix_table(:, 4), ...
|
|
fracture_table(:, 1), fracture_table(:, 3), fracture_table(:, 2), fracture_table(:, 4), ...
|
|
ow.density_o_sc, ow.ifpcow, r.rpt);
|
|
|
|
f.beta_non_Darcy_flow = ow.beta_non_darcy_flow;
|
|
f.Dwsi = ow.density_w_sc;
|
|
f.Dosi = ow.density_o_sc;
|
|
f.p_grad_threshold = ow.p_grad_threshold;
|
|
|
|
total_cells = r.nmc + r.nfc;
|
|
pressure0 = expand_initial_value(initial_config.pressure, total_cells);
|
|
sw0 = expand_initial_value(initial_config.sw, total_cells);
|
|
state0 = initialRS_oil_water_flow(pressure0, sw0);
|
|
end
|
|
|
|
function tf = has_complete_multicomponent_config(mc)
|
|
required_fields = { ...
|
|
'Ds', 'Db', 'c_ca_table', 'R', 'Vm', 'Temperature', ...
|
|
'cs_Nc', 'kr_nosurf', 'kr_surf', 'PC', ...
|
|
'cs_Nc_fracture', 'kr_nosurf_fracture', 'kr_surf_fracture', 'PC_fracture', ...
|
|
'density_g_sc', 'density_w_sc', 'gas_model', 'prw', 'Bwi', 'cw', 'vwi', 'cvw'};
|
|
|
|
tf = true;
|
|
for i = 1:numel(required_fields)
|
|
field_name = required_fields{i};
|
|
if ~isfield(mc, field_name) || isempty(mc.(field_name))
|
|
tf = false;
|
|
return;
|
|
end
|
|
end
|
|
|
|
if ~isfield(mc, 'gas_prop') || isempty(mc.gas_prop)
|
|
tf = false;
|
|
return;
|
|
end
|
|
|
|
if mc.gas_model == 1
|
|
tf = all(isfield(mc, {'prg', 'Bgi', 'cg', 'vgi', 'cvg'})) && ...
|
|
~isempty(mc.prg) && ~isempty(mc.Bgi) && ~isempty(mc.cg) && ...
|
|
~isempty(mc.vgi) && ~isempty(mc.cvg);
|
|
else
|
|
tf = all(isfield(mc, {'Ppr', 'BG', 'MUG'})) && ...
|
|
~isempty(mc.Ppr) && ~isempty(mc.BG) && ~isempty(mc.MUG);
|
|
end
|
|
end
|
|
|
|
function [f, state0] = build_multi_component_flow_from_config(mc, initial_config, r)
|
|
cs_data = mc.c_ca_table(:, 1);
|
|
cb_data = mc.c_ca_table(:, 1);
|
|
csa_data = mc.c_ca_table(:, 2) * 1e-3;
|
|
cba_data = mc.c_ca_table(:, 3) * 1e-3;
|
|
Nc_nosurf = min(mc.cs_Nc(:, 2));
|
|
Nc_surf = max(mc.cs_Nc(:, 2));
|
|
|
|
if mc.gas_model == 1
|
|
[Ppr, BG, MUG] = cal_gas_prop(mc.prg, mc.Bgi, mc.cg, mc.vgi, mc.cvg);
|
|
else
|
|
Ppr = mc.Ppr;
|
|
BG = mc.BG;
|
|
MUG = mc.MUG;
|
|
end
|
|
|
|
f = fluidPVT_new( ...
|
|
Ppr, BG, MUG, mc.Bwi, mc.prw, mc.cw, mc.vwi, mc.cvw, ...
|
|
mc.ifpcgl, r.rpt, cs_data, csa_data, cb_data, cba_data, ...
|
|
mc.cs_Nc, Nc_nosurf, Nc_surf, mc.kr_nosurf, mc.kr_surf, mc.PC, ...
|
|
mc.cs_Nc_fracture, mc.kr_nosurf_fracture, mc.kr_surf_fracture, mc.PC_fracture);
|
|
|
|
gas_prop = mc.gas_prop;
|
|
c1 = 1;
|
|
c2 = 8 / 9;
|
|
gas_prop.Kn_modified_factor = 1 + 8 * c1 * gas_prop.Kn + 16 * c2 * gas_prop.Kn^2;
|
|
if isfield(gas_prop, 'beta_non_darcy_flow') && ~isfield(gas_prop, 'beta_non_Darcy_flow')
|
|
gas_prop.beta_non_Darcy_flow = gas_prop.beta_non_darcy_flow;
|
|
elseif ~isfield(gas_prop, 'beta_non_Darcy_flow')
|
|
gas_prop.beta_non_Darcy_flow = 0;
|
|
end
|
|
|
|
f.Dwsi = mc.density_w_sc;
|
|
f.Dgsi = mc.density_g_sc;
|
|
f.gas_prop = gas_prop;
|
|
f.Ds = mc.Ds;
|
|
f.Db = mc.Db;
|
|
f.chemistry_cof = mc.R * mc.Temperature / mc.Vm * 1e-3 / 10;
|
|
f.p_grad_threshold = mc.p_grad_threshold;
|
|
f.number_state_variables = mc.number_state_variables;
|
|
|
|
total_cells = r.nmc + r.nfc;
|
|
matrix_x = default_if_empty(mc, 'x_matrix', 0.9);
|
|
fracture_x = default_if_empty(mc, 'x_fracture', 1.0);
|
|
f.x_total = [matrix_x * ones(r.nmc, 1); fracture_x * ones(r.nfc, 1)];
|
|
f.chemistry_potential = f.chemistry_cof * log(f.x_total);
|
|
|
|
pressure0 = expand_initial_value(initial_config.pressure, total_cells);
|
|
sw0 = expand_initial_value(initial_config.sw, total_cells);
|
|
cs0 = expand_initial_value(initial_config.cs, total_cells);
|
|
cb0 = expand_initial_value(sanitize_positive_initial_value(initial_config.cb, 50.0), total_cells);
|
|
state0 = initialRS(pressure0, sw0, cs0, cb0);
|
|
end
|
|
|
|
function value = expand_initial_value(raw_value, total_cells)
|
|
if isscalar(raw_value)
|
|
value = raw_value * ones(total_cells, 1);
|
|
else
|
|
value = raw_value;
|
|
end
|
|
end
|
|
|
|
function value = sanitize_positive_initial_value(raw_value, fallback)
|
|
if isempty(raw_value)
|
|
value = fallback;
|
|
return;
|
|
end
|
|
|
|
value = raw_value;
|
|
invalid_mask = ~isfinite(value) | value <= 0;
|
|
if all(invalid_mask(:))
|
|
value = fallback;
|
|
elseif any(invalid_mask(:))
|
|
value(invalid_mask) = fallback;
|
|
end
|
|
end
|
|
|
|
function value = default_if_empty(s, field_name, fallback)
|
|
if isfield(s, field_name) && ~isempty(s.(field_name))
|
|
value = s.(field_name);
|
|
else
|
|
value = fallback;
|
|
end
|
|
end
|
|
|
|
function [well2, welloc] = build_fracture_wells(wells_config, r)
|
|
num_fracture_wells = wells_config.num_fracture_wells;
|
|
welloc = wells_config.welloc;
|
|
perfnum = cell(num_fracture_wells, 1);
|
|
|
|
for i = 1:num_fracture_wells
|
|
perfnum{i, 1} = findWelloc(r, welloc{i, 1});
|
|
end
|
|
|
|
well2 = cell(num_fracture_wells, 6);
|
|
for i = 1:num_fracture_wells
|
|
if numel(wells_config.well2) >= i && ~isempty(wells_config.well2{i, 1})
|
|
well2(i, :) = wells_config.well2(i, :);
|
|
well2{i, 2} = length(perfnum{i, 1});
|
|
well2{i, 3} = perfnum{i, 1};
|
|
else
|
|
well2(i, :) = {sprintf('wf%d', i), length(perfnum{i, 1}), perfnum{i, 1}, 0.178/2, 0, 4};
|
|
end
|
|
end
|
|
end
|