function [r, Times, OutputRs, Wellpara, trun] = run_case(config)
%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
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);

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);
    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