Function solvefi.m
tutorial main page computeComposition.m computeWaterComp.m equationCompositional.m flash_calculation.m getHenryCoef.m
getR.m initStateBravo.m omega_l.m quadraticRelPerm.m runBravo.m setNonlinearSolverParameters.m setupControls.m
setupGeometry.m setupSystem.m solvefi.m vaporPressure.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);

We start with the Newton iterations

   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, ...
      converged = all(residuals <= system.nonlinear.tol);

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

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

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

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

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};
   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});
   state.sL = min(1, max(0, state.sL + step*dsL));


Published October 8, 2014