Threshold initial oil saturation to 1e-6 minimum.

[atgeirr]

This is to avoid spurious oil saturation in the gas or water zones that puts us erroneously at the boundary of saturated/unsaturated behaviour, when we should be well away from that boundary.

This can occur when using endpoint scaling and minimum oil saturation is set to 1 - sgu + swl, and this becomes a very small number instead of zero (as it should have been).

This has follow-on consequences: being on the saturated/unsaturated boundary, the jacobian of the storage term may become singular, causing trouble for the "trueimpes" variant for calculating CPR weights which is currently the default.

[atgeirr]

jenkins build this please

[michal-toth]

Could you please point to the place where the follow-up calculation takes place?
It puzzles me that the commit sets to zero not just negative numbers but also small positive numbers. Are small positive values worse than 0? If yes, is it because we treat 0 differently (in CPR weights calculation)?

[atgeirr]

Could you please point to the place where the follow-up calculation takes place?

After equilibration we will have stored states in FlowProblemBlackoil::initialFluidStates_. This will be overdetermined, in the sense that it contains all saturations, as well as rs, rv etc. Making up some values, say we have sw = 0.1, sg = 0.899999999 and so = 1e-9, rv = 2e-7 (from RVVD or PDVD). That rv value is below the saturated rv value, rvsat = 7e-7.

The critical step is when we first construct primary variables by calling BlackoilPrimaryVariables::assignNaive(). That function has this logic:

        const bool oilPresent =
            fluidState.fluidSystem().phaseIsActive(oilPhaseIdx)
                ? fluidState.saturation(oilPhaseIdx) > 0.0
                : false;

The primary variables stored will therefore for this cell be (po, sw, sg) instead of (po, sw, rv). At that point the rv value we had is lost. Instead of as intended being well into the undersaturated state, we are almost exactly at the boundary. So small positive values are worse than zero, since that puts us at the place where the CPR trueimpes weights for the cell become the solution of a singular system. This is not because the weights treat a zero in any special way, but because of the change of variables that is wrong: rv should be the variable here, not sg.

Later in the simulation, assignNaive() is not used, instead adaptPrimaryVariables() is used to decide what the correct primary variable should be, which has different logic. So the "exact zero or not" only matters in the initialization. At least I think so! Also, ending up exactly at the saturated/unsaturated boundary is not likely to occur, unless extremely unlucky. In it happens during simulation the worst that happens is failure to converge the linear system and a timestep cut. Whereas with the initial state on that sat/unsat boundary for some cell deep into the gas zone, the oil is very likely immobile, and may cause the same trouble over and over again, at every step throughout the simulation.

I also tried modifying the logic in assignNaive() as an alternative, but I prefer the current solution.

[michal-toth]

Thank you for the explanation. If I understand it correctly, the key problem is that the assignNaive does not properly handle the saturation close to zero like it is done in adaptPrimaryVariables.

I also tried modifying the logic in assignNaive() as an alternative, but I prefer the current solution.

My preference is the opposite - I would rather change how variables are interpreted (and put threshold in assignNaive) than change the variables. I also do not like breaking the condition that the sum of saturations is 1. But it does not matter which approach is used if only the initialization is affected, I am just worried whether the changed saturation is used elsewhere.

Is deriveOilSat or assignNaive used when the computation is restarted from a report step? If yes, then the values at the boundary of saturated/unsaturated region can be reached naturally.

[atgeirr]

I also do not like breaking the condition that the sum of saturations is 1

We are not actually! After assignNaive() with the case I gave, the sg will not be a primary variable and its value will be implicitly 1-sw (instead of 0.899999999). The so is never a primary variable anyway, and its value will implicitly be 0 since we are in the undersaturated gas state.

If yes, then the values at the boundary of saturated/unsaturated region can be reached naturally.

I did not understand what you mean here, please explain.

[michal-toth]

The so is never a primary variable anyway, and its value will implicitly be 0 since we are in the undersaturated gas state.

That is a relief. Changes in assignNaive and deriveOilSat will be then equivalent.

the values at the boundary of saturated/unsaturated region can be reached naturally.

I was thinking about a problem where a domain has a region with zero oil saturation and the rest has positive saturation. The boundary of the region can move but not disappear (unless it shrinks to nothing..), and there will be cells with very low saturation.

[atgeirr]

This is undergoing some more testing to ensure it does not have any adverse effects. I hope to get it merged in time for the release, but it should not be merged until the testing is done.

[bska]

Tiny nit: You probably meant to say sgu+swl here rather than sgu+sgl. If we're going to add special case treatment for this I think we'd better be specific about which case we're treating.

[atgeirr]

Thanks, will fix!

[atgeirr]

jenkins build this please

[atgeirr]

Testing on this was a bit less conclusive than I had hoped, although I still believe this is a move in the right direction. I will therefore continue to investigate it a bit, and also look at deciding saturation status by comparing rv to saturated rv values instead of looking at oil saturation. (And similar with rs.)