Darcy flow test cases and example

Awareness of whether the Darcy velocity or interstitial velocity is solved for and how/if the other is calculated needs to be clarified and included into test cases and documentation.

The same applies to the transport of tracer or tracer*porosity fields.

Also it appears the test case porous_medium_p1p1_2d uses the porous media option.

Update: actually it appears the porous_media options model requires a viscosity field under the material phase velocity field to calculate the absorption (internal algorithm). This implies that the stress terms are also added which is incorrect for darcy flow. Also the above test case doesnt remove the time term for p1p1 - maybe this isnt possible. Basically the porous_media model and test case has two extra terms (time+stress) that shouldnt be there, so probably best not to use it but persist with the dg velocity method using python diagnostics.

Further Applications: A darcy flow model may also be useful for severe accident applications so this should be taken into consideration in development/testing.

More Again: There will be issues with shock changes between elements in material properties. A solution may be to use dg for pressure or to represent the material properties with a continuous mesh (such as porosity).

Weak BC: It appears that for the square domain test case with pressure BC generating a gradient left to right with no y flow out the top and bottom BC doesnt run correctly if instead a aligned with surface no normal component velocity is used. The pressure is fine, the velocity is fubar. This is an issue.

Correction to last point: The pressure is also fubar'd if the aligned with surface vel BC are applied weakly, if strongly then pressure is fine. This was checked in both 2d and 3d. I wander if a different element type would have this issue.

The no_normal_flow BC is possible if using CV for pressure of CG for pressure with cty integrated by parts.

However using P0P1CV with a pressure gradient doesnt work. Applying a weak inlet velocity BC instead has the flow velocity in the wrong direction (strong is fine).

Using P0P1CG with cty integrated by parts with strong BC pressure gradient is fine (schema suggest this shouldnt work if vel dg). Also fine if weak pressure BC. Not integrating cty by parts makes no difference and P1dgP2 seems fine also with no_normal_flow (schema suggests this isnt the case)

Using p1dgp2cg (or p0dgp1cg) for expanding channel case with no_normal_flow BC generates strange velocities at the corners where the channel expands. This needs checking. Using CV for pressure has same issues as stated above.

Update: It is also possible to use generic scalar and vector fields to solve a pressure poisson equation and calculate a velocity field from this after using python diagnostics (although it appears shape_dshape doesnt work very well in 1d)

Comparing solving a pressure poisson to mixed projection method there appears to be little difference in the solutions for the expanding channel test case.

Another model showed that it also had the issue at the corner of the flow. However on inspection of the velocity vector they looked fine. This may not be such a big issue.

With regard to tracer advection it is now considered best to advect tracer rather than tracer*porosity. This implies solving for velocity*porosity and including the porosity in the time term of the tracer equation. This could be included generically as scalar_field(TimeTermCoeff}. It may also be useful to have vector_field(AdvectionVelocity) and scalar_field(DiscretisedSource).
With regard to the TimeTermCoeff options should be included for whether the time derivative of the coefficient should be considered, and if yes the theta values to use. Note with regard to CV this is happening for density for ConservationMass and InternalEnergy equations already. The TimeTermCoeff will need to be on the same mesh as the tracer solution so for CV if this is porosity (which could be element wise) it needs a mapping/projecting/averaging to the CV mesh. Or just represent porosity in the model on the CV mesh initially in a similar way as other volume fraction fields.

It would be preferable to use CV for pressure as this implies a CV weighting of the cty equation which will induce a divergence free velocity field which is important for tracer advection conservation. The bug associated with CV pressure with a pressure BC should not in principle be an issue with weakly applied version.

On another note the prototype code will solve pressure with CG (more accurate) but weight the cty with CV to remain consistent.

darcy_p0p1cv_pressBCinlet test case added - fails due to pressure BC.

darcy_p0p1cv_velBCinlet test case added - this found that there is a bug for CV pressure with weak velocity BC in that the sign is wrong in the results (flows wrong way) This doesnt happen for CG or for strong velocity BC.

On inspection of the momentum code if using CV for pressure then one cannot use pressure dirichlet BC.

A test case with a second phase representing the solid has been added. This has the issue that a zero on the diagonal is generated when the solid volf goes to 1.0

A change to the code has now occurred. A TimeTermCoefficient field is available for any scalar field. This is currently only applied in the CV module. This will scale the time term by this value at n+theta. Currently the time derivative of it is ignored, but it could later be added. This permits one phase to be use for single phase darcy flow, solving for pressure, dary velocity then advecting the actual tracer concentration with the new time term coefficient being used to apply the porosity (which may be mapped from a DGP0 version). Three test cases for this have been added. Two have a step change in porosity such as to cause a spatial acceleration. With mesh adaptivity the mesh was found to be left behind at the interface. A solution to this may be to include the time term coeff for the metric advection, advect the metric with the interstitial velocity, iterate between adapt and solve within a time step or take smaller time step.

Future work should consider the two channels cases, changes in permeability as well as porosity, expanding channel, multiphase. The pressure CV BC bug if fixed would be useful. Weighting the CTY by CV when pressure is CG would also perhaps be useful.

UPDATE 30th OCTOBER 2011 +++++++++++++++++++++++++++++++++

A bug associated with the sign of weak velocity BC using CV pressure has been fixed and included in a darcy test case.

Weak pressure CV boundary conditions has been added. This is tested with a darcy test case.

The ability to have CG pressure but with a CV dual weighting of the continuity has been added, with an option. This is currently only for incompressible flow and is not an optimal implementation as ct and ct_rhs are formed to often or twice. This is tested with some darcy test cases plus a column of fluid under gravity, 2 identical fluid column under gravity and a divergence free velocity test case. This appears to be useful for p1dgp2 and p2p1 element pairs. For p1dgp2 it was found that a GMRES solver was required rather than CG for pressure - even though the CG solver diverged the answer was correct, which is strange.

All darcy test cases are now active apart from multiphase. This includes the expanding channel, step change channel, two channel and two channel with step change cases. Representing step changes in porosity could do with more investigation. Perhaps a project to CV dual is required.

 Using CV for pressure for multiphase was discovered in the code to not be possible yet as volf terms are missing from the ct matrix. This could be added as happens for CG pressure.

Future: different permiabilities, compressible, CV pressure multiphase, CG-CV pressure multiphase, darcy momentum equation, check absorption is latest, use non conservative form with explicit time for field/volf advection to be conservative and bounded, overlapping (one phase then two with special upwinding of CV surface), block matrix solves for tightly coupled multiphase systems.

31st October +++++++++++++++++++++++++++++

After discussion with CPain regarding CG pressure with CV testing of continuity he said that the pressure matrix may well require a GMRES as it isnt symmetric, CG solver may work for some problems

CV multiphase pressure has an initial implementation.

3rd Nov ++++++++++++++++++++++++++++++++

Some simple 2 phase problems have been added in 2d and 3d with constant saturations and momentum absorptions. Also I realised that p1dgp2 in 1d is not LBB stable so have decided to go straight for 2d 2 phase BL. This however has issues with skewed meshes and will run for a while then blow up - hence the need to check simpler test cases.

Also it was found that there may be an issue in 3d with the diagnostic SumVelocityDivergence when testing with CV as it oscillates between e-06 and e-012 when the problem is actually steady state.

Comparing with the prototype for incompressible it was found that the proto will solve all saturation prognostically then include a sum(Sat_{i}} into the cty_rhs as it may not be 1 and a corresponding term in cty_rhs associated with the end of time level sum - which is enforced to 1.0. Fluidity will calc the last saturation as a diagnostic. Which is better is unknown. The prototype solves all saturations as they are solved as a block with potential coupling terms.

3rd Feb 2012 ++++++++++++++++++++++++++++++++++++++++++

Blueprint achieved!