Better represent equations and how they are coupled/solved

One of the reasons why StGermain / StgNumerics works, is that it has very good computational mathematics abstractions for parallel scalable discretisations. In addition to this, there is a clear map from the linear algebra Ax=b system to the desired discretisation. However, the layer of abstractions above this is horrible... something along the lines of one big solver.

At the very beginning I (Steve) knew what was going on with these abstractions, but nowadays I only have an understanding of the model, and can see its limitations. What I'd would like to achieve (as is now a deliverable for AuScope?), is to both:

• Create templated PDEs... such that:
  • a PDE relates to one (or possibly more Ax=b)
  • PDE variables automatically inject into the system as expected
• Have it such that another solver can be applied to the single or coupled equations.

What we do now

Disclaimer... some of what I say below may be wrong... needs some verification!

The root component to the code or solver is a System of Linear Equations which in this case is the coupled Stokes PDEs. It states that the solver that's attached to it is the Uzawa one. The matrix representation takes the form of [K G G^T C] = [u p], which is mapped to the momentum and pressure Ax=bs.

    <struct name="stokesEqn">
      <param name="Type">Stokes_SLE</param>
      <param name="SLE_Solver">uzawa</param>
      <param name="StressTensorMatrix">k_matrix</param>
      <param name="GradientMatrix">g_matrix</param>
      <param name="CompressibilityMatrix">c_matrix</param>
      <param name="VelocityVector">solutionVelocity</param>
      <param name="PressureVector">solutionPressure</param>
      <param name="ForceVector">mom_force</param>
      <param name="ContinuityForceVector">cont_force</param>
      ...
    </struct>

Our Uzawa its actually doing a conjugate gradient of the Shur's complement (Mayhem!?!). In any case, a linear solve is performed on momentum equation to obtain a velocity solution that is substituted into the pressure equation to obtain the pressure solution, and this is iterated until things are nice. In this example the linear solve (the solve for velocity) is done using multigrid.

    <struct name="uzawa">
      <param name="Type">Stokes_SLE_UzawaSolver</param>
      <param name="velocitySolver">matrixSolver</param>
      ...
    </struct>
    <struct name="matrixSolver">
      <param name="Type">PETScMGSolver</param>
      <param name="levels">mgLevels</param>
      ...
    </struct>

The momentum equation... k_matrix is A, solutionVelocity is x and mom_force is b. The matrix is told that it can't be empty.

    <struct name="k_matrix">
      <param name="Type">StiffnessMatrix</param>
      <param name="RowVariable">VelocityField</param>
      <param name="ColumnVariable">VelocityField</param>
      <param name="RHS">mom_force</param>
      <param name="allowZeroElementContributions">False</param>
    </struct>
    <struct name="solutionVelocity">
      <param name="Type">SolutionVector</param>
      <param name="FeVariable">VelocityField</param>
    </struct>
    <struct name="mom_force">
      <param name="Type">ForceVector</param>
      <param name="FeVariable">VelocityField</param>
      <param name="ExtraInfo">context</param>
    </struct>

The pressure equation... c_matrix is A, solutionPressure is x and cont_force is b. The matrix is told it can be empty... which is important, as its not actually needed because of the solver method, and hence the matrix memory allocation shouldn't actually be created (I suppose I need to verify this really happens).

    <struct name="c_matrix">
      <param name="Type">StiffnessMatrix</param>
      <param name="RowVariable">PressureField</param>
      <param name="ColumnVariable">PressureField</param>
      <param name="RHS">cont_force</param>
      <param name="allowZeroElementContributions">True</param>
    </struct>
    <struct name="solutionPressure">
      <param name="Type">SolutionVector</param>
      <param name="FeVariable">PressureField</param>
    </struct>
    <struct name="cont_force">
      <param name="Type">ForceVector</param>
      <param name="FeVariable">PressureField</param>
      <param name="ExtraInfo">context</param>
    </struct>

This lot is really specific to the the matrix representation [K G G^T C] = [u p], which is used by Uzawa method, and to the Stokes PDEs. constitutiveMatrix is how the Lagrangian material constitutive behaviour is injected onto the K (so I suppose K is k_matrix + constitutiveMatrix). Similarly, G is g_matrix + gradientStiffnessMatrixTerm. Both constitutiveMatrix and gradientStiffnessMatrixTerm are parts of the ??? PDE.

    <struct name="constitutiveMatrix">
      <param name="Type">ConstitutiveMatrixCartesian</param>
      <param name="Swarm">picIntegrationPoints</param>
      <param name="StiffnessMatrix">k_matrix</param>
    </struct>

    <struct name="g_matrix">
      <param name="Type">StiffnessMatrix</param>
      <param name="RowVariable">VelocityField</param>
      <param name="ColumnVariable">PressureField</param>
      <param name="RHS">cont_force</param>
      <param name="allowZeroElementContributions">False</param>
    </struct>
    <struct name="gradientStiffnessMatrixTerm">
      <param name="Type">GradientStiffnessMatrixTerm</param>
      <param name="Swarm">gaussSwarm</param>
      <param name="StiffnessMatrix">g_matrix</param>
    </struct>

Steve's idea

TODO.