mrstModule add diagnostics gridtools mrst-gui grdecl = fullfile(ROOTDIR, 'examples', 'data', 'SAIGUP', 'SAIGUP.GRDECL'); if ~exist(grdecl, 'file'), error('SAIGUP model data is not available.') end grdecl = readGRDECL(grdecl); actnum = grdecl.ACTNUM; grdecl.ACTNUM = ones(prod(grdecl.cartDims),1); G = processGRDECL(grdecl, 'checkgrid', false); G = computeGeometry(G(1)); rock = grdecl2Rock(grdecl, G.cells.indexMap); is_pos = rock.perm(:, 3) > 0; rock.perm(~is_pos, 3) = min(rock.perm(is_pos, 3)); rock.perm = rock.perm.*milli*darcy;
The wells are simply placed by a double for loop over the logical indices. The resulting wells give producers along the edges of the domain and injectors near the center.
pv = poreVolume(G, rock); ijk = gridLogicalIndices(G); W = []; gcz = G.cells.centroids; [pi,ii] = deal(1); for i = 5:12:G.cartDims(1) for j = 5:25:G.cartDims(2) c = ijk{1} == i & ijk{2} == j; if any(c) c = find(c); x = gcz(c(1), 1); y = gcz(c(1), 2); if x > 500 && x < 2000 && y < 7000 && y > 1000 % Set up rate controlled injectors for cells in the middle % of the domain val = sum(pv)/(1000*day); type = 'rate'; name = ['I' num2str(ii)]; ii = ii + 1; else % Set up BHP controlled producers at the boundary of the % domain val = 250*barsa; type = 'bhp'; name = ['P' num2str(pi)]; pi = pi + 1; end W = addWell(W, G, rock, c, 'Type', type, 'Val', val, 'Name', name); end end end % Plot the resulting well setup clf plotCellData(G, rock.poro), plotWell(G, W) view(-100, 25)
This is not strictly needed for the diagnostics application, but it allows us to see the reservoir component distribution around producers grouped by time of flight. The distribution of phases is done based on a simple hydrostatic approximation by using cell centroids. Some phase mixing is present.
state = initResSol(G, 200*barsa, [0 0 0]); gcz = G.cells.centroids(:, 3); height = max(gcz) - min(gcz); pos = 1 - (gcz - min(gcz))./height; oil = pos > 0.4 & pos < 0.8; gas = pos > 0.7; wat = pos < 0.45; state.s(wat, 1) = 1; state.s(oil, 2) = 1; state.s(gas, 3) = 1; state.s = state.s ./ repmat(sum(state.s, 2), 1, 3); clf; rgb = @(x) x(:, [2 3 1]); plotCellData(G, rgb(state.s), 'EdgeColor', 'k', 'EdgeAlpha', 0.1) view(-100, 25) axis tight
We print the help page to show functionality and then invoke the interactive diagnostics tool.
help interactiveDiagnostics clf; interactiveDiagnostics(G, rock, W, 'state', state); view(-100, 25) axis tight
Published April 18, 2013