﻿ Function solvefi.m

## MRST - MATLAB Reservoir Simulation Toolbox

Function solvefi.m

## Single step non-linear solver

Compute the next state of the system for given previous state state0 and time step dt.

The discretized residual equation are given by equation which is assembled in the function equationCompositional.

```function [state, convergence] = solvefi(state0, dt, bc, system, equation, param,  varargin)
```
```   opt = struct('verbose', false);
opt = merge_options(opt, varargin{:});

fluid = system.fluid;

meta = struct('converged'  , false                       , ...
'stopped'    , false                       , ...
'relax'      , system.nonlinear.relaxation , ...
'stagnate'   , false                       , ...
'iteration'  , 0                           , ...
'res_history', []                          );

timer = tic;

converged = false;
state = state0;

fprintf('%13s%-26s%-36s\n', '', 'CNV (oil, water)', 'MB (oil, water)');

equation = @(state) equation(state0, state, dt, bc, system);
```

```   while ~ (converged || meta.stopped),
```
```      % Save iteration number in meta info
meta.iteration = meta.iteration + 1;
```

Initial saturation solve

At each Newton step, we start by solving the flash equations and update the liquid saturation variable.

```      [C, p] = deal(state.C, state.pressure);
nc = numel(p);
init_sL = state.sL;
Cvec = cell2mat(C);
[cg, cl, cw, s, Cw] = flash_calculation(Cvec, p, system, init_sL);
state.sL = s(:, 2);
```

The residual equations for the whole system (pressure, total concentrations, liquid saturation) are assembled.

```      eqs = equation(state);
```

We call a standard linear solve to compute the Newton step dx.

```      dx = SolveEqsADI(eqs, []);
```

We update state.

```      state      = updateState(state, dx, system);
```

We compute the residual values by calling getResiduals. This function detects oscillation and stagnation

```      [meta, residuals] = getResiduals(meta, eqs, system, false);
```

We test for convergence. Here we use a very stringent test given by a max norm

```      fprintf('Newton iteration %d, max residual: %5g\n', meta.iteration, ...
max(abs(residuals)));
converged = all(residuals <= system.nonlinear.tol);

if meta.stagnate,
warning('newt:stagnate', 'Non-linear solver stagnated...')
end

if meta.iteration > system.nonlinear.maxIterations
meta.stopped = true;
end
```
```   end

if meta.stopped,
warning('newt:maxit', ...
['Non-linear solver did not converge, stopped ', ...
'by max iterations...']);
convergence = false;
else
convergence = true;
end

dispif(mrstVerbose, 'Completed %d iterations in %1.2f s\n', ...
meta.iteration, toc(timer));
```
```end
```

Update function

Given a Newton step, the state variables are updated. We cap the saturation variable.

```function state = updateState(state, dx, system)

nComp = system.nComp;
dp = dx{1};
for ic = 1 : nComp
dC{ic} = dx{ic + 1};
end
dsL = dx{nComp + 2};

step = 1;

state.pressure = state.pressure + step*dp;
for ic = 1 : nComp
state.C{ic} = max(0, state.C{ic} + step*dC{ic});
end
state.sL = min(1, max(0, state.sL + step*dsL));

end
```

Published October 8, 2014