Please help.

I'm having trouble getting my Matlab code to run:

clc; clear all; close all;

%%Construct simple Cartesian test case nx = 40; ny = 40; nz = 1;

G = meshgrid([nx ny nz]); G = computeGeometry(G); rock.perm = repmat(100*milli*darcy, [G.cells.num, 1]); rock.poro = repmat(0.3 , [G.cells.num, 1]);

x = linspace(0, 1, 11) .'; y = linspace(1, 0, 11) .'; %%Define fluid and capillary pressure curve %We define the relative permeability and the capillary pressure in form of tables, and let the relative permeability curves be quadratic and the capillary function linear. The strength of the capillary pressure is decided by cap_scale. The capillary pressure is defined in the non-wetting phase, p_c = p_nw - p_w.

pc_form = 'nonwetting'; cap_scale = 10; [kr, pc] = tabulatedSatFunc([x, x.^2, y.^2, y.*cap_scale*barsa]); %%Define constant properties for viscosity and density

props = constantProperties([ 1, 10] .* centi*poise, ... [1000, 700] .* kilogram/meter^3); %%Here we put together a valid fluid object from the above defined functions. To read more about the fluid structure write help fluid_structure in MRST. First make a fluid without capillary pressure

fluid = struct('properties', props , ... 'saturation', @(x, varargin) x.s , ... 'relperm' , kr); %%Then make another fluid object identical to the one above except for the capillary pressure term 'pc'.

fluid_pc = struct('properties', props , ... 'saturation', @(x, varargin) x.s , ... 'relperm' , kr , ... 'pc' , @(x, varargin) pc(x.s)); %%Plot the pc-curve %Make a dummy state/solution structure to plot the pc curve since 'fluid.pc' demands state as an input

xDummy = initState(G, [], [0, 1]); xDummy.s = linspace(0, 1, numel(xDummy.s))'; ... pc = convertTo(fluid_pc.pc(xDummy), barsa);

clf plot(xDummy.s, pc); xlabel('s_w'); ylabel('pc [bar]'); title('Capillary pressure curve')

%%Set wells rate = 0.5*meter^3/day; bhp = 1*barsa;

W = verticalWell([], G, rock, 1, 1, 1:nz, ... 'Type', 'rate', 'Val', rate, ... 'Radius', .1, 'Name', 'I', 'Comp_i', [1 0]); W = verticalWell(W, G, rock, nx, ny, 1:nz, ... 'Type','bhp', 'Val', bhp, ... 'Radius', .1, 'Dir', 'x', 'Name', 'P', 'Comp_i', [0 1]); %%Set up solution structures and assemble linear system rSol = initState(G, W, 0, [0.2, 0.8]); rSol_pc = initState(G, W, 0, [0.2, 0.8]);

gravity off verbose = false;

S = computeMimeticIP(G, rock, 'Verbose', verbose,'InnerProduct','ip_tpf'); %%Set up pressure and transport solvers %This example uses an implicit transport solver, an explicit solver can be used if the time step restriction for the parabolic term is less than for the hyperbolic term. This is the case if 'cap_scale' is small. We let 'fluid' be a parameter in 'psolve' and 'tsolve' so that we can use the solvers for simulation both with and without capillary pressure by supplying different fluid objects. For this case we use the verbose = false for the transport solver. If more information about the convergence of the method is required; use verbose = true.

psolve = @(state, fluid) solveIncompFlow(state, G, S, fluid, 'wells', W); tsolve = @(state, dT, fluid) implicitTransport(state, G, dT, rock, ... fluid, 'wells', W, ... 'verbose', verbose); %%Alternatively we could have defined an explicit transport solver by

% tsolve = @(state, dT) explicitTransport(state, G, dT, rock, fluid, ... % 'wells', W, 'verbose', verbose); %%Solve initial pressure in reservoir %Observe that we supply different fluid objects for the two solutions, one with capillary pressure and one without.

rSol = psolve(rSol, fluid); rSol_pc = psolve(rSol_pc, fluid_pc); %%Transport loop %We solve the two-phase system using a sequential splitting in which the pressure and fluxes are computed by solving the flow equation and then held fixed as the saturation is advanced according to the transport equation.

T = 300*day(); dT = T/15; dTplot = 100*day(); % plot only every 100th day N = fix(T/dTplot); pv = poreVolume(G,rock); Start the main loop t = 0; plotNo = 1; h1 = 'No pc - '; H2 = 'Linear pc - '; e = []; p_org = []; p_pc = []; figure;

while t

% Check for inconsistent saturations s = [rSol.s(:,1); rSol_pc.s(:,1)]; assert(max(s) -eps);

% Update solution of pressure equation. rSol = psolve(rSol, fluid); rSol_pc = psolve(rSol_pc, fluid_pc);

% Measure water saturation in production cells in saturation e = [e; sum(abs(rSol.s(:,1) - rSol_pc.s(:,1)).*pv)/sum(pv)]; %#ok p_org = [p_org; rSol.s(W(2).cells,1)' ]; %#ok p_pc = [p_pc; rSol_pc.s(W(2).cells,1)']; %#ok

% Increase time and continue if we do not want to plot saturations t = t + dT; if ( t

% Plot saturation heading = [num2str(convertTo(t,day)), ' days']; r = 0.01; subplot('position',[(plotNo-1)/N+r, 0.50, 1/N-2*r, 0.48]), cla plotCellData(G, rSol.s(:,1)); caxis([0 1]), view(60,50), axis equal off, title([h1 heading])

subplot('position',[(plotNo-1)/N+r, 0.02, 1/N-2*r, 0.48]), cla plotCellData(G, rSol_pc.s(:,1)); caxis([0 1]), view(60,50), axis equal off, title([H2 heading])

plotNo = plotNo+1;

end

%%Plot water breakthrough at heel %As we clearly see from the plots in the figure, the simulation with capillary pressure has much more diffusion than the simulation without capillary pressure. This is confirmed by the water breakthrough curve.

clf n = numel(p_org(:,1)); plot(1:n,p_org(:,1),'-o',1:n,p_pc(:,1),'--*') legend('No capillary pressure','Linear capillary pressure','Location','Best'); title('Water breakthrough at heel');

Output graphs should look like this:

Plot the pc-curve Capilary prescue 10 91 0 01 02 04 05 OT 03 09 1 Transport loop Nope - 100 days No pc-200 days No pt - 300 days Linear pc - 100 days Linear pc-200 days Linear pe - 300 days Plot water breakthrough at heel 0.32 Water breakthrough at heel No capillary pressure Linear capillary pressure 0.31 0.28 0.26 0.24 0.22 0.2 0 10 15 Plot the pc-curve Capilary prescue 10 91 0 01 02 04 05 OT 03 09 1 Transport loop Nope - 100 days No pc-200 days No pt - 300 days Linear pc - 100 days Linear pc-200 days Linear pe - 300 days Plot water breakthrough at heel 0.32 Water breakthrough at heel No capillary pressure Linear capillary pressure 0.31 0.28 0.26 0.24 0.22 0.2 0 10 15