~skimura04/fluidity/fluidityDGSlope

scalar_boundary_conditions =
  (
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               ## Apply the dirichlet bc weakly.  Available
               ## automatically with discontinuous_galerkin and
               ## control_volume
               ## spatial_discretisations.
               ## If not selected boundary conditions are applied strongly.
               element apply_weakly {
                  ## If the initial condition and boundary conditions
                  ## differ, setting this option will cause the initial
                  ## condition on the boundary to be overwritten with
                  ## the boundary condition. Since you are applying the
                  ## boundary condition weakly, you probably do *not*
                  ## want this.
                  element boundary_overwrites_initial_condition {
                     empty
                  }?
               }?,
               input_choice_real
            }|
            element type {
               attribute name { "neumann" },
               input_choice_real
            }|
            robin_bc_scalar|
            ## Prevent the field from fluxing out of the boundary.
            ## Only applicable to control volume spatial discretisations.
            element type {
              attribute name { "zero_flux" },
              empty
            }|
            ## Add a bulk formulae flux to the scale field. Should
            ## really be added to Temperature, Salinity, or PhotosyntheticRadation
            ## fields, or nothing will happen. Do not add another Neumann boundary
            ## onto the same surface, or you will get an error.
            element type {
              attribute name { "bulk_formulae" },
              empty
            }|
            ## Sediment reentrainment boundary. Any sediment in the correct
            ## sediment class will be rentrained back into the flow
            ## depending on the bed shear stress and the parameters of the 
            ## sediment grain
            element type {
               attribute name { "sediment_reentrainment" },
               empty
            }|
            ## [UNDER DEVELOPMENT - DOES NOT WORK YET]
            ## Special type of Dirichlet boundary condition for the k-epsilon
            ## turbulence model. Can be used on TurbulentKineticEnergy (k)
            ## and/or TurbulentDissipation (epsilon) fields.
            ## e.g. use a Dirichlet BC on inlets and k_epsilon on walls.
            element type {
               attribute name { "k_epsilon" },
               ## Select high/low Reynolds number wall functions for k and epsilon fields.
               element string_value {
               "low_Re"|"high_Re"
               }
            }|
            ## Flux boundary condition. Currently only works with a control_volume 
            ## discretisation.
            ##
            ## This weakly enforces the total (advective plus diffusive) flux across 
            ## a boundary when solving the advection-diffusion equation.
            element type {
               attribute name { "flux" },
               input_choice_real
            }
         )
      }
)


robin_bc_scalar = 
   (
      ##  A robin boundary condition of the form
      ##  C1*T + n.(k*grad(T)) = C0
      ##  where k is the diffusivity tensor,
      ##  n the outward normal vector to the surface,
      ##  T the scalar field value on the surface,
      ##  C0 is the input order zero coefficient and
      ##  C1 is the input order one coefficient.
      ##  THIS WILL ONLY WORK FOR CONTINUOUS GALERKIN SPATIAL DISCRETISATION
      element type {
         attribute name { "robin" },
         ##  The order zero coefficient represented as C0 in
         ##  C1*T + n.(k*grad(T)) = C0
         element order_zero_coefficient {
            input_choice_real
         },
         ##  The order one coefficient represented as C1 in
         ##  C1*T + n.(k*grad(T)) = C0
         element order_one_coefficient {
            input_choice_real
         }
      }      
   )


prognostic_scalar_field =
   (
      scalar_equation_choice,
      spatial_discretisation_options,
      temporal_discretisation_options,
      (
         ## Solver
         element solver {
            linear_solver_options_asym
         }|
         ## Assume this field is being solved explicitly and skip the solver.
         ##
         ## ONLY AVAILABLE FOR PURE CONTROL VOLUME SPATIAL DISCRETISATIONS.
         ##
         ## Assumes lhs matrix only has diagonal lumped mass (times
         ## density if appropriate for equation)
         ## and divides the rhs by this.
         element explicit {
            empty
         }
      ),
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_real
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_real
         }
      )+,
      scalar_boundary_conditions*,
      ## Choice of subgridscale model to apply to this field. 
      ## Note that the selected parameterisation must be switched 
      ## on for this material phase.
      ##
      ## At this time, do not switch on this option if you are using Mellor-Yamada.
      (
         element subgridscale_parameterisation {
            attribute name { "Gent_McWilliams"}
         }|
         element subgridscale_parameterisation {
            attribute name { "GLS"}
         }|
         ## [UNDER DEVELOPMENT - DOES NOT WORK YET]
         ## Set the diffusivity to the eddy diffusivity calculated by the
         ## k-epsilon turbulence model. If using this, set the diffusivity field to
         ## diagnostic/algorithm(internal).
         element subgridscale_parameterisation {
            attribute name { "k-epsilon"},
            ## Option to set the Prandtl number (ratio of momentum diffusivity to field diffusivity)
            ##  for this field. Default = 1.0 if not specified.
            element Prandtl_number {
               input_choice_real
            }?
         }|
         element subgridscale_parameterisation {
            attribute name { "prescribed_diffusivity"}
         }
      )?,
      ## Buoyancy adjustment (vertical stabilization) schemes.
      element buoyancy_adjustment {
         ## Vertical mixing by diffusion.
         ## Stabilises unstable stratifications through mixing by diffusion.
         ##
         ## Note this scheme currently only applies to tracer fields represented in discontinuous spaces.
         element by_vertical_diffusion{
            ## Scaling amplitude of the applied diffusion. 
            ## If left unspecified, the default value of 1.0 is applied.
            element amplitude  {
                real
            }?,
            ## Calculates the rate of change of buoyancy in the direction of gravity
            ## on the field projected to continuous space - i.e. the buoyancy field is first projected to continous space.
            ## Without this option, this calculation is performed directly on the discontinuous field.
            ##
            ## Note this only make a difference when the buoyancy field is represented in a discontinuous space, and in fact,
            ## makes no difference to the routines called when this field is represented in a continuous space.
            ##
            ## This is one approach to allow the scheme to account for differences with respect to neighbouring
            ## elements when buoyancy is a discontinuous field.
            element project_buoyancy_to_continuous_space {
                empty
            }?
         }?
      }?,
      ## Diffusivity for field
      element tensor_field {
         attribute name { "Diffusivity" },
         attribute rank { "2" },
         (
            element prescribed {
                mesh_choice?,
                prescribed_tensor_field_no_adapt
            }|
            element diagnostic {
               mesh_choice?,
               (
                  tensor_python_diagnostic_algorithm |
                  internal_algorithm
               ),
               diagnostic_tensor_field
            }
         )
      }?,
      ## Source term. If the scalar field to be 
      ## solved for is subcycled then variations 
      ## in this source field across the subcycled 
      ## time steps are not considered.
      element scalar_field {
         attribute name { "Source" },
         attribute rank { "0" },
         (
            element prescribed {
               prescribed_scalar_field_no_adapt
            }|
            element diagnostic {
               (
                  scalar_python_diagnostic_algorithm |
                  internal_algorithm
               ),
               ## Add the source field directly to the 
               ## matrix system right hand side. To be 
               ## consistent this term should be tested 
               ## using the function space that 
               ## is used for the spatial discretisation.
               element add_directly_to_rhs {
                  comment
               }?,
               diagnostic_scalar_field_no_adapt
            } 
         )
      }?,
      ## Absorption term
      element scalar_field {
         attribute name { "Absorption" },
         attribute rank { "0" },
         (
            element prescribed {
               prescribed_scalar_field_no_adapt
            }|
            element diagnostic {
               (
                  scalar_python_diagnostic_algorithm |
                  internal_algorithm
               ),
               diagnostic_scalar_field_no_adapt
            }
         )
      }?,
      ## Velocity at which this substance sinks through the water column.
      ## 
      ## This velocity is in the direction of gravity so if the substance
      ## floats or swims upwards, this field should be negative.
      element scalar_field {
         attribute name { "SinkingVelocity" },
         attribute rank { "0" },
         (
            element prescribed {
               prescribed_scalar_field_no_adapt
            }|
            element diagnostic {
               (
                  scalar_python_diagnostic_algorithm |
                  internal_algorithm
               ),
               diagnostic_scalar_field_no_adapt
            }
         )
      }?,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar,
      discrete_properties_algorithm_scalar?,
      ## Set the priority of this field
      ## This determines the order in which scalar_fields are solved for:
      ##  - higher numbers have the highest priority
      ##  - lower numbers (including negative) have the lowest priority
      ##  - default if not set is 0
      element priority {
         integer
      }?
   )

prognostic_velocity_field =
   (
      velocity_equation_choice,
      ## Spatial discretisation options
      element spatial_discretisation {
         (
            ## A new version of continuous galerkin assembly.
            element continuous_galerkin {
               advection_stabilisation_options,
               ## Discretisation options for the mass terms in the velocity equation.
               element mass_terms{
                  ## Lump the mass matrix - currently required if solving for pressure
                  ## and not using the schur complement scheme under Pressure
                  element lump_mass_matrix {
                     ## Lump on the submesh.
                     ## This only works for simplex meshes and is only
                     ## strictly valid on 2d meshes.
                     element use_submesh {
                       empty
                     }?
                  }?,
                  ## Remove the mass terms from the equation.
                  element exclude_mass_terms {
                     empty
                  }?
               },
               ## Discretisation options for the advection terms in the velocity equation.
               element advection_terms {
                  ## Integrate the advection terms of the momentum equation by parts.
                  ## This allows for the imposition of weak boundary conditions.
                  ## If activated the element advection matrix takes the form:
                  ##    /                                            /
                  ##  - | (grad N_A dot nu) N_B rho dV - (1. - beta) | N_A ( div nu ) N_B rho dV
                  ##    /                                            /
                  ## otherwise it takes the standard form:
                  ##    /                                     /
                  ##    | N_A (nu dot grad N_B) rho dV + beta | N_A ( div nu ) N_B rho dV
                  ##    /                                     /
                  ## where beta is set in conservative_advection, N is
                  ## a shape function and nu is the relative nonlinear
                  ## velocity.
                  element integrate_advection_by_parts {
                     empty
                  }?,
                  ## Remove the advection terms (u.grad u rho + beta
                  ## div u rho u) from the equation.
                  ## This overrides any other advection term options
                  ## (including conservative_advection below).
                  element exclude_advection_terms {
                     empty
                  }?
               },
               ## Discretisation options for the stress terms in the velocity equation.
               element stress_terms {
                  (
                     ## Use tensor form of the stress terms:
                     ##
                     ## mu u_{i,jj}
                     ##
                     ## This is only valid for incompressible
                     ## simulations as it is basically a simplication
                     ## of full stress form when divergent elements can
                     ## be cancelled out.
                     ##
                     ## Only diagonal components of viscosity
                     ## (i.e. either isotropic or
                     ## diagonal) are physical for isotropic materials.
                     ##
                     ## If components differ from each other
                     ## this must be for numerical reasons (i.e. not
                     ## physical variations in viscosity otherwise
                     ## simplification is not valid).
                     ##
                     ## If activated, the dim x dim (in this example
                     ## 3d) discrete stress matrix takes the form:
                     ##
                     ##  /  mu_xx*N_a,x*N_b,x + mu_yy*N_a,y*N_b,y + mu_zz*N_a,z*N_b,z
                     ##  |   0                                             ...
                     ##  \   0
                     ##
                     ##      0
                     ##  ... mu_xx*N_a,x*N_b,x + mu_yy*N_a,y*N_b,y + mu_zz*N_a,z*N_b,z   ...
                     ##      0
                     ##
                     ##      0                                                           \
                     ##  ... 0                                                           |
                     ##     mu_xx*N_a,x*N_b,x + mu_yy*N_a,y*N_b,y + mu_zz*N_a,z*N_b,z   /
                     ##
                     ## which is derived from b_a^T c b_b, where:
                     ##
                     ##  b_a = / N_a,x  \   c = /  mu_xx    0      0    \
                     ##        | N_a,y  |       |    0    mu_yy    0    |
                     ##        \ N_a,z  /       \    0      0    mu_zz  /
                     ##
                     ## where N_a and N_b are shape functions of the
                     ## ath and bth node respectively and mu are the
                     ## components of the viscosity tensor.
                    element tensor_form {
                      empty
                    }|
                    ## Use partial stress form of the stress tensor:
                    ##
                    ## 2*mu (u_{(i,j)})
                    ##
                    ## This couples the velocity components together.
                    ##
                    ## If using a viscosity ALL COMPONENTS OF
                    ## VISCOSITY MUST BE SET (i.e. either
                    ## anisotropic_symmetric or
                    ## anisotropic_asymmetric tensors).
                    ##
                    ## If components differ form each other this must
                    ## be for numerical reasons (i.e. not physical
                    ## variations in viscosity as the material is isotropic).
                    ##
                    ## If activated, the dim x dim (in this example
                    ## 3d) discrete stress matrix takes the form:
                    ##
                    ##  /   2*N_a,x*N_b,x*mu_xx + N_a,y*N_b,y*mu_xy + N_a,z*N_b,z*mu_xz
                    ##  |   N_a,x*N_b,y*mu_xy                                             ...
                    ##  \   N_a,x*N_b,z*mu_xz 
                    ##
                    ##      N_a,y*N_b,x*mu_xy 
                    ##  ... N_a,x*N_b,x*mu_xy + 2*N_a,y*N_b,y*mu_yy + N_a,z*N_b,z*mu_yz   ...
                    ##      N_a,y*N_b,z*mu_yz 
                    ##
                    ##      N_a,z*N_b,x*mu_xz                                             \
                    ##  ... N_a,z*N_b,y*mu_yz                                             |
                    ##      N_a,x*N_b,x*mu_xz + N_a,y*N_b,y*mu_yz + 2*N_a,z*N_b,z*mu_zz   /
                    ##
                    ## which is derived from b_a^T c b_b, where:
                    ##
                    ##  b_a = / N_a,x   0     0   \   c = /  2*mu_xx        0          0      0    0    0   \
                    ##        |  0    N_a,y   0   |       |     0        2*mu_yy       0      0    0    0   |
                    ##        |   0     0   N_a,z |       |     0           0       2*mu_zz   0    0    0   |
                    ##        | N_a,y N_a,x   0   |       |     0           0          0    mu_xy  0    0   |
                    ##        | N_a,z   0   N_a,x |       |     0           0          0      0  mu_xz  0   |
                    ##        \   0   N_a,z N_a,y /       \     0           0          0      0    0  mu_yz /
                    ##
                    ## where N_a and N_b are shape functions of the ath and bth node respectively and mu are the components of the viscosity tensor.
                    element partial_stress_form {
                      empty
                    }|
                    ## Use full stress form of the stress tensor:
                    ##
                    ## 2*mu (u_{(i,j)} - 1/3*\delta_{ij}u_{k,k})
                    ##
                    ## This couples the velocity components together.
                    ## It is required if performing a compressible simulation.
                    ##
                    ## If using a viscosity ALL COMPONENTS OF
                    ## VISCOSITY MUST BE SET (i.e. either
                    ## anisotropic_symmetric or
                    ## anisotropic_asymmetric tensors).
                    ##
                    ## If components differ form each other this must
                    ## be for numerical reasons (i.e. not physical
                    ## variations in viscosity as the material is isotropic).
                    ##
                    ## If activated, the dim x dim (in this example
                    ## 3d) discrete stress matrix takes the form:
                    ##
                    ##  /   4/3*N_a,x*N_b,x*mu_xx + N_a,y*N_b,y*mu_xy + N_a,z*N_b,z*mu_xz
                    ##  |   N_a,x*N_b,y*mu_xy - 2/3*N_a,y*N_b,x*mu_yx                       ...
                    ##  \   N_a,x*N_b,z*mu_xz - 2/3*N_a,z*N_b,x*mu_zx
                    ##
                    ##      N_a,y*N_b,x*mu_xy - 2/3*N_a,x*N_b,y*mu_xy
                    ##  ... N_a,x*N_b,x*mu_xy + 4/3*N_a,y*N_b,y*mu_yy + N_a,z*N_b,z*mu_yz   ...
                    ##      N_a,y*N_b,z*mu_yz - 2/3*N_a,z*N_b,y*mu_zy
                    ##
                    ##      N_a,z*N_b,x*mu_xz - 2/3*N_a,x*N_b,z*mu_xz                       \
                    ##  ... N_a,z*N_b,y*mu_yz - 2/3*N_a,y*N_b,z*mu_yz                       |
                    ##      N_a,x*N_b,x*mu_xz + N_a,y*N_b,y*mu_yz + 4/3*N_a,z*N_b,z*mu_zz   /
                    ##
                    ## which is derived from b_a^T c b_b, where:
                    ##
                    ##  b_a = / N_a,x   0     0   \   c = /  4/3*mu_xx  -2/3*mu_xy -2/3*mu_xz  0    0    0   \
                    ##        |  0    N_a,y   0   |       | -2/3*mu_yx   4/3*mu_yy -2/3*mu_yz  0    0    0   |
                    ##        |   0     0   N_a,z |       | -2/3*mu_zx  -2/3*mu_zy  4/3*mu_zz  0    0    0   |
                    ##        | N_a,y N_a,x   0   |       |     0           0          0     mu_xy  0    0   |
                    ##        | N_a,z   0   N_a,x |       |     0           0          0       0  mu_xz  0   |
                    ##        \   0   N_a,z N_a,y /       \     0           0          0       0    0  mu_yz /
                    ##
                    ## where N_a and N_b are shape functions of the ath and bth node respectively and mu are the components of the viscosity tensor.
                    element stress_form {
                      empty
                    }
                  )
               },
               ## large-eddy simulation models (currently restricted to use in incompressible flow cases, where the discrete velocity is divergence-free and the eddy viscosity tensor is traceless)
               element les_model {
                  ## similar to the original Smag model
                  element second_order {
                     ## literature symbol: Cs
                     ##
                     ## suggested value 0.1
                     element smagorinsky_coefficient {
                        real
                     },
                     ## eddy viscosity (turbulent diffusion of velocity).
                     ## This is an anisotropic symmetric viscosity, added to normal viscosity field, that
                     ## carries the influence of turbulence onto the velocity field.
                     element tensor_field {
                        attribute rank { "2" },
                        attribute name { "EddyViscosity" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_tensor_field
                           }
                        )
                     }?
                  }
               |
                  ## this adds a fourth order operator to tackle the
                  ## dissipation issues of the original Smag model
                  ##
                  ## use this if your turbulence dissipates faster than it should
                  ##
                  ## requires a fine mesh to perform well
                  element fourth_order {
                     ## literature symbol: Cs
                     ##
                     ## suggested value 0.1
                     element smagorinsky_coefficient {
                        real
                     }
                  }
               |
                  ## similar to the original WALE model
                  element wale {
                     ## literature symbol: Cs
                     ##
                     ## suggested value 0.1
                     element smagorinsky_coefficient {
                        real
                     }
                  }
               |
                  ## [UNDER DEVELOPMENT] The Germano dynamic method (See Germano 1991, Lilly 1991, Pope)
                  ## accounts for turbulent flow anisotropy by implementing the Smagorinsky model
                  ## with a spatially varying coefficient. At its heart is the difference between turbulent
                  ##  stress tensors calculated from the velocity field filtered at 2 filter widths, G_t and G_m.
                  ## See the manual for further details.
                  element dynamic_les {
                     ## Filtering scale factor alpha: the ratio between the filter widths.
                     ## The test filter G_t = G_m*alpha, where G_m is the filter
                     ## width associated with the local mesh size. Filtering is achieved by applying
                     ## the inverse Helmholtz operator: u_filter_2 = (1-alpha*M)^-1 * u_filter_1.
                     ##
                     ## default/recommended value: 2.0
                     element alpha {
                        real
                     },
                     ## Solver options are required for calculation of Helmholtz-smoothed velocity.
                     element solver {
                        linear_solver_options_sym
                     },
                     ## Use the Lilly modification to the model (see Lilly 1991).
                     ## 
                     element enable_lilly {
                        empty
                     }?,
                     ## Allow negative eddy viscosity (backscatter).
                     ## WARNING: this may result in the simulation blowing up.
                     ## 
                     element enable_backscatter {
                        empty
                     }?,
                     ## Filtered velocity (diagnostic). This is the Velocity field filtered
                     ## with the test filter - see manual.
                     element vector_field {
                        attribute rank { "2" },
                        attribute name { "FilteredVelocity" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_vector_field
                           }
                        )
                     }?,
                     ## Filter width tensor (diagnostic). This is the square of the mesh size.
                     element tensor_field {
                        attribute rank { "2" },
                        attribute name { "FilterWidth" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_tensor_field
                           }
                        )
                     }?,
                     ## Strain rate S1: the test filter applied to the strain rate of the mesh-filtered velocity.
                     element tensor_field {
                        attribute rank { "2" },
                        attribute name { "StrainRate" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_tensor_field
                           }
                        )
                     }?,
                     ## Strain rate S2: the strain rate of the test-filtered velocity.
                     element tensor_field {
                        attribute rank { "2" },
                        attribute name { "FilteredStrainRate" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_tensor_field
                           }
                        )
                     }?,
                     ## Dynamic eddy viscosity (turbulent diffusion of velocity).
                     ## This is an anisotropic viscosity, added to normal viscosity field, that
                     ## carries the influence of turbulence onto the velocity field.
                     element tensor_field {
                        attribute rank { "2" },
                        attribute name { "EddyViscosity" },
                        (
                           element diagnostic {
                              internal_algorithm,
                              velocity_mesh_choice,
                              diagnostic_tensor_field
                           }
                        )
                     }?
                  }
               }?,
               ## A gauss point viscosity, calculated in such a way that it depends
               ## upon temperature. This is only valid when you have a temperature field.
               element temperature_dependent_viscosity {
                   element reference_viscosity {
                    real
                   },
                   element activation_energy {
                    real
                   }
               }?
            }|
            ## Discontinuous galerkin formulation. Confusingly it is not necessary to provide
            ## a discontinuous velocity field for this to work!
            element discontinuous_galerkin {
               ## Discretisation options for the mass terms in the velocity equation.
               element mass_terms{
                  ## Lump the mass matrix
                  element lump_mass_matrix {
                    empty
                  }?,
                  ## Remove the mass terms from the equation.
                  element exclude_mass_terms {
                     empty
                  }?
               }?,
               element viscosity_scheme {
                  (
                     ## Compact discontinuous Galerkin scheme.
                     ## (Peraire and Persson SIAM J. Sci. Comput. 30, p1806)
                     element compact_discontinuous_galerkin {
                        ## Penalty_parameter
                        ## Add penalty term Int [u][v] dS on element boundaries
                        ## scaled by C_0
                        element penalty_parameter {
                           real
                        }?
                     }|
                     ## Classical scheme from Bassi and Rebay 
                     ## (JCP 131 267-179 1997)
                     element bassi_rebay {
                        empty
                     }|
                     ## Scheme in which upwinding is applied in
                     ## alternating directions. Devised by C.Pain.
                     element arbitrary_upwind {
                        empty
                     }|
                     ## Classical interior penalty scheme
                     ## see, e.g., SIAM Journal on Numerical Analysis
                     ## Vol. 39, No. 5 (2002), pp. 1749-1779 
                     element interior_penalty {
                        ## Penalty_parameter
                        ## The penalty term Int [u][v] dS on element boundaries
                        ## is scaled by C = C_0 h**p
                        ## This option specifies the C_0
                        ## There is a theoretical lower bound for 
                        ## stability and hence convergence
                        element penalty_parameter {
                           real
                        },
                        ## Penalty_parameter
                        ## The penalty term Int [u][v] dS on element boundaries
                        ## is scaled by C = C_0 h**p
                        ## This option specifies p
                        ## Theoretically p=-1 is required for linear elements
                        element edge_length_power {
                           real
                        },
                        ## Option for how to compute the edge length h
                        element edge_length_option {
                           ## Use face integral (take sqrt in 3D)
                           element use_face_integral {
                              empty
                           }|
                           ## Use difference between element centre 
                           ## and neighbour centre
                           ## Use 2x distance to face centre on boundaries
                           element use_element_centres {
                              empty
                           }
                        },
                        ## Switch on debugging output
                        element debug {
                           ## Bound for testing element gradient matrix
                           element gradient_test_bound {
                              real
                           },
                           ## Remove the elemental integral:
                           ## Int grad u.kappa.grad v dV
                           element remove_element_integral {
                              empty
                           }?,
                           ## Remove the primal fluxes
                           element remove_primal_fluxes {
                              empty
                           }?,
                           ## Remove the penalty fluxes
                           element remove_penalty_fluxes {
                              empty
                           }?
                        }?
                     }
                  )
               },
               element advection_scheme {
                  (
                     ## Straightforward upwinding of the nonlinear velocity.
                     element upwind {
                        empty
                     }|
                     ## Disable advection
                     element none {
                        empty
                     }
                  ),
                  ## Project the advecting velocity to continuous
                  ## space. This is useful for obtaining bounded
                  ## advection schemes.
                  element project_velocity_to_continuous {
                     ## The mesh to which the projection should occur.
                     element mesh {
                        attribute name { "CoordinateMesh" }
                     }|
                     element mesh {
                        attribute name { xsd:string }
                     }
                  }?,
                  ## Integrate the advection terms of the momentum equation by parts.
                  ##
                  ## Integrating the advection term by parts is
                  ## necessary for a discontinuous
                  ## galerkin discretisation however it is possible to
                  ## select how many times the
                  ## integration by parts is performed.
                  ## Twice is the norm.
                  element integrate_advection_by_parts {
                    (
                      ## If activated the element advection matrix takes the form:
                      ##    /                                 /
                      ##    | N_A (nu dot grad N_B) dV + beta | N_A ( div nu ) N_B dV
                      ##    /                                 /
                      ##      /                                         /
                      ##  + I | N_A_i (nu dot n) N_B_o ds + [(1-I) - 1] | N_A_i (nu dot n) N_B_i ds
                      ##      /                                         /
                      ## where beta is set in conservative_advection,
                      ## N is a shape function (uppercase
                      ## subscripts indicate nodes A or B while
                      ## lowercase subscripts indicate inner or outer
                      ## faces i and o respectively), nu is the
                      ## nonlinear velocity and n is the outward
                      ## pointing normal from the element.
                      element twice {
                        empty
                      }|
                      ## If activated the element advection matrix takes the form:
                      ##    /                                        /
                      ##  - | (grad N_A dot nu) N_B dV - (1. - beta) | N_A ( div nu ) N_B dV
                      ##    /                                        /
                      ##      /                                   /
                      ##  + I | N_A_i (nu dot n) N_B_o ds + (1-I) | N_A_i (nu dot n) N_B_i ds
                      ##      /                                   /
                      ## where beta is set in conservative_advection,
                      ## N is a shape function (uppercase
                      ## subscripts indicate nodes A or B while
                      ## lowercase subscripts indicate inner or outer
                      ## faces i and o respectively), nu is the
                      ## nonlinear velocity and n is the outward
                      ## pointing normal from the element.
                      element once {
                        empty
                      }
                    )
                  },
                  ## If activated the conservation term:
                  ##  /
                  ##  | N_A ( div nu ) N_B dV
                  ##  /
                  ## is integrated_by_parts such that the element
                  ## advection matrix becomes:
                  ##         /                                        /
                  ##  - beta | (grad N_A dot nu) N_B dV + (1. - beta) | N_A (nu dot grad N_B) dV
                  ##         /                                        /
                  ##      /                                                /
                  ##  + I | N_A_i (nu dot n) N_B_o ds + [(1-I) - (1-beta)] | N_A_i (nu dot n) N_B_i ds
                  ##      /                                                /                  
                  ## where beta is set in conservative_advection, N is
                  ## a shape function (uppercase
                  ## subscripts indicate nodes A or B while lowercase
                  ## subscripts indicate inner or outer
                  ## faces i and o respectively), nu is the nonlinear
                  ## velocity and n is the outward pointing normal
                  ## from the element.
                  ## This is invariant regardless of whether the main
                  ## advection term is integrated by parts once or
                  ## twice.
                  element integrate_conservation_term_by_parts {
                    empty
                  }?
               }
            }
         ),
         ## Conservative discretisation of momentum equations
         ##  BETA=1. -- conservative (divergence form)
         ##  BETA=0. -- non-conservative
         ##  0. < BETA < 1.
         element conservative_advection {
            real
         }
      },
      ## Temporal discretisation options
      element temporal_discretisation {
         ## Implicit/explicit control (THETA)
         ##  =0.  -- explicit
         ##  =0.5 -- Crank-Nicolson
         ##  =1.  -- implicit
         element theta {
            real
         },
         ## Non-linear relaxation term
         ##  = 0  -- previous timestep velocity solution used in non-linear terms of momentum equations
         ##  = 1  -- previous iteration velocity solution used in non-linear terms of momentum equations
         ##  0 <= ITHETA <= 1
         element relaxation {
            real
         },
         ## Implicit/explicit control of velocity-divergence temporal
         ##  discretisation, theta_divergence
         ##  =0.  -- explicit
         ##  =0.5 -- Crank-Nicolson
         ##  =1.  -- implicit
         ##  Specification of this parameter will only affect
         ##  solution procedure when a free-surface is present or 
         ##  compressible projection is used. If left unspecified
         ##  the code defaults to theta_divergence = theta.
         element theta_divergence {
            real
         }?,
         element discontinuous_galerkin {
            (
               ## Use timestep subcycling to solve this equation.
               ## Specify the maximum courant number per subcycle.
               element maximum_courant_number_per_subcycle {
                  real
               }
            )
         }?
      },
      ## <b>ONLY MAKES SENSE WITH STOKES</b>
      ##
      ## Reference node (Node at which all components of velocity = 0.)
      ##
      ## Must be less than the total number of nodes.
      ## If parallel must be less than the total number of nodes of the first processor.
      ##
      ## WARNING: This only makes sense if you are solving a Stokes
      ## equation (i.e. excluding mass and advection terms in the 
      ## spatial_discretisation above) with all Neumann or periodic boundaries.
      ## 
      ##
      ## Note: it is also an option to remove the null-space of the residual vector. This
      ## option is available under solvers but only removes a single null space so will
      ## only work if your Stokes velocities are coupled in some way (not generally the 
      ## case unless using stress form viscosity and/or coriolis).
      element reference_node {
         integer
      }?,
      ## Solver
      element solver {
         linear_solver_options_asym
      },
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_vector
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_vector
         }
      )+,
      (
         ## Prescribed different regions of the field
            element prescribed_region {
            attribute name { string },
            region_ids,
            input_choice_initial_condition_vector
         }
      )*,
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         velocity_boundary_conditions
      }*,
      ## For a Newtonian fluid this is the shear viscosity.
      ##
      ## For continuous_galerkin see stress_terms to see how the
      ## viscosity tensor is dealt with in the momentum equation.
      element tensor_field {
         attribute name { "Viscosity" },
         attribute rank { "2" },
         (
            element prescribed {
               mesh_choice?,
               prescribed_tensor_field_no_adapt
            }|
            element diagnostic {
               mesh_choice?,
               (
                  internal_algorithm |
                  bulk_viscosity_algorithm |
                  tensor_python_diagnostic_algorithm
               ),
               diagnostic_tensor_field
            }
         )
      }?,
      ## Source
      element vector_field {
         attribute name { "Source" },
         attribute rank { "1" },
         (
            element prescribed {
               mesh_choice?,
               prescribed_vector_field_no_adapt
            }|
            element diagnostic {
               mesh_choice?,
               (
                  vector_python_diagnostic_algorithm |
                  imposed_material_velocity_source_algorithm |
                  internal_velocity_source_algorithm
               ),
               diagnostic_vector_field
            }
         ),
         element lump_source {
            empty
         }?
      }?,
      ## Absorption
      ##
      ## Note: When in spherical geometry the absorption is now automatically rotated.
      ## The input values below correspond to setting the diagonal of the absorption matrix
      ## in the rotated frame of reference. The columns correspond to phi, theta and r 
      ## respectively.
      element vector_field {
         attribute name { "Absorption" },
         attribute rank { "1" },
         (
            element prescribed {
               mesh_choice?,
               prescribed_vector_field_no_adapt
            }|
            element diagnostic { 
               mesh_choice?,
               (
                  vector_python_diagnostic_algorithm |
                  imposed_material_velocity_absorption_algorithm |
                  internal_velocity_absorption_algorithm
               ),
               diagnostic_vector_field
            }
         ),
         (
           ## Default absorption: no lumping, is fully evaluated before the
           ## the pressure correction.
           element default_absorption {
              empty
           }|
           ## Lump the inclusion of absorbtion terms.
           element lump_absorption {
              ## Lump on the submesh.
              ## This only works for simplex meshes and is only
              ## strictly valid on 2d meshes.
              element use_submesh {
                empty
              }?
           }|
           ## Includes the pressure correction to the velocity in the
           ## absorption term (for theta>0). This makes the absorption
           ## term more implicit. The absorption term is lumped if and
           ## only if the mass matrix is lumped (lump_mass_matrix).
           element include_pressure_correction {
              empty
           }
         )
      }?,
      ## Vertical stabilization options
      ##
      ## Note: Switching on one of these options will result in all
      ## absorption terms being included in the pressure correction
      ## if running a planar simulation.
      element vertical_stabilization {
         ## A depth dependent absorption term to
         ## increase stability in shallow regions 
         ## given by n_grav*g*dt*rho/d
         element vertical_velocity_relaxation{
            ## This option scales the vertical velocity
            ## relaxation absorption by the chosen factor.
            ## (Note: Typical value for a time step of 100s
            ##  is ~ 0.0001. i.e. it probably needs to be small.)
            element scale_factor{
              real
            }
         }?,
         ## Implicit buoyancy, calculated via
         ## theta*g*dt*drho/dr_grav
         element implicit_buoyancy{
            ## Use this option to set a 'minimum' drho/dr_grav
            ## to be used in the implicit buoyancy calculation.
            ## Note that the default value is zero. 
            element min_gradient{
              real
            }?
         }?
      }?,
      ## SurfaceTension
      element tensor_field {
         attribute name { "SurfaceTension" },
         attribute rank { "2" },
         (
            element diagnostic {
                internal_algorithm,
                attribute field_name { "MaterialVolumeFraction" },
                ## Choose whether the mass matrix is lumped or not for the calculation of the gradient
                element lump_mass_matrix {
                  empty
                }?,
                ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg
                element solver {
                  linear_solver_options_sym
                }?,
                ## Choose whether the surface tension term in the momentum equation is integrated by parts or not
                element integrate_by_parts {
                  empty
                }?,
                diagnostic_tensor_field
            }
         )
      }?,
      ## DrainageK1
      ## Used for calculating drainage for foams.
      ## It is based on the properties of the liquid within the Plateau borders.
      ## K1 = (0.0, -density*gravity/3*drag_coefficient*viscosity, 0.0)
      element vector_field {
         attribute name { "DrainageK1" },
         attribute rank { "1" },
         (
            element prescribed {
               mesh_choice?,
               prescribed_vector_field_no_adapt
            }|
            element diagnostic {
               diagnostic_vector_field
            }
         )
      }?,
      ## DrainageK2
      ## Used for calculating drainage for foams
      ## It is based on the properties of the liquid within the Plateau borders.
      ## K2 = sqrt(sqrt(3)-pi/2)*surface_tension/(6*drag_coefficient*viscosity)
      element scalar_field {
         attribute name { "DrainageK2" },
         attribute rank { "1" },
         (
            element prescribed {
               mesh_choice?,
               prescribed_scalar_field_no_adapt
            }|
            element diagnostic { 
               diagnostic_scalar_field
            }
         )
      }?,
      prognostic_vector_output_options,
      prognostic_velocity_stat_options,
      vector_convergence_options,
      prognostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_prognostic_vector_field,
      interpolation_algorithm_vector,
      discrete_properties_algorithm_vector?
   )

prognostic_pressure_field =
   (
      element spatial_discretisation {
         (
            element continuous_galerkin {
               ## remove the  fourth order pressure stabilisation term KCMC
               ## must be removed for multimaterial and free surface calculations
               element remove_stabilisation_term {
                  empty
               }?,
               ## Integrate the continuity equation by parts.
               ##
               ## This allows for the imposition of weak velocity boundary conditions with continuous_galerkin.
               ## If activated this means that the pressure gradient operator is not integrated by parts.
               element integrate_continuity_by_parts {
                  empty
               }?,
               ## ** UNDER DEVELOPMENT **
               ##
               ## Fix for low reynolds numbers by adding the viscous terms to the inversed lumped mass matrix for the pressure correction.
               ## 
               ## ** ONLY WORKS with continuous_galerkin Velocity and lumping the mass matrix and not with lumping on submesh.
               element low_re_p_correction_fix {
                  empty
               }?,
               ## Test the continuity equation with the control volume dual mesh of the finite element pressure mesh.
               ## Only works for incompressible flow. This will make the pressure matrix non symmetric which must be 
               ## considered when selecting the pressure solver options. This will not work with the free surface and 
               ## wetting and drying models. Also the mesh assoicated with the pressure field must use the lagrangian 
               ## element type (the default). This will not work with the implicit solids model using two way coupling. 
               element test_continuity_with_cv_dual {
                  comment
               }?
            }|
            ## Select this discretisation if Pressure is defined on a discontinous mesh.
            ## Only works if your Velocity (or Momentum) is continuous. With this option
            ## the continuity equation is never integrated by parts.
            element discontinuous_galerkin {
               empty
            }|
            element control_volumes {
               empty
            }
         )
      },
      (
         ## Reference node (Node at which pressure = 0.) Note that the node number must be less than the total number of nodes.
         ## If running in parallel, the node number must be less than the total number of nodes of the first processor.
         ## ** Note - it is also an option to remove the null-space of the residual vector. This
         ## option is available under solvers
         element reference_node {
            integer
         }|
         ## Input coordinates of desired reference node. If a node does not exist at these
         ## coordinates, the nearest vertex will be selected.
         element reference_coordinates {
            real_dim_vector
         }
      )?,
      ## **UNDER DEVELOPMENT**
      ## This searches the CMC matrix diagonal looking for nodes that are less than the maximum value time epsilon(0.0) (i.e. nodes that are effectively zero).
      ## It then zeros that row and column and places a one on the diagonal and a zero on the rhs.
      ## At a debug level of 2 it also prints out the value and the sum of the row values.
      ## This is useful as a debugging tool if PETSc complains about zeros on the diagonal (i.e. if you have a stiff node in your mesh) but doesn't necessary produce nice answers at the end.
      element repair_stiff_nodes {
         empty
      }?,
      ## Atmospheric pressure
      ##
      ## Manual suggests 1.01325E+6
      element atmospheric_pressure {
         real
      }?,
      ## scheme
      element scheme {
         ## Use a poisson pressure equation to calculate a first guess at pressure.
         ## This does not necessarily satisfy continuity.
         ##   = 1 -- use a poisson guess at every timestep
         ##   = 0 -- never use a poisson guess
         ##   =-1 -- use a poisson guess at the first timestep only
         ## Manual suggests -1
         element poisson_pressure_solution {
            (
               element string_value{
                  # Lines is a hint to the gui about the size of the text box.
                  # It is not an enforced limit on string length.
                  attribute lines { "1" },
                  ( "never" | "only first timestep")
               },
               comment
            )
         },
         (
            ## Use the incompressible projection method to determine
            ## the pressure and satisfy continuity
            element use_projection_method {
               ## Assemble and use the full schur complement.
               ## This allows you to not lump the mass matrix if you're using
               ## cg and to use the full momentum matrix in the projection if
               ## you so desire.
               element full_schur_complement {
                 ( 
                   ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). 
                   ## Use the full mass matrix.
                   ##
                   ## Make sure you've not lumped your mass in the velocity spatial_discretisation if you want to be consistent!
                   element inner_matrix {
                     attribute name { "FullMassMatrix" },
                     element solver {
                       linear_solver_options_sym
                     }
                   }|
                   ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). 
                   ## Use the full momentum matrix.
                   ##
                   ## Doesn't really matter if you've lumped your mass or not but why would you if you're doing a full inner solve anyway?
                   element inner_matrix {
                     attribute name { "FullMomentumMatrix" },
                     element solver {
                       linear_solver_options_asym
                     }
                   }
                 ),
                 (
                   ## Specify the preconditioner matrix to use on the schur complement.
                   ##
                   ## For DG, the LumpedSchurComplement is our best approximation to CMC.
                   element preconditioner_matrix {
                     attribute name { "LumpedSchurComplement" },
                     element lump_on_submesh {
                       empty
                     }?
                   }|
                   ## Specify the preconditioner matrix to use on the schur complement.
                   ##
                   ## DiagonalSchurComplement = C_P^T * [(Big_m)_diagonal]^-1 * C
                   element preconditioner_matrix {
                     attribute name { "DiagonalSchurComplement" },
                     empty
                   }|
                   ## Specify the preconditioner matrix to use on the schur complement.
                   ##
                   ## Pressure Mass Matrix, scaled with the inverse of viscosity. This is
                   ## shown to be spectrally equivalent to the Schur complement.
                   ## Note that this currently only works with isoviscous and/or isotropic
                   ## viscosity tensors.
                   element preconditioner_matrix {
                     attribute name { "ScaledPressureMassMatrix" },
                     empty
                   }|
                   ## Specify the preconditioner matrix to use on the schur complement.
                   element preconditioner_matrix {
                     attribute name { "NoPreconditionerMatrix" },
                     empty
                   }
                 )
               }?
            }|
            ## Use the compressible projection method to determine the
            ## pressure and satisfy continuity and the eos.
            ## This is only currently compatible with control volume
            ## pressure spatial discretisations and requires a
            ## multimaterial eos.
            element use_compressible_projection_method {
               (
                  ## Variable (normally a density) used to normalise
                  ## each materials contribution
                  ## to the C_P^T matrix.  Leave unselected for no normalisation.
                  ## Selects the MaterialDensity field.
                  element normalisation {
                     attribute name{ "MaterialDensity" },
                     empty
                  }|
                  ## Variable (normally a density) used to normalise
                  ## each materials contribution
                  ## to the C_P^T matrix.  Leave unselected for no normalisation.
                  ## Selects the bulk Density field.
                  element normalisation {
                     attribute name{ "Density" },
                     empty
                  }|
                  ## Variable (normally a density) used to normalise
                  ## each materials contribution
                  ## to the C_P^T matrix.  Leave unselected for no normalisation.
                  ## Allows the selection of an arbitrary field.
                  element normalisation {
                     attribute name{ string },
                     empty
                  }
               )?

            }
         ),
         ## rediscretise the equations at every timestep and iteration
         ## (this is useful as a debugging tool but shouldn't be necessary for any application runs)
         element update_discretised_equation {
            empty
         }?
      },
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      (
         ## Initial condition for WholeMesh
         ##asc
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_pressure
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_pressure
         }
      )*,
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         element type {
            attribute name { "dirichlet" },
            ## Apply the dirichlet bc weakly.  Available
            ## automatically with discontinuous_galerkin and
            ## control_volume spatial_discretisations.
            ## If not selected boundary conditions are applied strongly.
            element apply_weakly {
              ## If the initial condition and boundary conditions
              ## differ, setting this option will cause the initial
              ## condition on the boundary to be overwritten with
              ## the boundary condition. Since you are applying the
              ## boundary condition weakly, you probably do *not*
              ## want this.
              element boundary_overwrites_initial_condition {
                 empty
              }?
            }?,
            input_choice_real_plus_boundary_forcing
         }
      }*,
      pressure_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      detector_options_disabled_default,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )

prognostic_geostrophic_pressure_field =
   (
      element spatial_discretisation {
         ## Enables / disables RHS terms in the geopressure solver:
         ##
         ##   include_buoyancy - Include both the buoyancy and Coriolis terms on the RHS
         ##   exclude_buoyancy - Include only the Coriolis term on the RHS
         ##   exclude_coriolis - Include only the buoyancy term on the RHS
         element geostrophic_pressure_option {
            element string_value {
               "include_buoyancy" | "exclude_buoyancy" | "exclude_coriolis"
            }
         }
      },
      (
         ## Sets node 1 in the mesh as a reference node
         element reference_node {
            attribute name { "node_1" },
            element integer_value {
              attribute rank { "0" },
              attribute shape { "1" },
              "1"
            },
            comment
         }|
         ## Sets a custom node in the mesh as a reference node
         element reference_node {
            attribute name { "custom"},
            integer
         }|
         ## Sets the value of the field to zero at a supplied coordinate.
         ## This is a post-processing step after the solve, and hence should
         ## be used with the solver/remove_null_space option.
         element zero_coord {
           real_dim_vector
         }
      )?,
      (
         ## Solver
         element solver {
            linear_solver_options_sym
         }
      ),
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_real
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_real
         }
      )*,
      ## Apply a strong dirichlet boundary condition to GeostrophicPressure.
      ## If applied, this would normally be a homogeneous bc on the top but
      ## this only makes sense when excluding coriolis.
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         element type {
            attribute name { "dirichlet" },
            input_choice_real_plus_boundary_forcing
         }
      }*,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )
   
# Vertical balance pressure field, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_vertical_balance_pressure_field =
   (
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      ## Apply a strong dirichlet boundary condition to VerticalBalancePressure.
      ## This is normally be a homogeneous bc on the top surface.
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         element type {
            attribute name { "dirichlet" },
            input_choice_real_plus_boundary_forcing
         }
      }+,      
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )
   
# Hydrostatic pressure field, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_hydrostatic_pressure_field =
   (
      ## Spatial discretisation options
      element spatial_discretisation {
        (
          ## Uses an advancing front technique to integrate
          ## downwards through the mesh.
          ##
          ## Requires a discontinuous mesh!
          element discontinuous_galerkin {
            empty
          }|
          ## Solves a continuous steady state equation.
          ##
          ## Requires a continuous mesh and solver options below.
          element continuous_galerkin {
            advection_stabilisation_options,
            ## By default when the gradient of the HydrostaticPressure is
            ## subtracted from the rhs of the momentum equation, it is 
            ## integrated by parts.  This is the most general case as the
            ## HydrostaticPressure can be discontinuous.
            ## Use this option to turn off this behaviour, which will
            ## be valid for a continuous HydrostaticPressure.
            element do_not_integrate_gradient_by_parts {
              empty
            }?
          }
        )
      },
      ## Solver
      ## Only required for continuous_galerkin spatial_discretisations!
      element solver {
         linear_solver_options_asym
      }?,
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )
# Hydrostatic pressure gradient field
prognostic_hydrostatic_pressure_gradient_field =
   (
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      prognostic_vector_output_options,
      prognostic_velocity_stat_options,
      vector_convergence_options,
      prognostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_prognostic_vector_field,
      interpolation_algorithm_vector_full
   )

# Foam Velocity Potential field.
# Used to calculate the velocity of flowing foams.
prognostic_foam_velocity_potential_field =
   (
      ## Spatial discretisation options
      element spatial_discretisation {
         element foam_velocity_option {
            empty
         }
      },
     (
         ## Solver
         element solver {
            linear_solver_options_sym
         }
      ),
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_real
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids,
            input_choice_initial_condition_real
         }
      )+,
      ## Boundary conditions
      element boundary_conditions {
         attribute replaces { "boundary, TTPER1 TTPER2 TTPERI" },
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               ## Apply the dirichlet bc weakly.  Available
               ## automatically with discontinuous_galerkin,
               ## control_volume, and mixed_cv_cg
               ## spatial_discretisations.
               ## If not selected boundary conditions are applied strongly.
               element apply_weakly {
                  ## If the initial condition and boundary conditions
                  ## differ, setting this option will cause the initial
                  ## condition on the boundary to be overwritten with
                  ## the boundary condition. Since you are applying the
                  ## boundary condition weakly, you probably do *not*
                  ## want this.
                  element boundary_overwrites_initial_condition {
                     empty
                  }?
               }?,
               input_choice_real
            }|
            element type {
               attribute name { "neumann" },
               input_choice_real
            }|
            robin_bc_scalar|
            ## Prevent the field from fluxing out of the boundary.
            ## Only applicable to control volume spatial discretisations.
            element type {
              attribute name { "zero_flux" },
              empty
            }
         )
      }*,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )
   
# free surface field, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_free_surface_field =
   (
      (
         ## Spatial discretisation options
         element spatial_discretisation {
            ## Form a full 3D system for the free surface
            element free_surface_3D {
               empty
            }?,
            element fourth_order_dissipation {
               empty
            }?,
            ## low order (linear) free surface
            element low_order_free_surface {
               empty
            }?,
            (
               ## Select free surface filter
               ##
               ## With PN-PN we need some filter to supress spurious modes.
               element default_free_surface_filter {
                  empty
               }|
               ## Select free surface filter
               ##
               ## With PN-PN we need some filter to supress spurious modes.
               element user_specified_free_surface_filter {
                  ## Default is to apply 0.01 and for wetting and drying cases 1.0
                  element non_linear_filter_coefficient {
                     real
                  }
               }|
               ## Switch off free surface filter, this is more efficient than setting the coefficient to 0.
               element switch_off_free_surface_filter {
                  empty
               }
            ),
            ## Apply wetting and drying routines
            element wetting_drying {
               empty
            }?,
            ## Tidal forcing options 
            element tidal_forcing {
               ## M2
               element M2 {
                 empty
               }?,
               ## S2
               element S2 {
                  empty
               }?,
               ## N2
               element N2 {
                  empty
               }?,
               ## K2
               element K2 {
                  empty
               }?,
               ## K1
               element K1 {
                  empty
               }?,
               ## O1
               element O1 {
                  empty
               }?,
               ## P1
               element P1 {
                  empty
               }?,
               ## Q1
               element Q1 {
                  empty
               }?,             
               ## Switch on all tidal components
               element all_tidal_components {
                  empty
               }?,
               ## Switches on a Love number of 0.3
               element love_number {
                  empty
               }?,
               ## Use static tidal force for testing
               element static_tidal_force {
                  empty
               }?
            }?
         }
      ),
      # atheta, ctheta and fstheta (absorption, coriolis and free surface)
      # need to go in temporal discretisation
      # they are currently hard-coded however
      element temporal_discretisation {
         ## Implicit/explicitness for the free surface.
         ##
         ## Suggested value 1.0 (should be at least bigger than 0.5).
         ##  =0.  -- explicit
         ##  =0.5 -- Crank-Nicolson
         ##  =1.  -- implicit
         element theta {
            real
         },
         # Maybe this should go under a proper absorption field under free surface?
         ## Implicit/explicitness for absorption
         ## =0.  -- explicit (default)
         ## =0.5 -- Crank-Nicolson
         ## =1.  -- implicit
         element absorption_theta {
            real
         }?
      },
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      (
         ## Initial condition for WholeMesh
         ##
         ## Only specify one condition if not using mesh regions.
         ## Otherwise select other initial_condition option, specify region_ids
         ## and distinct names.  Then add extra intial conditions for other regions.
         element initial_condition {
            attribute name { "WholeMesh" },
            input_choice_initial_condition_real
         }|
         ## Multiple initial_conditions are allowed if specifying
         ## different values in different
         ## regions of the mesh (defined by region_ids).  In this case
         ## each initial_condition
         ## requires a distinct name for the options dictionary.
         element initial_condition {
            attribute name { string },
            region_ids?,
            input_choice_initial_condition_real
         }
      )+,
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               (
                  input_choice_real_contents|
                  element from_file {
                     element tidal {
                        attribute file_name { string },
                        attribute variable_name_amplitude { string },
                        attribute variable_name_phase { string },
                        ## See E.W. Schwiderski - Rev. Geophys. Space
                        ## Phys. Vol. 18 No. 1 pp. 243--268, 1980
                        ## for details of these constituent.
                        attribute name {"M2"|"S2"|"N2"|"K2"|"K1"|"O1"|"P1"|"Q1"|"Mf"|"Mm"|"Ssa"}
                     }+
                  }|
                  ## Set the boundary free-surface height from NEMO data.
                  ## A prescribed NEMO pressure scalar field must be set to use this option.
                  ## Set the name of the prescribed NEMO pressure scalar field below.
                  element NEMO_data {
                     attribute field_name { string }
                  }
               )
            }|
            element type {
               attribute name { "neumann" },
               input_choice_real
            }|
            robin_bc_scalar
         )
      }*,
      # no Diffusivity for field
      # no source term
      # no Absorption term
      # no Adaptive timestepping option
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar,
      discrete_properties_algorithm_scalar?
   )

# stream function, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_stream_function_field =
   (
      ## Solver
      element solver {
         linear_solver_options_asym
      },
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      # no Diffusivity for field
      # no source term
      # no Absorption term
      # no Adaptive timestepping option
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      prognostic_detector_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )


# stream function, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_multipath_stream_function_field =
   (
      ## Solver
      element solver {
         linear_solver_options_asym
      },
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## The streamfunction will be zero on the primary boundary. There must be exactly one primary boundary.
          element primary_boundary{
            empty
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               element internally_calculated {
                  empty
               }
            }
         )
      },
      ## Boundary conditions
      element boundary_conditions {
         attribute name { string },
         ## Surface id:
         element surface_ids {
            integer_vector
         },
         ## Secondary boundaries have a value given by the flux between this boundary and the primary boundary.  
         element secondary_boundary{
            ## A point on or behind the primary boundary *from* which the flux line should extend. 
            ## Note: Path should not go through periodic boundary
            element primary_point {
               real_dim_vector
            },
            ## A point on or behind the secondary boundary *to* which the flux line should extend. 
            ## Note: Path should not go through periodic boundary
            element secondary_point {
               real_dim_vector
            }
         },
         ## Type
         (
            element type {
               attribute name { "dirichlet" },
               element internally_calculated {
                  empty
               }
            }
         )
      }*,
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      # no Diffusivity for field
      # no source term
      # no Absorption term
      # no Adaptive timestepping option
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_prognostic_scalar_field,
      interpolation_algorithm_scalar_full,
      discrete_properties_algorithm_scalar?
   )

scalar_equation_choice =
   (
      (
         ## Select the equation used to solve for this field.
         ## Advection Diffusion is the norm for scalar fields.
         ## Works for all discretisation types.
         element equation { 
            attribute name { "AdvectionDiffusion" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the equation used to solve for this field.
         ## Heat Transfer equation - requires the selection of a Density field.
         ##
         ## ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS.
         ##
         ## This equation is very similar to a standard advection of temperature equation
         ## except that a coefficient density field may be spatially and/or temporally
         ## varying.
         element equation {
            attribute name { "HeatTransfer" },
            equation_coefficients
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the equation used to solve for this field.
         ## Conservation of Mass equation - requires the selection of a Density field.
         ## ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS
         element equation {
            attribute name { "ConservationOfMass" },
            equation_coefficients
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the equation used to solve for this field.
         ## Reduced Conservation of Mass equation - requires the selection of a Density field.
         ##
         ## ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS
         ##
         ## This equation is very similar to a standard conservation of mass equation
         ## except that the time discretisation uses only a single time level of density.
         ## This enables consistency between the
         ## MaterialVolumeFraction (ReducedConservationOfMass) and
         ## MaterialDensity (Advection) equations in compressible multimaterial simulations.
         element equation {
            attribute name { "ReducedConservationOfMass" },
            equation_coefficients
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the equation used to solve for this field.
         ## Internal Energy equation - requires the selection of a Density field.
         ## ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS
         ## Solve the internal energy equation for this field.
         ## Requires pressure and velocity fields to be present.
         ## Uses a nonconservative time discretisation.
         element equation {
            attribute name { "InternalEnergy" },
            equation_coefficients
         }|
         ## Option to solve for electrical potential from
         ## electrokinetic, electrochemical or electrothermal sources 
         element equation { 
            attribute name { "ElectricalPotential" }
         }
      )
   )

equation_coefficients =
   (
      ## Select density to use in the equation
      ## Use the MaterialDensity - useful for multimaterial simulations
      ## Clearly this requires a MaterialDensity field to be present
      ## Whatever field is selected must be present.
      element density {
        attribute name { "MaterialDensity" },
        coefficient_discretisation_options?
      }|
      ## Select density to use in the equation
      ## Use the bulk Density
      ## Clearly this requires a Density field to be present
      ## Whatever field is selected must be present.
      element density {
        attribute name { "Density" },
        coefficient_discretisation_options?
      }|
      ## Select density to use in the equation
      ## Whatever field is selected must be present.
      element density {
        attribute name { string },
        coefficient_discretisation_options?
      }
   )

coefficient_discretisation_options =
  (
    ## Provide discretisation options for the coefficient density field.
    ##
    ## If not provided then the discretisation options will default to those
    ## under the field that is named (hence it will generally have to be a prognostic
    ## field itself).
    element discretisation_options {
      element spatial_discretisation {
        element control_volumes {
          spatial_control_volume_options_excluding_none
        }
      },
      element temporal_discretisation {
        ## Implicit/explicit control (TTHETA)
        ##  =0.  -- explicit
        ##  =0.5 -- Crank-Nicolson
        ##  =1.  -- implicit
        element theta {
          real
        },
        element control_volumes {
          ## Only works if a control volume or coupled_cv spatial discretisation is selected.
          ## If not active then the theta specified above will be used.
          ## Otherwise use variable limited theta on individual faces.
          element limit_theta {
              empty
          }?
        }
      }
    }
  )

velocity_equation_choice =
   (
      ## Select the equation used to solve for velocity.
      ## LinearMomentum is the norm and works for all discretisation types.
      element equation {
         attribute name { "LinearMomentum" }
      }|
      ## Select the equation used to solve for velocity.
      ## Boussinesq only works for continuous_galerkin and discontinuous_galerkin.
      element equation {
         attribute name { "Boussinesq" }
      }|
      ## Select the equation used to solve for velocity.
      ## Drainage only works for continuous_galerkin and discontinuous_galerkin.
      element equation {
         attribute name { "Drainage" }
      }
   )