Function solvefi.m
Single step non-linear solverCompute 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, ... 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 |