fix: Add HU simplified 2-phase version, refactor IHU implementation and fix bugs where possible, change phase/component flux logic

[paveltomin]

1. Another implementation of for hybrid upwinding, limited to 2 phases, tries to do total mobility treatment closer to published schemes
2. Refactored IHU implementation, added few helper functions to minimize code duplication, fixed/attempted to fix the following bugs:

• densMean * dGravD_dP should be outside of for( localIndex j = 0; j < numFluxSupportPoints; ++j ) subloop (not very important since dTrans_dPres is often zero), similar fix in capillary flux
• 'totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp]' and derivatives should not be in the for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke ) subloop - but i could not change this one without getting bad convergence for barrier test case - results with wrong treatment are definitely off but i am not sure what to do with that, @jafranc any thoughts?
• derivatives for total mobility here: dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob ) should be taken from ke not k_up_ppu (same in other similar places)

3. Change phase/component flux logic - assemble full phase flux first and then do x_cp upwinding once to compute component flux. The danger of old treatment is that phaseFlux being used in subsequent computations and it might be tricky to separate parts of phase flux to compute component flux

[paveltomin]

Waiting on code owner review from @CusiniM, @castelletto1, @rrsettgast, and/or @ryar9534.

[paveltomin]

@CusiniM, @castelletto1, @rrsettgast, and/or @ryar9534 please review