Nonlinear residuals should not be mesh-dependent

[klevzoff]

Problem In single-phase flow solvers (FVM, HybridFVM, and by extension, single-phase poromechanics), the nonlinear residual is dependent on mesh resolution. E.g. for a particular problem (fixed-size cube with left-to-right pressure drop prescribed) we observe the following initial residuals:

mesh	R
35x35x35	5.32e+02
70x70x70	4.26e+03
140x140x140	3.40e+04
280x280x280	2.72e+05

As you can see, each one is exactly 8 times larger that the previous. This poses a significant challenge when designing mesh refinement studies, as nonlinear tolerance has to be adapted to mesh resolution.

Proposed solution We should make flow residuals independent of mesh resolution (as is already the case for mechanics and multiphase flow). Dependence comes from dividing the normalizer by the total cell count: https://github.com/GEOSX/GEOSX/blob/11cbe7d403087bd341f0986598269f6bb5d67bb4/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp#L161 which effectively ends up multiplying the residual. It't not clear to me what the purpose of that division is, at least in the context of normal porous media flow. We end up dividing a mesh-total quantity (L2-norm of real mass balance residual) by a cell-average quantity (average mass per cell), which doesn't make sense to me. Is it something that was added for hydrofracturing problems? What is the rationale? I recall us having a discussion about this residual evaluation a long time ago, but I can't find a written record of it.

My proposal is to remove globalResidualNorm[2] (cell count) altogether.

@rrsettgast @francoishamon @CusiniM @castelletto1 and anyone else, objections?

[CusiniM]

No, I do not recall a reason for this. I do remember that we added the m_fluxEstimate to help the HydrofractureSolver in the very first time-step after an element has been created but I cannot entirely remember the reason for multiplying by the number of elements. I think it may have to do with the fact that newly created elements may have some quantities equal to zero.

In other words, I have no objections.

[castelletto1]

Assuming the nonlinear residual is evaluated using the $L^{\infty}$-norm, normalizing each residual component by the corresponding cell mass ($\phi_n \rho_n V$) seems reasonable. Doing this for the $L^2$-norm doesn't make a lot of sense. Dividing by the total mass only (globalResidualNorm[1]), as suggested by @klevzoff, seems more appropriate. I'm not sure I understand the role played by m_fluxEstimate.

Additional comments:

• Irrespective of whether we think it is correct or not, why is the number of cells incremented by 1 (globalResidualNorm[2]+1)?
• Is m_fluxEstimate used anywhere? It seems to me that it is initialized using the default constructor (here) but never updated.

[klevzoff]

• Is m_fluxEstimate used anywhere? It seems to me that it is initialized using the default constructor (here) but never updated.

It's an input parameter used in hydrofracturing tests. Basically it's meant to be an estimate of the total fluid mass (flux*dt) injected in the first time step, to avoid dividing by zero fluid mass present in the fracture at the start of the simulation. I don't think it plays a significant role beyond the first few time steps?

Assuming the nonlinear residual is evaluated using the $L^\infty$-norm, normalizing each residual component by the corresponding cell mass ($\phi_n \rho_n V$) seems reasonable. Doing this for the $L^2$-norm doesn't make a lot of sense. Dividing by the total mass only (globalResidualNorm[1]), as suggested by @klevzoff, seems more appropriate.

We're not using $L^\infty$-norm unfortunately. But combined with cell-wise normalization, that would be my preferred choice actually, as it forces the solver to maintain stricter mass balance in areas of the reservoir with local refinement. It's more important for transport problems than pure flow though. It's likely to be somewhat problematic when the mesh contains degenerate cells with really low volumes, and I'm not sure about fractures either.

[rrsettgast]

• The intent of the normalization is to give some the residual norm some sort of indication of numerical precision. Normalizing the residual error in mass for a cell by the mass in a cell seemingly seems like it achieves this to me. Is this incorrect?
• What is driving the initial residual in the stated mesh convergence case? I suspect there is a source BC that is being dumped into a single cell of decreasing size as we refine? This would lead to sort of increases in initial residual shown.
• Why not use the L-inf norm? I don't remember why we choose the L2 over the L-inf.

[rrsettgast]

@klevzoff Is this resolved?

[klevzoff]

@rrsettgast this will be auto closed by @francoishamon's PR (which we need to review)