refactor: detangle and unify flow initialization

[paveltomin]

also make sure non-public functions are under protected

[paveltomin]

this double update is still confusing but i don't really want to touch the constitutives...

[paveltomin]

one case requires rebaseline:

Rank 0 is comparing /tmp/tmp.MtQrIy5KpC/GEOS_integratedTests_working/poromechanics/PoroElastic_hybridHexPrism_co2_sequential_3d_01/0to20_restart_000000020/rank_0000000.hdf5 with /tmp/tmp.MtQrIy5KpC/GEOS_integratedTests_baselines/poromechanics/PoroElastic_hybridHexPrism_co2_sequential_3d_01/0to20_restart_000000020/rank_0000000.hdf5 
********************************************************************************
Error: /Problem/Solvers/multiphasePoroelasticity/SolverStatistics/numSuccessfulNonlinearIterations/__values__
	Scalar values of types int32 and int32 differ: 58, 60.
********************************************************************************
********************************************************************************
Error: /Problem/Solvers/linearElasticity/SolverStatistics/numSuccessfulLinearIterations/__values__
	Scalar values of types int32 and int32 differ: 58, 60.
********************************************************************************
********************************************************************************
Error: /Problem/Solvers/linearElasticity/SolverStatistics/numSuccessfulNonlinearIterations/__values__
	Scalar values of types int32 and int32 differ: 58, 60.
********************************************************************************
********************************************************************************
Error: /Problem/Solvers/linearElasticity/maxForce/__values__
	Scalar values of types float64 and float64 differ: 59056644.21638392, 59057351.806163654.
********************************************************************************
********************************************************************************
Error: /Problem/Solvers/twoPhaseFlow/SolverStatistics/numSuccessfulLinearIterations/__values__
	Scalar values of types int32 and int32 differ: 122, 130.
********************************************************************************
********************************************************************************
Error: /Problem/Solvers/twoPhaseFlow/SolverStatistics/numSuccessfulNonlinearIterations/__values__
	Scalar values of types int32 and int32 differ: 122, 130.
********************************************************************************
********************************************************************************
Error: /Problem/domain/MeshBodies/mesh/meshLevels/Level0/nodeManager/incrementalDisplacement
	Arrays of types float64 and float64 have 10164 values of which 308 fail both the relative and absolute tests.
		Max absolute difference is at index (np.int64(3367), np.int64(2)): value = 3.329290118136317e-05, base_value = 2.0853878413888033e-05
		Max relative difference is at index (np.int64(31), np.int64(1)): value = 2.834348701588979e-18, base_value = -1.1883888138405466e-20
	Statistics of the q values greater than 1.0 defined by absolute tolerance: N = 308
		max = 1.2439022767475134, mean = 0.03769400838628123, std = 0.2132295115203364
		max is at index (np.int64(3367), np.int64(2)), value = 3.329290118136317e-05, base_value = 2.0853878413888033e-05
	Statistics of the q values greater than 1.0 defined by relative tolerance: N = 0
********************************************************************************
********************************************************************************
Error: /Problem/domain/MeshBodies/mesh/meshLevels/Level0/nodeManager/totalDisplacement
	Arrays of types float64 and float64 have 10164 values of which 1232 fail both the relative and absolute tests.
		Max absolute difference is at index (np.int64(3337), np.int64(2)): value = 0.052993050777275695, base_value = 0.05295400749413587
		Max relative difference is at index (np.int64(2643), np.int64(0)): value = 3.63995352012823e-18, base_value = 3.467236708214421e-21
	Statistics of the q values greater than 1.0 defined by absolute tolerance: N = 1232
		max = 3.9043283139826577, mean = 0.31018273308305666, std = 0.8999556883939195
		max is at index (np.int64(3337), np.int64(2)), value = 0.052993050777275695, base_value = 0.05295400749413587
	Statistics of the q values greater than 1.0 defined by relative tolerance: N = 0
********************************************************************************
********************************************************************************
Error: /Problem/domain/MeshBodies/mesh/meshLevels/Level0/ElementRegions/elementRegionsGroup/reservoir/elementSubRegions/hexahedra/globalCompDensity
	Arrays of types float64 and float64 have 300 values of which 300 fail both the relative and absolute tests.
		Max absolute difference is at index (np.int64(127), np.int64(1)): value = 990.4965609336606, base_value = 1033.6158879891204
		Max relative difference is at index (np.int64(145), np.int64(0)): value = 9.830636800993778, base_value = 0.00028436378927075074
	Statistics of the q values greater than 1.0 defined by absolute tolerance: N = 45
		max = 2434156.209658291, mean = 235406.16614701573, std = 610313.7952750626
		max is at index (np.int64(127), np.int64(0)), value = 37.21953027905491, base_value = 12.877968182471994
	Statistics of the q values greater than 1.0 defined by relative tolerance: N = 255
		max = 761325.7271238628, mean = 129221.33479600628, std = 174294.33208988397
		max is at index (np.int64(97), np.int64(0)), value = 62.131763000015084, base_value = 53.92141144485521
********************************************************************************

the reason is a bug that appears here:

GEOS/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

Line 886 in e601fad

updateCompAmount( subRegion );

where compAmount will not include the effect of netToGross (don't ask me why).

The issue is addressed here

GEOS/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

Line 523 in 1313a67

void FlowSolverBase::initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames )