Skip to content

Geothermal updates #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 14 commits into
base: main
Choose a base branch
from
Open

Geothermal updates #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
63db403
Remove empty file
strene Aug 5, 2025
c3c76ae
Better handling of cases where well has been added from trajectory
strene Aug 5, 2025
6fc5747
Ensure wellSpace is a row vector
strene Aug 5, 2025
30cd4b2
No functional changes
strene Aug 6, 2025
59be544
Support cases with missing fields (schedule, name)
strene Aug 6, 2025
1b0e780
Update wellbore mass flux eq for improved convergence
strene Aug 8, 2025
de0b09f
Add support for composite models
strene Sep 23, 2025
48f51ef
Expose more options through varargin + set dummy "number of cells" to…
strene Sep 23, 2025
0f7f960
Support cell input to plotToolbar
strene Sep 23, 2025
8ed643c
Add boundary face centroid depths to equilibriation
strene Sep 23, 2025
eedb56f
Update example to use new linear solver setup
strene Sep 24, 2025
6c409ab
Attempt at improved handling of very low wellbore mass rates
strene Sep 24, 2025
c496c62
Remove debug disp
strene Sep 24, 2025
6fd80f4
Slightly lower perm
strene Sep 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
117 changes: 74 additions & 43 deletions autodiff/test-suite/TestCase.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,8 +101,8 @@
test.toolbarOptions([1:2:end, 2:2:end]) = horzcat(names, values);
% Set plotting function if not given
if isempty(test.plotFn)
test.plotFn = @(G, v, varargin) ...
plotToolbar(G, v, test.toolbarOptions{:}, varargin{:});
test.plotFn = @(G, v, cells, varargin) ...
plotToolbar(G, v, cells, test.toolbarOptions{:}, varargin{:});
end
% Compute test case ID
if opt.computeID
Expand Down Expand Up @@ -278,12 +278,20 @@
require mrst-gui

if nargin == 1 || ischar(varargin{1})
v = test.model.rock;
if isa(test.model, 'CompositeModel')
assert(isfield(test.model.submodels, 'Reservoir'), ...
['Did not find submodel called ''Reservoir'', ', ...
'please provide values for plotting.']);
v = test.model.submodels.Reservoir.rock;
else
v = test.model.rock;
end
else
v = varargin{1};
varargin = varargin(2:end);
end
opt = struct('Name' , '' , ...
'cells' , [] , ...
'plotWells', true, ...
'wellOpts' , {{}}, ...
'camlight' , true);
Expand All @@ -294,7 +302,10 @@
end
h = test.figure('Name', Name);
G = test.getVisualizationGrid();
test.plotFn(G, v, varargin{:});
if isempty(opt.cells)
opt.cells = true(G.cells.num, 1);
end
test.plotFn(G, v, opt.cells, varargin{:});
if opt.plotWells
test.plotWells(opt.wellOpts{:});
end
Expand Down Expand Up @@ -322,47 +333,55 @@
function varargout = plotWells(test, varargin)
% Plot wells in the model

if ~isfield(test.schedule.control(1), 'W')
W = [];
if isa(test.model, 'CompositeModel')
if isfield(test.model.submodels, 'Wellbore')
W = test.model.submodels.Wellbore.wells;
else
warning(['Did not find submodel called ''Wellbore'', ', ...
'cannot plot wells.']);
end
elseif isfield(test.schedule.control(1), 'W')
W = test.schedule.control(1).W;
end
if isempty(W)
return;
end
W = test.schedule.control(1).W;
if ~isempty(W)
G = test.getVisualizationGrid();
if G.griddim == 3
ax = gca; dz = ax.ZLim(2) - ax.ZLim(1);
varargout = cell(nargout, 1);
[varargout{:}] = plotWell(G, W, 'color' , 'k' , ...
'height', 0.15*dz, ...
varargin{:} );
else
opt = struct('fontsize', 12 , ...
'color' , 'k', ...
'color2' , 'k');
[opt, varargin] = merge_options(opt, varargin{:});
nw = numel(W);
[hw, htext] = deal(zeros(nw,1));
gray = [1,1,1]*0.8;
hold on
for w = 1:numel(W)
color = opt.color;
if W(w).sign < 0, color = opt.color2; end
x = G.cells.centroids(vertcat(W(w).cells), :);
hw(w) = plot(x(1), x(2) , 'o' , ...
'Color' , color, ...
'markerFaceColor', gray , ...
'markerSize' , 8 , ...
varargin{:} );
if opt.fontsize > 0
htext(w) = text(x(1), x(2), W(w).name, ...
'FontSize' , opt.fontsize, ...
'VerticalAlignment', 'bottom' , ...
'interp' , 'latex' , ...
'Color' , color );
end
G = test.getVisualizationGrid();
if G.griddim == 3
ax = gca; dz = ax.ZLim(2) - ax.ZLim(1);
varargout = cell(nargout, 1);
[varargout{:}] = plotWell(G, W, 'color' , 'k' , ...
'height', 0.15*dz, ...
varargin{:} );
else
opt = struct('fontsize', 12 , ...
'color' , 'k', ...
'color2' , 'k');
[opt, varargin] = merge_options(opt, varargin{:});
nw = numel(W);
[hw, htext] = deal(zeros(nw,1));
gray = [1,1,1]*0.8;
hold on
for w = 1:numel(W)
color = opt.color;
if W(w).sign < 0, color = opt.color2; end
x = G.cells.centroids(vertcat(W(w).cells), :);
hw(w) = plot(x(1), x(2) , 'o' , ...
'Color' , color, ...
'markerFaceColor', gray , ...
'markerSize' , 8 , ...
varargin{:} );
if opt.fontsize > 0
htext(w) = text(x(1), x(2), W(w).name, ...
'FontSize' , opt.fontsize, ...
'VerticalAlignment', 'bottom' , ...
'interp' , 'latex' , ...
'Color' , color );
end
hold off
varargout = {hw, htext};
end
hold off
varargout = {hw, htext};
end

end
Expand All @@ -371,8 +390,20 @@
function ok = hasDrivingForce(test, force)
% Check if the test case has a given driving force

ok = isfield(test.schedule.control(1), force) ...
&& ~isempty(test.schedule.control(1).(force));
if isa(test.model, 'CompositeModel')
switch force
case 'W'
mdl = 'Wellbore';
case {'bc', 'src'}
mdl = 'Reservoir';
end
ok = isfield(test.model.submodels, mdl) ...
&& isfield(test.schedule.control(1).(mdl), force) ...
&& ~isempty(test.schedule.control(1).(mdl).(force));
else
ok = isfield(test.schedule.control(1), force) ...
&& ~isempty(test.schedule.control(1).(force));
end

end

Expand Down
1 change: 0 additions & 1 deletion autodiff/test-suite/setup-functions/layered_nls_test_wog.m
View file
Open in desktop

This file was deleted.

5 changes: 1 addition & 4 deletions modules/geothermal/examples/geothermalExampleWellbore.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,7 @@
test.schedule.control(2).Wellbore = ctrl;

%% Set up linear solver
lsol = setUpGeothermalWellboreLinearSolver(test.model);
test.model.G.cells.num = test.model.submodels.Reservoir.G.cells.num + ...
nnz(test.model.submodels.Wellbore.G.cells.type == 0 & ...
test.model.submodels.Wellbore.G.cells.hybrid == 0);
[lsol, test.model] = setUpGeothermalWellboreLinearSolver(test.model);

%% Simulate
problem = test.getPackedSimulationProblem('LinearSolver', lsol);
Expand Down
20 changes: 15 additions & 5 deletions modules/geothermal/models/wellbore/WellboreModel.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -489,14 +489,14 @@
function [eqs, names, types, state] = getMassFluxEquations(model, state)
% Equations realting segment mass fluxes to the pressure by means
% of a wellbore friction model


dp = model.getProp(state, 'FrictionLoss');
[p, s, rho] = model.parentModel.getProps(state, ...
'pressure' , ...
's' , ...
'Density' ...
);
v = model.getProp(state, 'massFlux');

nph = model.getNumberOfPhases();

Expand All @@ -511,10 +511,20 @@
g = norm(model.parentModel.gravity);
dz = model.parentModel.operators.Grad(model.G.cells.centroids(:,3));
dpw = model.parentModel.operators.Grad(p);
pot = dpw - rhoMix.*g.*dz;

% Assemble equation. We add the mass flux so that this reduces to a
% Darcy-type equation for very low friction loss to avoid singular
% systems.
% Convert mass flux to velocity
[Di, Do] = deal(0, model.G.faces.radius.*2);
if size(Do, 2) == 2
Di = Do(:,1);
Do = Do(:,2);
end
v = v./(pi*rhoMix.*((Do/2).^2 - (Di/2).^2));
eqs = (pot - dp) + 10*v;
Comment on lines +525 to +526
Copy link

Copilot AI Oct 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The magic number '10' in the equation should be defined as a named constant to improve code maintainability and make the scaling factor's purpose clearer.

Copilot uses AI. Check for mistakes.

pot = dpw - rhoMix.*g.*dz;

eqs = (pot - dp)./(1*atm);
eqs = {eqs};
names = {'flux'};
types = {'face'};
Expand Down Expand Up @@ -676,7 +686,7 @@
% re-initialize the values and change the controls so that the next
% step keeps within the prescribed ranges.


if ~isfield(state, 'wells'), state.wells = model.wells; end
if ~isfield(state, 'groups'), state.groups = model.groups; end
% Get driving forces and check that it's non-empty
Expand Down
16 changes: 13 additions & 3 deletions modules/geothermal/utils/addHybridCellsToWell.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,10 @@
%
% SYNOPSIS:
% W = addHybridCellsToWell(W, G, rock)
% setup = spe10_wo(W, G, rock, 'pn1', pv1, ...)
%
% DESCRIPTION:
% For each well in W, this function finds all hybrid cells adjacent ot
% the cells of the well and adds them to the well. Hybrid are grid faces
% For each well in W, this function finds all hybrid cells adjacent to the
% cells of the well and adds them to the well. Hybrid cells are grid faces
% that are assigned a volume in the computational grid. This function is
% useful if the wells are added using a function that does not recognize
% hybrid cells, e.g., addWellFromTrajectory.
Expand Down Expand Up @@ -79,6 +78,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
function W = addWellCells(W, G, rock, cells, aperture, varargin)
% Add cells

cells0 = W.cells;
nc0 = numel(W.cells);
W.cells = [W.cells; cells];
nc = numel(W.cells);
Expand All @@ -103,6 +103,16 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
W.cells = W.cells(ix);
W.WI = W.WI(ix);
W.dZ = W.dZ(ix);

if isfield(W, 'trajectory')
assert(all(W.trajectory.cell == cells0), ...
'W.trajectory.cell must match W.cells');
W.trajectory.cell = W.cells;
vec = [W.trajectory.vec; [zeros(numel(cells), 2), dz]];
W.trajectory.vec = vec(ix,:);
weight = [W.trajectory.weight; ones(numel(cells), 1)];
W.trajectory.weight = weight(ix);
end

% Add remaining fields by repeating W.(fn)(1)
for fn = fieldnames(W)'
Expand Down
64 changes: 35 additions & 29 deletions modules/geothermal/utils/convertToReservoirWellboreModel.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,15 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
% Optional input arguments
opt = struct('name', []);
[opt, varargin] = merge_options(opt, varargin{:});
if isempty(opt.name), opt.name = [setup.name, '_wb']; end
if isfield(setup, 'name') && isempty(opt.name), opt.name = [setup.name, '_wb']; end

% Get reservoir model and construct wellbore model
reservoirModel = setup.model;
W = setup.schedule.control(1).W;
if isfield(setup, 'W')
W = setup.W;
else
W = setup.schedule.control(1).W;
end
wellboreModel = WellboreModel(reservoirModel, W, varargin{:});
% Make composite model from the two
model = CompositeModel({reservoirModel, wellboreModel}, ...
Expand Down Expand Up @@ -70,36 +74,38 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
state0.Wellbore.s = state0.Wellbore.s(wc);
state0.Wellbore.T = state0.Wellbore.T(wc);

% Convert schedule
schedule = setup.schedule;
control = schedule.control;
for i = 1:numel(control)
[rctrl, wctrl] = deal(control(i));
[rctrl.W, rctrl.groups] = deal([]);
if ~isfield(rctrl, 'bc' ), rctrl.bc = []; end
if ~isfield(rctrl, 'src'), rctrl.src = []; end
[wctrl.bc, wctrl.src] = deal([]);
schedule.control(i).Reservoir = rctrl;
schedule.control(i).Wellbore = wctrl;
end
if isfield(schedule.control, 'W')
schedule.control = rmfield(schedule.control, 'W');
end
if isfield(schedule.control, 'src')
schedule.control = rmfield(schedule.control, 'src');
end
if isfield(schedule.control, 'bc')
schedule.control = rmfield(schedule.control, 'bc');
end
if isfield(schedule.control, 'groups')
schedule.control = rmfield(schedule.control, 'groups');
end

% Update setup
setup.state0 = state0;
setup.model = model;
setup.schedule = schedule;
setup.name = opt.name;
if ~isempty(opt.name), setup.name = opt.name; end

if isfield(setup, 'schedule') || isa(setup, 'TestCase')
% Convert schedule
schedule = setup.schedule;
control = schedule.control;
for i = 1:numel(control)
[rctrl, wctrl] = deal(control(i));
[rctrl.W, rctrl.groups] = deal([]);
if ~isfield(rctrl, 'bc' ), rctrl.bc = []; end
if ~isfield(rctrl, 'src'), rctrl.src = []; end
[wctrl.bc, wctrl.src] = deal([]);
schedule.control(i).Reservoir = rctrl;
schedule.control(i).Wellbore = wctrl;
end
if isfield(schedule.control, 'W')
schedule.control = rmfield(schedule.control, 'W');
end
if isfield(schedule.control, 'src')
schedule.control = rmfield(schedule.control, 'src');
end
if isfield(schedule.control, 'bc')
schedule.control = rmfield(schedule.control, 'bc');
end
if isfield(schedule.control, 'groups')
schedule.control = rmfield(schedule.control, 'groups');
end
setup.schedule = schedule;
end

% Set visualization grid
setup.visualizationGrid = model.submodels.Reservoir.G;
Expand Down
Loading