Skip to content

Commit 900ce30

Browse files
moynermoyner
committed
Update subtree
2 parents 8bb0961 + 9ee0698 commit 900ce30

File tree

5 files changed

+103
-77
lines changed

5 files changed

+103
-77
lines changed

‎modules/geothermal/models/wellbore/WellboreModel.m‎

Lines changed: 91 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -469,6 +469,7 @@
469469
%-----------------------------------------------------------------%
470470
function [eqs, names, types, state] = getModelEquations(model, state0, state, dt, drivingForces)
471471

472+
state = model.applyLimits(state, drivingForces);
472473
% Get parent model equations
473474
[eqs, names, types, state] ...
474475
= model.parentModel.getModelEquations(state0, state, dt, drivingForces);
@@ -526,7 +527,7 @@
526527
% Equations imposing well control for each well
527528

528529
% Get control types, targets, and mask for open wells
529-
[type, target, compi, targetTemp, is_open] = model.getControls(drivingForces);
530+
[type, target, compi, targetTemp, is_open] = model.getControls(state, drivingForces);
530531

531532
% Convenient well and group indices
532533
nw = model.numWells();
@@ -633,14 +634,14 @@
633634
%-----------------------------------------------------------------%
634635

635636
%-----------------------------------------------------------------%
636-
function [type, target, compi, targetTemp, is_open] = getControls(model, drivingForces)
637+
function [type, target, compi, targetTemp, is_open] = getControls(model, state, drivingForces)
637638

638639
% Get driving forces and check that it's non-empty
639640
[wCtrl, gCtrl] = deal([]);
640641
hasWellCtrl = isfield(drivingForces, 'W') && ~isempty(drivingForces.W);
641642
hasGroupCtrl = isfield(drivingForces, 'groups') && ~isempty(drivingForces.groups);
642-
if hasWellCtrl , wCtrl = model.wells; end
643-
if hasGroupCtrl, gCtrl = model.groups; end
643+
if hasWellCtrl , wCtrl = state.wells; end
644+
if hasGroupCtrl, gCtrl = state.groups; end
644645
if isempty(wCtrl) && isempty(gCtrl), return; end
645646

646647
nw = model.numWells();
@@ -650,7 +651,7 @@
650651
compi = ones(1, nph)./nph;
651652
ctrl0 = struct('type', 'none', 'val' , nan, ...
652653
'compi', compi, 'T', nan, 'status', true);
653-
if ~hasWellCtrl , wCtrl = repmat(ctrl0, nw, 1) ; end
654+
if ~hasWellCtrl , wCtrl = repmat(ctrl0, nw, 1); end
654655
if ~hasGroupCtrl, gCtrl = repmat(ctrl0, ng, 1); end
655656

656657
target = [vertcat(wCtrl.val); vertcat(gCtrl.val)];
@@ -669,81 +670,94 @@
669670
%-----------------------------------------------------------------%
670671

671672
%-----------------------------------------------------------------%
672-
function model = applyLimits(model, state)
673+
function state = applyLimits(model, state, drivingForces)
673674
% Update solution variables and wellSol based on the well limits.
674675
% If limits have been reached, this function will attempt to
675676
% re-initialize the values and change the controls so that the next
676677
% step keeps within the prescribed ranges.
677678

678-
withinLimits = true;
679-
is_open = model.getActiveWells();
680-
% Check if we can change controls, and return if not
681-
if ~any(model.allowCtrlSwitching), return; end
682-
% Check if any wells are open, and return if not
683-
if ~any(is_open), return; end
679+
680+
if ~isfield(state, 'wells'), state.wells = model.wells; end
681+
if ~isfield(state, 'groups'), state.groups = model.groups; end
682+
% Get driving forces and check that it's non-empty
683+
[wCtrl, gCtrl] = deal([]);
684+
hasWellCtrl = isfield(drivingForces, 'W') && ~isempty(drivingForces.W);
685+
hasGroupCtrl = isfield(drivingForces, 'groups') && ~isempty(drivingForces.groups);
686+
if hasWellCtrl , wCtrl = state.wells; end
687+
if hasGroupCtrl, gCtrl = state.groups; end
688+
if isempty(wCtrl) && isempty(gCtrl), return; end
689+
690+
nw = model.numWells();
691+
ng = model.numGroups();
692+
693+
nph = model.getNumberOfPhases;
694+
compi = ones(1, nph)./nph;
695+
ctrl0 = struct('type', 'none', 'val' , nan, ...
696+
'compi', compi, 'T', nan, 'status', true);
697+
if ~hasWellCtrl , wCtrl = repmat(ctrl0, nw, 1) ; end
698+
if ~hasGroupCtrl, gCtrl = repmat(ctrl0, ng, 1); end
699+
700+
target = [vertcat(wCtrl.val); vertcat(gCtrl.val)];
701+
type = vertcat({wCtrl.type}', {gCtrl.type}');
702+
compi = [vertcat(wCtrl.compi); vertcat(gCtrl.compi)];
703+
704+
targetTemp = [];
705+
if model.thermal
706+
targetTemp = [vertcat(wCtrl.T); vertcat(gCtrl.T)];
707+
end
708+
709+
% Wells that are open to flow
710+
is_open = vertcat(wCtrl.status, gCtrl.status);
684711

685712
limits = model.getLimits();
686713

687-
[bhp, qs] = model.getProps(state, 'BottomholePressure', 'SurfaceRate');
688-
qsTot = sum(cell2mat(qs), 2);
689-
is_inj = qsTot > 0;
690-
691-
692-
% Injectors have three possible limits:
693-
% bhp: Upper limit on pressure.
694-
% rate: Upper limit on total surface rate.
695-
% vrat: Lower limit on total surface rate.
696-
flags_inj = [bhp > limits.bhp, ...
697-
qsTot > limits.rate, ...
698-
qsTot < limits.vrat];
699-
700-
% Producers have several possible limits:
701-
% bhp: Lower limit on pressure.
702-
% orat: Lower limit on surface oil rate
703-
% lrat: Lower limit on surface liquid (water + oil) rate
704-
% grat: Lower limit on surface gas rate
705-
% wrat: Lower limit on surface water rate
706-
% vrat: Upper limit on total volumetric surface rate
707-
qsm = zeros(model.numWells() + model.numGroups(), 3);
708-
qsm(:, model.getActivePhases()) = cell2mat(qs);
709-
flags_prod = [bhp < limits.bhp , ...
710-
qsm(:,2) < limits.orat, ...
711-
sum(qsm(:,1:2),2) < limits.lrat, ...
712-
qsm(:,3) < limits.grat, ...
713-
qsm(:,1) < limits.wrat, ...
714-
sum(qsm, 2) > limits.vrat];
714+
% Convenient well and group indices
715+
nw = model.numWells();
716+
ng = model.numGroups();
717+
[iic, iif] = model.getInletSegments();
718+
719+
[bhp, rhoS, Qt] = model.getProps(state, ...
720+
'BottomholePressure', 'SurfaceDensity', 'qt');
721+
qs = Qt./rhoS{1}(iic);
722+
is_inj = vertcat(model.wells.sign) > 0;
723+
724+
lims = [limits.bhp, limits.rate];
725+
vals = [value(bhp), value(qs) ];
726+
d = vals - lims;
727+
728+
flags = false(nw + ng, 2);
729+
flags( is_inj,:) = d( is_inj,:) > 0;
730+
flags(~is_inj,:) = d(~is_inj,:) < 0;
715731

716732
% limits we need to check (all others than w.type):
717-
types = [{model.wells.type}, {model.groups.type}]';
733+
718734
modes = fieldnames(limits)';
719-
chkInx = cellfun(@(t) strcmpi(t,modes), types, 'UniformOutput', false);
720-
vltInx = find(flags(chkInx), 1);
721-
if ~isempty(vltInx)
722-
withinLimits = false;
723-
modes = modes(chkInx);
724-
switchMode = modes{vltInx};
725-
fprintf('Well %s: Control mode changed from %s to %s.\n', wellSol.name, wellSol.type, switchMode);
726-
wellSol.type = switchMode;
727-
wellSol.val = limits.(switchMode);
735+
check = cellfun(@(t) ~strcmpi(t,modes), type, 'UniformOutput', false);
736+
check = cell2mat(check);
737+
738+
violate = flags & check;
739+
has_violation = any(violate, 2);
740+
741+
modes = repmat(modes, nw + ng, 1);
742+
type(has_violation) = modes(violate);
743+
target(has_violation) = lims(violate);
744+
745+
for i = 1:nw
746+
if has_violation(i)
747+
fprintf('Well %s: Control mode changed from %s to %s.\n', ...
748+
model.wells(i).name, state.wells(i).type, type{i});
749+
end
750+
state.wells(i).type = type{i};
751+
state.wells(i).val = target(i);
752+
state.wells(i).T = targetTemp(i);
728753
end
729-
730-
if ~withinLimits
731-
v = wellSol.val;
732-
switch wellSol.type
733-
case 'bhp'
734-
wellSol.bhp = bhp;
735-
case 'rate'
736-
for ix = 1:numel(phases)
737-
wellSol.(['q', phases(ix), 's']) = v*well.W.compi(ix);
738-
end
739-
case 'orat'
740-
wellSol.qOs = v;
741-
case 'wrat'
742-
wellSol.qWs = v;
743-
case 'grat'
744-
wellSol.qGs = v;
745-
end % No good guess for qOs, etc...
754+
755+
for i = nw + (1:ng)
756+
state.groups(i-nw).type = type{i};
757+
state.groups(i-nw).val = target(i);
758+
state.groups(i-nw).T = targetTemp(i);
746759
end
760+
747761
end
748762

749763
%-----------------------------------------------------------------%
@@ -758,15 +772,15 @@
758772
limits.(type{:}) = facility.val;
759773
end
760774

761-
modes = {'bhp', 'rate', 'orat', 'lrat', 'grat', 'wrat', 'vrat'};
775+
modes = {'bhp', 'rate'};
762776
modes = setdiff(modes, facility.type);
763777
missing = modes(~isfield(limits, modes));
764778

765779
sgn = sign(facility.sign);
766780
if sgn == 0, sgn = 1; end
767781
val = sgn.*inf;
768782
for f = missing
769-
v = val; if strcmpi(f, 'vrat'), v = -v; end
783+
v = val; %if strcmpi(f, 'vrat'), v = -v; end
770784
limits = setfield(limits, f{:}, v);
771785
end
772786
facility.lims = limits;
@@ -825,6 +839,8 @@
825839
qh = model.parentModel.getProp(state, 'HeatFlux');
826840
state.effect = model.parentModel.operators.groupSum(qh(iif));
827841
end
842+
843+
state = rmfield(state, {'wells', 'groups'});
828844

829845
if ~model.thermal, return; end
830846

@@ -999,6 +1015,12 @@
9991015
end
10001016
%-----------------------------------------------------------------%
10011017

1018+
%-----------------------------------------------------------------%
1019+
function types = getTypes(model)
1020+
types = {model.wells.types};
1021+
end
1022+
%-----------------------------------------------------------------%
1023+
10021024
%-----------------------------------------------------------------%
10031025
% Conveient parent model calls
10041026
%-----------------------------------------------------------------%

‎modules/geothermal/utils/getForceProperties.m‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -172,11 +172,12 @@
172172
Tbc = model.AutoDiffBackend.convertToAD(force.T, propsRes.pressure);
173173
is_Hflux = ~isnan(force.Hflux);
174174
if any(is_Hflux)
175-
ThR = propsRes.Thr;
176-
ThF = propsRes.Thf;
177-
T = propsRes.T;
178-
dT = force.Hflux(is_Hflux)./(ThR(is_Hflux) + ThF(is_Hflux));
179-
Tbc(is_Hflux) = T(is_Hflux) + dT;
175+
ThR = propsRes.Thr(is_Hflux);
176+
ThF = propsRes.Thf(is_Hflux);
177+
Q = force.Hflux(is_Hflux);
178+
T = propsRes.T(is_Hflux);
179+
dT = Q./(ThR + ThF);
180+
Tbc(is_Hflux) = T + dT;
180181
end
181182
end
182183

‎modules/geothermal/utils/getWellSolsFromWellboreState.m‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,10 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
5656

5757
wellSols = cell(n,1);
5858
wellNames = {model.submodels.(wbname).wells.name};
59-
groupNames = {model.submodels.(wbname).groups.name};
59+
groupNames = {};
60+
if model.submodels.(wbname).numGroups > 0
61+
groupNames = {model.submodels.(wbname).groups.name};
62+
end
6063

6164
nw = model.submodels.(wbname).numWells();
6265
ng = model.submodels.(wbname).numGroups();

‎modules/geothermal/utils/initializeGeothermalEquilibrium.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
127127
assert(~(isempty(solUp) && isempty(solDown)), 'Something went wrong');
128128
% Evaluate in cells
129129
[p,T] = evaluateFn(solUp, solDown, z0, z);
130-
sol = @(z) evaluateFn(solUp, solDown, z);
130+
sol = @(z) evaluateFn(solUp, solDown, z0, z);
131131

132132
end
133133

‎modules/geothermal/utils/processWellboreTrajectories.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -203,7 +203,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
203203

204204
endpoints = find(sum(A,2) == 1);
205205

206-
[~, ix] = min(G.cells.centroids(endpoints,3));
206+
[~, ix] = min(G.cells.centroids(cells(endpoints),3));
207207
first = endpoints(ix);
208208

209209
cells = nan(n,1);

0 commit comments

Comments
 (0)