dg:

Contents

MODELS

Files
TransportModelDG -
Contents

UTILS

Files
computeBoundaryFluxesDG - Undocumented Utility Function
computeBoundaryFluxesDG(model, state, bc)

Undocumented Utility Function

Contents

AD

Files
SpatialVector - Class for supporting ad as spatial vectors, for use in e.g.,
class SpatialVector(varargin)

Class for supporting ad as spatial vectors, for use in e.g., Discontinuous Galerkin discretizations

subsref(v, s)

Subscripted reference. Called for h = u(v).

times(u, v)

assert(all(size(u) == size(v)));

Contents

CUBATURES

Files
Cubature - Cubature class for cubatures on MRST grids LineCubature - 1D line cubatureature for faces in 2D MRST grids MomentFitting2DCubature - Cubature based on moment-fitting for MRST grids MomentFitting3DCubature - Cubature based on moment-fitting for MRST grids TetrahedronCubature - TriangleCubature - Triangle cubature class for cubatures on MRST grids
class Cubature(G, prescision, dim, varargin)

Cubature class for cubatures on MRST grids

Cubature(G, prescision, dim, varargin)

Setup up cubature properties. Acutual construction of cubature is handled by children class.

assignCubaturePointsAndWeights(cubature, x, w, n)

Assign points and weights

getCubature(cubature, elements, type, varargin)

Get cubature for given elements (either a set of cells or a set of faces) of the grid.

Parameters:
  • elements – Vector of cells or faces we want cubature for
  • type

    String equal to either * ‘volume’ : Volume integral, elements

    interpreted as cells
    • ’face’ : Surface integral, elements
      interpreted as faces
    • ’surface’: Surface integral over all faces of
      elements, which are interpreted as cells
Keyword Arguments:
 
  • excludeBoundary – Exclude boundary faces in cell surface integrals
  • internalConn – Boolean indicating what faces are internal connections, and therefore included in surface integrals. If empty, calculated from G.faces.neighbors
  • outwardNormal – Surface integrals over cells are often on the form int(vdot n). outwardNormal = true changes signs based on face orientation so that normal points outwards.
Returns:
  • W – Integration matrix - integral = W*integrand(x)
  • x – Integration points
  • cellNo, faceNo – cell/face number telling what cubature point goes where
makeCubature(cubature)

#ok

transformCoords(cub, x, cells, inverse)

Transfor coordinates from physical to reference coordinates.

Parameters:
  • x – Coordinates in physical space
  • cells – Cells we want reference coordinates for, cells(ix) are used when transforming x(ix,:)
  • inverse – Boolean indicatiing if we are mapping to (inverse = false) or from (inverse = true) reference coordiantes. Default = false.
  • useParent – Boolean indicating if we are working on the full grid (G.parent) or a subgrid.
Returns:
  • xhat – Transformed coordinates
  • translation – Translation applied to x
  • scaling – Scaling applied to x
G = None

Grid the cubature is defined on

dim = None

Dimension of cubature

internalConn = None

Boolean indicating if face is internal connection

numPoints = None

Number of cubature points

points = None

Cubature points

pos = None

Cubature points/weights for element e can be found

prescision = None

Prescision of cubature

weights = None

Cubature weights

class LineCubature(G, prescision)

1D line cubatureature for faces in 2D MRST grids

LineCubature(G, prescision)

Set up line cubatureature

makeCubature(cubature)

Make cubatureature

mapCoords(cubature, x, xR)

Map cubatureature points from reference to physical coordinates

class MomentFitting2DCubature(G, prescision, varargin)

Cubature based on moment-fitting for MRST grids

MomentFitting2DCubature(G, prescision, varargin)

Set up cubature

makeCubature(cubature)

Dimension of cubature

class MomentFitting3DCubature(G, prescision)

Cubature based on moment-fitting for MRST grids

MomentFitting3DCubature(G, prescision)

Set up cubatureature

makeCubature(cubature)

Dimension of cubature

class TriangleCubature(G, prescision)

Triangle cubature class for cubatures on MRST grids

TriangleCubature(G, prescision)

Set up triangle cubature

getTriAreas(cubature)

Get triangle areas

getTriangulation(cubature, G)

#ok Triangulate all cells or faces of G

makeCubature(cubature)

Make cubature

mapCoords(cubature, x, xR)

Map cubature points from reference to physical coordinates

mapCoordsBary(cubature, x)

Map cubature points from Barycentric to physical coordinates

triangulation = None

Struct containing triangulation information

Contents

DISCRETIZATIONS

Files
DGDiscretization - SpatialDiscretization -
Contents

STATEFUNCTIONS

Files
FluxDiscretizationDG -
Contents

FLOWPROPS

Files
MultipliedPoreVolumeDG - Effective pore-volume after pressure-multiplier
class MultipliedPoreVolumeDG(model, varargin)

Effective pore-volume after pressure-multiplier

evaluateOnDomain(prop, model, state)

Get effective pore-volume, accounting for possible rock-compressibility

Contents

FLUX

Files
ComponentPhaseVelocityFractionalFlowDG - ComponentTotalVelocityDG - FixedTotalFluxDG - FixedTotalVelocityDG - GravityPermeabilityGradientDG - GravityPotentialDifferenceDG - PhasePotentialUpwindFlagDG -

Examples

Buckley-Leverett displacement

Generated from dgExampleBLDisplacement.m

In this example, we consider an incompressible Buckley-Leverett displacement in a horizontal 1D channel aligned with the x-axis, and compare higher-order discontinuous Galerkin methods with the standard finite-volume discretization.

mrstModule add dg ad-core ad-props ad-blackoil blackoil-sequential
saveeps = @(a,b) disp(b);  % Dummy function

Set up problem

Construct grid, compute geometry and cell dimensions, and set petrophysical properties

G    = computeGeometry(cartGrid([30,1]));
G    = computeCellDimensions(G);
rock = makeRock(G, 1, 1);

% We consider Buckley-Leverett-type displacement with quadratic relperms
fluid = initSimpleADIFluid('phases', 'WO' , 'n', [2,2], ...
                           'mu', [1,1],  'rho', [1,1]);

% The base model is a generic black-oil model with oil and water
model  = GenericBlackOilModel(G, rock, fluid, 'gas', false);

% Initial state: filled with oil and unit volumetric flow rate
state0 = initResSol(G, 1, [0,1]);
state0.flux(1:G.cells.num+1) = 1;

% Boundary conditions
bc = fluxside([], G, 'left' ,  1, 'sat', [1,0]); % Inflow
bc = fluxside(bc, G, 'right', -1, 'sat', [1,0]); % Outflow

% Define schedule: unit time steps, simulate almost to breakthrough
schedule = simpleSchedule(ones(30,1), 'bc', bc);

Finite volume simulation

Use standard discretization in MRST: single-point upwind (SPU)

tmodel           = TransportModel(model);
[~, stFV, repFV] = simulateScheduleAD(state0, tmodel, schedule);
Solving timestep 01/30:            -> 1 Second
Solving timestep 02/30: 1 Second   -> 2 Seconds
Solving timestep 03/30: 2 Seconds  -> 3 Seconds
Solving timestep 04/30: 3 Seconds  -> 4 Seconds
Solving timestep 05/30: 4 Seconds  -> 5 Seconds
Solving timestep 06/30: 5 Seconds  -> 6 Seconds
Solving timestep 07/30: 6 Seconds  -> 7 Seconds
Solving timestep 08/30: 7 Seconds  -> 8 Seconds
...

Lowest-order dG simulation

Simulate using dG(0). This should be equivalent to SPU. The model is built as a wrapper around the base model, so that all physics are consistently handled by this. The base model is stored as the property ‘’parentModel’’

tmodelDG0 = TransportModelDG(model, 'degree', 0);
disp(tmodelDG0)
Getting supports.
Elapsed time is 0.078643 seconds.
Inverting systems...
  TransportModelDG with properties:

          discretization: [1×1 DGDiscretization]
             dgVariables: {'s'  'sW'  'sO'  'sT'}
                limiters: [2×1 struct]
...

Simulate dG(0)

[~, stDG0, repDG0] = simulateScheduleAD(state0, tmodelDG0, schedule);
Solving timestep 01/30:            -> 1 Second
Solving timestep 02/30: 1 Second   -> 2 Seconds
Solving timestep 03/30: 2 Seconds  -> 3 Seconds
Solving timestep 04/30: 3 Seconds  -> 4 Seconds
Solving timestep 05/30: 4 Seconds  -> 5 Seconds
Solving timestep 06/30: 5 Seconds  -> 6 Seconds
Solving timestep 07/30: 6 Seconds  -> 7 Seconds
Solving timestep 08/30: 7 Seconds  -> 8 Seconds
...

First-order dG simulation

TransportModelDG also takes optional arguments to DGDiscretization, such as degree. DGDiscretization supports using different polynomial order in each coordinate direction. We are simulating a 1D problem on a 2D grid, so we use higher order in the x-direction only.

tmodelDG1 = TransportModelDG(model, 'degree', [1,0]);
Getting supports.
Elapsed time is 0.041544 seconds.
Inverting systems...

Limiters

Higher-order dG methods typically exhibit spurious oscillations near discontinuities. We countermand this with a limiter that effectively introduces enough numerical diffusion to dampen oscillations. We use two limiters by default: a TVB (Total Variation Bounded) limiter that adjusts the gradient whenever the jump across an interface is greater than 0, and a physical limiter that scales the solution so that the minimum and maximum in each cell is within physical limits. To see the limiters in action, we use ‘’afterStepFn’’ to plot the unlimited and limited solution after each timestep

tmodelDG1.storeUnlimited = true; % Store the unlimited state in each step
fn = plotLimiter(tmodelDG1, 'plot1d', true, 'n', 500); % afterStepFn
_images/dgExampleBLDisplacement_01.png

Simulate dG(1)

[~, stDG1, repDG1] = ...
    simulateScheduleAD(state0, tmodelDG1, schedule, 'afterStepFn', fn);
Solving timestep 01/30:            -> 1 Second
Solving timestep 02/30: 1 Second   -> 2 Seconds
Solving timestep 03/30: 2 Seconds  -> 3 Seconds
Solving timestep 04/30: 3 Seconds  -> 4 Seconds
Solving timestep 05/30: 4 Seconds  -> 5 Seconds
Solving timestep 06/30: 5 Seconds  -> 6 Seconds
Solving timestep 07/30: 6 Seconds  -> 7 Seconds
Solving timestep 08/30: 7 Seconds  -> 8 Seconds
...
_images/dgExampleBLDisplacement_02.png

Higher-order dG simulations

In principle, the dg module supports arbitrary order, as long as you can provide a sufficiently accurate quadrature rule. These are stored in e.g., ‘’getELEMENTCubaturePointsAndWeights’’, ELEMENT in {Triangle, Square, Tetrahedron, Cube}. We simulate the problem with dG(2)

tmodelDG2 = TransportModelDG(model, 'degree', [2,0]);
tmodelDG2.storeUnlimited = true;
fn = plotLimiter(tmodelDG2, 'plot1d', true, 'n', 500);
Getting supports.
Elapsed time is 0.038878 seconds.
Inverting systems...
_images/dgExampleBLDisplacement_03.png

Simulate dG(2)

[~, stDG2, repDG2] = ...
    simulateScheduleAD(state0, tmodelDG2, schedule, 'afterStepFn', fn);
Solving timestep 01/30:            -> 1 Second
Solving timestep 02/30: 1 Second   -> 2 Seconds
Solving timestep 03/30: 2 Seconds  -> 3 Seconds
Solving timestep 04/30: 3 Seconds  -> 4 Seconds
Solving timestep 05/30: 4 Seconds  -> 5 Seconds
Solving timestep 06/30: 5 Seconds  -> 6 Seconds
Solving timestep 07/30: 6 Seconds  -> 7 Seconds
Solving timestep 08/30: 7 Seconds  -> 8 Seconds
...
_images/dgExampleBLDisplacement_04.png

Compare results

Finally, we plot the results, and compare them to the exact solution

close all
% Plotting options
coords = getPlotCoordinates(G, 'n', 500, 'plot1d', true);    % Plot coords
opt    = {'linewidth', 1.5, 'coords', coords, 'plot1d', true}; % Plot options
s      = buckleyLeverettProfile('nW', 2, 'nN', 2);       % Exact sol
tsteps = [5,10,20]; % Timesteps

for t = tsteps
    figure('Position', [100, 100, 800, 400])
    % Tweak FV solution so we can plot it as dG
    stFV_plot      = stDG0{t};
    stFV_plot.sdof = stFV{t}.s;
    hold on
    plotSaturationDG(tmodelDG0.discretization, stDG0{t}, opt{:});
    plotSaturationDG(tmodelDG1.discretization, stDG1{t}, opt{:});
    plotSaturationDG(tmodelDG2.discretization, stDG2{t}, opt{:});
    plotSaturationDG(tmodelDG0.discretization, stFV_plot, ...
        'plot1d', true, 'coords', coords, 'linewidth', 3, 'lineStyle', '--');
    plot(coords.points(:,1), s(coords.points(:,1), t), 'linewidth', 1.5, 'color', 'k')
    hold off
    legend({'dG(0)', 'dG(1)', 'dG(2)', 'SPU', 'exact'});
    box on
    ylim([-0.2, 1.2])
end
_images/dgExampleBLDisplacement_05.png
_images/dgExampleBLDisplacement_06.png
_images/dgExampleBLDisplacement_07.png

Plot dG(1) solution

close all
pos = [0,0,800,400];

figure('Position', pos)
t = 15;
hold on
    plotSaturationDG(tmodelDG1.discretization, stDG1{t}.ul, opt{:}, 'linewidth', 1.5, 'color', 'k');
    plotSaturationDG(tmodelDG1.discretization, stDG1{t}, opt{:}, 'linewidth', 4, 'color', 'k', 'linestyle', '--');
hold off
box on;
ylim([-0.2, 1.2])
ax = gca;
ax.FontSize = 14;
bx = [19.5, 22.5, -0.1, 0.25];
rectangle('Position', [bx([1,3]), bx([2,4]) - bx([1,3])], 'linestyle', '--');
drawnow(); pause(0.5), saveeps('buckley-leverett', 'limiter');

figure('Position', [0,0,350,400])
hold on
    plotSaturationDG(tmodelDG1.discretization, stDG1{t}.ul, opt{:}, 'linewidth', 1.5, 'color', 'k');
    plotSaturationDG(tmodelDG1.discretization, stDG1{t}, opt{:}, 'linewidth', 4, 'color', 'k', 'linestyle', '--');
hold off
box on;
axis(bx)
ax = gca;
ax.FontSize = 14;
[ax.XTick, ax.YTick] = deal([]);
legend({'Unlimited', 'Limited'}, 'fontSize', 15, 'Location', 'northeast')
drawnow(); pause(0.5), saveeps('buckley-leverett', 'limiter-zoom');
limiter
limiter-zoom
_images/dgExampleBLDisplacement_08.png
_images/dgExampleBLDisplacement_09.png

Compare FV and dG(0) - dG(2)

t = 15;
figure('Position', [100, 100, 800, 400])
% Tweak FV solution so we can plot it as dG
stFV_plot      = stDG0{t};
stFV_plot.sdof = stFV{t}.s;
hold on
plotSaturationDG(tmodelDG0.discretization, stDG0{t}, opt{:});
plotSaturationDG(tmodelDG1.discretization, stDG1{t}, opt{:});
plotSaturationDG(tmodelDG2.discretization, stDG2{t}, opt{:});
plotSaturationDG(tmodelDG0.discretization, stFV_plot, ...
    'plot1d', true, 'coords', coords, 'linewidth', 4, 'lineStyle', '--');
plot(coords.points(:,1), s(coords.points(:,1), t), 'linewidth', 1.5, 'color', 'k')
hold off
legend({'dG(0)', 'dG(1)', 'dG(2)', 'SPU', 'exact'}, 'fontsize', 14, 'location', 'southwest');
box on
ylim([-0.2, 1.2])
ax = gca;
ax.FontSize = 14;
bx = [16, 20, 0.35, 0.75];
rectangle('Position', [bx([1,3]), bx([2,4]) - bx([1,3])], 'linestyle', '--');
drawnow(); pause(0.5), saveeps('buckley-leverett', 'solutions')

figure('Position', [100, 100, 400, 400])
% Tweak FV solution so we can plot it as dG
stFV_plot      = stDG0{t};
stFV_plot.sdof = stFV{t}.s;
hold on
plotSaturationDG(tmodelDG0.discretization, stDG0{t}, opt{:});
plotSaturationDG(tmodelDG1.discretization, stDG1{t}, opt{:});
plotSaturationDG(tmodelDG2.discretization, stDG2{t}, opt{:});
plotSaturationDG(tmodelDG0.discretization, stFV_plot, ...
    'plot1d', true, 'coords', coords, 'linewidth', 4, 'lineStyle', '--');
plot(coords.points(:,1), s(coords.points(:,1), t), 'linewidth', 1.5, 'color', 'k')
hold off
box on
axis(bx)
ax = gca;
[ax.XTick, ax.YTick] = deal([]);
saveeps('buckley-leverett', 'solutions-zoom')
solutions
solutions-zoom
_images/dgExampleBLDisplacement_10.png
_images/dgExampleBLDisplacement_11.png

Components of a Discontinuous Galerkin Discretization in MRST

Generated from dgExampleDiscretization.m

This educational example goes through the main components of a discontinuous Galerkin (dG) discretization in MRST

mrstModule add dg
saveeps = @(a,b) disp(b);  % Dummy function
savepng = @(a,b) disp(b);  % Dummy function

Make grid

We construct a small PEBI grid

mrstModule add upr            % For generating PEBI grids
G = pebiGrid(1/8, [1,1]);     % PEBI grid
G = computeGeometry(G);       % Compute geometry
G = computeCellDimensions(G); % Compute cell dimensions

gray    = [1,1,1]*0.8; % For plotting
addText = @(x, string) text(x(:,1), x(:,2), string, 'fontSize', 12, 'HorizontalAlignment', 'left');

Plot the grid

We plot the grid and look at some geometric properties of a single cell

c = 22; % Cell number 17
% Plot cell 17 and its topological neighbors
figure();
plotGrid(G, 'facecolor', 'none');
cn = G.faces.neighbors(any(G.faces.neighbors == c,2),:);
cn = cn(cn>0);
plotGrid(G, cn(:), 'facecolor', 'none', 'linewidth', 2); axis equal off
plotGrid(G, c, 'facecolor', gray, 'linewidth', 2);
drawnow(), pause(0.2), saveeps('discretization', 'grid');
% By calling omputeCellDimensions(G), we computed the dimensions of a the
% smallest bounding box aligned with the coordinate axes so that the
% centroid of the box coincides with the cell centroid. We will use this to
% define our basis functions
xc = G.cells.centroids(c,:); % Centroid of cell 10
dx = G.cells.dx(c,:);        % Bounding box dimensions of cell 10
figure()
hold on
plotGrid(G, c, 'facecolor', gray);
rectangle('Position', [xc - dx/2, dx]);
plot(xc(1), xc(2), '.k', 'markerSize', 30);
plot([xc(1); xc(1)], [xc(2) - dx(2)/2; xc(2) + dx(2)/2], '--k');
plot([xc(1) - dx(1)/2; xc(1) + dx(1)/2], [xc(2); xc(2)], '--k');
axis equal off
drawnow(), pause(0.2), saveeps('discretization', 'cell');
grid
cell
_images/dgExampleDiscretization_01.png
_images/dgExampleDiscretization_02.png

Construct discretization

mrstModule add dg                        % Load module
G    = computeCellDimensions(G);         % Compute cell dimensions
disc = DGDiscretization(G, 'degree', 2); % Construct dG(2) discretization

Inspect basis

The basis functions are conveniently constructed by a Polynomial class. A polynomial is then described as a set of exponents k and weights w, and the class has overloaded operators for elementary arithmetic operations (addition, subtraction, multiplication, etc.), as well as tensor products, derivatives and gradients. We have construct a Legendre-type basis of polynomials of order <= 2.

disp(disc.basis)        % Inspect basis
disp(disc.basis.psi{3}) % The third basis function is simply y
psi: {6×1 cell}
    gradPsi: {6×1 cell}
     degree: 2
          k: [6×2 double]
       nDof: 6
       type: 'legendre'

  Polynomial with properties:
...

Plot polynomial basis

The basis functions are evaluated in local coordinates in each cell: $psi(hat{x}) = psi(frac{x - x_c}{Delta_x/2})

h = plotDGBasis(G, disc.basis, c, 'edgealpha', .05);
for i = 1:numel(h)
    set(0, 'CurrentFigure', h(i));
    ax = gca;
    [ax.XTickLabel, ax.YTickLabel] = deal({});
    ax.FontSize = 20; caxis([-1 1.2])
    savepng('discretization', ['basis-', num2str(i)]);
end
Warning: DistMesh did not converge in maximum number of iterations.
basis-1
basis-2
basis-3
basis-4
basis-5
basis-6
_images/dgExampleDiscretization_03.png
_images/dgExampleDiscretization_04.png
_images/dgExampleDiscretization_05.png
_images/dgExampleDiscretization_06.png
_images/dgExampleDiscretization_07.png
_images/dgExampleDiscretization_08.png

Construct triangle cubature

DG involves excessive evaluation of integrals. This is handeled by cubatures, which are implemented in the Cubature class. The simplest possible approach for polygonal cells is to subdivide the cells into simplices, and use a cubature rule for each simplex. This is implemented in the TriangleCUbature/TetrahedronCubature classes, which is the default cubature of the DGDiscretization

triCubature = disc.cellCubature;
figure(), plotCubature(G, triCubature, c, 'faceColor', gray); axis equal tight
savepng('discretization', 'cubature-tri');
cubature-tri
_images/dgExampleDiscretization_09.png

Test the triangle cubature

[W, x, ~, cells] = triCubature.getCubature(c, 'cell'); % Get points/weights
x = triCubature.transformCoords(x, cells); % Transform to reference coords
% The integral is easily evalauted using the matrix W, which has the
% cubature weights w placed so that int(psi)dx = W*psi(x). The weights sum
% to one, so that by construction, the integral of the first basis function
% (which is constant) should equal one, whereas the linear basis functions
% (1 and 2) should be zero
cellfun(@(psi) W*psi(x), disc.basis.psi)
ans =

    1.0000
    0.0000
    0.0000
   -0.1178
   -0.0006
...

Construct moment fitting cubature

The cubature above uses much more points than strictly needed. We can construct a more efficient cubature by moment fitting. In this approach, we use a known cubature to compute the integrals for all functions in the basis for the polyomials we want our rule to be exact for. Then, given an initial set of cubature points, we compute corresponding weights by linear least squares, successively removing points until it is no longer possible to take away more points. In this case, the result is a cubature with six points, which equals the number of basis functions Moment-fitting without redcution

mrstModule add vem vemmech
degree = 2;
mfCubature_nr = MomentFitting2DCubature(G, 2*degree, 'chunkSize', 1, 'reduce', false);
figure(), plotCubature(G, mfCubature_nr, c, 'faceColor', gray); axis equal tight
savepng('discretization', 'cubature-mf-nr');
% Moment-fitting with redcution
mfCubature = MomentFitting2DCubature(G, 2*degree, 'chunkSize', 1);
figure(), plotCubature(G, mfCubature, c, 'faceColor', gray); axis equal tight
savepng('discretization', 'cubature-mf');
Compressing quadrature for element 1 to 1 of 87 ... Compressed from 28 to 28 points in 0.010448 second
Compressing quadrature for element 2 to 2 of 87 ... Compressed from 28 to 28 points in 0.004636 second
Compressing quadrature for element 3 to 3 of 87 ... Compressed from 28 to 28 points in 0.003819 second
Compressing quadrature for element 4 to 4 of 87 ... Compressed from 28 to 28 points in 0.005281 second
Compressing quadrature for element 5 to 5 of 87 ... Compressed from 28 to 28 points in 0.003920 second
Compressing quadrature for element 6 to 6 of 87 ... Compressed from 28 to 28 points in 0.003870 second
Compressing quadrature for element 7 to 7 of 87 ... Compressed from 28 to 28 points in 0.003774 second
Compressing quadrature for element 8 to 8 of 87 ... Compressed from 28 to 28 points in 0.003800 second
...
_images/dgExampleDiscretization_10.png
_images/dgExampleDiscretization_11.png

Test the moment fitting cubature

[W, x, w, cells] = mfCubature.getCubature(c, 'cell');
x = mfCubature.transformCoords(x, cells);
cellfun(@(psi) W*psi(x), disc.basis.psi)
ans =

    1.0000
    0.0000
    0.0000
   -0.1178
   -0.0006
...

Construct line cubature

We also need a cubature for integrals over cell faces. The cubature for a an interface is used for both cells sharing the interface, but the functions are evaluated from different sides depending on which cell we are considering

lineCubature = disc.faceCubature;
faces     = G.cells.faces(G.cells.facePos(c):G.cells.facePos(c+1)-1);
figure(), plotCubature(G, lineCubature, faces(3), 'smax', 10, ...
            'faceColor', gray, 'plotBoundingBox', false); axis equal tight
savepng('discretization', 'cubature-line');
cubature-line
_images/dgExampleDiscretization_12.png

Simulation of an inverted five-spot pattern with dG

Generated from dgExampleIFS.m

In this example, we simulate water injection in an inverted five-spot pattern posed on a PEBI mesh with dG(0) and dG(1) and visualize the results

Add modules

mrstModule add dg upr vem vemmech spe10 ad-props ad-core ad-blackoil ...
    blackoil-sequential

Set up fine-scale model

We extract the lower half of layer 13 of SPE10 2

[state0Ref, modelRef, scheduleRef] = setupSPE10_AD('layers', 13           , ...
                                                   'J'     , (1:110) + 30);
GRef    = modelRef.G;
rockRef = modelRef.rock;
WRef    = scheduleRef.control(1).W;
fluid = modelRef.fluid;

Make PEBI grid

We use the upr module to construct a PEBI grid with refinement around the wells

rng(2019) % For reproducibility
n = 20;   % Approx number of cells in each direction
% Get well coordinates
l = max(GRef.nodes.coords(:,1:2));
wellLines = mat2cell(GRef.cells.centroids(vertcat(WRef.cells),1:2), ...
                                                    ones(numel(WRef),1), 2)';
% Construct PEBI grid
G = pebiGrid(max(l)/n, l, 'wellLines'     , wellLines, ...
% Well coords
                          'wellRefinement', true     , ...
% Refine
                          'wellGridFactor', 0.4      );
G = computeGeometry(G);       % Compute geometry
G = computeCellDimensions(G); % Compute cell dimensions
Warning: DistMesh did not converge in maximum number of iterations.

Sample rock properties

We assign rock properties in the PEBI grid cells by sampling from the fine grid using sampleFromBox

poro = sampleFromBox(G, reshape(rockRef.poro, GRef.cartDims));
perm = zeros(G.cells.num,G.griddim);
for i = 1:G.griddim
    perm(:,i) = sampleFromBox(G, reshape(rockRef.perm(:,i), GRef.cartDims));
end
rock = makeRock(G, perm, poro);

Set up schedule and initial state

W        = WRef;
schedule = scheduleRef;
x  = G.cells.centroids;
xwR = GRef.cells.centroids(vertcat(WRef.cells),1:2);
% Slick oneliner to find corresponding cells in the new grid
[~, c] = min(sum(bsxfun(@minus, reshape(xwR, [], 1 , G.griddim), ...
                             reshape(x , 1 , [], G.griddim)).^2,3), [], 2);
for i = 1:numel(W)
    W(i).cells = c(i);
end
xw = G.cells.centroids(vertcat(W.cells),1:2);
schedule.control.W = W;
state0 = initResSol(G, state0Ref.pressure(1), state0Ref.s(1,:));

Inspect the geological model

figure('Position', [0,0,800,410]), subplot(1,2,1)
Kref = convertTo(rockRef.perm(:,1),milli*darcy);
plotCellData(GRef,log10(Kref),'edgeAlpha',.1); axis tight
set(gca,'FontSize',12)
mrstColorbar(Kref,'South',true); cx = caxis();
hold on; plot(xwR(:,1),xwR(:,2),'.r','MarkerSize',18); hold off

subplot(1,2,2);
K = convertTo(rock.perm(:,1),milli*darcy);
plotCellData(G,log10(K),'edgeAlpha',.1); axis tight; caxis(cx);
set(gca,'FontSize',12)
mrstColorbar(K,'South',true); axis tight
hold on; plot(xw(:,1),xw(:,2),'.r','MarkerSize',18); hold off
_images/dgExampleIFS_01.png

Set base model

The base model is a two-phase oil-water model

model        = GenericBlackOilModel(G, rock, fluid, 'gas', false);
pmodel       = PressureModel(model); % Pressure model
makeSeqModel = @(tmodel) ...
    SequentialPressureTransportModel(pmodel, tmodel,'parentModel', model);

Simulate with dG(0)

tmodelDG0 = TransportModelDG(model, 'degree', 0);
modelDG0  = makeSeqModel(tmodelDG0);
[wsDG0, stDG0, repDG0] = simulateScheduleAD(state0, modelDG0, schedule);
Getting supports.
100 of 521...
200 of 521...
300 of 521...
400 of 521...
500 of 521...
Elapsed time is 0.422063 seconds.
Inverting systems...
...

Simulate with dG(1)

tmodelDG1 = TransportModelDG(model, 'degree', 1);
modelDG1  = makeSeqModel(tmodelDG1);
[wsDG1, stDG1, repDG1] = simulateScheduleAD(state0, modelDG1, schedule);
Getting supports.
100 of 521...
200 of 521...
300 of 521...
400 of 521...
500 of 521...
Elapsed time is 0.377109 seconds.
Inverting systems...
...

Visualize results

We plot the evolving saturation front computed using the two discuretizations as surface plots. To this end, we use plotSaturationDG, which supports surface plots of cell-wise discontinuous data using the patch function. Axis properties

setAxProps = @(ax) set(ax, 'Projection'        , 'Perspective', ...
                           'PlotBoxAspectRatio', [1,1,0.3]    , ...
                           'View'              , [-75,52]     , ...
                           'XLim'              , [0,l(1)]     , ...
                           'YLim'              , [0,l(2)]     , ...
                           'ZLim'              , [0,1]        , ...
                           'Box'               , 'on'         );
% For plotting wells
pw = @() plotWell(GRef, WRef, 'height', -1, 'color', 'k');
% Get coordinates for plotting (to be used in patch)
coords = getPlotCoordinates(G);
close all; figure('Position', [0,0,1000,500]);
subplot(1,2,1); % Plot dG(0)
[hs0, satDG0] = plotSaturationDG(tmodelDG0.discretization, stDG0{1}, ...
                                                         'coords', coords);
pw(); setAxProps(gca); caxis([0.2,0.8]); camlight; % Set axis properties
subplot(1,2,2); % Plot dG(1)
[h1, satDG1] = plotSaturationDG(tmodelDG1.discretization, stDG1{1}, ...
                                                         'coords', coords);
pw(); setAxProps(gca); caxis([0.2,0.8]); camlight; % Set axis properties

for i = 1:numel(stDG1)
    % Update dG(0) patch
    s0 = satDG0(stDG0{i});
    hs0.Vertices(:,3)   = s0;
    hs0.FaceVertexCData = s0;
    % Update dG(1) patch
    s1 = satDG1(stDG1{i});
    h1.Vertices(:,3)   = s1;
    h1.FaceVertexCData = s1;
    % Pause
    pause(0.2),
end
_images/dgExampleIFS_02.png

Plot the difference between the two solutions

figure
plotCellData(G,stDG1{end}.s(:,1)-stDG0{end}.s(:,1),'EdgeColor','none'); axis tight
view(-90,90), box on
colormap(interp1([0; 0.5; 1], [1, 0, 0; 1, 1, 1; 0, 0, 1], 0:0.01:1))
caxis([-.151 .151]);
colorbar('EastOutside');
_images/dgExampleIFS_03.png

Simulation of an inverted five-spot pattern with dG

Generated from dgExampleNess.m

In this example, we simulate water injection in an inverted five-spot pattern posed on a subset of SPE10 using dG(0) and dG(1) visualize the results

Add modules

mrstModule add dg spe10 ad-props ad-core ad-blackoil blackoil-sequential

Set up fine-scale model

We extract part of layer 51 of SPE10 2

[state0, imodel, schedule] = ...
    setupSPE10_AD('layers', 51, 'J',1:110,'make2D', true, 'T', 3*year, 'dt', 20*day);
G  = computeCellDimensions(imodel.G);
schedule.control.W(3).cells = 6550;
schedule.control.W(4).cells = 6578;
schedule.control.W(5).cells = 3811;

Set base model

model   = GenericBlackOilModel(G, imodel.rock, imodel.fluid, 'gas', false);
pmodel  = PressureModel(model); % Pressure model
tmodel0 = TransportModelDG(model, 'degree', 0);
model0  = SequentialPressureTransportModel(pmodel, tmodel0,'parentModel', model);
model0.AutoDiffBackend = DiagonalAutoDiffBackend('useMex', true);
[ws0, states0, reports0] = simulateScheduleAD(state0, model0, schedule);
Getting supports.
100 of 6600...
200 of 6600...
300 of 6600...
400 of 6600...
500 of 6600...
600 of 6600...
700 of 6600...
...

Simulate with dG(1)

tmodel1 = TransportModelDG(model, 'degree', 1, 'dsMaxAbsDiv', 2);
model1  = SequentialPressureTransportModel(pmodel, tmodel1,'parentModel', model);
model1.AutoDiffBackend = DiagonalAutoDiffBackend('useMex', true);
model1.transportNonLinearSolver = NonLinearSolver('useLinesearch', true);
[ws1, states1, reports1] = simulateScheduleAD(state0, model1, schedule);
Getting supports.
100 of 6600...
200 of 6600...
300 of 6600...
400 of 6600...
500 of 6600...
600 of 6600...
700 of 6600...
...

Make animation of the test case

To produce nice plots, we use three axes on top of each other. The first plots the permeability using a pink colormap, the second plots cells with water saturation exceeding 0.2 using a green-blue colormap, whereas the third only shows the well positions.

figure('Position', [0,0,900,400])
xw = model.G.cells.centroids(vertcat(schedule.control.W.cells),1:2);
bx1 = subplot(1,2,1);
plotCellData(G,log10(model.rock.perm(:,1)),'EdgeColor','none'); axis tight off
colormap(bx1,flipud(pink));
axlim = axis(bx1);
ax1 = axes('Position',bx1.Position);
colormap(ax1,flipud(winter));
cx1 = axes('Position',bx1.Position);
plot(xw(:,1),xw(:,2),'.r','MarkerSize',18);
axis(cx1,axlim); axis off

bx2 = subplot(1,2,2);
plotCellData(G,log10(model.rock.perm(:,1)),'EdgeColor','none'); axis tight off
colormap(bx2,flipud(pink));
ax2 = axes('Position',bx2.Position);
colormap(ax2,flipud(winter));
cx2 = axes('Position',bx2.Position);
plot(xw(:,1),xw(:,2),'.r','MarkerSize',18);
axis(cx2,axlim); axis off


[h1,h2]=deal([]);
for i=1:numel(states0)
    delete([h1 h2]);
    set(gcf,'CurrentAxes',ax1);
    h1 = plotCellData(G,states0{i}.s(:,1),states0{i}.s(:,1)>.205,'EdgeColor','none');
    axis(ax1,axlim), axis off
    set(gcf,'CurrentAxes',ax2);
    h2 = plotCellData(G,states1{i}.s(:,1),states1{i}.s(:,1)>.205,'EdgeColor','none');
    axis(ax2,axlim), axis off
    drawnow;
    pause(.2)
end
_images/dgExampleNess_01.png

Plot production curves

plotWellSols({ws0, ws1},reports0.ReservoirTime,...
    'datasetnames',{'dG(0)','dG(1)'}, 'field', 'qWs');
_images/dgExampleNess_02.png

Show Sparsity Pattern

Generated from dgShowSparsity.m

This example compares the sparsity patterns of dG(0)/SPU and dG(1) for a small PEBI grid and shows how one can use potential ordering to permute the systems to (block)triangular form.

mrstModule add upr dg ad-core ad-props ad-blackoil blackoil-sequential

Set up the model and compute a pressure solution

We construct a PEBI grid using a function from the UPR module, specify a quarter five-spot type well pattern, and set up an incompressible oil-water system with unit properties and slight rock compressibility. We then split the generic black-oil model into a pressure and a transport model (with dG) and use a standalone pressure solver to compute pressure. This is used as input to compute the linearized transport equations for the dG transport solver.

G        = pebiGrid(1/4, [1,1]);
G        = computeGeometry(G);
G        = computeCellDimensions(G);
rock     = makeRock(G, 1, 1);
W        = addWell([], G, rock, 1, 'type', 'rate', 'Radius', 1e-3, ...
                    'val',  1, 'compi', [1 0], 'Name', 'I');
W        = addWell(W, G, rock, G.cells.num, 'type', 'rate', 'Radius', 1e-3,...
                    'val', -1, 'compi', [1 0], 'Name', 'P');
fluid    = initSimpleADIFluid('phases', 'WO', 'mu', [1, 1], 'rho', [1, 1], 'cR', 1e-6/barsa);
schedule = simpleSchedule(1, 'W', W);

model    = GenericBlackOilModel(G, rock, fluid, 'gas', false);
pmodel   = PressureModel(model); % Pressure model
tmodel0  = TransportModelDG(model, 'degree', 0);
model0   = SequentialPressureTransportModel(pmodel, tmodel0,'parentModel', model);

forces   = schedule.control(1);
model0   = model0.validateModel(forces);
state0   = model0.validateState(initResSol(G, 1*barsa, [0, 1]));
state    = standaloneSolveAD(state0, model0, schedule.step.val(1), 'W', forces.W);
problem  = model0.transportModel.getEquations(state0, state, 1, forces);
Getting supports.
Elapsed time is 0.067339 seconds.
Inverting systems...
Solving timestep 1/1:          -> 1 Second
*** Simulation complete. Solved 1 control steps in 3 Seconds, 17 Milliseconds ***

Sparsity pattern for dG(0)/SPU

Show sparsity patter before and after applying potential ordering

figure('Position',[520 360 985 420]);
subplot(1,2,1)
A0 = problem.equations{1}.jac{1};
spy(A0);
set(gca,'XTick',[],'YTick',[]); xlabel('');
axes('Position',[.31 .555 .15 .35]);
plotGrid(G,'FaceColor',[1 1 .9]);
h=text(G.cells.centroids(:,1),G.cells.centroids(:,2),num2str((1:G.cells.num)'));
set(h,'FontSize',8,'HorizontalAlignment','center','VerticalAlignment','middle')
axis off

subplot(1,2,2)
[~,q] = sort(state.pressure,'descend');
spy(A0(q,q));
set(gca,'XTick',[],'YTick',[]); xlabel('');
axes('Position',[.75 .555 .15 .35]);
plotGrid(G,'FaceColor',[1 1 .9]);
h=text(G.cells.centroids(q,1),G.cells.centroids(q,2),num2str((1:G.cells.num)'));
set(h,'FontSize',8,'HorizontalAlignment','center','VerticalAlignment','middle')
axis off
_images/dgShowSparsity_01.png

Get the discretization matrix for the dG(1) scheme

As for the dG(0)/SPU scheme, we first construct a transport model and a corresponding sequential model, before we compute one pressure step and use the result as input to construct the linear transport system

tmodel1  = TransportModelDG(model, 'degree', 1);
model1   = SequentialPressureTransportModel(pmodel, tmodel1,'parentModel', model);

model1   = model1.validateModel(forces);
state0   = model1.validateState(initResSol(G, 1*barsa, [0, 1]));
state    = standaloneSolveAD(state0, model1, schedule.step.val(1), 'W', forces.W);
problem  = model1.transportModel.getEquations(state0, state, 1, forces);
Getting supports.
Elapsed time is 0.048801 seconds.
Inverting systems...
Solving timestep 1/1:          -> 1 Second
*** Simulation complete. Solved 1 control steps in 3 Seconds, 648 Milliseconds ***

Sparsity pattern for dG(1)

Show the sparsity system before and after potential ordering

figure('Position',[520 360 985 420]);
subplot(1,2,1)
A1 = problem.equations{1}.jac{1};
spy(A1,8);
set(gca,'XTick',[],'YTick',[]); xlabel('');
axes('Position',[.31 .555 .15 .35]);
plotGrid(G,'FaceColor',[1 1 .9]);
h=text(G.cells.centroids(:,1),G.cells.centroids(:,2),num2str((1:G.cells.num)'));
set(h,'FontSize',8,'HorizontalAlignment','center','VerticalAlignment','middle')
axis off

subplot(1,2,2)
qq = 3*ones(G.cells.num,1); qq([1 end])=1; qq=cumsum([0; qq]);
map = cell(G.cells.num,1);
for i=1:G.cells.num
    map{i} = (qq(i)+1:qq(i+1))';
end
p=vertcat(map{q});
spy(A1(p,p),8);
set(gca,'XTick',[],'YTick',[]); xlabel('');
axes('Position',[.75 .555 .15 .35]);
plotGrid(G,'FaceColor',[1 1 .9]);
h=text(G.cells.centroids(q,1),G.cells.centroids(q,2),num2str((1:G.cells.num)'));
set(h,'FontSize',8,'HorizontalAlignment','center','VerticalAlignment','middle')
axis off
_images/dgShowSparsity_02.png