Changing state's step PVT props

[Angelo Kennedy]

Hello, 

I have this situation where I need to modify phase PVT properties, for example, inside the FlowDiscretization function; however, since the current property values, such as in density and viscosity, are received with respective derivatives (ADI objects), I end up with nonconvergence when I simply replace the said properties in the state.PVTProps.

So, could you clarify to me if: (a) modifying and replacing specific properties in PVTProps handlestruct is a good approach, or (b) the alternative I use in assigning the new values inside respective density, viscosity, and shrinkage factors new functions is the better option. 

The main issue here is how you could modify PVTProps for a given timestep. I will provide the new CO2FlowDiscretization function below for reference.

Thank you.

Best regards,

Angelo

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

classdef CO2FlowDiscretizationCW < FlowDiscretization
   methods
    function [acc, v, names, types, statex] = componentConservationEquations(fd, model, state, state0, dt)
                % Compute discretized conservation equations in the interior of the domain.
                % REQUIRED PARAMETERS:
                %   model  - ReservoirModel derived class
                %   state  - State corresponding to the given model at the end
                %            of the current time-integration interval. f
                %            As a new update, the state's PVT props handle
                %            struct is modified upon the solution of the
                %            geochemical run.
                %   state0 - State at previous time-step.
                %   dt     - Time-step length.
                % RETURNS:
                %   acc    - Cell array of component accumulation terms
                %            discretized with a first order finite-difference
                %            scheme.
                %   v      - Cell array of component fluxes in the interior of
                %            the domain.
                %   names  - The names of each equation (corresponds to
                %            component names.
                %   types  - Types of each equation (defaults to 'cell', to
                %            indicate cell-wise equations on the whole domain)
                ncomp = model.getNumberOfComponents();
                [acc, types] = deal(cell(1, ncomp));
                names = model.getComponentNames();
                [types{:}] = deal('cell');
               
                % CSM Modification HU, June 2025...
                if state0.activateGeochem
                    state.passChemD = true;
                    state = curStateGC(state, model);
                    % Previous state geochemistry (just the very first
                    % state0)
                    if ~isfield(state0, 'PVTProps')
                        state0.passChemD = true;
                        state0.PVTProps  = state.PVTProps;
                        state0           = prevStateGC(state0,model);
                    end

                    % Check for any nan/inf etc............................
                    gcFields = {'rhoO_x', 'rhoG_x', 'bO_x', 'bG_x', ...
                                'muO_x',  'muG_x',  'rs_x'}; % Fields to check
                    [~] = gcfPhysicalCheck(gcFields, state);
                   
                                        % New adjustment trial
                    [~, ~] = model.getProps(state, 'Density', 'Viscosity');
                    state.PVTProps = value(state.PVTProps);
                   
                    state.PVTProps.Density = [state.rhoO_x, state.rhoG_x];             % Density
                    state.PVTProps.ShrinkageFactors = [state.bO_x, state.bG_x];        % SF
                    state.PVTProps.Viscosity = [state.muO_x, state.muG_x];             % Viscosity
                   
                    mass = model.getProps(state, 'ComponentTotalMass');
                    mass0 = model.getProps(state0, 'ComponentTotalMass');  
                    flowState = fd.buildFlowState(model, state, state0, dt);
                else
                    mass = model.getProps(state, 'ComponentTotalMass');
                    mass0 = model.getProps(state0, 'ComponentTotalMass');  
                    flowState = fd.buildFlowState(model, state, state0, dt);
                end
                statex = state;
                %...................End of modif...........................

                v = model.getProps(flowState, 'ComponentTotalFlux');
                if ~iscell(mass0)
                    % May have been reduced to a double array
                    mass0 = {mass0};
                end
                for c = 1:ncomp
                    acc{c} = (mass{c} - mass0{c})./dt;
                end
    end
   end
end

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

[Angelo Kennedy]

[Addition context]

So, the code above is an example of the approach (a) aforementioned, where I try replacing a property, e.g., shrinkage factor with the line: "state.PVTProps.ShrinkageFactors = [state.bO_x, state.bG_x];        % SF". And once all properties are replaced, we insert back the new state.PVTProps.

For approach (b), I replace the values once they are computed as follows in their functions, e.g., shrinkage factor: "b = state.bO_x; % replace with the modified SF value."  Full modified shrinkage factor code: 

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

classdef CO2BlackOilShrinkageFactors < ShrinkageFactors

% Shrinkage factors for black-oil

properties

useSaturatedFlag = true;

disgas = false;

vapoil = false;

end

methods

function b = CO2BlackOilShrinkageFactors(model, varargin)

b@ShrinkageFactors(model, varargin{:});

assert(isa(model, 'ThreePhaseBlackOilModel'), ...

['Model must be derived from the black-oil model. Did you want', ...

' the regular ShrinkageFactors class instead?'])

b.disgas = model.disgas;

b.vapoil = model.vapoil;

if b.disgas

b = b.dependsOn({'rs'}, 'state');

end

if b.vapoil

b = b.dependsOn({'rv'}, 'state');

end

if b.useSaturatedFlag

b = b.dependsOn({'s'}, 'state');

end

end

function b = evaluatePhaseShrinkageFactor(prop, model, state, name, p)

if prop.disgas && strcmp(name, 'O')

% Oileic phase with dissolved gas component

rs = model.getProp(state, 'rs');

if prop.useSaturatedFlag

sG = model.getProp(state, 'sg');

flag = sG > 0;

else

flag = false(numelValue(rs), 1);

end

b = prop.evaluateFluid(model, 'bO', p, rs, flag);

elseif prop.vapoil && strcmp(name, 'G')

% Gaseous phase with vaporized oil component

rv = model.getProp(state, 'rv');

if prop.useSaturatedFlag

sO = model.getProp(state, 'so');

flag = sO > 0;

else

flag = false(numelValue(rv), 1);

end

b = prop.evaluateFluid(model, 'bG', p, rv, flag);

else

% Can use base class directly

b = evaluatePhaseShrinkageFactor@ShrinkageFactors(prop, model, state, name, p);

end

% CSM Modification. HU., May, 2025...................

fluid = model.fluid;

if isfield(fluid, 'cox') && fluid.cox

if state.gcPVTProps && state.activateGeochem

stateGC = solverGeochem_CW2(statestrt, model);

state = stateGC;

end

if strcmp(name, 'O')

b = state.bO_x; % replace with the modified SF value

end

if strcmp(name, 'G')

b = b .* state.bG_x;

end

end

%..............................................................

end

end

end

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

[Olav Møyner]

Hi,

 

Apologies for the delay in response – July is usually holiday here in Norway.

 

The approach you have set up seems correct (i.e. inherit the class and make sure that the dimensions get correct). You will double evaluate some things, but that is not so important for getting things to work. The main issue you are running into is that you need to get ADI types out from solverGeochem_CW2. The ADI types wrap both the values (which are probably already what you want) and the derivatives, which need to be implemented. This can get fairly complicated, but you can see examples of how this is done in the ADI class (for primitive operations like multiplication and addition) and EquationOfStateModel in the compositional model (for some manual tricks for implicit differentiation). The free second MRST textbook also goes into some detail on this (how evaluations are structured, and how to extended them), as well as later parts of the first textbook (for the AD concept in general).

 

A first question may be if you have the option to generate derivatives in your solverGeochem (i.e. shrinkage factors or densities with respect to pressure, and any other values that are used). If so, the problem is then how to most easily get the values into the right AD format.

 

Best regards,

Olav

 

--
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 visit https://groups.google.com/d/msgid/sintef-mrst/85551f9e-efac-4a80-a6a1-abac4120c2b6n%40googlegroups.com.

[Angelo Kennedy]

Hi Olav,

 

Thank you for the helpful suggestions; never mind the delay.

 

I attempted to generate derivatives for the new PVT values in solverGeochem, as you suggested. In the function (assigngGeochemPVT), I generate derivatives only as a function of the ADI phase pressure (compute slopes and multiply with Jacobian: you may inspect that in the helper function makeADfromValSlope).

As usual, that’s succeeded by solving for the masses in the discretization function:

[state = assigngGeochemPVT(state, model);                      

mass = model.getProps(state, 'ComponentTotalMass');

mass0 = model.getProps(state0, 'ComponentTotalMass');  

flowState = fd.buildFlowState(model, state, state0, dt);].

 

 

Now, again, I am running into the error:

 

“Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  8.066229e-31.
> In BackslashSolverAD/solveLinearSystem (line 21)
In LinearSolverAD/solveLinearProblem (line 361)
In PhysicalModel/stepFunction (line 739)
In ReservoirModel/stepFunction (line 307)
In NonLinearSolver/solveMinistep (line 374)
In NonLinearSolver/solveTimestep (line 210)
In simulateScheduleAD (line 295)
In seq5_oneQgrid (line 281)”.

 

Thus, besides these fixes, I also tried inserting just one PVT parameter, such as density, and I still ran into long iterations. So, I really can’t tell if the issue is with the new PVT structure or values. I will note that not assigning any of the new PVTs, as shown in the screenshot, it just runs. Additionally, examining the derivatives, I noticed a difference in my oleic phase shrinkage factor with the normal, suggesting an inconsistency/inaccuracy in my derivative approach (Jacobians).

I would appreciate any comments on the assignment of parameters (generation of derivatives) in the function below or otherwise.

 

Thank you.

 

Best regards,

Angelo

----------------------------

function state = assigngGeochemPVT(state, model)

% Build ADI PVT properties from geochem arrays in state.*_x

% Uses phase pressures pO/pG when available; falls back to total pressure.

 

%% 0) Phase pressures (ADI) and numeric copies

p  = [];  pO = [];  pG = [];

PO = [];  PG = [];

 

if nargin >= 2 && ~isempty(model)

    % Try to get per-phase pressures

    try

        pp = model.getProps(state, 'PhasePressures');  % 1xNph cell of ADI

    catch

        pp = [];

    end

    if ~isempty(pp)

        % Phase indices (brine/"O" and gas/"G")

        try

            iO = model.getPhaseIndex('O');

            iG = model.getPhaseIndex('G');

        catch

            % Fallback: assume [O,G] ordering

            iO = 1; iG = 2;

        end

        pO = pp{iO};      PO = value(pO);

        pG = pp{iG};      PG = value(pG);

    else

        % Fallback: single pressure

        p  = model.getProp(state, 'pressure');

        PO = value(p);    pO = p;                 % use same for both

        PG = PO;          pG = p;

    end

else

    % No model passed in: use state.pressure

    p  = state.pressure;

    PO = value(p);       pO = p;

    PG = PO;             pG = p;

end

 

N = numel(PO);

 

%% 1) Pack geochem vectors (cellwise doubles)

names = {'bO_array','bG_array', ...

         'muO_array','muG_array', ...

         'rhoO_array','rhoG_array', ...

         'rs_array'};

 

vals  = {state.bO_x(:), state.bG_x(:), ...

         state.muO_x(:), state.muG_x(:), ...

         state.rhoO_x(:), state.rhoG_x(:), ...

         state.rs_x(:)};

 

pvt_val = cell2struct(vals, names, 2);

 

% Basic size sanity

for k = 1:numel(names)

    assert(numel(pvt_val.(names{k})) == N, ...

        'assigngGeochemPVT: %s must be length %d (got %d).', ...

        names{k}, N, numel(pvt_val.(names{k})));

end

 

%% 2) Slopes w.r.t. the matching phase pressure

% Brine/"O" props use pO; Gas/"G" props use pG; Rs uses pO.

d_bO_dPO   = dvaldP(pvt_val.bO_array  , PO);

d_bG_dPG   = dvaldP(pvt_val.bG_array  , PG);

d_muO_dPO  = dvaldP(pvt_val.muO_array , PO);

d_muG_dPG  = dvaldP(pvt_val.muG_array , PG);

d_rhoO_dPO = dvaldP(pvt_val.rhoO_array, PO);

d_rhoG_dPG = dvaldP(pvt_val.rhoG_array, PG);

d_rs_dPO   = dvaldP(pvt_val.rs_array  , PO);

 

%% 3) Build ADI variables using the *same* phase pressure blocks

bO   = makeADfromValSlope(pO, pvt_val.bO_array  , d_bO_dPO);

bG   = makeADfromValSlope(pG, pvt_val.bG_array  , d_bG_dPG);

 

muO  = makeADfromValSlope(pO, pvt_val.muO_array , d_muO_dPO);

muG  = makeADfromValSlope(pG, pvt_val.muG_array , d_muG_dPG);

 

rhoO = makeADfromValSlope(pO, pvt_val.rhoO_array, d_rhoO_dPO);

rhoG = makeADfromValSlope(pG, pvt_val.rhoG_array, d_rhoG_dPG);

 

rs   = makeADfromValSlope(pO, pvt_val.rs_array  , d_rs_dPO);   %#ok<NASGU>

% ^ keep as local unless you explicitly want to override state.rs (often

%   rs also depends on other primaries, so replacing it can change Jacobian

%   block structure). We persist the double copy below as before.

 

%% 4) Store in state (mirror MRST PVTProps layout you actually use)

% state.PVTProps.ShrinkageFactors = {bO, bG};

% state.PVTProps.Viscosity        = {muO, muG};

state.PVTProps.Density          = {rhoO, rhoG};

% state.rs                        = rs;           % ADI rs for downstream use

 

% Keep doubles as well if you want downstream diagnostics

state.rs_x = pvt_val.rs_array;

end

 

% ===== helpers ============================================================

 

function d = dvaldP(val, P)

% Safe ∂val/∂P for possibly non-uniform / constant P

if numel(unique(P)) <= 1

    d = zeros(numel(val),1);

else

    d = gradient(val(:), P(:));

end

d = sanitizeSlope(d);

end