add result button making plot

This commit is contained in:
xinxiao
2026-03-20 17:01:50 +08:00
parent 36f0592bff
commit 44cfe7c1e3
8 changed files with 1059 additions and 419 deletions
@@ -0,0 +1,46 @@
function actions = get_case_plot_actions(config)
%GET_CASE_PLOT_ACTIONS Return available legacy plot actions for a case.
case_id = lower(string(config.meta.case_id));
base_actions = struct( ...
'id', "plot_well_response", ...
'label', "Well Response", ...
'script', "plotWellResponse", ...
'description', "Plot the legacy well response figure.");
actions = base_actions;
switch case_id
case "case01"
actions(end + 1) = make_action("plot_2d_layer", "2D Layer", "plot_2D_layer", "Plot the default 2D layer pressure map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_3d_distribution", "3D Distribution", "plot_3D_dis", "Plot the default 3D distribution map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_perm", "Permeability", "plot_perm", "Plot the permeability layer map."); %#ok<AGROW>
case "case02"
actions(end + 1) = make_action("plot_2d_layer", "2D Layer", "plot_2D_layer", "Plot the default 2D layer pressure map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_3d_distribution", "3D Distribution", "plot_3D_dis", "Plot the default 3D distribution map."); %#ok<AGROW>
case "case03"
actions(end + 1) = make_action("plot_2d_layer", "2D Layer", "plot_2D_layer", "Plot the default 2D layer pressure map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_3d_distribution", "3D Distribution", "plot_3D_dis", "Plot the default 3D distribution map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_sp_dp", "SP/DP Map", "plot_SP_DP", "Plot the SP/DP shape-factor map."); %#ok<AGROW>
case "case04"
% well response only
case "case05"
% well response only
case "case06"
actions(end + 1) = make_action("plot_2d_layer", "2D Layer", "plot_2D_layer", "Plot the default 2D layer pressure map."); %#ok<AGROW>
actions(end + 1) = make_action("plot_3d_distribution", "3D Distribution", "plot_3D_dis", "Plot the default 3D distribution map."); %#ok<AGROW>
case "case07"
% well response only
otherwise
% no additional actions
end
end
function action = make_action(id, label, script, description)
action = struct( ...
'id', string(id), ...
'label', string(label), ...
'script', string(script), ...
'description', string(description));
end
@@ -0,0 +1,261 @@
function render_case_plot_in_axes(ax, config, results, action_id, options)
%RENDER_CASE_PLOT_IN_AXES Render a legacy plot directly into UIAxes.
arguments
ax
config (1,1) struct
results (1,1) struct
action_id
options struct = struct()
end
action_id = lower(string(action_id));
cla(ax);
switch action_id
case "plot_well_response"
render_well_response(ax, results, options);
case "plot_2d_layer"
render_2d_layer(ax, config, results, options);
case "plot_3d_distribution"
render_3d_distribution(ax, config, results, options);
case "plot_sp_dp"
render_sp_dp(ax, results, options);
case "plot_perm"
render_perm(ax, results, options);
otherwise
error('render_case_plot_in_axes:UnsupportedAction', ...
'Unsupported plot action: %s', string(action_id));
end
end
function render_well_response(ax, results, options)
times = results.Times(:);
wellpara = results.Wellpara;
well_index = get_option(options, 'well_index', 1);
metric = upper(string(get_option(options, 'metric', "BHP")));
set(ax, 'XScale', 'linear', 'YScale', 'linear');
if isempty(wellpara)
error('render_case_plot_in_axes:NoWellpara', 'Well response data is empty.');
end
n = numel(times);
data = zeros(n, 1);
switch metric
case "GPR"
field_name = 'qg';
y_label = 'Gas production rate, m^3/d';
for i = 1:n, data(i) = wellpara{i}{1, well_index}.(field_name); end
plot(ax, times, data, 'k^-', 'LineWidth', 1.2);
case "OPR"
field_name = 'qo';
y_label = 'Oil production rate, m^3/d';
for i = 1:n, data(i) = wellpara{i}{1, well_index}.(field_name); end
plot(ax, times, data, 'k^-', 'LineWidth', 1.2);
case "WPR"
field_name = 'qw';
y_label = 'Water production rate, m^3/d';
for i = 1:n, data(i) = wellpara{i}{1, well_index}.(field_name); end
plot(ax, times, data, 'k^-', 'LineWidth', 1.2);
case "DPWF"
y_label = 'Pressure derivative, MPa';
for i = 1:n
if i == 1
data(i) = (wellpara{i + 1}{1, well_index}.pwf - wellpara{i}{1, well_index}.pwf) / ...
(log(times(i + 1)) - log(times(i)));
elseif i == n
data(i) = (wellpara{i}{1, well_index}.pwf - wellpara{i - 1}{1, well_index}.pwf) / ...
(log(times(i)) - log(times(i - 1)));
else
data(i) = (wellpara{i + 1}{1, well_index}.pwf - wellpara{i - 1}{1, well_index}.pwf) / ...
(log(times(i + 1)) - log(times(i - 1)));
end
end
plot(ax, times, data, 'k^-', 'LineWidth', 1.2);
set(ax, 'XScale', 'log', 'YScale', 'log');
otherwise
field_name = 'pwf';
y_label = 'BHP, MPa';
for i = 1:n, data(i) = wellpara{i}{1, well_index}.(field_name); end
plot(ax, times, data, 'k^-', 'LineWidth', 1.2);
end
grid(ax, 'on');
xlabel(ax, 'Time, day');
ylabel(ax, y_label);
title(ax, sprintf('Well Response: well %d, %s', well_index, metric));
end
function render_2d_layer(ax, config, results, options)
r = results.r;
output = results.OutputRs{get_option(options, 'time_step', numel(results.OutputRs))};
layer_index = get_option(options, 'layer_index', 1);
property_id = string(get_option(options, 'property_id', "pressure"));
values = extract_plot_values(config, output, property_id);
coordinate = r.coordinates;
nodes = r.nodes;
nz = r.nz;
nx = r.nx;
ny = r.ny;
nel = 4;
dis = zeros(nel, nx * ny);
x = zeros(nel, nx * ny);
y = zeros(nel, nx * ny);
for i = 1:nx * ny
grid_numbering = i + (nz - layer_index) * nx * ny;
for j = 1:nel
x(j, i) = coordinate(nodes(grid_numbering, j), 1);
y(j, i) = coordinate(nodes(grid_numbering, j), 2);
if j == 3
x(j, i) = coordinate(nodes(grid_numbering, 4), 1);
y(j, i) = coordinate(nodes(grid_numbering, 4), 2);
elseif j == 4
x(j, i) = coordinate(nodes(grid_numbering, 3), 1);
y(j, i) = coordinate(nodes(grid_numbering, 3), 2);
end
dis(j, i) = values(grid_numbering);
end
end
fill(ax, x, y, dis, 'EdgeColor', 'interp');
axis(ax, 'equal');
axis(ax, 'tight');
colormap(ax, jet);
colorbar(ax);
xlabel(ax, 'x, m');
ylabel(ax, 'y, m');
title(ax, sprintf('2D Layer: layer %d, %s', layer_index, char(property_id)));
end
function render_3d_distribution(ax, config, results, options)
r = results.r;
output = results.OutputRs{get_option(options, 'time_step', numel(results.OutputRs))};
property_id = string(get_option(options, 'property_id', "pressure"));
v = extract_plot_values(config, output, property_id);
[downgrid, frontgrid, leftgrid, rightgrid, backgrid, upgrid] = preplot_dis(r);
coordinate = r.coordinates;
nodes = r.nodes;
hold(ax, 'on');
draw_surface(ax, downgrid, [1,2,4,3,1], coordinate, nodes, v);
draw_surface(ax, leftgrid, [1,5,7,3,1], coordinate, nodes, v);
draw_surface(ax, rightgrid, [2,4,8,6,2], coordinate, nodes, v);
draw_surface(ax, frontgrid, [1,2,6,5,1], coordinate, nodes, v);
draw_surface(ax, backgrid, [3,4,8,7,3], coordinate, nodes, v);
draw_surface(ax, upgrid, [5,6,8,7,5], coordinate, nodes, v);
hold(ax, 'off');
colormap(ax, jet);
colorbar(ax);
xlabel(ax, 'x, m');
ylabel(ax, 'y, m');
zlabel(ax, 'z, m');
title(ax, sprintf('3D Distribution: %s', char(property_id)));
view(ax, 3);
grid(ax, 'on');
axis(ax, 'tight');
end
function render_sp_dp(ax, results, options)
r = results.r;
layer_index = get_option(options, 'layer_index', 1);
render_static_layer_map(ax, r, r.sigma, layer_index, 'SP/DP Sigma');
end
function render_perm(ax, results, options)
r = results.r;
layer_index = get_option(options, 'layer_index', 1);
render_static_layer_map(ax, r, r.kx, layer_index, 'Permeability');
end
function render_static_layer_map(ax, r, values, layer_index, title_text)
coordinate = r.coordinates;
nodes = r.nodes;
nz = r.nz;
nx = r.nx;
ny = r.ny;
nel = 4;
dis = zeros(nel, nx * ny);
x = zeros(nel, nx * ny);
y = zeros(nel, nx * ny);
for i = 1:nx * ny
grid_numbering = i + (nz - layer_index) * nx * ny;
for j = 1:nel
x(j, i) = coordinate(nodes(grid_numbering, j), 1);
y(j, i) = coordinate(nodes(grid_numbering, j), 2);
if j == 3
x(j, i) = coordinate(nodes(grid_numbering, 4), 1);
y(j, i) = coordinate(nodes(grid_numbering, 4), 2);
elseif j == 4
x(j, i) = coordinate(nodes(grid_numbering, 3), 1);
y(j, i) = coordinate(nodes(grid_numbering, 3), 2);
end
dis(j, i) = values(grid_numbering);
end
end
fill(ax, x, y, dis, 'EdgeColor', 'interp');
axis(ax, 'equal');
axis(ax, 'tight');
colormap(ax, jet);
colorbar(ax);
xlabel(ax, 'x, m');
ylabel(ax, 'y, m');
title(ax, sprintf('%s: layer %d', title_text, layer_index));
end
function draw_surface(ax, grid_indices, order, coordinate, nodes, values)
nel = numel(order);
nmc = size(grid_indices, 2);
dis = zeros(nel, nmc);
x = zeros(nel, nmc);
y = zeros(nel, nmc);
z = zeros(nel, nmc);
z_shift = max(coordinate(:, 3)) + 2000;
for i = 1:nmc
for j = 1:nel
x(j, i) = coordinate(nodes(grid_indices(i), order(j)), 1);
y(j, i) = coordinate(nodes(grid_indices(i), order(j)), 2);
z(j, i) = coordinate(nodes(grid_indices(i), order(j)), 3) - z_shift;
dis(j, i) = values(grid_indices(i));
end
end
fill3(ax, x, y, z, dis, 'EdgeColor', 'interp');
end
function values = extract_plot_values(config, output, property_id)
property_id = lower(string(property_id));
switch property_id
case "pressure"
values = output.p;
case "water_saturation"
values = output.sw;
case "oil_saturation"
values = 1 - output.sw;
case "gas_saturation"
if isfield(output, 'sg')
values = output.sg;
else
values = 1 - output.sw;
end
otherwise
error('render_case_plot_in_axes:UnsupportedProperty', ...
'Unsupported plot property: %s', char(property_id));
end
end
function value = get_option(options, field_name, fallback)
if isfield(options, field_name) && ~isempty(options.(field_name))
value = options.(field_name);
else
value = fallback;
end
end
+104
View File
@@ -0,0 +1,104 @@
function run_case_plot(config, results, action_id)
%RUN_CASE_PLOT Run a legacy plot script for the current case/results.
arguments
config (1,1) struct
results (1,1) struct
action_id
end
if isempty(fieldnames(results))
error('run_case_plot:NoResults', 'No results are available for plotting.');
end
project_root = fileparts(fileparts(fileparts(mfilename('fullpath'))));
source_case_folder = resolve_source_case_folder(project_root, config.meta);
if isempty(source_case_folder) || ~isfolder(source_case_folder)
error('run_case_plot:MissingCaseFolder', ...
'Cannot resolve source case folder for case %s.', string(config.meta.case_id));
end
old_dir = pwd;
cleanup_obj = onCleanup(@() cd(old_dir)); %#ok<NASGU>
addpath(genpath(project_root));
cd(source_case_folder);
action_id = lower(string(action_id));
switch action_id
case "plot_well_response"
plotWellResponse(results.Times, results.Wellpara);
case "plot_2d_layer"
assert_plot_function_exists('plot_2D_layer');
plot_2D_layer(default_time_step(results), default_layer(results.r), results.r, results.OutputRs, default_plot_type(config));
case "plot_3d_distribution"
assert_plot_function_exists('plot_3D_dis');
plot_3D_dis(default_time_step(results), results.r, results.OutputRs, default_plot_type(config));
case "plot_sp_dp"
assert_plot_function_exists('plot_SP_DP');
plot_SP_DP(default_layer(results.r), results.r);
case "plot_perm"
assert_plot_function_exists('plot_perm');
plot_perm(default_layer(results.r), results.r);
otherwise
error('run_case_plot:UnsupportedAction', ...
'Unsupported plot action: %s', string(action_id));
end
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));
case_num = sscanf(char(case_id), 'case%d');
if isempty(case_num)
source_case_folder = '';
return;
end
folder_candidates = dir(fullfile(project_root, sprintf('*%d-*', case_num)));
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 idx = default_time_step(results)
if isfield(results, 'OutputRs') && ~isempty(results.OutputRs)
idx = numel(results.OutputRs);
else
idx = 1;
end
end
function k = default_layer(r)
k = 1;
if isfield(r, 'nz') && ~isempty(r.nz)
k = min(1, r.nz);
end
end
function type = default_plot_type(config)
switch config.model.flow_model
case 1
type = 1;
case 2
type = 1;
case 3
type = 1;
otherwise
type = 1;
end
end
function assert_plot_function_exists(function_name)
if exist(function_name, 'file') ~= 2
error('run_case_plot:MissingPlotFunction', ...
'Plot function not found in current case folder: %s', function_name);
end
end