-
Notifications
You must be signed in to change notification settings - Fork 0
/
computeTOFAndFirstArrival.m
385 lines (341 loc) · 13.6 KB
/
computeTOFAndFirstArrival.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
%%%Copyright 2018 SINTEF AS
function [T, A, q] = computeTOFAndFirstArrival(state, G, rock, varargin)
%Compute time of flight using finite-volume scheme.
%
% SYNOPSIS:
% T = computeTimeOfFlight(state, G, rock)
% T = computeTimeOfFlight(state, G, rock, 'pn1', pv1, ...)
% [T, A] = computeTimeOfFlight(...)
% [T, A, q] = computeTimeOfFlight(...)
%
% DESCRIPTION:
% Compute time-of-flight by solving
%
% v·\nabla T = \phi
%
% using a first-order finite-volume method with upwind flux. Here, 'v' is
% the Darcy velocity, '\phi' is the porosity, and T is time-of-flight.
% The time-of-flight T(x) is the time it takes an inert particles that is
% passively advected with the fluid to travel from the nearest inflow
% point to the point 'x' inside the reservoir. For the computation to
% make sense, inflow points must be specified in terms of inflow
% boundaries, source terms, wells, or combinations of these.
%
% REQUIRED PARAMETERS:
% G - Grid structure.
%
% rock - Rock data structure.
% Must contain a valid porosity field, 'rock.poro'.
%
% state - Reservoir state structure, must contain valid cell interface
% fluxes, 'state.flux'. Typically, 'state' will contain the
% solution produced by a flow solver like 'incompTPFA'.
%
% OPTIONAL PARAMETERS (supplied in 'key'/value pairs ('pn'/pv ...)):
% wells - Well structure as defined by function 'addWell'. May be empty
% (i.e., wells = []) which is interpreted as a model without any
% wells.
%
% src - Explicit source contributions as defined by function
% 'addSource'. May be empty (i.e., src = []) which is
% interpreted as a reservoir model without explicit sources.
%
% bc - Boundary condition structure as defined by function 'addBC'.
% This structure accounts for all external boundary conditions
% to the reservoir flow. May be empty (i.e., bc = []) which is
% interpreted as all external no-flow (homogeneous Neumann)
% conditions.
%
% reverse - Boolean variable. If true, the reverse time-of-flight will be
% computed, i.e., the travel time from a cell inside the
% reservoir to the nearest outflow point (well, source, or
% outflow boundary). Default: FALSE, i.e, compute forward tof.
%
% tracer - Cell-array of cell-index vectors for which to solve a
% stationary tracer equation:
%
% \nabla·(v C) = 0, C(inflow,t)=1
%
% One equation is solved for each vector with tracer injected in
% the cells with the given indices. Each vector adds one
% additional right-hand side to the original tof-system. Output
% given as additional columns in T.
%
% allowInf - Switch to turn off (true) or on (false) maximal TOF
% thresholding. Default is false.
%
% maxTOF - Maximal TOF thresholding to avoid singular/ill-conditioned
% systems. Default (empty) is 50*PVI (pore-volumes-injected).
% Only takes effect if 'allowInf' is set to false.
%
% solver - Function handle to solver for use in TOF/tracer equations.
% Default (empty) is matlab mldivide (i.e., \)
%
% processCycles - Extend TOF thresholding to strongly connected
% components in flux graph by considering Dulmage-Mendelsohn
% decomposition (dmperm). Recommended for highly cyclic flux
% fields. Only takes effect if 'allowInf' is set to false.
%
%
%
% RETURNS:
% T - Cell values of a piecewise constant approximation to time-of-flight
% computed as the solution of the boundary-value problem
%
% (*) v · \nabla T = \phi
%
% using a finite-volume scheme with single-point upwind approximation
% to the flux.
%
% A - Discrete left-hand side of (*), a G.cells.num-by-G.cells.num matrix
% whose entries are
%
% A_ij = min(F_ij, 0), and
% A_ii = sum_j max(F_ij, 0) + max(2*q_i, 0),
%
% where F_ij = -F_ji is the flux from cell i to cell j
%
% F_ij = A_ij n_ij · v_ij.
%
% and n_ij is the outward-pointing normal of cell i for grid face ij.
% The discretization uses a simple model for cells containing inflow.
% If q_i denotes the rate and V_i the volume of cell i, then T_i is
% set to half the time it takes to fill the cell, T_i = V_i/(2q_i).
%
% OPTIONAL. Only returned if specifically requested.
%
% q - Aggregate source term contributions (per grid cell) from wells,
% explicit sources and boundary conditions. These are the
% contributions referred to as 'q_i' in the definition of the matrix
% elements, 'A_ii'. Measured in units of m^3/s.
%
% OPTIONAL. Only returned if specifically requested.
%
% SEE ALSO:
% simpleTimeOfFlight, solveIncompFlow.
%{
Copyright 2009-2016 SINTEF ICT, Applied Mathematics.
This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).
MRST is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MRST is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with MRST. If not, see <http://www.gnu.org/licenses/>.
%}
opt = struct('bc', [], ...
'src', [], ...
'wells', [], ...
'reverse', false,...
'allowInf', false, ...
'maxTOF', [], ...
'tracer', {{}}, ...
'solver', [], ...
'processCycles', false, ...
'computeWellTOFs', false);
opt = merge_options(opt, varargin{:});
checkInput(G, rock, opt)
tr = opt.tracer;
if ~iscell(tr), tr = {tr}; end
% Find external sources of inflow: q contains the contribution from from
% sources (src) and wells (W), while qb contains the cell-wise sum of
% fluxes from boundary conditions.
[q,qb] = computeSourceTerm(state, G, opt.wells, opt.src, opt.bc);
pv = poreVolume(G, rock);
if opt.reverse,
q = -q;
qb = -qb;
state.flux = -state.flux;
end
% Build upwind flux matrix in which we define v_ji = max(flux_ij, 0) and
% v_ij = -min(flux_ij, 0). Then the diagonal of the discretization matrix
% is obtained by summing rows in the upwind flux matrix. This will give the
% correct diagonal in all cells except for those with a positive fluid
% source. In these cells, the average time-of-flight is set to half the
% time it takes to fill the cell, which means that the diagonal entry
% should be equal twice the fluid rate inside these cells.
i = ~any(G.faces.neighbors==0, 2);
n = double(G.faces.neighbors(i,:));
nc = G.cells.num;
qp = max(q+qb, 0);
out = min(state.flux(i), 0);
in = max(state.flux(i), 0);
% Cell wise total inflow
inflow = accumarray([n(:, 2); n(:, 1)], [in; -out]);
% The diagonal entries are equal to the sum of outfluxes minus divergence
% which equals the influx plus the positive source terms.
d = inflow + qp;
%--------------------------------------------------------------------------
% Handling of t -> inf (if opt.allowInf == false):
% We locate cells that have sufficiently small influx to reach maxTOF
% locally, i.e., t-t_up >= maxTOF <=> pv/influx >= maxTOF. System is
% modified such that these cells are set to maxTOF and upstream connections
% are removed.
% If no maxTOF is given, it is set to the time it takes to fill 50
% pore volumes (subject to change)
if ~opt.allowInf
if isempty(opt.maxTOF)
if ~opt.reverse, str = 'Forward '; else str = 'Backward'; end
opt.maxTOF = 50*sum(pv)/sum(qp);
dispif(mrstVerbose, ...
'%s maximal TOF set to %5.2f years.\n', str, opt.maxTOF/year)
end
% Find cells that reach max TOF locally
maxIn = max(d);
aboveMax = full(pv./ (max(d, eps*maxIn)) ) > opt.maxTOF;
% Set aboveMax-cells to maxTOF
d(aboveMax) = maxIn;
pv(aboveMax) = maxIn*opt.maxTOF;
% remove upstream connections
in(aboveMax(n(:,2))) = 0;
out(aboveMax(n(:,1))) = 0;
end
% Inflow flux matrix
A = sparse(n(:,2), n(:,1), in, nc, nc)...
+ sparse(n(:,1), n(:,2), -out, nc, nc);
A = -A + spdiags(d, 0, nc, nc);
if ~opt.allowInf && opt.processCycles
[A, pv] = thresholdConnectedComponents(A, pv, maxIn, opt);
end
% Build right-hand sides for tracer equations. Since we have doubled the
% rate in any cells with a positive source, we need to also double the rate
% on the right-hand side here.
numTrRHS = numel(tr);
TrRHS = zeros(nc,numTrRHS);
for i=1:numTrRHS,
TrRHS(tr{i},i) = qp(tr{i});
end
% Time of flight for velocity field.
if isempty(opt.solver)
T = A \ [pv TrRHS];
else % if other solver, iterate over RHSs
T = zeros(size(TrRHS)+[0, 1]);
T(:, 1) = opt.solver(A, pv);
for k = 2:size(T, 2)
T(:, k) = opt.solver(A, TrRHS(:, k-1));
end
end
% reset all tof > maxTOF to maxTOF
if ~opt.allowInf
T(T>opt.maxTOF) = opt.maxTOF;
end
% compute individual well-tofs A*(c_i.*tau_i) = c_i.*pv
if opt.computeWellTOFs
%if false
C = T(:, 2:end);
pvi = bsxfun(@times, C, pv);
pvi(pvi<0) = 0;
if isempty(opt.solver)
X = A \ pvi;
else % if other solver, iterate over RHSs
X = zeros(size(pvi));
for k = 1:size(X, 2)
X(:, k) = opt.solver(A, pvi(:, k));
end
end
X(X<0) = 0;
% disregard tracer values < sqrt(eps)
ix = and(pvi*opt.maxTOF > X, C > sqrt(eps));
X(~ix) = opt.maxTOF;
X(ix) = X(ix)./C(ix);
T = [T, X];
%else
n(out<0, [1 2]) = n(out<0, [2 1]);
% add edges for well-to-connection cell
nconn = cellfun(@numel, tr);
wnum = rldecode((1:numel(tr))', nconn(:));
wc = vertcat(tr{:});
n1 = [n(:,1); wnum+nc];
n2 = [n(:,2); wc];
fl = [inflow(n(:,2)); qp(wc)];
fl = max(fl, max(fl)*eps);
AM = sparse(n1, n2, pv(n2)./fl, nc+max(wnum), nc+max(wnum));
gr = digraph(AM);
% gr = digraph(n1, n2, pv(n2)./fl);
X = T(:, 2:(numTrRHS+1));
for k = 1:size(X,2)
[~, t] = shortestpathtree(gr, nc+k, 'OutputForm', 'vector');
X(:,k) = t(1:nc);
end
X = min(X, opt.maxTOF);
%end
%X = min(X, opt.maxTOF);
T = [T, X];
end
end
function [q, qb] = computeSourceTerm(state, G, W, src, bc)
qi = []; % Cells to which sources are connected
qs = []; % Actual strength of source term (in m^3/s).
% Contribution from wells
if ~isempty(W),
qi = [qi; vertcat(W.cells)];
qs = [qs; vertcat(state.wellSol.flux)];
end
% Contribution from sources
if ~isempty(src),
qi = [qi; src.cell];
qs = [qs; src.rate];
end
% Assemble all source and sink contributions to each affected cell.
q = sparse(qi, 1, qs, G.cells.num, 1);
% Contribution from boundary conditions
if ~isempty(bc),
ff = zeros(G.faces.num, 1);
isDir = strcmp('pressure', bc.type);
i = bc.face(isDir);
if ~isempty(i)
ff(i) = state.flux(i) .* (2*(G.faces.neighbors(i,1)==0) - 1);
end
isNeu = strcmp('flux', bc.type);
ff(bc.face(isNeu)) = bc.value(isNeu);
is_outer = ~all(double(G.faces.neighbors) > 0, 2);
qb = sparse(sum(G.faces.neighbors(is_outer,:), 2), 1, ...
ff(is_outer), G.cells.num, 1);
else
qb = sparse(G.cells.num,1);
end
end
function [A, pv] = thresholdConnectedComponents(A, pv, maxIn, opt)
% Find strongly connected components in flux-matrix:
[p,r,r]= dmperm(A); %#ok
% Pick components containing more than a single cell
ix = find(diff(r)>1);
if ~isempty(ix)
nc = numel(p);
% Retrieve cell-indices to components, and construct sparse index-mapping
c = arrayfun(@(b,e)p(b:e)', r(ix), r(ix+1)-1, 'UniformOutput', false);
rc = rldecode( (1:numel(c))', cellfun(@numel, c)');
C = sparse(vertcat(c{:}), rc, 1, nc, numel(c));
% Compute influx to each component
q_in = full(diag(C'*A*C));
% Threshold
compAboveMax = full((C'*pv)./ (max(q_in, eps*maxIn)) ) > opt.maxTOF;
if any(compAboveMax)
badCells = (vertcat(c{compAboveMax}));
dispif(mrstVerbose, 'Found %d strongly connected components, ', nnz(compAboveMax));
dispif(mrstVerbose, 'total of %d cells, with influx below threshold.\n', numel(badCells));
% Modify system setting tof to maxTOF and remove upstream connections
A(badCells,:) = sparse((1:numel(badCells))', badCells, maxIn, numel(badCells), nc);
pv(badCells) = maxIn*opt.maxTOF;
end
end
end
function checkInput(G, rock, opt)
assert (~all([isempty(opt.src), isempty(opt.bc), isempty(opt.wells)]), ...
'Must have inflow described as boundary conditions, sources, or wells');
assert (isfield(rock, 'poro') && ...
numel(rock.poro)==G.cells.num, ...
['The rock input must have a field poro with porosity ',...
'for each cell in the grid.']);
assert(min(rock.poro) > 0, 'Rock porosities must be positive numbers.');
if ~isempty(opt.maxTOF) && opt.allowInf
warning('Input value for ''maxTOF'' ignored since option ''allowInf'' has value ''true''.')
end
if opt.processCycles && opt.allowInf
warning('Input request to process cycles ignored since option ''allowInf'' has value ''true''.')
end
end