test_advection_diffusion_options.rnc

include "spud_base.rnc"

include "diagnostic_algorithms.rnc"
include "stabilisation.rnc"

start =
   (
      ## The root node of the options dictionary.
      element fluidity_options {
         comment,
         ## Model output files are named according to the simulation
         ## name, e.g. [simulation_name]_0.vtu. Non-standard
         ## characters in the simulation name should be avoided.
         element simulation_name {
            anystring
         },
         ## Options dealing with the specification of geometry
         element geometry {
            ## Dimension of the problem.
            ## <b>This can only be set once</b>
            element dimension {
               element integer_value {
                  attribute rank {"0"},
                  ("3"|"2"|"1")
               }
            },
            ## The position mesh
            element mesh {
               attribute name { "CoordinateMesh" },
               mesh_info
            },
            ## The velocity mesh
            element mesh {
               attribute name { "VelocityMesh" },
               mesh_info
            },
            element mesh {
               attribute name { xsd:string },
               mesh_info,
               element exclude_from_mesh_adaptivity{empty}?
            }*,
            ## Quadrature
            element quadrature {
               ## Quadrature degree
               ## 
               ## note: this specifies the degree of quadrature,
               ## not the number of gauss points
               element degree {
                  integer
               },
               ## Surface quadrature degree
               ## 
               ## note: this specifies the degree of surface
               ## quadrature not the number of surface gauss points
               element surface_degree {
                  integer
               }?,
               ## Select which family of quadrature rules to use.
               ## The default is family_cools.
               ## family_wandzura allows for degree up to 30
               ## on triangular meshes.
               ## family_grundmann_moeller allows for degree up to
               ## 29 on simplicial meshes in arbitrary dimension.
               element quadrature_family {
                 ( "family_cools" | "family_grundmann_moeller" | "family_wandzura" ) 
               }?
            }
         },
         ## Input/output options
         element io {
            ## Format for dump files. Only vtk for now.
            element dump_format {
               element string_value{
                  "vtk"
               }
            },
            (
               ## Period between dumps in time units.
               ##
               ## Specifies the period between each dump of the solution to disk.
               ## A value of 0.0 indicates that there would be a dump at every timestep.
               element dump_period {
               (
                 element constant {
                   real
                 }|
                ## Python function prescribing real input. Functions should be of the form:
                ##
                ##  def val(t):
                ##     # Function code
                ##     return # Return value
                ##
                ##
                element python {
                  python_code
                }
              )
              }|
              ## Dump period, in timesteps.
              ## 
              ## Specifies the number of timesteps between each dump of the solution to disk.
              ## A value of 0 indicates a dump at every timestep.
              element dump_period_in_timesteps {
              (
                element constant {
                  integer
                }|
                ## Python function prescribing real input. Functions should be of the form:
                ##
                ##  def val(t):
                ##     # Function code
                ##     return # Return value
                ##
                ##
                element python {
                  python_code
                }
              )
              }
            ),
            # every CPUDUM seconds write results to disc.
            ## This is usually disabled.
            element cpu_dump_period {
               real
            }?,
            ## The period between dumps in walltime seconds. This is usually disabled.
            element wall_time_dump_period {
               real
            }?,
            (
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "VelocityMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "PressureMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { "CoordinateMesh" }
               }|
               ## The mesh on to which all the fields will be
               ## interpolated for VTK output.
               element output_mesh {
                  attribute name { xsd:string }
               }
            )
         },
         ## Options dealing with time discretisation
         element timestepping {
            ## Current simulation time. At the start of the simulation this
            ## is the start time.
            element current_time {
               real,
               ## The following excerpt from the Udunits
               ## documentation explains the time unit encoding by
               ## example:
               ##
               ## The specification:
               ##
               ## seconds since 1992-10-8 15:15:42.5 -6:00
               ##
               ## indicates seconds since October 8th, 1992 at 3
               ## hours, 15 minutes and 42.5 seconds in the afternoon
               ## in the time zone which is six hours to the west of
               ## Coordinated Universal Time (i.e.  Mountain Daylight
               ## Time). The time zone specification can also be
               ## written without a colon using one or two-digits
               ## (indicating hours) or three or four digits
               ## (indicating hours and minutes).
               ##
               ## Time units are particularly required in situations
               ## where the problem (time-varying) boundary conditions
               ## and/ initial conditions are a function of time as
               ## defined by a calendar.  Examples include atmospheric
               ## forcing and climatology. The current time, specified
               ## above, is zero at the reference data/time.
               element time_units{attribute date { xsd:string }}?
            },
            ## The time step size. If adaptive time stepping is used
            ## then this is the initial time step size.
            element timestep {
               real
            },
            ## Simulation time at which the simulation should end.
            element finish_time {
               real
            },
            ## Timestep after which the simulation should end.
            element final_timestep {
               integer
            }?,
            ## Maximum CPU time (secs) taken up before
            ## simulation terminates writing results to disc.
            ## 
            ## Manual suggests 1.E+20
            element cpu_time_limit {
               real
            }?,
            ## Maximum wall time (secs) taken up before
            ## simulation terminates writing results to disc.
            ## 
            ## This is usually disabled.
            element wall_time_limit {
               real
            }?
         },
         ## The material or phase options
         element material_phase {
            attribute name { "Fluid" },
            (
               ## Velocity vector and momentum options
               element vector_field {
                  attribute rank { "1" },
                  attribute name { "Velocity" },
                  ## Field type
                  element prescribed {
                     element mesh {
                        attribute name { "VelocityMesh" }
                     },
                     prescribed_vector_field
                  }
               }?,
               ## Passive Tracer
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "Tracer" },
                  element prognostic {
                     velocity_mesh_choice,
                     prognostic_scalar_field
                  }
               },
               ## CFLNumber
               ##
               ## See http://amcg.ese.ic.ac.uk/index.php?title=Local:Diagnostics#CFL_Number
               ##
               ## Adapting to this field is not recommended
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "CFLNumber" },
                  (
                     element diagnostic {
                        internal_algorithm,
                        velocity_mesh_choice,
                        diagnostic_scalar_field
                     }
                  )
               }?
               #scalar_field_choice*,
               #vector_field_choice*,
               #tensor_field_choice*
            )
         }
      }
   )      

# Choice of input method, e.g. for boundary conditions
input_choice_real =
   (
      input_choice_real_contents
   )

input_choice_real_plus_boundary_forcing =
   (
      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"}
           }+
         }
   )

input_choice_real_plus_file =
   (
      input_choice_real_contents|
      ## Initialise the field from an existing file (indended primarily for picking up prescribed fields from previously run prognostic simulations). The file mesh must match the mesh of this field (except for piecewise constant fields which will be remapped back from the discontinuous nodal values).
      ##
      ## THIS WILL NOT WORK FOR PRESCRIBED FIELDS NOT DIRECTLY UNDERNEATH /material_phase
      element from_file {
         attribute file_name { xsd:string },
         ## The format of the input file containing field data.
         element format {
            element string_value {
               "vtu"
            }
         },
         comment
      }
   )

input_choice_real_contents =
   ## Constant value
   element constant {
      real
   }|
   ## Python function prescribing real input. Functions should be of the form:
   ##
   ##  def val(X, t):
   ##     # Function code
   ##     return # Return value
   ##
   ## where X is a tuple of length geometry dimension.
   element python {
      python_code
   }
   
# Choice of input method for initial conditions
# Note: combine = "choice" should be used here to combine with input_choice_real, but Diamond doesn't support it
input_choice_initial_condition_real =
   (
      ## Constant value
      element constant {
         real
      }|
      ## Python function prescribing real input. Functions should be of the form:
      ##
      ##  def val(X, t):
      ##     # Function code
      ##     return # Return value
      ##
      ## where X is a tuple of length geometry dimension.
      element python {
         python_code
      }|
      ## Initialise the field from an existing file (indended
      ## primarily for use in checkpointing). The file mesh must match
      ## the mesh of this field  (except for piecewise constant fields 
      ## which will be remapped back from the discontinuous nodal values). 
      ## In parallel the process number is
      ## appended to the filename, e.g. if the file_name attribute is
      ## set to "input.vtu", process 0 reads from "input-0.vtu".
      element from_file {
         attribute file_name { xsd:string },
         ## The format of the input file containing field
         ## data. Supported formats include: NetCDF CF 1.4
         ## (http://cf-pcmdi.llnl.gov/)
         element format {
            element string_value {
               "vtu"|"NetCDF - CF 1.x"
            }
         },
         comment
      }
   )

# Choice of input method, e.g. for boundary conditions
input_choice_real_dim_vector =
   (
      input_choice_real_dim_vector_contents
   )

# Choice of input method, e.g. for prescribed fields
input_choice_real_dim_vector_plus_file =
   (
      input_choice_real_dim_vector_contents|
      ## Initialise the field from an existing file (indended primarily for picking up prescribed fields from previously run prognostic simulations). The file mesh must match the mesh of this field (except for piecewise constant fields which will be remapped back from the discontinuous nodal values).
      ##
      ## THIS WILL NOT WORK FOR PRESCRIBED FIELDS NOT DIRECTLY UNDERNEATH /material_phase
      element from_file {
         attribute file_name { xsd:string },
         ## The format of the input file containing field data.
         element format {
            element string_value {
               "vtu"
            }
         },
         comment
      }
   )

input_choice_real_dim_vector_contents =
   ## Constant value
   element constant {
      real_dim_vector
   }|
   ## Python function prescribing dimensional vector input. Functions should be of the form:
   ##
   ##  def val(X, t):
   ##     # Function code
   ##     return # Return value
   ##
   ## where X and the return value are tuples of length geometry dimension.
   element python {
      python_code
   }

# Choice of input method, e.g. for boundary conditions
# this one specifies a vector field of dim minus one
input_choice_real_dim_minus_one_vector =
   (
      ## Constant value
      element constant {
         real_dim_minus_one_vector
      }|
      ## Python function prescribing dimensional vector input. Functions should be of the form:
      ##
      ##  def val(X, t):
      ##     # Function code
      ##     return # Return value
      ##
      ## where X and the return value are tuples of length geometry dimension.
      element python {
         python_code
      }
   )

## Import data from NetCDF CF-1.x file.
input_choice_netcdf =
   (
      element from_file {
         ## The format of this file should conform to NetCDF CF 1.x
         ## (http://cf-pcmdi.llnl.gov/)
         attribute file_name { xsd:string },
         comment
      }
   )

# Choice of input method for initial conditions
# Note: combine = "choice" should be used here to combine with input_choice_real, but Diamond doesn't support it
input_choice_initial_condition_vector =
   (
      ## Constant value
      element constant {
         real_dim_vector
      }|
      ## Python function prescribing dimensional vector input. Functions should be of the form:
      ##
      ##  def val(X, t):
      ##     # Function code
      ##     return # Return value
      ##
      ## where X and the return value are tuples of length geometry dimension.
      element python {
         python_code
      }|
      ## Initialise the field from an existing file (indended primarily for use in checkpointing). The file mesh must match the mesh of this field (except for piecewise constant fields which will be remapped back from the discontinuous nodal values).
      element from_file {
         attribute file_name { xsd:string },
         ## The format of the input file containing field data.
         element format {
            element string_value {
               "vtu"
            }
         },
         comment
      }
   )  

# Choice of input method for initial/boundary conditions
# version for real symmetric tensor
input_choice_real_dim_symmetric_tensor =
   (
      ## Constant symmetric tensor
      element constant {
         real_dim_symmetric_tensor
      }|
      ## Python command prescribing symmetric tensor input.
      ##
      ## Note that it is for the python function to determine 
      ## that the results it produces are, in fact, symmetric.
      ##
      ## An example that returns the three-dimensional identity:
      ##
      ##  def val(X, t):
      ##    return [[1, 0, 0],
      ##                [0, 1, 0],
      ##                [0, 0, 1]]
      element python {
         python_code
      }
   )

# Choice of input method for initial/boundary conditions
# version for real tensor
input_choice_real_dim_tensor =
   (
      ## Constant tensor
      element constant {
         real_dim_tensor
      }|
      ## Python command prescribing tensor input.
      ##
      ## An example that returns the three-dimensional identity:
      ##
      ##  def val(X, t):
      ##    return [[1, 0, 0],
      ##                [0, 1, 0],
      ##                [0, 0, 1]]
      element python {
         python_code
      }
   )

prognostic_velocity_field =
   (
      velocity_equation_choice,
      ## Spatial discretisation options
      element spatial_discretisation {
         (
            ## A new version of continuous galerkin assembly.
            element continuous_galerkin {
               ## Stabilisation options for the galerkin discretisation
               element stabilisation{
                  (
                     no_stabilisation|
                     su_stabilisation|
                     supg_stabilisation
                  )
               },
               ## Discretisation options for the mass terms in the velocity equation.
               element mass_terms{
                  ## Lump the mass matrix - currently required if solving for 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.
                     ##
                     ## 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 CAN BE
                     ## SET (i.e. either isotropic or
                     ## anistropic_symmetric with zero off diagonals
                     ## tensors).
                     ##
                     ## If diagonal 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) 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 full stress form of the stress tensor.
                     ##
                     ## This 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).
                     ##
                     ## If activated, the dim x dim (in this example
                     ## 3d) stress matrix takes the form:
                     ##
                     ##  /   2*N_a,x*N_b,x*mu_xx + N_a,y*N_b,y*mu_yy + N_a,z*N_b,z*mu_zz - 2/3*N_a,x*N_b,x*mu_xx
                     ##  |   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_xx + 2*N_a,y*N_b,y*mu_yy + N_a,z*N_b,z*mu_zz - 2/3*N_a,y*N_b,y*mu_yy   ...
                     ##      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_xx + N_a,y*N_b,y*mu_yy + 2*N_a,z*N_b,z*mu_zz - 2/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
                    }
                  )
               },
               element les_model {
                  ## suggested value 0.1
                  element smagorinsky_coefficient {
                     real
                  },
                  element form {
                     element tensor_form{empty}
                   |
                     element stress_form{empty}
                  },
                  element order {  
                     element second_order{empty}
                   |
                     element fourth_order{empty}
                  }
               }?
            }|
            ## Discontinuous galerkin formulation. This causes Momentum_DG to be
            ## called instead of diff3d. 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
                  }?
               }?,
               element viscosity_scheme {
                  (
                     ## 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
                        },
                        ## 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
                     }
                  ),
                  ## 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
                  }?
               }
            }|
            ## Use the legacy finite element discretisation
            element legacy_continuous_galerkin {
               (
                  ## balancing diffusion based on (x,y) space.
                  element balancing_diffusion_x {
                     empty
                  }|
                  ## Laxwendrof balancing diffusion.
                  element laxwendroff_balancing_diffusion {
                     empty
                  }|
                  ## (x,y,t) -balancing diffusion.
                  element balancing_diffusion_xt {
                     empty
                  }|
                  ## No balancing diffusion.
                  element no_balancing_diffusion {
                     empty
                  }|
                  ## nonlinear streamline and cross stream diffusion.
                  element nonlinear_streamline_w_crossstream_diffusion {
                     empty
                  }|
                  ## nonlinear upwind in steapest direction.
                  element nonlinear_upwind_steepest {
                     empty
                  }|
                  ## nonlinear streamline+ cross stream diffusion(but restricted)
                  element nonlinear_streamline_w_restricted_crossstream_diffusion {
                     empty
                  }|
                  ## LES option using constant length scale.
                  element les_constant_length_scale {
                     empty
                  }|
                  ## LES option using isotropic length scale.
                  element les_isotropic_length_scale {
                     empty
                  }|
                  ## LES option which uses no balancing diffusion.
                  element les_no_balancing_diffusion {
                     empty
                  }|
                  ## LES option which uses no balancing diffusion.
                  element les_no_balancing_diffusion_2 {
                     empty
                  }|
                  ## same as 45 but with 4th order dissipation.
                  element les_no_balancing_diffusion_fourth_order_dissipation {
                     empty
                  }|
                  ## LES but in tensor form like hart3d
                  element les_tensor_form {
                     empty
                  }|
                  ## LES 4th order version of 47
                  element les_fourth_order {
                     empty
                  }|
                  ## NO balancing diffusion(DISOPT=4)and take out non-linear terms.
                  element no_balancing_diffusion_remove_nonlinear_terms {
                     empty
                  }
               ),
               ## Lump the mass matrix in the momentum equation
               element lump_mass_matrix {
                  empty
               }?
            }|
            element legacy_discretisation {
               ## Legacy discretisation option (DISOPT)
               ##
               ##  From diff3d comments (other possibilities are known to exist!):
               ##  ==============================================================
               ##  DISOPT=1 - balancing diffusion based on (x,y) space.
               ##  DISOPT=2 - Laxwendrof balancing diffusion.
               ##  DISOPT=3 - (x,y,t) -balancing diffusion.
               ##  DISOPT=4 - No balancing diffusion.
               ##  DISOPT=5 - nonlinear streamline and cross stream diffusion.
               ##  DISOPT=6 - nonlinear upwind in steapest direction.
               ##  DISOPT=7 - nonlinear streamline+ cross stream diffusion(but restricted)
               ##  
               ##  DISOPT=42- LES option using constant length scale.
               ##  DISOPT=43- LES option using isotropic length scale.
               ##  DISOPT=44- LES option which uses no balancing diffusion.
               ##  DISOPT=45- LES option which uses no balancing diffusion.
               ##  DISOPT=46- same as 45 but with 4th order dissipation.
               ##  DISOPT=47 -LES but in tensor form like hart3d.
               ##  DISOPT=48 -LES 4th order version of 47.
               ##  
               ##  DISOPT=125 - NO balancing diffusion(DISOPT=4)and take out non-linear terms.
               element legacy_disopt {
                  integer
               },
               ## Lump the mass matrix in the momentum equation
               element legacy_mlump {
                  empty
               }?,
               ## Legacy discretisation option for control volume advection of momentum (DISOPT)
               ## Set to 0 if not activated
               ##            Method for face-value est.   Time-stepping     Limiting
               ##  ------------------------------------------------------------------
               ##    =0      1st order in space          Theta=specified    UNIVERSAL
               ##    =1      1st order in space          Theta=non-linear   UNIVERSAL
               ##    =2      Trapazoidal rule in space   Theta=specified    UNIVERSAL
               ##    =3      Trapazoidal rule in space   Theta=non-linear   UNIVERSAL
               ##    =4      Finite elements in space    Theta=specified    UNIVERSAL
               ##    =5      Finite elements in space    Theta=non-linear   UNIVERSAL
               ##    =6      Finite elements in space    Theta=specified    NONE
               ##    =7      Finite elements in space    Theta=non-linear   NONE
               ##    =8      Finite elements in space    Theta=specified    DOWNWIND+
               ##    =9      Finite elements in space    Theta=non-linear   DOWNWIND+
               element legacy_ndisop {
                  integer
               }?
            }
         ),
         ## Conservative discretisation of momentum equations
         ##  BETA=1. -- conservative (divergence form)
         ##  BETA=0. -- non-conservative
         ##  0. < BETA < 1.
         element conservative_advection {
            real
         },
         inner_element_velocity?
      },
      ## 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
         },
         element discontinuous_galerkin {
            (
               ## Use timestep subcycling to solve this equation.
               ## Specify the number of subcycles.
               ## This only works for pure control volume discretisations.
               element number_advection_subcycles {
                  integer
               }
            )?
         }?
      },
      ## Solver
      element solver {
         linear_solver_options_asym
      },
      constitutive_laws,
      (
         ## 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
         }
      )+,
      ## 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 {
               prescribed_values_tensor_field
            }|
            element diagnostic {
               internal_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 {
               internal_algorithm,
               diagnostic_vector_field
            }
         ),
         element lump_source {
            empty
         }?
      }?,
      ## Absorption
      element vector_field {
         attribute name { "Absorption" },
         attribute rank { "1" },
         (
            element prescribed {
               prescribed_vector_field_no_adapt
            }|
            element diagnostic { 
               internal_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 {
              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
           }
         )
      }?,
      ## 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
            }
         )
      }?,
      prognostic_vector_output_options,
      prognostic_vector_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_scalar_field =
   (
      ## 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
         }
      ),
      ## Diffusivity for field
      element tensor_field {
         attribute name { "Diffusivity" },
         attribute rank { "2" },
         element prescribed {
            mesh_choice?,
            prescribed_values_tensor_field
         }
      }?,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options
   )

# Default child of diagnostic scalar field
diagnostic_scalar_field =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?
   )
   
# Default child of diagnostic scalar field without adaptivity options
diagnostic_scalar_field_no_adapt =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      diagnostic_detector_options
   )

# Default child of diagnostic scalar field
diagnostic_scalar_field_tidal_range =
   (
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      diagnostic_detector_options,
      adaptivity_options_scalar_field,
      (
          element spin_up_time {
             real
          }
      )       
   )

# Default child of prescribed scalar field
# This is a choice of ways of inputing the prescribed field
prescribed_scalar_field =
   (
      prescribed_scalar_field_no_adapt,
      adaptivity_options_scalar_field,
      interpolation_algorithm_scalar?,
      discrete_properties_algorithm_scalar?,
      recalculation_options?
   )

# Default child of prescribed scalar field without adaptivity options
# This is a choice of ways of inputing the prescribed field
prescribed_scalar_field_no_adapt =
   (
      prescribed_values_scalar_field,
      prescribed_output_options,
      prescribed_scalar_stat_options,
      prescribed_detector_options
   )
   
prescribed_values_scalar_field =
   (
      (
         ## Value for WholeMesh
         ## Only specify one value if not using mesh regions.
         ## Otherwise select other value option, specify region_ids
         ## and distinct names.  Then add extra values for other regions.
         element value {
            attribute name { "WholeMesh" },
            input_choice_real_plus_file
         }|
         ## Multiple values are now allowed if using different value assignments
         ## in different regions of the mesh (specified by region_ids).
         ## In this case each value requires a distinct name for the options dictionary.
         element value {
            attribute name { string },
            region_ids,
            input_choice_real_plus_file
         }
      )+
   )

# Default child of diagnostic vector field
# Currently, this is empty, but in future this might include
# options that are general to all diagnostic vector fields
diagnostic_vector_field =
   (
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,      
      diagnostic_detector_options,
      vector_steady_state_options,      
      adaptivity_options_vector_field,
      recalculation_options?
   )


diagnostic_vector_field_bed_shear_stress =
   (
      (
          element density {
             real
          }
      ),       
      (
          element drag_coefficient {
             real
          }
      ),   
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      diagnostic_detector_options,
      adaptivity_options_vector_field
   )



# Default child of prescribed vector field
# This is a choice of ways of inputing the prescribed field
prescribed_vector_field =
   (
      prescribed_vector_field_no_adapt,
      adaptivity_options_vector_field,
      interpolation_algorithm_vector?,
      discrete_properties_algorithm_vector?,
      recalculation_options?
   )
   
# Default child of prescribed vector field without adaptivity options
# This is a choice of ways of inputing the prescribed field
prescribed_vector_field_no_adapt =
   (
      prescribed_values_vector_field,
      prescribed_output_options,
      prescribed_vector_stat_options,
      prescribed_detector_options
   )

prescribed_values_vector_field =
   (
      (
         ## Value for WholeMesh
         ##
         ## Only specify one value if not using mesh regions.
         ## Otherwise select other value option, specify region_ids
         ## and distinct names.  Then add extra values for other regions.
         element value {
            attribute name { "WholeMesh" },
            input_choice_real_dim_vector_plus_file
         }|
         ## Multiple values are now allowed if using different value assignments
         ## in different regions of the mesh (specified by region_ids).
         ## In this case each value requires a distinct name for the options dictionary.
         element value {
            attribute name { string },
            region_ids,
            input_choice_real_dim_vector_plus_file
         }
      )+
   )

# Default child of diagnostic tensor field
# Currently, this is empty, but in future this might include
# options that are general to all diagnostic tensor fields
diagnostic_tensor_field =
   (
      diagnostic_output_options,
      adaptivity_options_tensor_field
   )

# Default child of prescribed vector field
# This is a choice of ways of inputing the prescribed tensor field
# If the field is constant then a symmetric, or asymmetric tensor may be entered
prescribed_tensor_field =
   (
      prescribed_values_tensor_field,
      adaptivity_options_tensor_field
   )

prescribed_values_tensor_field =
   (
      (
         ## Value for WholeMesh
         ##
         ## Only specify one value if not using mesh regions.
         ## Otherwise select other value option, specify region_ids
         ## and distinct names.  Then add extra values for other regions.
         element value {
            attribute name { "WholeMesh" },
            input_choice_tensor_field
         }|
         ## Multiple values are now allowed if using different value assignments
         ## in different regions of the mesh (specified by region_ids).
         ## In this case each value requires a distinct name for the options dictionary.
         element value {
            attribute name { string },
            region_ids,
            input_choice_tensor_field
         }
      )+
   )

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
               }?
            }|
            element control_volumes {
               empty
            }
         )
      },
      ## Reference node (Node at which pressure = 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.
      ##
      ## 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
      }?,
      ## **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+5 Pa
      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" | "every timestep" | "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.
                element preconditioner_matrix {
                  attribute name { "NoPreconditionerMatrix" },
                  empty
                }
              )
            }?
         },
         ## rediscretise the equations at every timestep and iteration
         ## (for instance if using a compressible formulation
         ## or if density varies a lot or if not using a Boussinesque approximation)
         element update_discretised_equation {
            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
         }
      )*,
      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,
            ## control_volume, and legacy_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_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 {
         ## Select whether on not to include the buoyancy term "g" in the RHS:
         ##
         ##   f = - grad p_gp + g
         element geostrophic_pressure_option {
            element string_value {
               "include_buoyancy" | "exclude_buoyancy"
            }
         }
      },
      (
         ## 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
         }
      )*,
      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?
   )
   
# Balance pressure field, this is a copy of prognostic_scalar_field above
# removing all options that don't apply (mainly advection related)
prognostic_balance_pressure_field =
   (
      ## Spatial discretisation options
      element spatial_discretisation {
         # note I removed disott here as it
         # will switch on free-surface options in geoeli1p
         # should be hard-coded to 0 in comsca therefore
         # tlump is irrelevant
         # suftem should be hard-coded to .false. if nlevel is set
         ## Geostrophic pressure option
         element geostrophic_pressure_option {
            element string_value {
               "include_buoyancy"|"exclude_buoyancy"
            }
         }
         # unfortunately tbeta doesn`t make sense here
         # so we have to code an exception for not having it in comsca
      },
      # Temporal_discretisation doesn`t apply to balance pressure
      # (there`s no time derivative). Exception again
      # Solver block is the same as prognostic_scalar_field
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      # Alas, no initial_condition either, so we'd better not checkpointing it...
      ## Disables checkpointing of this field
      element exclude_from_checkpointing {
        comment
      },
      # There are boundary conditions, but nothing you can set
      # (all derived from velocity b.c.s)
      # no Diffusivity tensor_field
      # no Source scalar_field
      # no Absorption scalar_field
      # no adaptive time-stepping
      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 =
   (
      ## Spatial discretisation options
      element spatial_discretisation {
        # we don't have any yet
        empty
      },
      # Temporal_discretisation doesn`t apply to balance pressure
      # (there`s no time derivative). Exception again
      # no solver block as we don't do a PETSc solve
      # Alas, no initial_condition either...
      # boundary conditions are fixed (p=0 on top)
      # no Diffusivity tensor_field
      # no Source scalar_field
      # no Absorption scalar_field
      # no adaptive time-stepping
      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,
      discrete_properties_algorithm_scalar?
   )
   
# Richardson number field. This is a normal diagnostic scalar field, but with
# Richardson number metric options added
adaptivity_options_richardson_number_field.adaptivity_options =
   (
      ## Do not use an interpolation error driven metric for this field
      element no_interpolation_measure {
        comment
      }
   )
adaptivity_options_richardson_number_field.adaptivity_options |= adaptivity_options_scalar_field.adaptivity_options
adaptivity_options_richardson_number_field =
   (
      element adaptivity_options {
         adaptivity_options_richardson_number_field.adaptivity_options,
         ## An isotropic metric formulation based on the Richardson number. Uses
         ## the logic that wherever the Richardson number is small, we expect
         ## to need resolution. Generates edge lengths using:
         ##
         ##   Edge length = min_edge_length if Ri <= min_ri
         ##                 max_edge_length if Ri >= max_ri
         ##                 a linear fit between min_edge_length and max_edge_length otherwise
         element richardson_number_metric {
            ## Richardson number at which we have minimum edge length (default 0.0)
            element min_ri {
              real
            }?,
            ## Richardson number at which we have maximum edge length
            element max_ri {
              real
            },
            ## Minimum edge length that can be requested by the Richardson number
            ## metric
            element min_edge_length {
              real
            },
            ## Maximum edge length that can be requested by the Richardson number
            ## metric
            element max_edge_length {
              real
            },
            ## Enable to preserve anisotropy when merging with other metric
            ## formulations
            element anisotropy_preserving_merge {
               comment
            }?,
            comment
         }?,
         adaptivity_preprocessing
      }?
   )
diagnostic_richardson_number_field = diagnostic_scalar_field_no_adapt
diagnostic_richardson_number_field &= adaptivity_options_richardson_number_field

diagnostic_cv_gradient_vector_field =
   (
      ## Choose whether the mass matrix is lumped or not
      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
      }?,
      ## Normalise the gradient by its magnitude
      element normalise {
        empty
      }?,
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,
      diagnostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_vector_field,
      recalculation_options?
   )

diagnostic_gradient_vector_field =
   (
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      diagnostic_output_options,
      diagnostic_vector_stat_options,
      vector_convergence_options,
      diagnostic_detector_options,
      vector_steady_state_options,
      adaptivity_options_vector_field,
      recalculation_options?
   )

diagnostic_cv_divergence_scalar_field =
   (
      # No solver options because it can be solved directly!
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?
   )

diagnostic_fe_divergence_scalar_field =
   (
      ## Solver
      element solver {
         linear_solver_options_sym
      },
      diagnostic_output_options,
      diagnostic_scalar_stat_options,
      scalar_convergence_options,
      diagnostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      recalculation_options?
   )

# three optional input vectors for user-specified rotation matrix
rotation_matrix_components =
   (
      ## Select if you want to specify the normal direction
      ## of the rotation matrix.
      ## If off then fluidity computes the normal
      ## If on the tangents vectors must also be specified.
      element normal_direction {
         input_choice_real_dim_vector
      }?,
      ## specify first unit tangent vector to boundary
      element tangent_direction_1 {
         input_choice_real_dim_vector
      }?,
      ## specify second (if exists, i.e. if 3d) unit tangent vector to boundary
      element tangent_direction_2 {
         input_choice_real_dim_vector
      }?
   )

velocity_components_choice =
   (
      (
# rotated bcs are not implemented... this is where they should go when they are
#         element align_bc_with_surface {
#            element normal_component {
#               input_choice_real
#            }?,
#            element tangent_component_1 {
#               input_choice_real
#            }?,
#            element tangent_component_2 {
#               input_choice_real
#            }?,
#            rotation_matrix_components
#         }|
         element align_bc_with_cartesian {
            element x_component {
               input_choice_real_bc_component
            }?,
            element y_component {
               input_choice_real_bc_component
            }?,
            element z_component {
               input_choice_real_bc_component
            }?
         }
      )
   )
   
input_choice_real_bc_component = 
   (
      input_choice_real|
      element synthetic_eddy_method {
         ## use a large number to ensure Gaussian 
         ## behaviour of the fluctuating component
         element number_of_eddies {
           integer
         },
         element turbulence_lengthscale {
            input_choice_real
         },
         ## mean profile
         ##
         ## usually a function of height,
         ## for ABL simulations use a log profile
         element mean_profile {
            input_choice_real
         },
         ## Reynolds stresses profile
         ##
         ## usually a function of height,
         ## assumes that the remaining stresses are negligible 
         element Re_stresses_profile {
           input_choice_real
         }
      }
   )

# and again for robin b.c.s
robin_velocity_components_choice =
   (
      (
#         element align_bc_with_surface {
#            element normal_component {
#               element order_zero_coefficient {
#                  input_choice_real
#               },
#               element order_one_coefficient {
#                  input_choice_real
#               }
#            }?,
#            element tangent_component_1 {
#               element order_zero_coefficient {
#                  input_choice_real
#               },
#               element order_one_coefficient {
#                  input_choice_real
#               }
#            }?,
#            element tangent_component_2 {
#               element order_zero_coefficient {
#                  input_choice_real
#               },
#               element order_one_coefficient {
#                  input_choice_real
#               }
#            }?,
#            rotation_matrix_components
#         }|
         element align_bc_with_cartesian {
            element x_component {
               element order_zero_coefficient {
                  input_choice_real
               },
               element order_one_coefficient {
                  input_choice_real
               }
            }?,
            element y_component {
               element order_zero_coefficient {
                  input_choice_real
               },
               element order_one_coefficient {
                  input_choice_real
               }
            }?,
            element z_component {
               element order_zero_coefficient {
                  input_choice_real
               },
               element order_one_coefficient {
                  input_choice_real
               }
            }?
         }
      )
   )

velocity_boundary_conditions =
   (
      (
         element type {
            attribute name { "dirichlet" },
            ## Apply the dirichlet bc weakly.  Available automatically
            ## with a discontinuous_galerkin Velocity
            ## spatial_discretisation.  Available if you
            ## integrate_continuity_by_parts with a
            ## continuous_galerkin Pressure or use a control_volume
            ## Pressure spatial_discretisation and/or
            ## integrate_advection_by_parts under Velocity
            ## spatial_discretisation with continuous_galerkin.
            ##
            ## 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
               }?
            }?,
            velocity_components_choice
         }|
         element type {
            attribute name { "neumann" },
            velocity_components_choice
         }|
         element type {
            attribute name { "robin" },
            robin_velocity_components_choice
         }|
         element type {
            attribute name { "free_surface" }
         }|
         ## Apply quadratic drag. Specify drag coefficient. If you
         ## want to exactly replicate results from using the OCEDRA
         ## option, set this to 0.003 and remember to apply to both
         ## bottom and sides.
         element type {
            attribute name { "drag" },
            input_choice_real,
            (
              ## Use a quadratic drag.
              ##
              ## This means that the drag coefficient is nondimensional.
              element quadratic_drag {
                empty
              }|
              ## Use a linear drag (basically just a surface absorption term).
              ##
              ## This means that the drag coefficient has units of momentum.
              element linear_drag {
                empty
              }
            )
         }|

         ## Apply wind forcing specified by stress or wind velocity.
         ## Replaces windy.dat and windy.py
         element type {
            attribute name { "wind_forcing" },
            (
               ## Wind forcing with user specified wind stress
               ##
               ## <b> Note that the stress needs to be specified
               ## using the same density units as the reference_density 
               ## under equation of state.</b>So if you use the recommended
               ## non-dimensional value of 1.0 for reference_density and
               ## your calculated stress is in kg m^-1s^-2 and the dimensional
               ## reference_density is 1000.0 kg m^-3, you need to divide
               ## the calculated stress in SI units by 1000.0.
               element wind_stress {
                  input_choice_real_dim_minus_one_vector|
                  element from_netcdf {
                     ## The format of this file should conform to NetCDF CF 1.x
                     ## (http://cf-pcmdi.llnl.gov/).
                     attribute file_name { xsd:string },
                     attribute east_west { xsd:string },
                     attribute north_south { xsd:string },
                     comment
                  }
               }|
               ## Wind forcing with user specified 10m wind velocity
               element wind_velocity {
                  ## Specify wind drag coefficient (dimensionless)
                  ## Suggested value: 4.0e-4
                  element wind_drag_coefficient {
                     input_choice_real
                  },
                  ## Density of air. 
                  ##
                  ## <b>Note that you have to specify
                  ## this density in the same units as the 
                  ## reference_density under equation of state.</b>
                  ## So with a typicial value of rho_air=1.3 kgm^-3
                  ## and rho_water=1000 kgm^-3, if you fill in the 
                  ## recommended (non-dimensional) value of 1.0 for 
                  ## reference_density, this field needs to be 1.3e-3.
                  element density_air {
                     real
                  },
                  ## Specify wind velocity
                  element wind_velocity {
                     input_choice_real_dim_minus_one_vector|
                     element from_netcdf {
                        ## The format of this file should conform to NetCDF CF 1.x
                        ## (http://cf-pcmdi.llnl.gov/)
                        attribute file_name { xsd:string },
                        attribute east_west { xsd:string },
                        attribute north_south { xsd:string },
                        comment
                     }
                  }
               }
            )
         }|

         ## When using control_volumes under Pressure
         ## spatial_discretisation or when using
         ## integrate_continuity_by_parts with continuous_galerkin
         ## Pressure and continuous_galerkin Velocity this
         ## boundary condition type imposes a weak no normal flow
         ## boundary condition on the surface_ids specified.
         element type {
            attribute name { "no_normal_flow" },
            empty
         }
      )
   )
   
# Output options for prognostic fields
prognostic_scalar_output_options =
   (
      ## Specify what is written to vtu dump files
      element output {
         ## By default each field in the options file is written to vtu.
         ## Select this option to exclude this field.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?,
         ## Output a file details the convergence (or otherwise) of
         ## this field with every advective nonlinear
         ## iteration.
         ## ONLY WORKS FOR PURE CONTROL VOLUME DISCRETISATIONS.
         element convergence_file {
            comment
         }?
      }
   )

# Output options for pressure (can't have a convergence file)
pressure_output_options =
   (
      ## Specify what is written to vtu dump files
      element output {
         ## By default each field in the options file is written to vtu.
         ## Select this option to exclude this field.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?,
         ## Write out some extra debugging vtu files that can be used
         ## to analyse what goes on in the pressure projection steps.
         ## WARNING: this may create a huge amount of vtu files, as 
         ## multiple files are written per nonlinear iteration.
         element debugging_vtus {
            empty
         }?
      }
   )

# Output options for prognostic fields
prognostic_vector_output_options =
   (
      ## Specify what is written to vtu dump files
      element output {
         ## By default each field in the options file is written to vtu.
         ## Select this option to exclude this field.
         element exclude_from_vtu {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## on the previous timestep.
         ## (included under the name: Old<field_name> )
         element include_previous_time_step {
            empty
         }?,
         ## Select this option to also write the values of this field
         ## used in the nonlinear iteration.
         ## (included under the name: Nonlinear<field_name> )
         element include_nonlinear_field {
            empty
         }?
      }
   )
   
# Output options for all other fields
field_output_options =
   (
      ## Specify what is written to vtu dump files
      element output {
         ## By default each field in the options file is written to vtu.
         ## Select this option to exclude this field.
         element exclude_from_vtu {
            empty
         }?
      }
   )
   
diagnostic_output_options = field_output_options
prescribed_output_options = field_output_options

# Options for inclusion/exclusion of standard field statistics from the .stat
# file
include_stat =
   (
      ## Include this field in the .stat file (magnitude and components)
      element include_in_stat {
         comment
      }
   )
exclude_components_from_stat =
   (
      ## Include just the magnitude of this field in the .stat file
      ## (excluding the components)
      element exclude_components_from_stat {
         comment
      }
   )
exclude_stat =
   (
      ## Exclude this field from the .stat file.
      element exclude_from_stat {
         comment
      }
   )

# Diagnostic statistics options for prognostic scalar fields
prognostic_scalar_stat_options = 
   (
      ## Specify what is added to .stat files
      element stat {
        prognostic_scalar_stat_options.stat
      }
   )
  
# Diagnostic statistics for all other scalar fields
prognostic_scalar_stat_no_old_or_nonlinear_options =
   (
      ## Specify what is added to .stat files
      element stat {
         prognostic_scalar_stat_no_old_or_nonlinear_options.stat
         
      }
   )

diagnostic_scalar_stat_options = prognostic_scalar_stat_no_old_or_nonlinear_options
prescribed_scalar_stat_options = prognostic_scalar_stat_no_old_or_nonlinear_options

# Combining of stat elements for scalar fields
prognostic_scalar_stat_options.stat = prognostic_scalar_stat_no_old_or_nonlinear_options.stat
prognostic_scalar_stat_options.stat &=
   (
      ## Enable to include the previous timestep value of this field in the .stat file.
      element include_previous_time_step {
         comment
      }?,
      ## Enable to include the values of this field in the nonlinear
      ## iteration in the .stat file.
      element include_nonlinear_field {
         comment
      }?
   )
prognostic_scalar_stat_no_old_or_nonlinear_options.stat = 
   (
      exclude_stat?,
      cv_stats?,
      surface_integral_stats_scalar*,
      mixing_stats*
   )   
   
# Diagnostic statistics options for vector fields, with enabled by default
vector_field_stat_options_enabled_default = include_stat
vector_field_stat_options_enabled_default |= exclude_components_from_stat
vector_field_stat_options_enabled_default |= exclude_stat

# Diagnostic statistics options for vector fields, with enabled by default
vector_field_stat_options_disabled_default = exclude_stat
vector_field_stat_options_disabled_default |= exclude_components_from_stat
vector_field_stat_options_disabled_default |= include_stat

# Diagnostic statistics for prognostic vector fields
prognostic_vector_stat_options =
   (
      ## Specify what is added to .stat files
      element stat {
         (
            prognostic_vector_stat_options.stat
         )
      }      
   )

# Diagnostic statistics for all other vector fields
prognostic_vector_stat_no_old_or_nonlinear_options =
   (
      ## Specify what is added to .stat files
      element stat {
         prognostic_vector_stat_no_old_or_nonlinear_options.stat
      }
   )
diagnostic_vector_stat_options = prognostic_vector_stat_no_old_or_nonlinear_options
prescribed_vector_stat_options = prognostic_vector_stat_no_old_or_nonlinear_options

# Combining of stat elements for vector fields
prognostic_vector_stat_options.stat = prognostic_vector_stat_no_old_or_nonlinear_options.stat
prognostic_vector_stat_options.stat &=
   (
      ## Specify how the previous timestep value of this field is added to the .stat file.
      element previous_time_step {
         vector_field_stat_options_disabled_default
      },
      ## Specify how the values of this field used in the nonlinear iteration are added to the .stat file.
      element nonlinear_field {
         vector_field_stat_options_disabled_default
      },
      ## What surface IDs do you want to do the calculation over?
      element compute_body_forces_on_surfaces {
         integer_vector
      }?
   )
prognostic_vector_stat_no_old_or_nonlinear_options.stat =
   (
      vector_field_stat_options_enabled_default,
      surface_integral_stats_vector*
   )

# Convergence options for prognostic scalar fields
scalar_convergence_options =
   (
      ## Decide whether this field is tested for convergence
      ## during nonlinear iterations
      ## (if /timestepping/nonlinear_iterations and
      ## /timestepping/nonlinear_iterations/tolerance are
      ## enabled).
      ## Also specifies whether the field is added to the 
      ## convergence file (if /io/convergence_file is enabled).
      element convergence {
         (
            ## Include this field in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            element include_in_convergence {
               comment
            }|
            ## Exclude this field from convergence testing and file 
            element exclude_from_convergence {
               comment
            }
         )
      }
   )

# Convergence statistics options for prognostic vector fields (velocity)
vector_convergence_options =
   (
      ## Decide whether this field is tested for convergence
      ## during nonlinear iterations
      ## (if /timestepping/nonlinear_iterations and
      ## /timestepping/nonlinear_iterations/tolerance are
      ## enabled).
      ## Also specifies whether the field is added to the 
      ## convergence file (if /io/convergence_file is enabled).
      element convergence {
         (
            ## Include this field (magnitude and components)
            ## in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            element include_in_convergence {
               comment
            }|
            ## Include just the magnitude of this field 
            ## in convergence testing
            ## (if /timestepping/nonlinear_iterations and
            ## /timestepping/nonlinear_iterations/tolerance are
            ## enabled) and file (if /io/convergence_file is enabled)
            ## i.e. excluding the components
            element exclude_components_from_convergence {
               comment
            }|
            ## Exclude this field entirely from convergence testing and file 
            element exclude_from_convergence {
               comment
            }
         )
      }
   )

# Steady state options for prognostic scalar fields
scalar_steady_state_options =
   (
      ## Decide whether this field is tested for a steady state
      ## between timesteps
      ## (if /timestepping/steady_state is
      ## enabled).
      element steady_state {
         (
            ## Include this field in steady state testing
            ## (if /timestepping/steady_state is
            ## enabled)
            element include_in_steady_state {
               comment
            }|
            ## Exclude this field from steady state testing
            element exclude_from_steady_state {
               comment
            }
         )
      }
   )

# Steady state statistics options for prognostic vector fields (velocity)
vector_steady_state_options =
   (
      ## Decide whether this field is tested for a steady state
      ## between timesteps
      ## (if /timestepping/steady_state is
      ## enabled).
      element steady_state {
         (
            ## Include this field (magnitude and components)
            ## in steady state testing
            ## (if /timestepping/steady_state is enabled)
            element include_in_steady_state {
               comment
            }|
            ## Include just the magnitude of this field 
            ## in steady state testing
            ## (if /timestepping/steady_state is
            ## enabled)
            ## i.e. excluding the components
            element exclude_components_from_steady_state {
               comment
            }|
            ## Exclude this field entirely from convergence testing and file 
            element exclude_from_steady_state {
               comment
            }
         )
      }
   )

# Options for whether a field is to be included in detector output.
detector_options_enabled_default = 
   (
      ## Specify what is added to detector files
      element detectors {
         (
            ## This field is output at each detector location.
            element include_in_detectors {
               comment
            }|
            ## This field is not output at detector locations.
            element exclude_from_detectors {
               comment
            }
         )
      }
   )

# Options for whether a field is to be included in detector output.
detector_options_disabled_default = 
   (
      ## Specify what is added to detector files
      element detectors {
         (
            ## This field is not output at detector locations.
            element exclude_from_detectors {
               comment
            }|
            ## This field is output at each detector location.
            element include_in_detectors {
               comment
            }
         )
      }
   )

# Detector output defaults on for prognostic and diagnostic fields, 
# off for prescribed.
prognostic_detector_options = detector_options_enabled_default
diagnostic_detector_options = detector_options_enabled_default
prescribed_detector_options = detector_options_disabled_default

adaptivity_preprocessing =
      ## Occasionally, it is desirable to apply operations or filters
      ## to fields before using them for the purposes of adaptivity.
      element preprocessing {
        (
            ## Invert a helmholtz operator to smooth out the field
            ## before using it to adapt. This can help with noisy
            ## fields.
            element helmholtz_smoother {
                 element smoothing_length_scale {
                     real_dim_symmetric_tensor
                 },
                 element solver {
                   linear_solver_options_sym
                 }
            }
        )
      }?

adaptivity_options_prognostic_scalar_field =
   (
      element adaptivity_options {
         (
            ## When specifying absolute measure
            ## one specifies the absolute interpolation 
            ## error in the units of the field that is 
            ## being adapted, e.g. you can specify
            ## the error to be 1.3 units 
            element absolute_measure {
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_scalar_field_no_adapt
                  }
               }
            }|
            ##When specifying relative measure
            ##one specifies the interpolation error
            ##relative to the field that is
            ##being adapted, e.g. you can specify
            ##the error to be 5% (i.e. 0.05)
            element relative_measure {
               element scalar_field {
                  attribute rank { "0" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_scalar_field_no_adapt
                  }
               },
               ## The relative Hessian is calculated according to:
               ##
               ##   Q = H / max{ |psi|, psi_min}
               ##
               ## where H is the Hessian, psi is the field value and
               ## psi_min is the tolerance. The tolerance prevents
               ## division by zero errors.
               ##
               ## Source: Fluidity/ICOM manual draft version 1.2
               element tolerance {
                  real
               }
            }|
            ## Adapt using the anisotropic strategy of 
            ## Formaggia, Perotto, Micheletti.
            ## Rather than taking two derivatives
            ## and deriving the anisotropic information,
            ## this approach computes an anisotropic Zienkiewicz-Zhu
            ## error estimator for each element. The approach then
            ## optimises the element orientation and length scales
            ## to equidistribute the estimated error.
            element anisotropic_zienkiewicz_zhu {
              ## Tau is an anisotropic estimate for the H1 seminorm of the
              ## error. This estimator is efficient and reliable, under the
              ## caveat that the initial mesh is sufficiently fine so as to
              ## prevent data oscillation. (Micheletti & Perotto, 2006)
              ## Typically, tau will be ~= 6-8 * |e|_H1.
              element tau {
                real
              }
            }
         ),
         adaptivity_preprocessing
      }?
   )

adaptivity_options_scalar_field.adaptivity_options =
   (
      ## When specifying absolute measure
      ## one specifies the absolute interpolation 
      ## error in the units of the field that is 
      ## being adapted, e.g. you can specify
      ## the error to be 1.3 units 
      element absolute_measure {
         element scalar_field {
            attribute rank { "0" },
            attribute name { "InterpolationErrorBound" },
            element prescribed {
               prescribed_scalar_field_no_adapt
            }
         }
      }|
      ## When specifying relative measure
      ## one specifies the interpolation error
      ## relative to the field that is
      ## being adapted, e.g. you can specify
      ## the error to be 5% (i.e. 0.05)
      element relative_measure {
         element scalar_field {
            attribute rank { "0" },
            attribute name { "InterpolationErrorBound" },
            element prescribed {
               prescribed_scalar_field_no_adapt
            }
         },
         ## The relative Hessian is calculated according to:
         ##
         ##   Q = H / max{ |psi|, psi_min}
         ##
         ## where H is the Hessian, psi is the field value and
         ## psi_min is the tolerance. The tolerance prevents
         ## division by zero errors.
         ##
         ## Source: Fluidity/ICOM manual draft version 1.2
         element tolerance {
            real
         }
      }|
      ## Adapt using the anisotropic strategy of 
      ## Formaggia, Perotto, Micheletti.
      ## Rather than taking two derivatives
      ## and deriving the anisotropic information,
      ## this approach computes an anisotropic Zienkiewicz-Zhu
      ## error estimator for each element. The approach then
      ## optimises the element orientation and length scales
      ## to equidistribute the estimated error.
      element anisotropic_zienkiewicz_zhu {
        ## Tau is an anisotropic estimate for the H1 seminorm of the
        ## error. This estimator is efficient and reliable, under the
        ## caveat that the initial mesh is sufficiently fine so as to
        ## prevent data oscillation. (Micheletti & Perotto, 2006)
        ## Typically, tau will be ~= 6-8 * |e|_H1.
        element tau {
          real
        }
      }
   )
adaptivity_options_scalar_field =
   (
      element adaptivity_options {
         adaptivity_options_scalar_field.adaptivity_options,
         adaptivity_preprocessing
      }?
   )

adaptivity_options_prognostic_vector_field =
   (
      ## Adaptivity weights
      element adaptivity_options {
         (
            ## When specifying absolute measure
            ## one specifies the absolute interpolation 
            ## error in the units of the field that is 
            ## being adapted, e.g. you can specify
            ## the error to be 1.3 units 
            element absolute_measure {
               element vector_field {
                  attribute rank { "1" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_vector_field_no_adapt
                  }
               }
            }|
            ## When specifying relative measure
            ## one specifies the interpolation error
            ## relative to the field that is
            ## being adapted, e.g. you can specify
            ## the error to be 5% (i.e. 0.05)
            element relative_measure {
               element vector_field {
                  attribute rank { "1" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_vector_field_no_adapt
                  }
               },
               ## The relative Hessian is calculated according to:
               ##
               ##   Q = H / max{ |psi|, psi_min}
               ##
               ## where H is the Hessian, psi is the field value and
               ## psi_min is the tolerance. The tolerance prevents
               ## division by zero errors.
               ##
               ## Source: Fluidity/ICOM manual draft version 1.2
               element tolerance {
                  real_dim_vector
               }
            }
         ),
         adaptivity_preprocessing
      }?
   )

adaptivity_options_vector_field =
   (
      ## Adaptivity weights
      element adaptivity_options {
         (
            ## When specifying absolute measure
            ## one specifies the absolute interpolation 
            ## error in the units of the field that is 
            ## being adapted, e.g. you can specify
            ## the error to be 1.3 units 
            element absolute_measure {
               element vector_field {
                  attribute rank { "1" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_vector_field_no_adapt
                  }
               }
            }|
            ## When specifying relative measure
            ## one specifies the interpolation error
            ## relative to the field that is
            ## being adapted, e.g. you can specify
            ## the error to be 5% (i.e. 0.05)
            element relative_measure {
               element vector_field {
                  attribute rank { "1" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_vector_field_no_adapt
                  }
               },
               ## The relative Hessian is calculated according to:
               ##
               ##   Q = H / max{ |psi|, psi_min}
               ##
               ## where H is the Hessian, psi is the field value and
               ## psi_min is the tolerance. The tolerance prevents
               ## division by zero errors.
               ##
               ## Source: Fluidity/ICOM manual draft version 1.2
               element tolerance {
                  real_dim_vector
               }
            }
         ),
         adaptivity_preprocessing
      }?
   )

adaptivity_options_prognostic_tensor_field =
   (
      ## Adaptivity weights
      element adaptivity_options {
         (
            ## When specifying absolute measure
            ## one specifies the absolute interpolation 
            ## error in the units of the field that is 
            ## being adapted, e.g. you can specify
            ## the error to be 1.3 units 
            element absolute_measure {
               element tensor_field {
                  attribute rank { "2" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_values_tensor_field
                  }
               }
            }|
            ## When specifying relative measure
            ## one specifies the interpolation error
            ## relative to the field that is
            ## being adapted, e.g. you can specify
            ## the error to be 5% (i.e. 0.05)
            element relative_measure {
               element tensor_field {
                  attribute rank { "2" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_values_tensor_field
                  }
               },
               ## The relative Hessian is calculated according to:
               ##
               ##   Q = H / max{ |psi|, psi_min}
               ##
               ## where H is the Hessian, psi is the field value and
               ## psi_min is the tolerance. The tolerance prevents
               ## division by zero errors.
               ##
               ## Source: Fluidity/ICOM manual draft version 1.2
               element tolerance {
                  real_dim_tensor
               }
            }
         ),
         adaptivity_preprocessing
      }?
   )

adaptivity_options_tensor_field =
   (
      ## Adaptivity weights
      element adaptivity_options {
         (
            ## When specifying absolute measure
            ## one specifies the absolute interpolation 
            ## error in the units of the field that is 
            ## being adapted, e.g. you can specify
            ## the error to be 1.3 units 
            element absolute_measure {
               element tensor_field {
                  attribute rank { "2" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_values_tensor_field
                  }
               }
            }|
            ## When specifying relative measure
            ## one specifies the interpolation error
            ## relative to the field that is
            ## being adapted, e.g. you can specify
            ## the error to be 5% (i.e. 0.05)
            element relative_measure {
               element tensor_field {
                  attribute rank { "2" },
                  attribute name { "InterpolationErrorBound" },
                  element prescribed {
                     prescribed_values_tensor_field
                  }
               },
               ## The relative Hessian is calculated according to:
               ##
               ##   Q = H / max{ |psi|, psi_min}
               ##
               ## where H is the Hessian, psi is the field value and
               ## psi_min is the tolerance. The tolerance prevents
               ## division by zero errors.
               ##
               ## Source: Fluidity/ICOM manual draft version 1.2
               element tolerance {
                  real_dim_tensor
               }
            }
         ),
         adaptivity_preprocessing
      }?
   )
generic_aliased_field =
   (
      attribute material_phase_name { xsd:string },
      attribute field_name { xsd:string }
   )

# Most common mesh choices
mesh_choice = 
   (
      element mesh {
         attribute name { xsd:string }
      }|
      element mesh {
         attribute name { "CoordinateMesh" }
      }|
      element mesh {
         attribute name { "VelocityMesh" }
      }|
      element mesh {
         attribute name { "PressureMesh" }
      }
   )

# Not really a choice, for fields that have to be on the velocity mesh
# currently that's all scalar fields, except pressure
# and of course velocity itself
# If you want to implement scalar fields on other meshes, feel free to do so
# but bare in mind you need to make sure the field stays outside RMEM.
# Currently all scalar fields are packed in RMEM with length nonods
velocity_mesh_choice =
   (
      (
         element mesh {
            attribute name { "VelocityMesh" }
         }|
         element mesh {
            attribute name { "PressureMesh" }
         }|
         element mesh {
            attribute name { "CoordinateMesh" }
         }|
         element mesh {
            attribute name { string }
         }
      )
   )

pressure_mesh_choice =
   (
      (
         element mesh {
            attribute name { "PressureMesh" }
         }|
         element mesh {
            attribute name { "VelocityMesh" }
         }|
         element mesh {
            attribute name { "CoordinateMesh" }
         }|
         element mesh {
            attribute name { string }
         }
      )
   )

# This is the choice of additional scalar field to be solved for
scalar_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
      (
         element scalar_field {
            attribute rank { "0" },
            attribute name { xsd:string },
            ## Field type
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element diagnostic {
                  scalar_diagnostic_algorithms,
                  velocity_mesh_choice,
                  python_diagnostic_field_code?,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Prognostic scalar fields below this
         element ___Prognostic_Fields_Below___ {
            empty
         }|

# This is the long list of fields that FLUIDITY knows about
# -- First is a list of fields that are primarily prognostic,
#    but can be set to prescribed, or aliased...
# -- The list is in order of most frequently used.

         ## Salinity
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Salinity" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Temperature
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Temperature" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Background Temperature
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BackgroundTemperature" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Passive Tracer
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Tracer" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Free Surface
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FreeSurface" },
               ## Free Surface
               ## NOTE: the diagnostic FreeSurface field only works in combination
               ## with the free_surface boundary condition applied to the Velocity
               ## field. It gives you a 3D field (constant over the vertical)
               ## of the free surface elevation.
               element diagnostic {
                  internal_algorithm,
                  # this is hard-coded on the PressureMesh as long as the Pressure is
                  # if this is no longer true, it should be option-checked to be on the
                  # same mesh as Pressure
                  ## Must be on the same mesh as Pressure
                  element mesh {
                     attribute name { "PressureMesh" }
                  },
                  diagnostic_scalar_field
               }
         }|
         ## Second Fluid
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SecondFluid" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Diffuse Interface
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DiffuseInterface" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BalancePressure" },
            (
               element prognostic {
                  # Note the assumptions about quadratic/linear below and hard-coding
                  # of the mesh, this is because of restrictions of the
                  # current code, will change in the near future.
                  ## Note that this is not the quadratic mesh balance pressure is
                  ## actually calculated on, but the linear mesh it is projected back
                  ## on for output purposes.
                  element mesh {
                     attribute name {  "VelocityMesh" }
                  },
                  prognostic_balance_pressure_field
               }
            )
         }|
         ## If enabled, decomposes Pressure by solving for the balanced part of 
         ## Pressure using a "geopressure" solver:
         ## 
         ##   f = - grad p_gp + g
         ##
         ## By choosing an appropriate mesh (typically velocity mesh order + 1)
         ## for the balanced part of pressure, physical balance can be
         ## represented to a higher degree of accuracy.
         ##
         ## Expanded BalancePressure field. As BalancePressure, but with
         ## additional options, including the ability to choose a general
         ## mesh for the geopressure solver.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GeostrophicPressure" },
            (
               element prognostic {
                  ## The GeostrophicPressure mesh
                  ## 
                  ## <b>WARNING: It is usual for this to be a higher degree
                  ## mesh than the velocity mesh</b>
                  element mesh {
                     attribute name { xsd:string },
                     comment
                  },
                  prognostic_geostrophic_pressure_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "VerticalBalancePressure" },
            (
               element prognostic {
                  ## This needs to be a quadratic DG mesh
                  mesh_choice,
                  prognostic_vertical_balance_pressure_field
               }
            )
         }|
         ## MaterialVolumeFraction field:
         ##
         ## Volume fraction of this material.
         ## Required in multimaterial simulations.
         ##  - if prognostic solves for the volume fraction
         ##  - if prescribed uses a specified volume fraction
         ##  - if diagnostic solves for the final material volume fraction
         ## Only 1 diagnostic MaterialVolumeFraction field allowed per
         ## simulation or solves for all the volume fractions based on
         ## the SumMaterialVolumeFractions field.
         ## 
         ## A diagnostic MaterialVolumeFraction field is currently required for
         ## compressible multimaterial simulations (even if only 1 material).
         ## Generally also requires a MaterialDensity field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialVolumeFraction" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field,
                  cap_option?,
                  surface_tension_option?
               }|
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field,
                  cap_option?
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field,
                  cap_option?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialDensity field:
         ##
         ## Field for the density of this material.
         ## Required in multimaterial simulations.
         ##  - prescribed if an incompressible simulation
         ##  - diagnostic if using a linear equation of state
         ##  - prognostic if a compressible simulation
         ## (note that if you set a multimaterial
         ## equation of state and this field is
         ## prognostic then its initial condition
         ## will be overwritten by the density that
         ## satisfies the initial pressure and
         ## the equation of state)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialDensity" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialInternalEnergy field:
         ##
         ## Field for the internal energy of this material.
         ## Required in multimaterial compressible simulations
         ## with full stiffened_gas (perfect gas) eos.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialInternalEnergy" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## SumMaterialVolumeFractions field:
         ##
         ## Sums the prognostic MaterialVolumeFraction fields.
         ## - diagnostic: sums all the volume fractions in the other
         ##   material phases
         element scalar_field {
            attribute rank { "0" },
            attribute name { "SumMaterialVolumeFractions" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CopiedField - This field copies the previous timesteps
         ## values from another (specified) field at every iteration
         ## and then solves the field using different (again, specified)
         ## scheme and solution options.
         ## For instance, this field can be used to create a diffused
         ## field to adapt to.
         ## Unless someone requests otherwise this is only currently possible
         ## for fields within the same material_phase.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CopiedField" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  attribute copy_from_field { string },
                  prognostic_scalar_field
               }
            )
         }|
         ## Calculate the stream function of 2D incompressible flow. Note 
         ## that this *only* makes sense for proper 2D (not pseudo-2D) simulations.
         ## Requires a continuous mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "StreamFunction" },
            (
               element prognostic {
                  mesh_choice,
                  prognostic_stream_function_field
               }
            )
         }|
         ## Phytoplankton
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Phytoplankton" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Zooplankton
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Zooplankton" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Nutrient
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Nutrient" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|
         ## Detritus
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Detritus" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|

         ## PhaseVolumeFraction:
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PhaseVolumeFraction" },
            (
               element prognostic {
                  velocity_mesh_choice,
                  prognostic_scalar_field
               }|
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }
            )
         }|

         # Insert new prognostic scalar fields here using the template:
         #        element scalar_field {
         #            attribute rank { "0" },
         #            attribute name { "NewFieldName" },
         #            (
         #               element prognostic {
         #                  velocity_mesh_choice,
         #                  prognostic_scalar_field
         #               }|
         #               element prescribed {
         #                  velocity_mesh_choice,
         #                  prescribed_scalar_field
         #               }|
         #               element aliased {
         #                  generic_aliased_field
         #               }
         #            )
         #        }
         
# -- Second is a list of fields that are primarily prescribed,
#    but can be aliased. An example is wind velocity.
# -- The list is in order of most frequently used.

         ## Prescribed scalar fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DistanceToSideBoundaries" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialFrictionAngle for multimaterial
         ## plasticity problems
         ##
         ## Requires a diagnostic bulk FrictionAngle field
         ## - not tested yet
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialFrictionAngle" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## MaterialCohesion for multimaterial
         ## plasticity problems
         ##
         ## Requires a diagnostic bulk Cohesion field
         ## - not tested yet
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialCohesion" },
            (
               element prescribed {
                  velocity_mesh_choice,
                  prescribed_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
#
# Insert new prescribed scalar fields here using the template:
#        element scalar_field {
#            attribute rank { "0" },
#            attribute name { "NewFieldName" },
#            (
#               element prescribed {
#                  velocity_mesh_choice,
#                  prescribed_scalar_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
#
# -- Last is a list of fields that are primarily diagnostic,
#    but can be aliased. An example is Tidal Range.
# -- The list is in order of most frequently used.
#
         ## Diagnostic scalar fields below this
         element ___Diagnostic_Fields_Below___ {
            empty
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## ControlVolumeDivergence:
         ##
         ## div field
         ##
         ## Divergence of the velocity field where
         ## the divergence operator is defined using
         ## the control volume C^T matrix.
         ## This assumes that the test space is discontinuous
         ## control volumes.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ControlVolumeDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  diagnostic_cv_divergence_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## FiniteElementDivergence:
         ##
         ## div field
         ##
         ## Divergence of the velocity field where
         ## the divergence operator is defined using
         ## the finite element C^T matrix.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "FiniteElementDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  element integrate_divergence_by_parts {
                     empty
                  }?,
                  diagnostic_fe_divergence_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Diffusive dissipation
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DiffusiveDissipation" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Richardson Number:
         ##
         ##  Ri = \frac{N^2}{(\frac{\partial u}{\partial z})^2 + (\frac{\partial u}{\partial z})^2}
         ## with 
         ##  N^2 = -\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
         ##
         ## Limitations:
         ##  - Gravity must be constant.
         ##  - Assumes gravity is in -ve final coordinate direction.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "RichardsonNumber" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_richardson_number_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CFLNumber
         ##
         ## See http://amcg.ese.ic.ac.uk/index.php?title=Local:Diagnostics#CFL_Number
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## ControlVolumeCFLNumber
         ##
         ## Courant Number as defined on a control volume mesh
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ControlVolumeCFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## DG_CourantNumber
         ##
         ## Courant Number as defined on a DG mesh
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "DG_CourantNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## CVMaterialDensityCFLNumber
         ##
         ## Courant Number as defined on a control volume mesh and
         ## incorporating the MaterialDensity.
         ## Requires a MaterialDensity field!
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CVMaterialDensityCFLNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## This scalar field is meant to replace DENTRAF.
         ## Basically, if you use new options, DENTRAF is no longer needed
         ## No repointing is done from this field to DENTRAF.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "CopyofDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }
            )
         }|
         ## add a MaterialVolume scalar_field to calculate the spatially varying 
         ## volume of a material (requires a MaterialVolumeFraction)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialVolume" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## add a MaterialMass scalar_field to calculate the spatially varying 
         ## mass of a material (requires a MaterialVolumeFraction and a MaterialDensity)
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialMass" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the MaterialDensity based on the bulk Pressure
         ## (and MaterialInternalEnergy if appropriate) for the equation
         ## of state of this material.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialEOSDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the MaterialPressure based on the MaterialDensity
         ## (and MaterialInternalEnergy if appropriate) for the equation
         ## of state of this material.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "MaterialPressure" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculates the BulkMaterialPressure based on the MaterialDensity
         ## and MaterialVolumeFraction (and MaterialInternalEnergy if appropriate) 
         ## for the equation of state of all materials.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BulkMaterialPressure" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Grid Reynolds number
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GridReynoldsNumber" },
            (element diagnostic {
               internal_algorithm,
               element mesh {
                  attribute name { "VelocityMesh" }
               },
               diagnostic_scalar_field
            }
            | element aliased { generic_aliased_field })
         }|
         ## GridPecletNumber
         ##
         ## Peclet Number Pe = U*dx/2*diffusivity
         ##
         ## Also see the test case 'grid_peclet_number'
         ## if you wish to see the effect of changing the 
         ## diffusivity on a 1D, cg-discretised tracer-field
         ##
         ## Adapting to this field is not recommended
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GridPecletNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  ## Mesh on which to calculate dx
                  mesh_choice,
                  ## This is the name of the scalar field
                  ## to calculate the Peclet number for
                  ## Note this field needs to have a diffusivity
                  element field_name { string },
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Horizontal velocity divergence:
         ##
         ## div_H velocity
         ##
         ## Uses the gravity field direction to determine the horizontal plane.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HorizontalVelocityDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Velocity divergence:
         ##
         ## div velocity
         ##
         element scalar_field {
            attribute rank { "0" },
            attribute name { "VelocityDivergence" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Kinetic energy density:
         ##
         ##  1/2 rho_0*|u|^2
         ##
         ## where rho_0 is the (reference) density 
         ##
         ## Limitations:
         ##  - The Density, PerturbationDensity, KineticEnergyDensity and Velocity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "KineticEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Gravitational potential energy density:
         ##
         ## rho_0*(1.0 + rho')*(g dot (r - r_0))
         ##
         ## where rho_0 is the (reference) density, rho' is the perturbation density and r_0 is the potential energy zero point.
         ##
         ## Limitations:
         ##  - Requires a constant gravity direction.
         ##  - The Density, PerturbationDensity and GravitationalPotentialEnergyDensity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GravitationalPotentialEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field,
                  ## Coordinate of a point with a potential energy of zero.
                  element zero_point {
                     real_dim_vector
                  }
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Isopycnal coordinate
         ##
         ##  z_star(x,t) = 1/A int_V' H(rho(x',t)-rho(x,t)) dV'
         ##
         ## where rho is the density, A is the width/area of the domain
         ##
         ## Limitations:
         ##  - You need to specify a (fine) mesh to redistribute the PerturbationDensity onto
         ##  - Requires a constant gravity direction.
         ##  - The Density, PerturbationDensity and GravitationalPotentialEnergyDensity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "IsopycnalCoordinate" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  ## This is the mesh onto which we redistribute the PerturbationDensity
                  element fine_mesh {
                     attribute name { string }
                  },
                  diagnostic_scalar_field                 
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Background potential energy density:
         ##
         ## PE_b = rho*z_star
         ##
         ## where rho is the density, z_star is the isopycnal coordinate
         ##
         ## Limitations:
         ##  - Requires a constant gravity direction.
         ##  - The Density, PerturbationDensity and
         ##  GravitationalPotentialEnergyDensity fields must be on the
         ##  same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "BackgroundPotentialEnergyDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field                 
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Ertel potential vorticity:
         ##
         ##  (f + curl u) dot grad rho'
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PotentialVorticity" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Relative potential vorticity:
         ##
         ##   curl u dot grad rho'
         element scalar_field {
            attribute rank { "0" },
            attribute name { "RelativePotentialVorticity" },
            attribute depends { "Velocity,PerturbationDensity" },
            (
              element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Calculate the horizontal stream function psi where:
         ##   \partial_x \psi = -v
         ##   \partial_y \psi = u
         ## where u and v are perpendicular to the gravity direction. Applies a
         ## strong Dirichlet boundary condition of 0 on all boundaries.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HorizontalStreamFunction" },
            attribute depends { "Velocity" },
            (
              element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  ## Solver
                  element solver {
                     linear_solver_options_sym
                  },
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Speed:
         ##
         ##  |u|
         ##
         ## Limitations:
         ##  - The Speed and Velocity fields must be on the same mesh.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "Speed" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Absolute Difference between two scalar fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "AbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_scalar_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?,
                  ## Ignore boundary nodes (i.e. zero them when calculating the difference)
                  element ignore_boundaries {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Absolute Difference between two scalar fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ScalarAbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_scalar_field,
                  ## Evaluate the absolute difference once the average difference has been removed?
                  element relative_to_average {
                    empty
                  }?,
                  ## Ignore boundary nodes (i.e. zero them when calculating the difference)
                  element ignore_boundaries {
                    empty
                  }?
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Galerkin projection of one field onto another mesh.
         ##
         ## The field must be in this material_phase.
         ## 
         ## NOTE: you need the solver options if the mesh
         ## of this field is continuous.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "GalerkinProjection" },
            (
               element diagnostic {
                  internal_algorithm,
                  element source_field_name { string },
                  mesh_choice,
                  ## Lump the mass matrix of the galerkin projection
                  ## less accurate but faster and might give smoother result.                  
                  element lump_mass {
                     empty
                  }?,
                  element solver {
                    linear_solver_options_sym
                  }?,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Primary production of Phytoplankton. This is calculated by
         ## the ocean biology module and will not be calculated unless
         ## ocean biology is being simulated.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PrimaryProduction" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Grazing of Phytoplankton by Zooplankton. This is calculated by
         ## the ocean biology module and will not be calculated unless
         ## ocean biology is being simulated.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "PhytoplanktonGrazing" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )

         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "TidalRange" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field_tidal_range
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeM2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseM2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeS2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseS2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeN2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseN2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeK2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseK2" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeK1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseK1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeO1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseO1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeP1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseP1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeQ1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseQ1" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeMf" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseMf" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeMm" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicPhaseMm" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         element scalar_field {
            attribute rank { "0" },
            attribute name { "HarmonicAmplitudeSSa" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Output the universal numbering of the mesh on which this field is based.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "UniversalNumber" },
            (
               element diagnostic {
                  internal_algorithm,
                  mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Output the processors which own the nodes of the mesh on which this field is based.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "NodeOwner" },
            (
               element diagnostic {
                  internal_algorithm,
                  mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Output the processors which own the elements of the mesh on which this field is based.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "ElementOwner" },
            (
               element diagnostic {
                  internal_algorithm,
                  mesh_choice,
                  diagnostic_scalar_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }
                  
# Insert new diagnostic scalar fields here using the template:
#        element scalar_field {
#            attribute rank { "0" },
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  velocity_mesh_choice,
#                  diagnostic_scalar_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )

# This is the choice of additional vector field to be solved for
vector_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
# Prognostic vector fields are not possible (other than velocity and those known fields below).
      (
         ## Generic field variable (vector)
         element vector_field {
            attribute rank { "1" },
            attribute name { xsd:string },
            ## Field type
            (
               element prescribed {
                  mesh_choice,
                  prescribed_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }|
               element diagnostic {
                  vector_diagnostic_algorithms,
                  velocity_mesh_choice,
                  python_diagnostic_field_code?,
                  diagnostic_vector_field
               }
            )
         }|
#
# -- List of fields that are primarily prognostic,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Prescribed vector fields below this
         element ___Prognostic_fields_below___ {
            empty
         }|

#
# -- List of fields that are primarily prescribed,
#    but can be aliased. An example is Maximum bed shear stress.
# -- The list is in order of most frequently used.
#
         ## Prescribed vector fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|

         ## Gradient of a scalar field evaluated using the C gradient
         ## matrix constructed using finite elements.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "FiniteElementGradient" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  mesh_choice,
                  element integrate_gradient_by_parts {
                     empty
                  }?,
                  diagnostic_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Gradient of a scalar field evaluated using the transpose
         ## of the C^T divergence matrix constructed using finite
         ## elements.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "FiniteElementDivergenceTransposed" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  mesh_choice,
                  element integrate_divergence_by_parts {
                     empty
                  }?,
                  diagnostic_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Relative vorticity field - curl of the velocity field
         element vector_field {
            attribute rank { "1" },
            attribute name { "Vorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  ### Relative vorticity
                  #element algorithm {
                  #   attribute name { "curl" },
                  #   attribute material_phase_support { "single" },
                  #   attribute source_field_name { "Velocity" }
                  #},
                  element mesh {
                     attribute name { "VelocityMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|        
         ## Planetary vorticity
         ##
         ## Limitations:
         ##  - Requires geometry dimension of 3.
         element vector_field {
            attribute rank { "1" },
            attribute name { "PlanetaryVorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Absolute vorticity:
         ##
         ##   f + curl u
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         element vector_field {
            attribute rank { "1" },
            attribute name { "AbsoluteVorticity" },
            attribute depends { "Velocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Gradient of a scalar field evaluated using the transpose
         ## of the C^T matrix constructed using control volumes.
         ## Field must be in this material_phase.
         element vector_field {
            attribute rank { "1" },
            attribute name { "ControlVolumeDivergenceTransposed" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name { string },
                  velocity_mesh_choice,
                  diagnostic_cv_gradient_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Full velocity in an
         ## inner element SGS treatment of momentum
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "InnerElementFullVelocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "InnerElementMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Vorticity of the full velocity in an
         ## inner element SGS treatment of momentum
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "InnerElementFullVorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "InnerElementMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Vorticity of the SGS velocity in an
         ## inner element SGS treatment of momentum
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "InnerElementVorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "InnerElementMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## The continuous solution mapped to a discontinuous mesh
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "DgMappedVelocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "InnerElementMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Vorticity of the DG mapped Velocity
         ## Note vorticity is actually calculated over a DG field
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "DgMappedVorticity" },
            (
               element diagnostic {
                  internal_algorithm,
                  element mesh {
                     attribute name {  "InnerElementMesh" }
                  },
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## LinearMomentum field.
         ##  p = \rho*u 
         ## (where p is the linear momentum, \rho the density and u the velocity)
         element vector_field {
             attribute rank { "1" },
             attribute name { "LinearMomentum" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice,
                   diagnostic_vector_field
                }|
                element aliased {
                   generic_aliased_field
                }
             )
         }|

         ## Absolute Difference between two vector fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element vector_field {
            attribute rank { "1" },
            attribute name { "AbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         ## Absolute Difference between two vector fields.
         ##
         ## Both fields must be in this material_phase.
         ## Assumes both fields are on the same mesh as the AbsoluteDifference field.
         element vector_field {
            attribute rank { "1" },
            attribute name { "VectorAbsoluteDifference" },
            (
               element diagnostic {
                  internal_algorithm,
                  attribute field_name_a { string },
                  attribute field_name_b { string },
                  mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         
         ## Bed Shear Stress = density*drag_coeff*|u|*u
         ##
         ## at the moment this assumes the density and drag coefficients are constants.
         ## This diagnostic vector field is only calculated over surface elements/nodes, 
         ## interior nodes will have zero value.
         element vector_field {
             attribute rank { "1" },
             attribute name { "BedShearStress" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice,
                   diagnostic_vector_field_bed_shear_stress
                }|
                element aliased {
                   generic_aliased_field
                }
             )
         }|
         ## Max Bed Shear Stress.
         ##
         ## Note that you need BedShearStress turned on for this to work.
         element vector_field {
             attribute rank { "1" },
             attribute name { "MaxBedShearStress" },
             (
                element diagnostic {
                   internal_algorithm,
                   velocity_mesh_choice
                   #diagnostic_vector_field_max_bed_shear_stress
                }|
                element aliased {
                   generic_aliased_field
                }
             ),
            ## This is the time after which the max operator is
            ## applied to the bed shear stress.
            element spin_up_time {
               real
            }
         }|
         ## Galerkin projection of one field onto another mesh.
         ##
         ## The field must be in this material_phase.
         ## 
         ## NOTE: you need the solver options if the mesh
         ## of this field is continuous.
         element vector_field {
            attribute rank { "1" },
            attribute name { "GalerkinProjection" },
            (
               element diagnostic {
                  internal_algorithm,
                  element source_field_name { string },
                  mesh_choice,
                  ## Lump the mass matrix of the galerkin projection
                  ## less accurate but faster and might give smoother result.                  
                  element lump_mass {
                     empty
                  }?,                  
                  element solver {
                  linear_solver_options_sym
                  }?,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Compute the imbalanced component of velocity,
         ## that is,
         ## u - u_bal
         ## where u_bal is the velocity that puts the state in geostrophic
         ## balance.
         ## Note: if your VelocityMesh is continuous, then the solver option
         ## is necessary.
         element vector_field {
            attribute rank { "1" },
            attribute name { "ImbalancedVelocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  element solver {
                  linear_solver_options_sym
                  }?,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|
         ## Compute the balanced component of velocity,
         ## that is, the velocity that puts the state in geostrophic
         ## balance.
         ## This diagnostic depends on ImbalancedVelocity.
         element vector_field {
            attribute rank { "1" },
            attribute name { "BalancedVelocity" },
            (
               element diagnostic {
                  internal_algorithm,
                  velocity_mesh_choice,
                  diagnostic_vector_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }

# Insert new diagnostic vector field here using the template:
#        element vector_field {
#            attribute rank { "1" },
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  mesh_choice,
#                  diagnostic_vector_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )  

# This is the choice of additional tensor fields
tensor_field_choice =
   (
# The first is a generic field, which may be used for any user-defined field
# that FLUIDITY knows nothing about, or a generic diagnostic
# Prognostic tensor fields are not possible.
      (
         ## Generic field variable (tensor)
         element tensor_field {
            attribute rank { "2" },
            attribute name { xsd:string },
            ## Field type
            (
               element prescribed {
                  mesh_choice,
                  prescribed_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }|
               element diagnostic {
                  tensor_diagnostic_algorithms,
                  velocity_mesh_choice,
                  python_diagnostic_field_code?,
                  diagnostic_tensor_field
               }
            )
         }|
#
# -- Second is a list of tensor fields that are primarily prescribed,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Prescribed scalar fields below this
         element ___Prescribed_fields_below___ {
            empty
         }|

         ## MaterialViscosity field:
         ##
         ## Field for the viscosity of this material.
         ## Required if using a diagnostic bulk viscosity
         ## in a multimaterial simulation.
         element tensor_field {
            attribute rank { "2" },
            attribute name { "MaterialViscosity" },
            (
               element prescribed {
                  mesh_choice,
                  prescribed_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

         element tensor_field {
            attribute rank { "2" },
            attribute name { "MaterialElasticity" },
            (
               element prescribed {
                  mesh_choice,
                  prescribed_tensor_field
               }|
               element aliased {
                  generic_aliased_field
               }
            )
         }|

#
# Insert new prescribed tensor fields here using the template:
#        element tensor_field {
#            attribute rank { "2" },
#            attribute name { "NewFieldName" },
#            (
#               element prescribed {
#                  mesh_choice,
#                  prescribed_tensor_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }|
#
# -- Last is a list of fields that are primarily diagnostic,
#    but can be aliased.
# -- The list is in order of most frequently used.
#
         ## Diagnostic tensor fields below this
         element ___Diagnostic_Fields_Below___ {
            empty
         }

# Insert new diagnostic tensor field here using the template:
#        element tensor_field {
#            attribute name { "NewFieldName" },
#            (
#               element diagnostic {
#                  internal_algorithm,
#                  mesh_choice,
#                  diagnostic_tensor_field
#               }|
#               element aliased {
#                  generic_aliased_field
#               }
#            )
#        }
      )
   )

mesh_info =
   (
      ## Read mesh from file.
      element from_file {
         (
            ## Triangle mesh format.
            ##
            ## Enter the base name without the .edge .ele, .face or
            ## .node extensions, and without process numbers.
            element format {
              attribute name { "triangle" },
              # string_value elements are used only for backwards compatibility - any new format choices should NOT have these
              element string_value {
                 "triangle"
              },
              comment
            }|
            ## Read the mesh from a vtu. Note that the mesh will have no surface
            ## or region IDs.
            element format {
              attribute name { "vtu" },
              comment
            }
         ),
         attribute file_name { xsd:string },
         from_file_mesh_stat_options,
         comment
      }|
      ## Make mesh from existing mesh. The existing mesh cannot itself
      ## be made from an existing mesh (i.e. it must be read from a
      ## file).
      element from_mesh {
         mesh_choice,
         element mesh_shape {
            element polynomial_degree {
               integer
            }
         }?,
         element mesh_continuity {
            element string_value{
               "continuous" | "discontinuous"
            }
         }?,
         ## Make mesh periodic
         element periodic_boundary_conditions {
            attribute name { xsd:string },
            ## List of boundary ids that are aliased to
            element physical_boundary_ids {
               integer_vector
            },
            ## List of boundary ids that are aliased
            element  aliased_boundary_ids {
               integer_vector
            },
            ## Python code which takes coordinate of an aliased
            ## boundary node and returns the coordinate of a physical
            ## boundary node
            element coordinate_map {
               python_code
            }
         }*,
         ## Extrude a horizontal (1D or 2D) mesh in the vertical
         element extrude {
            ## Depth over which to extrude
            ## top will be at z=0
            ## bottom will be at z=-bottom_depth
            element bottom_depth {
               real
            },
            ## Constant or function to specify the depth of the
            ## layers. The function is a function of all coordinates
            ## (so in 2+1D: x,y and z) to specify a layer depth that
            ## varies both in the horizontal as in the vertical.
            element sizing_function {
               input_choice_real
            },
            ## surface_id to assign to the top of the extruded mesh
            element top_surface_id {
               integer
            }?,
            ## surface_id to assign to the bottom of the extruded mesh
            element bottom_surface_id {
               integer
            }?
         }?,
         derived_mesh_stat_options,
         comment
      }
   )

# Options for inclusion/exclusion of mesh statistics from the .stat file
include_mesh_in_stat =
   (
      ## Include this mesh in the .stat file
      element include_in_stat {
         comment
      }
   )
exclude_mesh_from_stat =
   (
      ## Exclude this mesh from the .stat file
      element exclude_from_stat {
         comment
      }
   )

# Diagnostic statistics options for meshed, with enabled by default
mesh_stat_options_enabled_default = include_mesh_in_stat
mesh_stat_options_enabled_default |= exclude_mesh_from_stat

# Diagnostic statistics options for meshed, with disabled by default
mesh_stat_options_disabled_default = exclude_mesh_from_stat
mesh_stat_options_disabled_default |= include_mesh_in_stat

from_file_mesh_stat_options =
   (
      ## Specify what is added to .stat files
      element stat {
         mesh_stat_options_enabled_default
      }
   )
derived_mesh_stat_options =
   (
      ## Specify what is added to .stat files
      element stat {
         mesh_stat_options_disabled_default
      }
   )

linear_solver_options_asym =
   (
      ## Iterative (Krylov) method to solve the linear discretised equation
      ## Given are the most frequently used methods. The solution is done
      ## by the PETSc library. Many more methods are provided.
      ##
      (
         ## GMRES
         ##
         ## Your safest bet for non-symmetric systems.
         element iterative_method {
            attribute name { "gmres" },
            ## Restart value for gmres iteration
            ## Higher values give better convergence but require more memory.
            ## Suggested value: 30
            element restart {
               integer
            }
         }|
         ## Conjugate gradient method
         ##
         ## Only works for symmetric systems.
         element iterative_method {
            attribute name { "cg" }
         }|
         ## Direct method
         ##
         ## This is for non-iterative methods
         ## Only makes sense in combination with preconditioners that do a complete solve, e.g. lu.
         element iterative_method {
            attribute name { "preonly" }
         }|
         ## Richardson iteration
         ##
         ## Only apply preconditioner each iteration, no krylov acceleration
         element iterative_method {
            attribute name { "richardson" }
         }|
         ## Other methods
         ##
         ## Any method provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html
         ## (available methods may depend on the PETSc library installed on your system)
         element iterative_method {
            attribute name { xsd:string }
         }
      ),
      ## Preconditioner to be used in combination with the iterative method.
      (
         ## Succesive Over-Relaxation
         ##
         ## This includes SSOR (symmetric sor)
         element preconditioner {
            attribute name { "sor" }
         }|
         ## The Eisenstat method
         ##
         ## This preconditioner is equivalent to SOR but only uses
         ## half the number of flops,
         ## i.e. same convergence rate but twice as fast per
         ## iteration. Because it computes
         ## a different preconditioned residual the convergence in
         ## practice may be quite different though.
         ##
         ## It does not work in parallel!
         element preconditioner {
            attribute name { "eisenstat" }
         }|
         ## Incomplete LU decomposition
         element preconditioner {
            attribute name { "ilu" }
         }|
         ## LU direct solver
         ##
         ## This performs a complete, direct solve of the equation and should only be used in combination with preonly as iterative method.
         element preconditioner {
            attribute name { "lu" }
         }|
         ## Fluidity`s own multigrid method
         ##
         ## Especially suited for ill-conditioned, large aspect ratio problems.
         element preconditioner {
            attribute name { "mg" }
         }|
         ## Prometheus multigrid method
         element preconditioner {
            attribute name { "prometheus" }
         }|
         ## Hypre preconditioners (includes boomeramg)
         element preconditioner {
            attribute name { "hypre" },
            (
               ## BoomerAMG multigrid method
               element hypre_type {
                  attribute name { "boomeramg" }
               }|
               ## Other Hypre preconditioners
               element hypre_type {
                  attribute name { string }
               }
            )
         }|
         ## Complete solve
         ##
         ## This only makes sense for solves where a different approximated preconditioner
         ## matrix is provided. For instance when solving pressure with the 
         ## option full_schur_complement and using a masslumped schur complement
         ## as preconditioner matrix.
         ##
         ## NOTE: If you are using a krylov method (cg/gmres) for this preconditioner 
         ## solve you either need to set your tolerances much stricter for it
         ## than in the outer solve (so that the preconditioner is close to an
         ## exact matrix inversion), or use fgmres in the outer solve.
         element preconditioner {
            attribute name { "ksp" },
            (
                ## Solver options for the full solve done by this preconditioner
                element solver {
                   pc_ksp_solver_options
                }
            )
         }|
         ## Other preconditioners
         ##
         ## Any preconditioner provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/PC/PCType.html
         ## (available preconditiors may depend on the PETSc library installed on your system)
         element preconditioner {
            attribute name { string }
         }
      ),
      generic_solver_options
   )
   
linear_solver_options_sym =
   (
      ## Iterative (Krylov) method to solve the linear discretised equation
      ## Given are the most frequently used methods. The solution is done
      ## by the PETSc library. Many more methods are provided.
      ##
      (
         ## Conjugate gradient method
         ##
         ## Only works for symmetric systems.
         element iterative_method {
            attribute name { "cg" }
         }|
         ## GMRES
         ##
         ## Your safest bet for non-symmetric systems.
         element iterative_method {
            attribute name { "gmres" },
            ## Restart value for gmres iteration
            ## Higher values give better convergence but require more memory.
            ## Suggested value: 30
            element restart {
               integer
            }
         }|
         ## Direct method
         ##
         ## This is for non-iterative methods
         ## Only makes sense in combination with preconditioners that do a complete solve, e.g. lu.
         element iterative_method {
            attribute name { "preonly" }
         }|
         ## Richardson iteration
         ##
         ## Only apply preconditioner each iteration, no krylov acceleration
         element iterative_method {
            attribute name { "richardson" }
         }|
         ## Other methods
         ##
         ## Any method provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html
         ## (available methods may depend on the PETSc library installed on your system)
         element iterative_method {
            attribute name { xsd:string }
         }
      ),
      ## Preconditioner to be used in combination with the iterative method.
      (
         ## Succesive Over-Relaxation
         ##
         ## This includes SSOR (symmetric sor)
         element preconditioner {
            attribute name { "sor" }
         }|
         ## The Eisenstat method
         ##
         ## This preconditioner is equivalent to SOR but only uses
         ## half the number of flops,
         ## i.e. same convergence rate but twice as fast per
         ## iteration. Because it computes
         ## a different preconditioned residual the convergence in
         ## practice may be quite different though.
         ## It does not work in parallel!
         element preconditioner {
            attribute name { "eisenstat" }
         }|
         ## Incomplete LU decomposition
         element preconditioner {
            attribute name { "ilu" }
         }|
         ## Incomplete Cholesky decomposition (only works for symmetric matrices)
         element preconditioner {
            attribute name { "icc" }
         }|
         ## LU direct solver
         ##
         ## This performs a complete, direct solve of the equation and should only be used in combination with preonly as iterative method.
         element preconditioner {
            attribute name { "lu" }
        }|
         ## Fluidity`s own multigrid method
         ##
         ## Especially suited for ill-conditioned, large aspect ratio problems.
         element preconditioner {
            attribute name { "mg" },
            ## apply vertical lumping from the full mesh to the surface mesh
            ## as the first coarsening step instead of the default
            ## aggregation method.
            element vertical_lumping {
               ## Does additional smoothing by solving the equation but with
               ## a dirichilet boundary condition on top given by the last iteration
               ## of the multigrid cycle. May be quite expensive per iteration
               ## but improves the solution quite a lot for difficult meshes.
               element internal_smoother {
                  empty
               }?
            }?
         }|
         ## Prometheus multigrid method
         element preconditioner {
            attribute name { "prometheus" }
         }|
         ## Hypre preconditioners (includes boomeramg)
         element preconditioner {
            attribute name { "hypre" },
            (
               ## BoomerAMG multigrid method
               element hypre_type {
                  attribute name { "boomeramg" }
               }|
               ## Other Hypre preconditioners
               element hypre_type {
                  attribute name { string }
               }
            )
         }|
         ## Complete solve
         ##
         ## This only makes sense for solves where a different approximated preconditioner
         ## matrix is provided. For instance when solving pressure with the 
         ## option full_schur_complement and using a masslumped schur complement
         ## as preconditioner matrix.
         ##
         ## NOTE: If you are using a krylov method (cg/gmres) for this preconditioner 
         ## solve you either need to set your tolerances much stricter for it
         ## than in the outer solve (so that the preconditioner is close to an
         ## exact matrix inversion), or use fgmres in the outer solve.
         element preconditioner {
            attribute name { "ksp" },
            (
                ## Solver options for the full solve done by this preconditioner
                element solver {
                   pc_ksp_solver_options
                }
            )
         }|
         ## Other preconditioners
         ##
         ## Any preconditioner provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/PC/PCType.html
         ## (available preconditiors may depend on the PETSc library installed on your system)
         element preconditioner {
            attribute name { string }
         }
      ),
      generic_solver_options
   )

# this is a copy linear_solver_options_sym, but with preconditioner "ksp"
# removed to avoid infinite recursion
pc_ksp_solver_options =
   (
      ## Iterative (Krylov) method to solve the linear discretised equation
      ## Given are the most frequently used methods. The solution is done
      ## by the PETSc library. Many more methods are provided.
      ##
      (
         ## Conjugate gradient method
         ##
         ## Only works for symmetric systems.
         element iterative_method {
            attribute name { "cg" }
         }|
         ## GMRES
         ##
         ## Your safest bet for non-symmetric systems.
         element iterative_method {
            attribute name { "gmres" },
            ## Restart value for gmres iteration
            ## Higher values give better convergence but require more memory.
            ## Suggested value: 30
            element restart {
               integer
            }
         }|
         ## Direct method
         ##
         ## This is for non-iterative methods
         ## Only makes sense in combination with preconditioners that do a complete solve, e.g. lu.
         element iterative_method {
            attribute name { "preonly" }
         }|
         ## Richardson iteration
         ##
         ## Only apply preconditioner each iteration, no krylov acceleration
         element iterative_method {
            attribute name { "richardson" }
         }|
         ## Other methods
         ##
         ## Any method provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html
         ## (available methods may depend on the PETSc library installed on your system)
         element iterative_method {
            attribute name { xsd:string }
         }
      ),
      ## Preconditioner to be used in combination with the iterative method.
      (
         ## Succesive Over-Relaxation
         ##
         ## This includes SSOR (symmetric sor)
         element preconditioner {
            attribute name { "sor" }
         }|
         ## The Eisenstat method
         ##
         ## This preconditioner is equivalent to SOR but only uses
         ## half the number of flops,
         ## i.e. same convergence rate but twice as fast per
         ## iteration. Because it computes
         ## a different preconditioned residual the convergence in
         ## practice may be quite different though.
         ## It does not work in parallel!
         element preconditioner {
            attribute name { "eisenstat" }
         }|
         ## Incomplete LU decomposition
         element preconditioner {
            attribute name { "ilu" }
         }|
         ## Incomplete Cholesky decomposition (only works for symmetric matrices)
         element preconditioner {
            attribute name { "icc" }
         }|
         ## LU direct solver
         ##
         ## This performs a complete, direct solve of the equation and should only be used in combination with preonly as iterative method.
         element preconditioner {
            attribute name { "lu" }
        }|
         ## Fluidity`s own multigrid method
         ##
         ## Especially suited for ill-conditioned, large aspect ratio problems.
         element preconditioner {
            attribute name { "mg" },
            ## apply vertical lumping from the full mesh to the surface mesh
            ## as the first coarsening step instead of the default
            ## aggregation method.
            element vertical_lumping {
               ## Does additional smoothing by solving the equation but with
               ## a dirichilet boundary condition on top given by the last iteration
               ## of the multigrid cycle. May be quite expensive per iteration
               ## but improves the solution quite a lot for difficult meshes.
               element internal_smoother {
                  empty
               }?
            }?
         }|
         ## Prometheus multigrid method
         element preconditioner {
            attribute name { "prometheus" }
         }|
         ## Hypre preconditioners (includes boomeramg)
         element preconditioner {
            attribute name { "hypre" },
            (
               ## BoomerAMG multigrid method
               element hypre_type {
                  attribute name { "boomeramg" }
               }|
               ## Other Hypre preconditioners
               element hypre_type {
                  attribute name { string }
               }
            )
         }|
         ## Other preconditioners
         ##
         ## Any preconditioner provided by the PETSc library
         ## http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/PC/PCType.html
         ## (available preconditiors may depend on the PETSc library installed on your system)
         element preconditioner {
            attribute name { string }
         }
      ),
      generic_solver_options
   )

generic_solver_options =
   (
      ## Relative error
      ##
      ## The solver finishes if the preconditioned error becomes smaller than the original preconditioned error times this value.
      ## Suggested value: 1.0e-7
      element relative_error {
         real
      },
      ## Absolute error bound
      ##
      ## The solver finishes if the preconditioned error becomes smaller than this value.
      element absolute_error {
         real
      }?,
      ## Maximum number of iterations allowed in the linear solver
      ## before giving up.
      element max_iterations {
         integer
      },
      ## Switch on to not use an initial guess from a previous solve but
      ## start with a zero vector. Note that some of the solves always
      ## start at zero in which case this switch will have no effect (see the log output).
      element start_from_zero {
         empty
      }?,
      ## Remove Null-space from residual after applying preconditioner.
      ## This often leads to better convergence rates, when compared to
      ## imposing a reference_node to pin the solution.
      element remove_null_space {
         empty
      }?,
      ## Miscellaneous PETSc options
      ##
      ## Any options not included above can be set using the PETSc option syntax
      ## For example: -ksp_dtol 1e5 -ksp_monitor_draw
      element petsc_options {
         string
      }?,
      (
         ## Solver failures are always treated as fatal errors. The
         ## model stops at the end of the time step in order to allow
         ## for the latest output to be written.
         element never_ignore_solver_failures {
            empty
         }|
         ## Allow for an initial period in which solver failures
         ## caused by non-convergence in the maximum number of
         ## iterations are ignored.
         element ignore_non_convergence_during_spin_up {
            ## As long as current_time < spin_up_time, solver failures
            ## due to non-convergence in the maximum number of
            ## iterations are ignored. This might be used for spinning
            ## up the model. As there is no guarantee we're actually
            ## solving the flow equations to any accuracy, the results
            ## in this period should not be trusted.
            element spin_up_time {
               real
            }
         }|
         ## Ignore all solver failures. This is a dangerous option
         ## that should only be used in exceptional cases.
         element ignore_all_solver_failures {
            empty
         }
      ),
      ## Extra diagnostics to help debug solver problems
      element diagnostics {
         ## Print out the norm of vectors and matrices before the
         ## solve, and that of the solution vector afterwards.
         ## Norms are printed at verbosity level 2, so run fluidity with -v2 or -v3
         element print_norms {
            empty
         }?,
         ## Options to give extra information for each iteration of the
         ## the solve. Some of those may really slow down your computation!
         element monitors {
            ## Prints the preconditioned residual for each iteration of the solve.
            ## This is the error estimation PETSc uses during the solve.
            element preconditioned_residual {
               empty
            }?,
            ## Prints the "true" residual for each iteration of the solve,
            ## i.e. PETSc computes the L2-norm of r=A-bx. This may mean
            ## PETSc has to do extra computations.
            element true_residual {
               empty
            }?,
            ## Draws a graph over the convergence of the preconditioned residual
            ## during the solve. This option only works for systems where PETSc
            ## has been linked with the X library.
            element preconditioned_residual_graph {
               empty
            }?,
            ## Prints the error by computing the difference with the provided
            ## exact solution each time step.
            element true_error {
               ## Give the field name of the field that contains the exact
               ## solution to be compared with each iteration
               attribute exact_solution_field {
                  string
               }
            }?
         }
      }
   )


input_choice_tensor_field =
   (
      (
         ## An isotropic tensor, i.e.
         ## one with no directional variation.
         ## Can be represented as a scalar real.
         element isotropic {
            input_choice_real
         }|
         ## A symmetric tensor, i.e.
         ## A^T = A
         element anisotropic_symmetric {
            input_choice_real_dim_symmetric_tensor
         }|
         ## A general asymmetric tensor.
         element anisotropic_asymmetric {
            input_choice_real_dim_tensor
         }
      )
   )

constitutive_laws =
   (
      (
         ## Constitutive laws for fluids
         element constitutive_law {
            attribute name { "fluid" }
         }|
         ## Constitutive laws for solids
         element constitutive_law {
            attribute name { "solid" }
         }
      )
   )
   
region_ids = 
   ( 
      ## Optional region ids to associate different values
      ## to different regions of the mesh.
      ## Leave unselected if you`re not using multiple regions or
      ## region_ids.
      ## Currently only works with triangle files created by gmsh2triangle.
      element region_ids {
         integer_vector
      }?
   )

temporal_control_volume_options =
   (
      ## Temporal discretisation options that are only relevant if a control volume or mixed control volume - continuous galerkin spatial discretisation is selected for this field.
      element control_volumes {
         ## Number of iterations within an advection solve.
         ## This increases the accuracy of the face values and ensures that
         ## the pivoted solution is cancelled out.
         ## Defaults to 1 if unselected.
         element number_advection_iterations {
            integer,
            ## Cut short advection_iterations if the specified tolerance
            ## is reached.
            ## This only works for pure control volume discretisations.
            element tolerance {
               real,
               (
                  ## Select the norm with which you want the tolerance to be tested.
                  ##
                  ## The infinity norm.
                  element infinity_norm {
                    empty
                  }|
                  ## Select the norm with which you want the tolerance to be tested.
                  ##
                  ## The l2 norm.
                  element l2_norm {
                    empty
                  }|
                  ## Select the norm with which you want the tolerance to be tested.
                  ##
                  ## The l2 norm evaluated on a control volume mesh.
                  element cv_l2_norm {
                    empty
                  }
               )               
            }?
         }?,
         (
            ## Use timestep subcycling to solve this equation.
            ## Specify the maximum courant number per subcycle.
            ## This only works for pure control volume discretisations.
            element maximum_courant_number_per_subcycle {
               real,
               field_based_cfl_number_options
            }|
            ## Use timestep subcycling to solve this equation.
            ## Specify the number of subcycles.
            ## This only works for pure control volume discretisations.
            element number_advection_subcycles {
               integer
            }
         )?,
         ## Only works if a control volume or mixed control volume -
         ## continuous galerkin 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
         }?,
         ## Only works if a control volume or mixed control volume -
         ## continuous galerkin spatial discretisation is selected.
         ## Time discretisation of upwind discretisation off which the
         ## higher order solution is pivotted.
         ##  - pivot_theta = 1 - implicit pivot (default if not set and 
         ##                      best choice if not intentionally modifying
         ##                      scheme to be explicit)
         ##  - pivot_theta = 0 - explicit pivot
         element pivot_theta {
            real
         }?
      }
   )

temporal_discontinuous_galerkin_options = 
   (
      ## This enables DG-specific timestepping options, such as
      ## explicit advection subcycling. 
      element discontinuous_galerkin {
         (
            ## Use timestep subcycling to solve this equation.
            ## Specify the maximum courant number per subcycle.
            element maximum_courant_number_per_subcycle {
               real
            }|
            ## Use timestep subcycling to solve this equation.
            ## Specify the number of subcycles.
            element number_advection_subcycles {
               integer
            }
         )?
      }
   )

legacy_mixed_cv_cg_options = 
   (
      ## Use a mixed control volume - continuous galerkin
      ## discretisation.
      ##
      ## This uses a control volume discretisation of the advective
      ## terms and a continuous galerkin representation of the
      ## diffusion, source and absorption terms.
      ## Clearly, in terms of discretisation this does the same thing
      ## as a pure control volume discretisation if diffusion terms
      ## are zero.
      ## However it still follows a partially legacy code path.
      element legacy_mixed_cv_cg {
         spatial_control_volume_options,
         
         ## Allows legacy options to be set.
         ## BE WARNED: this probably means you're using inconsistent discretisations
         ## on different terms - make sure you're happy with what your code path is doing
         ## before selecting this.
         element legacy_disott {
            integer
         }?
      }
   )

pure_cv_options = 
   (
      ## Use a pure control volume discretisations.
      ## Follows a new control volume code path.
      element control_volumes {
         spatial_control_volume_options,
         (
            ## Use the gradient of the field constructed using the
            ## basis functions of the parent finite element mesh to
            ## form the divergence.
            ##
            ## DOES NOT CURRENTLY WORK WITH ROBIN OR WEAK DIRICHLET BOUNDARY CONDITIONS!
            ##
            ## Based on schemes in Handbook of Numerical Analysis,
            ## P.G. Ciarlet, J.L. Lions eds, vol 7, pp 713-1020
            element diffusion_scheme {
              attribute name{"ElementGradient"}
            }|
            ## Use an auxiliary gradient equation to find the gradient of the field.
            ##
            ## DOES NOT CURRENTLY WORK WITH ROBIN BOUNDARY CONDITIONS!
            ##
            ## Based on scheme proposed in Bassi, F. & Rebay, S., A
            ## high-order accurate discontinuous finite element method
            ## for the numerical solution of the compressible
            ## Navier-Stokes equations, Journal Of Computational
            ## Physics, 1997, 131, 267-279
            element diffusion_scheme {
              attribute name{"BassiRebay"}
            }
         )
      }
   )

coupled_cv_options = 
   (
      ## Use a pure control volume discretisations with face value
      ## restrictions between different fields in different material_phases.
      ##
      ## THIS DOES NOT WORK WITH DIFFUSION!
      ##
      ## Follows a new control volume code path.
      element coupled_cv {
         coupled_spatial_control_volume_options,
         ## Set the maximum and minimum bounds for the sum up to and including this field.
         ## This defines the limiter used to enforce boundedness on this field.
         element parent_sum {
            element target_maximum {
              real
            },
            element target_minimum {
              real
            }
         }
      }
   )

spatial_control_volume_options = standard_control_volume_options
spatial_control_volume_options |= compressive_control_volume_options

standard_control_volume_options = 
   (
      ## First Order Upwind face value discretisation
      ##  face_value = donor_value, 
      ## where
      ##  donor_value = income*val_1 + (1.-income)*val_2, 
      ## where val_i is the value on the ith node neighbouring the face and
      ## income = [0, 1] depending on whether the flow is coming from node 1 or 2
      ## First order upwinding is monotonic so no limiting is ever required
      element face_value {
        attribute name { "FirstOrderUpwind" },
        empty
      }|
      ## Trapezoidal face value discretisation
      ##  face_value = 0.5*(val_1 + val_2), 
      ## where
      ##  val_i is the value on the ith node neighbouring the face
      ##
      ## Trapezoidal discretisation is unbounded so limiting is compulsory
      element face_value {
        attribute name { "Trapezoidal" },
        limiter_options
      }|
      ## Finite Element face value discretisation
      ##  face_value = finite element interpolation from surrounding nodes
      ##
      ## Finite element discretisation may become unbounded so limiting is often necessary.
      element face_value {
        attribute name { "FiniteElement" },
        limiter_options?
      }
   )
   
coupled_spatial_control_volume_options = coupled_control_volume_options
coupled_spatial_control_volume_options |= compressive_control_volume_options

# coupled control volume options are the same as the standard ones (annoyingly copied and pasted)
# except that firstorderupwind gets limiter options
coupled_control_volume_options = 
   (
      ## First Order Upwind face value discretisation
      ##  face_value = donor_value, 
      ## where
      ##  donor_value = income*val_1 + (1.-income)*val_2, 
      ## where val_i is the value on the ith node neighbouring the face and
      ## income = [0, 1] depending on whether the flow is coming from node 1 or 2
      ## First order upwinding is monotonic so no limiting is ever required
      element face_value {
        attribute name { "FirstOrderUpwind" },
        limiter_options?
      }|
      ## Trapezoidal face value discretisation
      ##  face_value = 0.5*(val_1 + val_2), 
      ## where
      ##  val_i is the value on the ith node neighbouring the face
      ##
      ## Trapezoidal discretisation is unbounded so limiting is compulsory
      element face_value {
        attribute name { "Trapezoidal" },
        limiter_options
      }|
      ## Finite Element face value discretisation
      ##  face_value = finite element interpolation from surrounding nodes
      ##
      ## Finite element discretisation may become unbounded so limiting is often necessary.
      element face_value {
        attribute name { "FiniteElement" },
        limiter_options?
      }
   )
   
compressive_control_volume_options = 
   (
      ## HyperC face value discretisation
      ##
      ## face_value calculated from upper bound of explicit TVD zone of NVD diagram
      ## Normally used for MaterialVolumeFraction fields
      element face_value {
        attribute name { "HyperC" },
        upwind_value_options?,
        cv_face_cfl_number_options
      }|
      ## UltraC face value discretisation
      ##
      ## face_value calculated from extended upper bound of
      ## explicit TVD zone of NVD diagram assuming
      ## values bounded by target_maximum and target_minimum.
      element face_value {
        attribute name { "UltraC" },
        ## Specifiy the upper bound which UltraC will use to
        ## calculate the maximum flux.
        element target_maximum {
            real
        },
        ## Specifiy the lower bound which UltraC will use to
        ## calculate the minimum flux.
        element target_minimum {
            real
        },
        upwind_value_options?,
        cv_face_cfl_number_options
      }|
      ## **UNDER TESTING**
      ##
      ## PotentialUltraC face value discretisation
      ##
      ## face_value calculated from extended upper bound of
      ## explicit TVD zone of NVD diagram if potential
      ## value of field is sufficient (as specified by
      ## target_maximum) to ensure the correct front advection
      ## velocity.
      ##
      ## If not then either switch to HyperC or use a modified flux
      ## based on the potential function.
      element face_value {
        attribute name { "PotentialUltraC" },
        ## Specifiy the upper bound which PotentialUltraC will use
        ## to calculate the maximum flux if the potential function
        ## value is sufficient to maintain the correct front
        ## advection velocity.
        element target_maximum {
            real
        },
        ## Specifiy the lower bound which PotentialUltraC will use to calculate the minimum flux.
        element target_minimum {
            real
        },
        (
            ## Select what PotentialUltraC should do if the
            ## potential function value does not reach the required
            ## value specified by the target_maximum.
            ##
            ## Switch to using HyperC face values.  This ensures
            ## that the advection velocity is correct however may
            ## create isolated regions beneath the target_maximum.
            element switch_to_hyperc {
              empty
            }|
            ## Select what PotentialUltraC should do if the
            ## potential function value does not reach the required
            ## value specified by the target_maximum.
            ##
            ## Modify the maximum nodal values (both downwind and
            ## upwind) so that the fluxes are at their maximum
            ## possible without affecting the front advection
            ## velocity.
            element use_potential_flux {
              empty
            }
        ),
        upwind_value_options?,
        cv_face_cfl_number_options
      }
   )

cap_option =
   (
      ## Cap the min and max values of this field when using
      ## it as a volume fraction to work out bulk material
      ## properties.
      ## No capping used if not selected.
      element cap_values {
         ## Set the upper bound on the field.
         ## Defaults to huge(0.0)*epsilon(0.0) if not set.
         element upper_cap {
            real
         }?,
         ## Set the lower bound on the field.
         ## Defaults to -huge(0.0)*epsilon(0.0) if not set.
         element lower_cap {
            real
         }?
      }
   )

surface_tension_option =
   (
      element surface_tension {
        ## Surface tension coefficient
        element surface_tension_coefficient {
          real
        },
        ## The equilibrium contact angle (in radians) with the boundaries identified by the surface ids
        element equilibrium_contact_angle {
          real,
          ## Surface ids:
          element surface_ids {
              integer_vector
          }
        }?
      }
   )

limiter_options =
  (
      ## Limit the face value to satisfy a boundedness criterion.
      element limit_face_value{
        (
          sweby_limiter|
          ultimate_limiter
        )
      }
  )

sweby_limiter = 
  ## See "High-Resolution Schemes Using Flux Limiters for
  ## Hyperblic Conservation-Laws", P. K. Sweby, 1984, Siam
  ## Journal on Numerical Analysis, 21, 995-1011
  element limiter {
    attribute name {"Sweby"},
    slope_options?,
    upwind_value_options?
  }

ultimate_limiter =
  ## See "The Ultimate Conservative Difference Scheme Applied
  ## to Unsteady One-Dimensional Advection", B. P. Leonard,
  ## 1991, Computer Methods in Applied Mechanics and
  ## Engineering, 88, 17-74
  element limiter {
    attribute name {"Ultimate"},
    field_based_cfl_number_options,
    upwind_value_options?
  }

slope_options =
   (
      ## Control the upper and lower slopes of the NVD limiter
      element slopes {
         ## Defaults to Sweby, 1984 limiter (= 1.0) if unselected
         element lower {
            real
         }?,
         ## Defaults to Sweby, 1984 limiter (= 2.0) if unselected
         element upper {
            real
         }?
      }
   )

upwind_value_options =
   (
      (
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## This method projects the upwind value from a point in the element just
         ## upwind of the node pair straddling the face.  It is otherwise known as 
         ## anisotropic limiting.
         ## This is only available on simplex meshes as it involes a search around
         ## the donor node to find the upwind element.
         element project_upwind_value_from_point {
            ## When the donor node is on a domain boundary reflect the projection
            ## back into the mesh.
            element reflect_off_domain_boundaries {
               empty
            }?,
            ## Constrain the projected value to be between the min and max of the
            ## element values which it was found from.
            element bound_projected_value_locally {
               empty
            }?,
            ## Store the locations of the elements where the upwind values
            ## are projected from for each node pair.
            ## This inserts an integer csr matrix into state so is memory expensive but
            ## saves a significant amount of time (searching around the neighbouring elements).
            ## This is unsafe for moving meshes but should be ok for adaptive meshes.
            element store_upwind_elements {
               ## Store the quadrature locations within the elements
               ## where the upwind values
               ## are projected from for each node pair.
               ## This inserts a real block csr matrix into state so is even more memory
               ## expensive than just storing the upwind elements and
               ## only saves a comparitively
               ## marginal amount of time (as actually searching the
               ## neighbouring elements is the
               ## slowest bit, finding the quadrature is relatively easy).
               element store_upwind_quadrature {
                  empty
               }?
            }?
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Projects the value of the advected variable from the downwind or donor node
         ## using the interpolated gradient at the donor node in the
         ## direction of the vector
         ## connecting the node pair straddling the face.
         ## This is available on all meshes (except if bounding the values).
         element project_upwind_value_from_gradient {
            (
               ## Select which node to project from:
               ## Project from the downwind node (Jasak et al., 1999) so that:
               ## upwind_value = downwind_value - 2*gradient.vector
               element project_from_downwind_value {
                  comment
               }|
               ## Select which node to project from:
               ## Project from the donor node so that:
               ## upwind_value = donor_value - gradient.vector
               element project_from_donor_value {
                  comment
               }
            ),
            ## When the donor node is on a domain boundary reflect the projection
            ## back into the mesh.
            element reflect_off_domain_boundaries {
               empty
            }?,
            ## Constrain the projected value to be between the min and max of the
            ## element values which surround it.
            ## This is only available on simplex meshes as it involes a search around
            ## the donor node to find the upwind element.
            element bound_projected_value_locally {
               ## Store the locations of the elements closest to the project value.
               ## This inserts an integer csr matrix into state so is
               ## memory expensive but
               ## saves a significant amount of time (searching around
               ## the neighbouring elements).
               ## This is unsafe for moving meshes but should be ok for adaptive meshes.
               element store_upwind_elements {
                  comment
               }?
            }?
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Chooses an upwind value by selecting the maximum or minimum of the neighbouring
         ## nodes depending on the local slope of the donor and downwind values.
         ## Otherwise known as isotropic limiting.
         ## This is available on all meshes.
         element locally_bound_upwind_value {
            empty
         }|
         ## Select the method to be used for calculating the upwind value.
         ## If not selected will default to project_upwind_value_from_point for
         ## simplex element meshes and to a locally_bound_upwind_value for cube
         ## element meshes.
         ##
         ## Chooses an upwind value by selecting the value at the node most directy
         ## upwind from the vector connecting the donor and downwind nodes.
         ## This is available on all meshes.
         element pseudo_structured_upwind_value {
            empty
         }
      )
   )

field_based_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "CVMaterialDensityCFLNumber" },
            empty
         }|
         ## Select the Courant Number definition to be used.
         ##
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
            empty
         }|
         ## Select the Courant Number definition to be used.
         element courant_number {
            attribute name { string },
            empty
         }
      )
   )

cv_face_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses the finite difference definition of the CFL Number
         ## consistent with the 1D version of HyperC (Leonard, 1981).
         ## This is the default that reproduces old behaviour.
         ## All others are under testing or construction.
         element courant_number {
            attribute name { "FiniteDifferenceCFLNumber" },
            empty
         }|
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses a control volume definition of the CFL Number
         ## that incorporates the MaterialDensity.
         ## Requires a MaterialDensity field in this material_phase!
         element courant_number {
            attribute name { "CVMaterialDensityCFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
            empty
         }|
         ## ***UNDER TESTING***
         ##
         ## Select the Courant Number definition to be used in the slope of
         ## the NVD diagram upper bound.
         element courant_number {
            attribute name { string },
            empty
         }
      )
   )

timestep_cfl_number_options =
   (
      (
         ## Select the Courant Number definition to be used for adaptive timestepping.
         ## This uses the finite element approximation of the CFL Number.
         element courant_number {
            attribute name { "CFLNumber" },
              ## Select the mesh on which you wish to evaluate the CFLNumber.
              velocity_mesh_choice
         }|
         ## Select the Courant Number definition to be used for adaptive timestepping.
         ## This uses the control volume definition of the CFL Number.
         element courant_number {
            attribute name { "ControlVolumeCFLNumber" },
              ## Select the mesh on which you wish to evaluate the ControlVolumeCFLNumber.
              velocity_mesh_choice
         }
      )
   )

mixing_stats =
   (
      ## Enable to include in the .stat file the fractions of the
      ## scalar field contained in
      ## bins specified by the user. This allows mixing of the field to be quantified.
      ## Replaces and expands upon the old heaviside.dat file
      element include_mixing_stats{
         attribute name { xsd:string },
         (
            ## Select whether to evaluate the volume fraction over the finite element
            ## (continuous galerkin) or within the control volume (control_volumes).
            ##
            ## NOTE: continuous_galerkin only works with linear tets
            ##
            ## NOTE: continuous_galerkin is not fully validated yet
            element continuous_galerkin {
               ## if select normalise the volume fractions will be
               ## divided by the total volume of the domain
               element normalise {
                  empty
               }?
            }|
            ## Select whether to evaluate the volume fraction over the finite element
            ## (continuous galerkin) or within the control volume (control_volumes).
            element control_volumes {
               ## if select normalise the volume fractions will be divided by the total volume of the domain  
               element normalise {
                  empty
               }?
            }
         ),
         ## The values of the bounds of the bins 
         ## e.g. the values 0 1 2 3 will return 4 bins 
         ## and the fraction of the field in each bin with,
         ## 0<=field<1, 1<=field<2, 2<=field<3, 3<=field, 
         ## will be calculated. 
         element mixing_bin_bounds { 
            real_vector 
         },
         ## Define the tolerance beneath the specified bins that should be included.
         ## Defaults to zero at machine tolerance (epsilon(0.0)) if not selected.
         element tolerance {
            real
         }?
      }
   )

cv_stats =
   (
      ## Include statistics evaluated on the control volume mesh.
      element include_cv_stats {
         empty
      }
   )

# Options for inclusion of calculations of surface integrals in the .stat file   
surface_integral_stats_base.surface_integral =
   (
      attribute name { xsd:string },
      ## Surface IDs defining the surface over which to integrate. If disabled, integrates over the whole surface.
      element surface_ids {
         integer_vector
      }?,
      ## Enable to normalise the integral by dividing by the surface area
      element normalise {
         comment
      }?
   )
surface_integral_stats_scalar =
   (
      ## Surface integral calculations. The following integral types are available:
      ##  value: Integrates the field
      ##  gradient_normal: Integrates the normal component of the gradient of the field
      element surface_integral {
         surface_integral_stats_scalar.surface_integral
      }
   )
surface_integral_stats_scalar.surface_integral = surface_integral_stats_base.surface_integral
surface_integral_stats_scalar.surface_integral &=
   (
      attribute type { "value" | "gradient_normal" }
   )
surface_integral_stats_vector =
   (
      ## Surface integral calculations. The following integral types are available:
      ##  normal: Integrates the normal component of the field
      element surface_integral {
         surface_integral_stats_vector.surface_integral
      }
   )
surface_integral_stats_vector.surface_integral = surface_integral_stats_base.surface_integral
surface_integral_stats_vector.surface_integral &=
   (
      attribute type { "normal" }
   )

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" }
      }
   )

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" }
         }|
         ## ***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 WITHOUT A
         ## DIFFUSIVITY, SOURCE OR ABSORPTION.
         element equation { 
            attribute name { "ConservationOfMass" },
            (
               ## Select density to use in the Conservation of Mass Equation
               ## Use the MaterialDensity - useful for multimaterial simulations
               ## Clearly this requires a MaterialDensity field to be present
               element density {
                  attribute name { "MaterialDensity" }
               }|
               ## Select density to use in the Conservation of Mass Equation
               ## Use the bulk Density
               ## Clearly this requires a Density field to be present
               element density {
                  attribute name { "Density" }
               }|
               ## Select density to use in the Conservation of Mass Equation
               element density {
                  attribute name { string }
               }
            )
         }|
         ## ***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 WITHOUT A
         ## DIFFUSIVITY, SOURCE OR ABSORPTION.
         ##
         ## 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" },
            (
               ## Select density to use in the Reduced Conservation of Mass Equation
               ## Use the MaterialDensity - useful for multimaterial simulations
               ## Clearly this requires a MaterialDensity field to be present
               element density {
                  attribute name { "MaterialDensity" }
               }|
               ## Select density to use in the Reduced Conservation of Mass Equation
               ## Use the bulk Density
               ## Clearly this requires a Density field to be present
               element density {
                  attribute name { "Density" }
               }|
               ## Select density to use in the Reduced Conservation of Mass Equation
               element density {
                  attribute name { string }
               }
            )
         }|
         ## ***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 WITHOUT A
         ## DIFFUSIVITY, SOURCE OR ABSORPTION.
         ## 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" },
            (
               ## Select density to use in the Internal Energy 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" }
               }|
               ## Select density to use in the Internal Energy 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" }
               }|
               ## Select density to use in the Internal Energy Equation
               ## Whatever field is selected must be present.
               element density {
                  attribute name { string }
               }
            )
         }|
      )
   )

inner_element_scalar = 
   (
      ## Inner element sub-grid scale model (Candy and Pain)
      ## Requires continuous galerkin selected above.
      element inner_element {
         ## Inner element solution of the scalar field.
         element scalar_field {
            attribute rank { "0" },
            attribute name { "InnerElement" },
            element prognostic {
               element mesh {
                  attribute name {  "InnerElementMesh" }
               },
               prognostic_scalar_output_options,
               prognostic_scalar_stat_options,
               prognostic_detector_options
            }
         }
      }
   )
      
inner_element_velocity =
   (
      ## Inner element sub-grid scale model (Candy and Pain)
      ## Requires continuous galerkin selected above.
      element inner_element {
         ## SGS velocity in an inner element SGS treatment of momentum
         ##
         ## Limitations:
         ##  - Requires a geometry dimension of 3.
         ##  - Requires inner element active for momentum
         element vector_field {
            attribute rank { "1" },
            attribute name { "InnerElement" },
            element prognostic {
               element mesh {
                  attribute name {  "InnerElementMesh" }
               },
               prognostic_vector_output_options,
               prognostic_vector_stat_options,
               vector_convergence_options
            }
         },
         ## A filter for the sub-grid scale equations
         ## Add diffusion to matrix D of the Inner Element model
         element use_filter {
            ## Strength of the diffusion term
            ## Suggested value: 0.01
            element strength {
               real
            },
            empty
         }?,
         element use_quadratic_pressure {
            empty
         }?,
         element apply_full_discontinuous_Galerkin {
            empty
         }?
      }?
   )
forcing =
  (
    ## Add forcing from ocean data
    ## If you enable this you MUST enable the /geometry/ocean_boundaries option too
    element ocean_forcing{
        ## The netCDF data file downloaded from ERA-40 reanalysis
        element input_file {
           attribute file_name {xsd:string} 
        },
        element mesh_choice {
            velocity_mesh_choice
        },
        element surface_stress {
            empty
        }?,
        element temperature_flux {
            empty
        }?,
        element salinity_flux {
            empty
        }?,
        element solar_radiation {
            empty
        }?
    }
  )

biology =
   (
      ## Model of biological processes in the ocean.
      element ocean_biology{
         ## A simple model of phytoplankton, zooplankton, general nutrient and detritus. 
         element pznd {
            (
               ## Python code specifying the source and sink relationships 
               ## between the biological tracers. This is usually achieved by 
               ## importing fluidity.ocean_biology and calling a scheme from there. 
               element source_and_sink_algorithm {
                  python_code
               }|
               ## Do not calculate sources and sinks. 
               ## This option is generally only useful for testing. 
               element disable_sources_and_sinks {
                  empty
               }
             ),
#            ## Phytoplankton
#            element scalar_field {
#               attribute rank { "0" },
#               attribute name { "Phytoplankton" },
#               (
#                  element prognostic {
#                     velocity_mesh_choice,
#                     prognostic_scalar_field
#                  }|
#                  element prescribed {
#                     velocity_mesh_choice,
#                     prescribed_scalar_field
#                  }
#               )
#            },
#            ## Zooplankton
#            element scalar_field {
#               attribute rank { "0" },
#               attribute name { "Zooplankton" },
#               (
#                  element prognostic {
#                     velocity_mesh_choice,
#                     prognostic_scalar_field
#                  }|
#                  element prescribed {
#                     velocity_mesh_choice,
#                     prescribed_scalar_field
#                  }
#               )
#            },
#            ## Nutrient
#            element scalar_field {
#               attribute rank { "0" },
#               attribute name { "Nutrient" },
#               (
#                  element prognostic {
#                     velocity_mesh_choice,
#                     prognostic_scalar_field
#                  }|
#                  element prescribed {
#                     velocity_mesh_choice,
#                     prescribed_scalar_field
#                  }
#               )
#            },
#            ## Detritus
#            element scalar_field {
#               attribute rank { "0" },
#               attribute name { "Detritus" },
#               (
#                  element prognostic {
#                     velocity_mesh_choice,
#                     prognostic_scalar_field
#                  }|
#                  element prescribed {
#                     velocity_mesh_choice,
#                     prescribed_scalar_field
#                  }
#               )
#            },
            ## Photosynthetically Active Radiation (PAR)
            element scalar_field {
               attribute rank { "0" },
               attribute name { "PhotosyntheticRadiation" },
               (
                  element prognostic {
                     velocity_mesh_choice,
                     prognostic_photosynthetic_radiation
                  }|
                  element prescribed {
                     velocity_mesh_choice,
                     prescribed_scalar_field
                  }
               )
            }
         }
      }
   )


prognostic_photosynthetic_radiation =
   (
      ## PAR equation.
      element equation { 
         attribute name { "PhotosyntheticRadiation" }
      },
      ## Spatial discretisation options
      element spatial_discretisation {
         (
            ## Discontinuous galerkin formulation. You can also use this
            ## formulation with a continuous field in which case a simple
            ## galerkin formulation will result. 
            element discontinuous_galerkin {
               empty
            }
         )
      },
      (
         ## Solver
         element solver {
            linear_solver_options_asym
         }
      ),
      ## Coefficients of absorption of photosynthetically active
      ## radiation for water and phytoplankton.
      element absorption_coefficients {
         ## Photosynthetically active radiation absorption coefficient for water.
         element water {
            real
         },
         ## Photosynthetically active radiation absorption coefficient for water.
         element phytoplankton {
            real
         }
      },
      ## 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.  Only available with
               ## discontinuous_galerkin, control_volume and
               ## legacy_mixed_cv_cg spatial_discretisations.
               ##
               ## If not selected boundary conditions are applied strongly.
               element apply_weakly {
                 empty
               },
               input_choice_real
            }|
            element type {
               attribute name { "neumann" },
               input_choice_real
            }|
            element type {
               attribute name { "robin" },
               element order_zero_coefficient {
                  input_choice_real
               },
               element order_one_coefficient {
                  input_choice_real
               }
            }
         )
      }+,
      prognostic_scalar_output_options,
      prognostic_scalar_stat_options,
      scalar_convergence_options,
      prognostic_detector_options,
      scalar_steady_state_options,
      adaptivity_options_scalar_field,
      ## 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
      }?
   )

recalculation_options =
   (
      ## Prevent this field from being recalculated at every timestep.
      ## This is cheaper especially if you are enforcing discrete properties on the field.
      element do_not_recalculate {
        empty
      }
   )

discrete_properties_algorithm_scalar =
   (
      ## Select discrete properties to enforce on the field
      ## either after being prescribed or interpolated
      element enforce_discrete_properties {
        ## Update this field using the lagrangian multiplier
        ## calculated in the solenoidal projection of a
        ## scalar field.
        ##
        ## Note this field must be specified as the update field
        ## underneath that vector field too.
        ##
        ## Note also this only really makes sense for coupled
        ## fields like velocity and pressure.
        element solenoidal_lagrange_update {
          empty
        }?
      }
   )

discrete_properties_algorithm_vector =
   (
      ## Select discrete properties to enforce on the field
      ## either after being prescribed or interpolated
      element enforce_discrete_properties {
        solenoidal_options?
      }
   )

interpolation_algorithm_disabled =
   (
      ## Disable interpolation
      element no_interpolation {
        comment
      }
   )

interpolation_algorithm_scalar =
    (
      ## Basis function interpolation.
      ## The standard algorithm. It is quick
      ## and bounded, but non-conservative and dissipative.
      ## All other algorithms require construction of a supermesh.
      element consistent_interpolation {
        empty
      }|
      ## Galerkin projection. By default, conservative, non-dissipative and
      ## non-bounded. The most accurate choice, in the sense of minimising
      ## the L2 norm of the residual
      element galerkin_projection {
        galerkin_projection_scalar
      }
    )
    
interpolation_algorithm_scalar_full = interpolation_algorithm_scalar
interpolation_algorithm_scalar_full |= interpolation_algorithm_disabled
    
interpolation_algorithm_vector =
    (
      ## Basis function interpolation.
      ## The standard algorithm. It is quick
      ## and bounded, but non-conservative and dissipative.
      ## All other algorithms require construction of a supermesh.
      element consistent_interpolation {
        consistent_interpolation_vector
      }|
      ## Galerkin projection. By default, conservative, non-dissipative and
      ## non-bounded. The most accurate choice, in the sense of minimising
      ## the L2 norm of the residual
      element galerkin_projection {
        galerkin_projection_vector
      }
    )

consistent_interpolation_vector = 
    (
      balanced_interpolation?
    )
    
galerkin_projection_vector =
    (
      continuous_discontinuous_projection,
      balanced_interpolation?,
      supermesh_free?,
      supermesh_conservation?
    )

galerkin_projection_scalar =
    (
      continuous_discontinuous_projection,
      supermesh_free?,
      supermesh_conservation?
    )

continuous_discontinuous_projection =
    (
      ## Continuous field Galerkin projection.
      ## If the field you are interpolating is continuous, then
      ## a linear solver is required to invert the mass matrix.
      element continuous {
        (
          ## Use a bounded Galerkin projection. Conservative, bounded in the
          ## limit, and minimally dissipative. This algorithm starts with the
          ## Galerkin projection and dissipates it until it achieves
          ## boundedness.
          ## If it does not converge, it may not be exactly bounded.
          ## Note well: this only works for linear fields.
          element bounded {
            attribute name {"Diffuse"},
            ## The number of dissipation iterations attempted to bound the
            ## Galerkin projection.
            element boundedness_iterations {
              integer,
              ## Specify the tolerance to which boundedness is to be tested.
              ## Defaults to computer precision if unspecified.
              element tolerance {
                real
              }?
            },
            ## If the bounds on this field are known then they can be set here.
            ## These can either further constrain the limits worked out by the
            ## lumped version of the projection (i.e. to make sure that errors 
            ## don't accumulate with succesive interpolations) or if apply_globally
            ## is set they are just made to be bounded within the bounds globally
            ## (i.e. anything between those bounds is not smoothed).
            element bounds {
              element upper_bound {
                real,
                ## If this is set the upper_bound is used everywhere.
                ## If left unset the upper_bound is only used to constrain
                ## the smoothed bounds calculated by the code
                element apply_globally {
                  empty
                }?,
                ## This field is to be considered as being coupled to another field
                ## such that the sum of the two fields is constrained to be less than
                ## the upper_bound specified above.
                ## 
                ## The relationships between fields are worked out according to their
                ## priority ordering.
                ##
                ## This method is akin to the coupled_cv advection method.
                element coupled {
                  empty
                }?
              }?,
              element lower_bound {
                real,
                ## If this is set the upper_bound is used everywhere.
                ## If left unset the upper_bound is only used to constrain
                ## the smoothed bounds calculated by the code
                element apply_globally {
                  empty
                }?
              }?
            }?,
            ## If, after performing all the boundedness_iterations, the field
            ## is still not bounded then perform surgery to redistribute the
            ## deviations to nodes that have less than their bounds.
            element repair_deviations {
              empty
            }?
          }
        )?,
        (
          ## Solver options for the linear solve.
          ## This method requires the inversion of a mass matrix. Note that
          ## conservation properties are affected by the tolerance of the
          ## linear solve.
          element solver {
            linear_solver_options_sym
          }|
          ## Lump the mass matrix on the left hand side of the galerkin projection.
          ## Hence solver options aren't necessary.
          element lump_mass_matrix {
            empty
          }
        )
      }|
      ## Discontinuous field Galerkin projection
      ## In this case, no linear solver is required to invert the mass matrix.
      element discontinuous {
        empty
      }
    )

solenoidal_options =
    ## Constrained divergence-free projection.
    ## This adds an additional constraint that ensures that the field
    ## is solenoidal, i.e. divergence-free.
    ## This is equivalent in cost to a pressure solve.
    ## This is expensive, and thus best left until
    ## needed.
    ##
    ## Note well: this only makes sense for nondivergent
    ## vector fields, such as incompressible velocity!
    element solenoidal {
      ## Options for the mass matrix of the field being interpolated
      element interpolated_field {
        (
          element continuous {
            ## Lump the mass matrix for the assembly of the projection matrix
            ## (not for the initial galerkin projection)
            ##
            ## Required when using interpolating continuous fields
            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
              }?
            }
          }|
          element discontinuous {
            ## Lump the mass matrix for the assembly of the projection matrix
            ## (not for the initial galerkin projection)
            element lump_mass_matrix {
              empty
            }?
          }
        )
      },
      ## Options for the lagrange multiplier
      ##
      ## Must be on a continuous mesh!
      element lagrange_multiplier {
        pressure_mesh_choice,
        element spatial_discretisation {
          (
            element continuous_galerkin {
              ## Remove the stabilisation term from the projection operator.
              ##
              ## Automatic when not using P1P1.
              element remove_stabilisation_term {
                empty
              }?,
              ## Integrate the divergence operator by parts.
              ##
              ## Automatic when projecting a discontinuous field
              element integrate_divergence_by_parts {
                empty
              }?
            }|
            element control_volumes {
              empty
            }
          )
        },
        element reference_node {
          integer
        }?,
        (
          ## Update a scalar field using the lagrange multiplier from
          ## the divergence free projection of this field.  The selected
          ## scalar field must have solenoidal selected in its interpolation
          ## options too and it must be on the same mesh as used for the
          ## solenoidal projection above.
          ##
          ## Note well: this only really makes sense for scalar fields linked to nondivergent
          ## vector fields, such as pressure to incompressible velocity!                  
          element update_scalar_field {
            attribute name { "Pressure" },
            empty
          }|
          ## Update a scalar field using the lagrange multiplier from
          ## the divergence free projection of this field.  The selected
          ## scalar field must have solenoidal selected in its interpolation
          ## options too and it must be on the same mesh as used for the
          ## solenoidal projection above.
          ##
          ## Note well: this only really makes sense for scalar fields linked to nondivergent
          ## vector fields, such as pressure to incompressible velocity!                  
          element update_scalar_field {
            attribute name { string },
            empty
          }
        )?,
        ## Solver options for the linear solve.
        ## This method requires the inversion of a projection matrix.
        element solver {
          linear_solver_options_sym
        }
      }
    }


balanced_interpolation =
    ## Geostrophically-balanced interpolation.
    ## During the interpolation, the velocity is split into
    ## balanced and imbalanced parts. The imbalanced part is
    ## interpolated whereas the balanced part is recovered
    ## on the other side of the interpolation from the pressure.
    ## This means that if your velocity field is balanced before
    ## interpolation, it will be balanced afterwards also.
    ##
    ## Note well: this only makes sense for velocity.
    element balanced_interpolation {
      ## The solver options for computing the balanced velocity to subtract.
      element solver {
        linear_solver_options_sym
      }
    }
    
supermesh_free =
    ## Enables a supermesh free Galerkin projection. Uses incomplete
    ## quadrature, and hence is not conservative.
    element supermesh_free {
      empty
    }

represcribe_before_interpolation =
    ## Represcribe the field before interpolation.
    ##
    ## This means the interpolation will not be conservative from the previous mesh so be careful what you're trying to achieve!
    element represcribe_before_interpolation {
      empty
    }
    
supermesh_conservation =
      ## Options for checking the supermesh conservation properties
      element supermesh_conservation {
        ## Specify the fraction of the original elemental area/volume
        ## to be used to check the conservation of the supermesh.
        ##
        ## Since all fields are supermeshed together the minimum tolerance
        ## specified over all fields will be used.
        ##
        ## Defaults to 0.001 if unspecified.
        ## i.e. 0.1% of the area/volume of an element in the new mesh may
        ## be lost without warning or attempts to fix (if compiled with cgal)
        ## during the construction of the supermesh between the old
        ## and new meshes.
        element tolerance {
          real
        }?,
        ## Compute the field integral after the interpolation and print the relative
        ## mass loss to the logfile (level 2 verbosity).
        ##
        ## Note this is a post interpolation step and offers no chance of
        ## fixing the conservation error (unlike the tolerance above if compiled
        ## with cgal)
        element print_field_integral {
          ## Relative tolerance with which to test the conservation of the field
          ## integral.  If the conservation fails this tolerance a warning is issued
          ## (level 0 verbosity) and vtus containing the field are output.
          element tolerance {
            real
          }
        }?
      }