Problem with Simulator workflow example

[Dāvis Bīgestāns]

Hi MRST community,

I am trying to get running a modified SimulatorWorkFlowExample but I keep getting an error of incomplete simulation. Could You please help me spot the mistake?

Thanks You!

=========================================================

mrstModule add ad-core ad-blackoil ad-props mrst-gui

close all;

cdims = [29,15,2];

pdims = [3000,3000,81.25];

G = cartGrid(cdims, pdims);

G = computeGeometry(G);

E_permx = 1:1:435;

R_permx = 1:1:435;

E_poro = 1:1:435;

R_poro = 1:1:435;

for i=1:435

E_poro(i) = 0.270;

R_poro(i) = 0.250;

E_permx(i) = 3507;

R_permx(i) = 230;

end

E_permx = reshape(E_permx, [435,1]);

E_poro = reshape(E_poro, [435,1]);

R_permx = reshape(R_permx, [435,1]);

R_poro = reshape(R_poro, [435,1]);

permx = [E_permx, R_permx];

poro = [E_poro, R_poro];

rock = makeRock(G, reshape(permx, [870,1]), reshape(poro, [870,1]));

simTime = 10*year;

nstep   = 25;

refine  = 5;

pv = poreVolume(G,rock);

%Wells ===================

%cellInx gives the position of the well. Vector; each integer refers to a place in grid (must be part of G.cells.num)

%cellInx is follows x axis. if grid 29x15 (nx * ny) then to get a well into 2 row must use 30,31,32 because first row ends @ 29

%Type = 'bhp' or 'rate'

%Val = if bhp then Pa (unit) if 'rate' then (m^3/s as unit)

W = [];

U_prate = 9.073*10^(-3); %m^3/s... Upper Brent prod rate

L_prate = 0.013; %m^3/s... Lower Brent prod rate

injRate = 1*sum(pv)/simTime;

W = addWell(W, G, rock, [427,428,429,430,431,432,433,434,435], 'Name', 'P1', 'comp_i', [1,0], 'Val', L_prate, 'Type', 'rate');

W = addWell(W, G, rock, [2,3,4,5,6,7,8,9,10], 'Name', 'I1', 'comp_i', [1 0], 'Val', injRate, 'Type', 'rate');

plotCellData(G, rock.perm, 'FaceAlpha', 0.5, 'EdgeAlpha', 0.3, 'EdgeColor', 'k'); %setting transparent grid to see the wells

plotWell(G,W); view(3) %plotWell must come after plotCellData

startSteps = repmat((simTime/(nstep + 1))/refine, refine, 1);

restSteps =  repmat(simTime/(nstep + 1), nstep, 1);

timesteps = [startSteps; restSteps];

% Set up the schedule containing both the wells and the timesteps

schedule = simpleSchedule(timesteps, 'W', W);

%Setting up the simulation model =======================

% Three-phase template model

%arguements in order: water, oil, gas

%table of 3 or 4 column array where 1st col = Sw; 2nd col = k,rw; 3rd col = k,ro

table = [0.0733, 0.0000, 1;

0.1057, 0.0002, 0.8752;

0.1440, 0.0064, 0.7175;

0.1723, 0.0086, 0.6172;

0.2027, 0.0127, 0.5086;

0.2148, 0.0168, 0.4636;

0.3081, 0.0375, 0.2875;

0.3853, 0.0561, 0.2057;

0.4321, 0.0726, 0.1751;

0.4850, 0.0953, 0.1403;

0.5318, 0.1138, 0.1100;

0.5542, 0.1200, 0.1036;

0.5745, 0.1344, 0.0913;

0.5945, 0.1444, 0.0813;

0.6145, 0.1544, 0.0713;

0.6345, 0.1644, 0.0613;

0.6545, 0.1744, 0.0513;

0.6745, 0.1844, 0.0413;

0.6945, 0.1944, 0.0313;

0.7145, 0.2044, 0.0213;

0.7345, 0.2144, 0.0113;

0.7545, 0.2244, 0];

fluid = initSimpleADIFluid('mu',[0.36, 0.34,0]*centi*poise, 'rho', [1025, 833,0]*kilogram/meter^3, 'b', [1, 4.116,0], 'c', [3*10^(-6), 0.001,0], 'cR', 4.12*10^(-6),'pRef', 38300032.02);

kr = tabulatedSatFunc(table)

fluid.relperm = kr

gravity reset on

model = TwoPhaseOilWaterModel(G, rock, fluid);

sW = ones(G.cells.num, 1);

for i=1:870;

if sW(i) < 436;

sW(i) = 0.11;

else sW(i) = 1;

end

end

sat = [sW, 1-sW];

% Compute initial pressure

p_res = [];

for i=1:870;

if i<435;

p_res(i) = 300*barsa;

else p_res(i) = 305*barsa;

end

end

p_res = reshape(p_res, [870,1]); 

state0 = initResSol(G, p_res, sat);

clf

plotCellData(G, state0.s(:,1))

plotWell(G,W)

view(50, 50), axis tight

fn = getPlotAfterStep(state0, model, schedule,'view',[50 50], ...

                     'field','s:1','wells',W);

[wellSols, states, report] = ...

   simulateScheduleAD(state0, model, schedule,'afterStepFn',fn);

[Olav Møyner]

Hi,

I have not run your script, but you may want to set the producer rate to a negative value (otherwise it will inject).

Also, the relperm format you are using will not have any effect since it is made for the incomp solvers. Have a look at assignSWOF under ad-props.

Best regards

Olav

Olav

---

From: sinte...@googlegroups.com <sinte...@googlegroups.com> on behalf of Dāvis Bīgestāns <davisbi...@gmail.com>
Sent: Friday, March 5, 2021 5:35:41 PM
To: MRST-users: The Matlab Reservoir Simulation Toolbox User Group <sinte...@googlegroups.com>
Subject: [MRST Users] Problem with Simulator workflow example

 

--
You received this message because you are subscribed to the Google Groups "MRST-users: The Matlab Reservoir Simulation Toolbox User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sintef-mrst...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sintef-mrst/e25ae8d0-b8ab-42b9-8075-1201cd2ab752n%40googlegroups.com.

[Dāvis Bīgestāns]

Hi Olav,

Thank You very much for suggesting new ideas!

I looked at the assignSWOF function but I hardly understand how to implement it since it appears that it has to be linked with external Eclipse file. Could You please show me how to use it in the given example?

Thank You in advance!

Davis

[Dāvis Bīgestāns]

Hi Olav,

MRST truly looks like an interesting toolbox - very helpful to learn reservoir simulation on own. I looked at previous threads and saw that you had suggested to use the concept of assignSWOF in this way for relative permeability to water:

fluid.krW = @(S) interpTable(table(:, 1), table(:, 2), S)

%where the first column in table is water saturation, second column krW and third - krO

Serveh from the perviously mentioned thread wanted to apply the fluid.krW concept to fluid.krO as well but you suggested to look at assignSWOF and posted this:

 krW{i} = @(sw) reg.interp1d(SW, swof(:, 2), sw);

 krOW{i} = @(so) reg.interp1d(SW, swof(:, 3), 1-so);

I tried to apply this concept on this kind of a way:

krW{i} = @(sw) reg.interp1d(S_table, table(:, 2), sw);

krOW{i} = @(so) reg.interp1d([1-S_table], table(:, 3), 1-so);

but once again I got the same error

In terms of the errors I keep getting something really similar to Serveh:

"

Solver did not converge in 0 iterations for timestep of length 21 Hours, 258 Seconds, 57 Milliseconds. Cutting timestep.

Warning: Matrix is singular to working precision. 

> In BackslashSolverAD/solveLinearSystem (line 21)

In LinearSolverAD/solveLinearProblem (line 271)

In PhysicalModel/stepFunction (line 690)

In ReservoirModel/stepFunction (line 300)

In NonLinearSolver/solveMinistep (line 352)

In NonLinearSolver/solveTimestep (line 203)

In simulateScheduleAD (line 270)

 

Error using NonLinearSolver/solveTimestep (line 272)

Did not find a solution: Model step resulted in failure state. Reason: Linear solver produced non-finite values.

Error in simulateScheduleAD (line 270)

            [state, report] = solver.solveTimestep(state0, dt, model,…

"

Thank You very much Olav!

Best Regards,

Davis

On Friday, March 5, 2021 at 5:11:22 PM UTC Olav Møyner wrote:

[Olav Møyner]

Hi Davis,

Here is a modified script that runs. The main issue was units. MRST uses strict SI units for all quantities, and it is very important that you do not mix in e.g. field units when writing scripts because the missing conversion factors are often very large or small. I added my best guesses for units to the permeability (milli*darcy), compressibility factors (1/barsa) and switched the producer to a bhp control. I also added the rel. perms. from the table. Then, the script runs.

Best regards,

Olav

---

From: sinte...@googlegroups.com <sinte...@googlegroups.com> on behalf of Dāvis Bīgestāns <davisbi...@gmail.com>

Sent: Friday, March 5, 2021 23:47

To: MRST-users: The Matlab Reservoir Simulation Toolbox User Group <sinte...@googlegroups.com>

Subject: Re: [MRST Users] Problem with Simulator workflow example

 

To view this discussion on the web visit https://groups.google.com/d/msgid/sintef-mrst/7a53a3eb-8719-4257-86e6-3eb84c24b0e5n%40googlegroups.com.

[Dāvis Bīgestāns]

Hi Olav,

Thank You very much for the help! 

I know that MRST has strict unit specifications (SI), but I forgot to specify them for the compressibility factors. Thanks for pointing that out -  You have really helped me!

Best Regards,

Davis Bigestans