Merge lp:~ctjacobs-multiphase/fluidity/enforce-hydrostatic-balance into lp:fluidity

[Christian Jacobs]

This merge will add a pressure/density splitting technique (see Giraldo and Restelli (2008), doi: 10.1016/j.jcp.2007.12.009), which helps to prevent spurious oscillations in the pressure field when using unbalanced finite element pairs.

When the enforce_hydrostatic_balance option (under the Velocity field) is enabled, the Density and Pressure fields are split up into a hydrostatic component (') and a perturbed component (''). The hydrostatic components, denoted p' and rho', satisfy the balance: grad(p') = rho'*g. The hydrostatic components are then subtracted from the pressure and density used in the pressure gradient and buoyancy terms in the momentum equation. This works with both CG and DG Velocity discretisations.

This merge also includes: single-phase and multiphase tests for this functionality, a description in the manual and schema, and some extra sanity checks regarding the application of pressure boundary conditions and inlet velocity boundary conditions for the multiphase flow model.

[Christian Jacobs (ctjacobs)]

The branch queue has a green light on Buildbot: http://buildbot-ocean.ese.ic.ac.uk:8080/builders/enforce-hydrostatic-balance/builds/0

[James Robert Percival (j-percival)]

Code looks fine to me, and the test appears sensible. Modulo concerns about the name enforce_hydrostatic_balance I think this is good to go. Actually, it's probably worth checking with Stephan what name he thinks is sensible for this. If we ever get keen we could also try implementing the sound proof equations using the same technology.

[Stephan Kramer (s-kramer)]

First of all: great job - as always, well written, tested and documented. We had a little discussion off-line (also with input from CianW) that I'll summarize below. Main points were:
* The approach here deviates from all the other hydrostatic/geostrophic balance options in that the hydrostatic pressure is not actually subtracted from the Pressure field, i.e. Pressure is still the combined Pressure, it's only subtracted in the momentum equation. This is not what happens with the subtract_out_hydrostatic_level option under eos/fluids/linear or with HydrostaticPressure or GeostrophicPressure fields. In fact the Giraldo paper that you cite does exactly what these other options do, i.e. only the perturbed pressure occurs in the momentum equation (of course the hydrostatic part then reappears in the linearized continuity equation). In any case, I don't think this really matters in practice so we're happy to go with it as it is, as long as it's clearly explained.

* Since this is an option that is sort of orthogonal to subtract_out_hydrostatic and probably only applicable for compressible flows - we thought it'd be best to have it as and option under equation_of_state/compressible/. Granted, subtract_out_hydrostatic probably shouldn't have been under eos in the first place, but that situation stems from long fluidity history and it's probably better to remain consistent with it for now. (This is a thing we can easily fix in the next generation of Fluidity).

* we didn't really like enforce_hydrostatic as the name for the option as we're not really enforcing anything - we're merely helping it find a balance in cases there is one. Suggestions were subtract_out_reference_(pressure_)profile. It would be good to emphasize the difference with the fluids/linear subtract option (i.e. the hydrostatic reference component is *not* subtracted from the Pressure field).

* In the naming of the fields currently called HydrostaticBalancePressure and HydrostaticBalanceDensity it's probably better to make a clearer distinction with the other hydrostatic balance fields (HydrostaticPressure and GeostrophicPressure) by choosing names such as HydrostaticReferencePressure and HydrostaticReferenceDensity. This makes it more clear that these are meant to be set up as time-independent, user-provided reference profiles. Also reference density/pressure (profile) seems to be a term that's used more often for this sort of thing (see a.o. Vallis) as a natural extension of constant reference density \rho_0

* It would be good to emphasize in the documentation that these fields are meant to be set up as vertical profiles (i.e. horizontally constant).

Some notes, reading through the diffs:

Momentum_DG.F90, line 2365: the sign of ele_val(pressure_bc, face) seems to change - not sure that's correct?

bibliography.bib: you seem to break the convention for citation keys that's used here

[lp:~ctjacobs-multiphase/fluidity/enforce-hydrostatic-balance]

4227. By Christian Jacobs on 2013-07-04

- Changed the sign of ele_val(pressure_bc, face) in Momentum_DG.F90. This shouldn't affect any existing setups.
- Removed the reference to Giraldo and Restelli (2008).
- Renamed HydrostaticBalancePressure and HydrostaticBalanceDensity to HydrostaticReferencePressure and HydrostaticReferenceDensity, respectively.
- Added the sentence in the schema and manual about the differences between this approach and the other approaches used in Fluidity, as noted by Stephan.
- Renamed enforce_hydrostatic_balance to subtract_out_reference_profile under the compressible EoS options.

4228. By Christian Jacobs on 2013-07-04

Corrected option name.

4229. By Christian Jacobs on 2013-07-05

Corrected field names.

4230. By Christian Jacobs on 2013-07-05

Renamed the tests to subtract_out_reference_profile and mphase_subtract_out_reference_profile.

4231. By Christian Jacobs on 2013-07-05

Corrected pass_tests in the .xml.

[Christian Jacobs (ctjacobs)]

Changes have now been made to address Stephan's comments. The buildbot queue has a green light. Stephan/James: can you please take a look and let me know what you think.

[Stephan Kramer (s-kramer)]

Brilliant. Thanks a lot for that. I'm entirely happy now.

review: Approve

[lp:~ctjacobs-multiphase/fluidity/enforce-hydrostatic-balance]

4232. By Christian Jacobs on 2013-07-10

Merge in trunk.

=== modified file 'assemble/Momentum_CG.F90'
--- assemble/Momentum_CG.F90	2013-03-06 13:20:27 +0000
+++ assemble/Momentum_CG.F90	2013-07-10 19:09:26 +0000
@@ -156,6 +156,13 @@
     ! wetting and drying switch
     logical :: have_wd_abs
 
+    ! If .true., the pressure and density fields will be split up into hydrostatic
+    ! and perturbed components. The hydrostatic components will be subtracted 
+    ! from the pressure and density used in the pressure gradient and buoyancy terms
+    ! in the momentum equation. This helps to maintain hydrostatic balance and prevent
+    ! spurious oscillations in the pressure field when using unbalanced finite element pairs.
+    logical :: subtract_out_reference_profile
+
     ! scale factor for the absorption
     real :: vvr_sf 
     ! scale factor for the free surface stabilisation
@@ -255,6 +262,9 @@
       ! for temperature dependent viscosity :
       type(scalar_field), pointer :: temperature
 
+      ! Fields for the subtract_out_reference_profile option under the Velocity field
+      type(scalar_field), pointer :: hb_density, hb_pressure
+
       integer :: stat, dim, ele, sele
 
       ! Fields for vertical velocity relaxation
@@ -360,6 +370,21 @@
       end if
       ewrite_minmax(buoyancy)
 
+      ! Splits up the Density and Pressure fields into a hydrostatic component (') and a perturbed component (''). 
+      ! The hydrostatic components, denoted p' and rho', should satisfy the balance: grad(p') = rho'*g
+      ! We subtract the hydrostatic component from the density used in the buoyancy term of the momentum equation.
+      if (have_option(trim(state%option_path)//'/equation_of_state/compressible/subtract_out_reference_profile')) then
+         subtract_out_reference_profile = .true.
+         hb_density => extract_scalar_field(state, "HydrostaticReferenceDensity", stat)
+         if(stat /= 0) then
+            FLExit("When using the subtract_out_reference_profile option, please set a (prescribed) HydrostaticReferenceDensity field.")
+            ewrite(-1,*) 'The HydrostaticReferenceDensity field, defining the hydrostatic component of the density field, needs to be set.'
+         end if
+      else
+         subtract_out_reference_profile = .false.
+         hb_density => dummyscalar
+      end if
+
       viscosity=>extract_tensor_field(state, "Viscosity", stat)
       have_viscosity = stat == 0
       if(.not. have_viscosity) then
@@ -684,7 +709,7 @@
          call construct_momentum_element_cg(state, ele, big_m, rhs, ct_m, mass, inverse_masslump, &
               x, x_old, x_new, u, oldu, nu, ug, &
               density, ct_rhs, &
-              source, absorption, buoyancy, gravity, &
+              source, absorption, buoyancy, hb_density, gravity, &
               viscosity, grad_u, &
               fnu, tnu, leonard, strainprod, alpha, gamma, &
               gp, surfacetension, &
@@ -730,6 +755,14 @@
            fs_sf=get_surface_stab_scale_factor(u)
          end if
 
+         if(subtract_out_reference_profile.and.integrate_continuity_by_parts.and.(assemble_ct_matrix_here .or. include_pressure_and_continuity_bcs)) then
+            hb_pressure => extract_scalar_field(state, "HydrostaticReferencePressure", stat)
+            if(stat /= 0) then
+               FLExit("When using the subtract_out_reference_profile option, please set a (prescribed) HydrostaticReferencePressure field.")
+               ewrite(-1,*) 'The HydrostaticReferencePressure field, defining the hydrostatic component of the pressure field, needs to be set.'
+            end if
+         end if
+
          surface_element_loop: do sele=1, surface_element_count(u)
             
             ! if no_normal flow and no other condition in the tangential directions, or if periodic
@@ -744,7 +777,7 @@
             call construct_momentum_surface_element_cg(sele, big_m, rhs, ct_m, ct_rhs, &
                  inverse_masslump, x, u, nu, ug, density, gravity, &
                  velocity_bc, velocity_bc_type, &
-                 pressure_bc, pressure_bc_type, &
+                 pressure_bc, pressure_bc_type, hb_pressure, &
                  assemble_ct_matrix_here, include_pressure_and_continuity_bcs, oldu, nvfrac)
             
          end do surface_element_loop
@@ -898,7 +931,7 @@
     subroutine construct_momentum_surface_element_cg(sele, big_m, rhs, ct_m, ct_rhs, &
                                                      masslump, x, u, nu, ug, density, gravity, &
                                                      velocity_bc, velocity_bc_type, &
-                                                     pressure_bc, pressure_bc_type, &
+                                                     pressure_bc, pressure_bc_type, hb_pressure, &
                                                      assemble_ct_matrix_here, include_pressure_and_continuity_bcs,&
                                                      oldu, nvfrac)
 
@@ -923,6 +956,7 @@
 
       type(scalar_field), intent(in) :: pressure_bc
       integer, dimension(:), intent(in) :: pressure_bc_type
+      type(scalar_field), intent(in) :: hb_pressure
       
       logical, intent(in) :: assemble_ct_matrix_here, include_pressure_and_continuity_bcs
 
@@ -1028,8 +1062,15 @@
                 !      /
                 ! add -|  N_i M_j \vec n p_j, where p_j are the prescribed bc values
                 !      /
-                call addto(rhs, dim, u_nodes_bdy, -matmul( ele_val(pressure_bc, sele), &
-                                                            ct_mat_bdy(dim,:,:) ))
+                if (subtract_out_reference_profile) then
+                   ! Here we subtract the hydrostatic component from the pressure boundary condition used in the surface integral when
+                   ! assembling ct_m. Hopefully this will be the same as the pressure boundary condition itself.
+                   call addto(rhs, dim, u_nodes_bdy, -matmul(ele_val(pressure_bc, sele)-face_val(hb_pressure, sele), &
+                                                            ct_mat_bdy(dim,:,:) ))
+                else
+                   call addto(rhs, dim, u_nodes_bdy, -matmul( ele_val(pressure_bc, sele), &
+                                                            ct_mat_bdy(dim,:,:) ))
+                end if
              end if
           end do
         end if
@@ -1125,7 +1166,7 @@
                                             mass, masslump, &
                                             x, x_old, x_new, u, oldu, nu, ug, &
                                             density, ct_rhs, &
-                                            source, absorption, buoyancy, gravity, &
+                                            source, absorption, buoyancy, hb_density, gravity, &
                                             viscosity, grad_u, &
                                             fnu, tnu, leonard, strainprod, alpha, gamma, &
                                             gp, surfacetension, &
@@ -1172,6 +1213,8 @@
       ! Temperature dependent viscosity:
       type(scalar_field), intent(in) :: temperature
 
+      type(scalar_field), intent(in) :: hb_density
+
       ! Non-linear approximation of the volume fraction
       type(scalar_field), intent(in) :: nvfrac
       ! Pointer to the nvfrac field's shape function
@@ -1347,7 +1390,7 @@
       
       ! Buoyancy terms
       if(have_gravity) then
-        call add_buoyancy_element_cg(x, ele, test_function, u, buoyancy, gravity, nvfrac, detwei, rhs_addto)
+        call add_buoyancy_element_cg(x, ele, test_function, u, buoyancy, hb_density, gravity, nvfrac, detwei, rhs_addto)
       end if
       
       ! Surface tension
@@ -1673,12 +1716,12 @@
       
     end subroutine add_sources_element_cg
     
-    subroutine add_buoyancy_element_cg(positions, ele, test_function, u, buoyancy, gravity, nvfrac, detwei, rhs_addto)
+    subroutine add_buoyancy_element_cg(positions, ele, test_function, u, buoyancy, hb_density, gravity, nvfrac, detwei, rhs_addto)
       type(vector_field), intent(in) :: positions
       integer, intent(in) :: ele
       type(element_type), intent(in) :: test_function
       type(vector_field), intent(in) :: u
-      type(scalar_field), intent(in) :: buoyancy
+      type(scalar_field), intent(in) :: buoyancy, hb_density
       type(vector_field), intent(in) :: gravity
       type(scalar_field), intent(in) :: nvfrac
       real, dimension(ele_ngi(u, ele)), intent(in) :: detwei
@@ -1687,12 +1730,14 @@
       real, dimension(ele_ngi(u, ele)) :: nvfrac_gi
       real, dimension(ele_ngi(u, ele)) :: coefficient_detwei
       
-      if(multiphase) then
-         nvfrac_gi = ele_val_at_quad(nvfrac, ele)
+      if (subtract_out_reference_profile) then
+        coefficient_detwei = gravity_magnitude*(ele_val_at_quad(buoyancy, ele)-ele_val_at_quad(hb_density, ele))*detwei
+      else
+        coefficient_detwei = gravity_magnitude*ele_val_at_quad(buoyancy, ele)*detwei
       end if
 
-      coefficient_detwei = gravity_magnitude*ele_val_at_quad(buoyancy, ele)*detwei
       if(multiphase) then
+         nvfrac_gi = ele_val_at_quad(nvfrac, ele)
          coefficient_detwei = coefficient_detwei*nvfrac_gi
       end if
 

=== modified file 'assemble/Momentum_DG.F90'
--- assemble/Momentum_DG.F90	2013-05-30 09:33:45 +0000
+++ assemble/Momentum_DG.F90	2013-07-10 19:09:26 +0000
@@ -135,6 +135,7 @@
   logical :: have_advection
   logical :: move_mesh
   logical :: have_pressure_bc
+  logical :: subtract_out_reference_profile
   
   real :: gravity_magnitude
 
@@ -215,6 +216,9 @@
     !! Surface tension field
     type(tensor_field) :: surfacetension
 
+    ! Fields for the subtract_out_reference_profile option under the Velocity field
+    type(scalar_field), pointer :: hb_density, hb_pressure
+
     !! field over the entire surface mesh, giving bc values
     type(vector_field) :: velocity_bc
     type(scalar_field) :: pressure_bc
@@ -395,6 +399,20 @@
     end if
     ewrite_minmax(buoyancy)
 
+    ! Splits up the Density and Pressure fields into a hydrostatic component (') and a perturbed component (''). 
+    ! The hydrostatic components, denoted p' and rho', should satisfy the balance: grad(p') = rho'*g
+    ! We subtract the hydrostatic component from the density used in the buoyancy term of the momentum equation.
+    if (have_option(trim(state%option_path)//'/equation_of_state/compressible/subtract_out_reference_profile')) then
+       subtract_out_reference_profile = .true.
+       hb_density => extract_scalar_field(state, "HydrostaticReferenceDensity")
+
+       if(l_include_pressure_bcs) then
+          hb_pressure => extract_scalar_field(state, "HydrostaticReferencePressure")
+       end if
+    else
+       subtract_out_reference_profile = .false.
+    end if
+
     Viscosity=extract_tensor_field(state, "Viscosity", stat)
     have_viscosity = (stat==0)
     if (.not.have_viscosity) then
@@ -661,7 +679,7 @@
        ele = fetch(colours(clr), nnid)
        call construct_momentum_element_dg(ele, big_m, rhs, &
             & X, U, advecting_velocity, U_mesh, X_old, X_new, &
-            & Source, Buoyancy, gravity, Abs, Viscosity, &
+            & Source, Buoyancy, hb_density, hb_pressure, gravity, Abs, Viscosity, &
             & P, Rho, surfacetension, q_mesh, &
             & velocity_bc, velocity_bc_type, &
             & pressure_bc, pressure_bc_type, &
@@ -718,7 +736,7 @@
   end subroutine construct_momentum_dg
 
   subroutine construct_momentum_element_dg(ele, big_m, rhs, &
-       &X, U, U_nl, U_mesh, X_old, X_new, Source, Buoyancy, gravity, Abs, &
+       &X, U, U_nl, U_mesh, X_old, X_new, Source, Buoyancy, hb_density, hb_pressure, gravity, Abs, &
        &Viscosity, P, Rho, surfacetension, q_mesh, &
        &velocity_bc, velocity_bc_type, &
        &pressure_bc, pressure_bc_type, &
@@ -748,6 +766,7 @@
     !! Viscosity
     type(tensor_field) :: Viscosity
     type(scalar_field) :: P, Rho
+    type(scalar_field), intent(in) :: hb_density, hb_pressure
     !! surfacetension
     type(tensor_field) :: surfacetension
     !! field containing the bc values of velocity
@@ -821,7 +840,7 @@
     integer :: start, finish
     
     ! Variable transform times quadrature weights.
-    real, dimension(ele_ngi(U,ele)) :: detwei, detwei_old, detwei_new
+    real, dimension(ele_ngi(U,ele)) :: detwei, detwei_old, detwei_new, coefficient_detwei
     ! Transformed gradient function for velocity.
     real, dimension(ele_loc(U, ele), ele_ngi(U, ele), mesh_dim(U)) :: du_t
     ! Transformed gradient function for grid velocity.
@@ -1252,22 +1271,28 @@
 
     if(have_gravity.and.acceleration.and.assemble_element) then
       ! buoyancy
+      if(subtract_out_reference_profile) then
+         coefficient_detwei = detwei*gravity_magnitude*(ele_val_at_quad(buoyancy, ele)-ele_val_at_quad(hb_density, ele))
+      else
+         coefficient_detwei = detwei*gravity_magnitude*ele_val_at_quad(buoyancy, ele)
+      end if
+
       if (on_sphere) then
       ! If were on a spherical Earth evaluate the direction of the gravity vector
       ! exactly at quadrature points.
         rhs_addto(:, :loc) = rhs_addto(:, :loc) + shape_vector_rhs(u_shape, &
                                     sphere_inward_normal_at_quad_ele(X, ele), &
-                                    detwei*gravity_magnitude*ele_val_at_quad(buoyancy, ele))
+                                    coefficient_detwei)
       else
       
         if(multiphase) then
           rhs_addto(:, :loc) = rhs_addto(:, :loc) + shape_vector_rhs(u_shape, &
                                     ele_val_at_quad(gravity, ele), &
-                                    detwei*gravity_magnitude*ele_val_at_quad(buoyancy, ele)*nvfrac_gi)
+                                    coefficient_detwei*nvfrac_gi)
         else
           rhs_addto(:, :loc) = rhs_addto(:, :loc) + shape_vector_rhs(u_shape, &
                                     ele_val_at_quad(gravity, ele), &
-                                    detwei*gravity_magnitude*ele_val_at_quad(buoyancy, ele))
+                                    coefficient_detwei)
         end if
         
       end if
@@ -1691,7 +1716,7 @@
                         & rhs_addto, Grad_U_mat_q, Div_U_mat_q, X,&
                         & Rho, U, U_nl, U_mesh, P, q_mesh, surfacetension, &
                         & velocity_bc, velocity_bc_type, &
-                        & pressure_bc, pressure_bc_type, &
+                        & pressure_bc, pressure_bc_type, hb_pressure, &
                         & subcycle_m_tensor_addto, nvfrac, &
                         & ele2grad_mat=ele2grad_mat, kappa_mat=kappa_mat, &
                         & inverse_mass_mat=inverse_mass_mat, &
@@ -1704,7 +1729,7 @@
                         & rhs_addto, Grad_U_mat_q, Div_U_mat_q, X,&
                         & Rho, U, U_nl, U_mesh, P, q_mesh, surfacetension, &
                         & velocity_bc, velocity_bc_type, &
-                        & pressure_bc, pressure_bc_type, &
+                        & pressure_bc, pressure_bc_type, hb_pressure, &
                         & subcycle_m_tensor_addto, nvfrac)
             end if
         end if
@@ -1973,7 +1998,7 @@
        & rhs_addto, Grad_U_mat, Div_U_mat, X, Rho, U,&
        & U_nl, U_mesh, P, q_mesh, surfacetension, &
        & velocity_bc, velocity_bc_type, &
-       & pressure_bc, pressure_bc_type, &
+       & pressure_bc, pressure_bc_type, hb_pressure, &
        & subcycle_m_tensor_addto, nvfrac, &
        & ele2grad_mat, kappa_mat, inverse_mass_mat, &
        & viscosity, viscosity_mat)
@@ -2002,6 +2027,7 @@
     integer, dimension(:,:), intent(in) :: velocity_bc_type
     type(scalar_field), intent(in) :: pressure_bc
     integer, dimension(:), intent(in) :: pressure_bc_type
+    type(scalar_field), intent(in) :: hb_pressure
 
     !! Computation of primal fluxes and penalty fluxes
     real, intent(in), optional, dimension(:,:,:) :: ele2grad_mat
@@ -2334,8 +2360,13 @@
        ! add -|  N_i M_j \vec n p_j, where p_j are the prescribed bc values
        !      /
        do dim = 1, U%dim
-          rhs_addto(dim,u_face_l) = rhs_addto(dim,u_face_l) - &
-               matmul( ele_val(pressure_bc, face), mnCT(1,dim,:,:) )
+          if(subtract_out_reference_profile) then
+            rhs_addto(dim,u_face_l) = rhs_addto(dim,u_face_l) - &
+                 matmul( ele_val(pressure_bc, face) - face_val(hb_pressure, face), mnCT(1,dim,:,:) )
+          else
+            rhs_addto(dim,u_face_l) = rhs_addto(dim,u_face_l) - &
+                 matmul( ele_val(pressure_bc, face), mnCT(1,dim,:,:) )
+          end if
        end do
     end if
 

=== modified file 'assemble/Momentum_Equation.F90'
--- assemble/Momentum_Equation.F90	2013-03-03 01:36:16 +0000
+++ assemble/Momentum_Equation.F90	2013-07-10 19:09:26 +0000
@@ -1515,6 +1515,11 @@
          type(vector_field) :: delta_u
          type(vector_field), pointer :: positions
 
+         ! Fields for the subtract_out_reference_profile option under the Velocity field
+         type(scalar_field), pointer :: hb_pressure
+         type(scalar_field) :: combined_p
+         integer :: stat
+
          ewrite(1,*) 'Entering advance_velocity'
 
 
@@ -1539,6 +1544,20 @@
                &have_option('/ocean_forcing/shelf')) then
             ewrite(1,*) "shelf: Entering compute_pressure_and_tidal_gradient"
                call compute_pressure_and_tidal_gradient(state(istate), delta_u, ct_m(istate)%ptr, p_theta, x)
+            else if (have_option(trim(state(istate)%option_path)//'/equation_of_state/compressible/subtract_out_reference_profile')) then
+               ! Splits up the Density and Pressure fields into a hydrostatic component (') and a perturbed component (''). 
+               ! The hydrostatic components, denoted p' and rho', should satisfy the balance: grad(p') = rho'*g
+               ! We subtract the hydrostatic component from the pressure used in the pressure gradient term of the momentum equation.
+               hb_pressure => extract_scalar_field(state(istate), "HydrostaticReferencePressure", stat)
+               if(stat /= 0) then
+                  FLExit("When using the subtract_out_reference_profile option, please set a (prescribed) HydrostaticReferencePressure field.")
+                  ewrite(-1,*) 'The HydrostaticReferencePressure field, defining the hydrostatic component of the pressure field, needs to be set.'
+               end if
+               call allocate(combined_p,p_theta%mesh, "PressurePerturbation")
+               call set(combined_p, p_theta)
+               call addto(combined_p, hb_pressure, scale=-1.0)
+               call mult_T(delta_u, ct_m(istate)%ptr, combined_p)
+               call deallocate(combined_p)
             else
                call mult_T(delta_u, ct_m(istate)%ptr, p_theta)
             end if

=== modified file 'manual/configuring_fluidity.tex'
--- manual/configuring_fluidity.tex	2013-05-29 18:26:30 +0000
+++ manual/configuring_fluidity.tex	2013-07-10 19:09:26 +0000
@@ -2483,6 +2483,17 @@
 
 \emph{Important note when using 'stress form' and 'partial stress form' with a CG velocity field:} when using isotropic viscosities, and not the 'tensor form' of the viscosity term, all components of viscosity must be set to equal to the isotropic viscosity due to the method of implementation. If using an anisotropic viscosity with the 'partial stress form' or 'stress form' consult the schema and query the code to understand how the viscosity tensor must be defined. This does not apply when using DG discretisations. 
 
+\subsubsection{Hydrostatic balance}
+The option \option{subtract\_out\_reference\_profile} splits up the Density and Pressure fields into a hydrostatic (reference) component (') and a perturbed component (''). The hydrostatic (reference) components, denoted $p^{\prime}$ and $\rho^{\prime}$, should satisfy the balance:
+\begin{equation}
+  \nabla p^{\prime} = \rho^{\prime}\mathbf{g}
+\end{equation}
+Enabling this option will subtract the hydrostatic components, specified here, from the pressure and density used in the pressure gradient and buoyancy terms in the momentum equation. This helps to maintain hydrostatic balance and prevent spurious oscillations in the pressure field when using unbalanced finite element pairs.
+
+To use this option you will also need to create two prescribed scalar fields, called HydrostaticReferencePressure and HydrostaticReferenceDensity, which define $p^{\prime}$ and $\rho^{\prime}$. These must be on the same mesh as pressure and density, respectively, and are meant to be set up as vertical profiles (i.e. constant in the horizontal direction). 
+
+Note that unlike all the other hydrostatic/geostrophic balance options in Fluidity (i.e. \option{subtract\_out\_hydrostatic\_level} under the linear incompressible EoS option, or with HydrostaticPressure or GeostrophicPressure fields), the hydrostatic pressure is not subtracted from the Pressure field itself. In other words, the Pressure field that gets solved for (and output in the .vtu files) is still the combined Pressure ($p = p^{\prime} + \rho^{\prime}$), and the hydrostatic pressure $p^{\prime}$ is only subtracted in the momentum equation.
+
 \subsection{InternalEnergy options}
 \begin{itemize}
 \item The optional heat flux term can be included via the InternalEnergy's Diffusivity field. Users will need to specify an isotropic value for $\frac{k}{C_v}$, where $k$ is the effective conductivity and $C_v$ is the specific heat at constant volume.

=== modified file 'schemas/equation_of_state.rnc'
--- schemas/equation_of_state.rnc	2012-05-30 07:26:05 +0000
+++ schemas/equation_of_state.rnc	2013-07-10 19:09:26 +0000
@@ -91,62 +91,76 @@
    (
       ## Equations of state for compressible applications
       element compressible {
-         ## Stiffened Gas EoS
-         ##
-         ## Used with compressible simulations
-         element stiffened_gas {
-            ## reference uncompressed density
-            ##
-            ## if activated then either Liquid EoS or
-            ## full Stiffened Gas EoS
-            element reference_density {
-               real
-            }?,
-            ## Ratio of specific heats at constant
-            ## pressure to that at constant volume
-            ##
-            ## Requires an energy field.
-            ## If activated then a full Stiffened Gas EoS 
-            element ratio_specific_heats {
-               real
-            }?,
-            ## bulk_sound_speed_squared = isothermal_bulk_modulus/reference_density
-            ##
-            ## if activated then either full or partial Liquid EoS or full
+         (
             ## Stiffened Gas EoS
-            element bulk_sound_speed_squared {
-               real
-            }?
-         }|
-         ## Giraldo et. al. 2008, J. Comp. Phys.
-         ##
-         ## Used with compressible simulations
-         element giraldo {
-            ## reference pressure
-            ##
-            ## The pressure at which potential temperature
-            ## equals actual temperature. Normally taken to
-            ## be the surface pressure, 1.0e05.
-            element reference_pressure {
-               real
-            },
-            ## Heat capacity at constant pressure
-            element C_P {
-               real
-            },
-            ## Heat capacity at constant volume
-            element C_V {
-               real
+            ##
+            ## Used with compressible simulations
+            element stiffened_gas {
+               ## reference uncompressed density
+               ##
+               ## if activated then either Liquid EoS or
+               ## full Stiffened Gas EoS
+               element reference_density {
+                  real
+               }?,
+               ## Ratio of specific heats at constant
+               ## pressure to that at constant volume
+               ##
+               ## Requires an energy field.
+               ## If activated then a full Stiffened Gas EoS 
+               element ratio_specific_heats {
+                  real
+               }?,
+               ## bulk_sound_speed_squared = isothermal_bulk_modulus/reference_density
+               ##
+               ## if activated then either full or partial Liquid EoS or full
+               ## Stiffened Gas EoS
+               element bulk_sound_speed_squared {
+                  real
+               }?
+            }|
+            ## Giraldo et. al. 2008, J. Comp. Phys.
+            ##
+            ## Used with compressible simulations
+            element giraldo {
+               ## reference pressure
+               ##
+               ## The pressure at which potential temperature
+               ## equals actual temperature. Normally taken to
+               ## be the surface pressure, 1.0e05.
+               element reference_pressure {
+                  real
+               },
+               ## Heat capacity at constant pressure
+               element C_P {
+                  real
+               },
+               ## Heat capacity at constant volume
+               element C_V {
+                  real
+               }
+            }|
+            ## Foam EoS
+            ## Used with compressible simulations of liquid drainage in foams.
+            ## It describes the liquid content in the foam as the product of the 
+            ## Plateau border cross sectional area and the local Plateau 
+            ## border length per unit volume (lambda). 
+            element foam {
+               empty
             }
-         }|
-         ## Foam EoS
-         ## Used with compressible simulations of liquid drainage in foams.
-         ## It describes the liquid content in the foam as the product of the 
-         ## Plateau border cross sectional area and the local Plateau 
-         ## border length per unit volume (lambda). 
-         element foam {
+         ),
+         ## This splits up the Density and Pressure fields into a hydrostatic (reference) component (') and a perturbed component (''). 
+         ## The hydrostatic (reference) components, denoted p' and rho', should satisfy the balance:
+         ## grad(p') = rho'*g
+         ##
+         ## Enabling this option will subtract the hydrostatic components, specified here, from the pressure and density used in the pressure gradient and buoyancy terms in the momentum equation. This helps to maintain hydrostatic balance and prevent spurious oscillations in the pressure field when using unbalanced finite element pairs.
+         ##
+         ## Note 1: You will also need to create two prescribed scalar fields, called HydrostaticReferencePressure and HydrostaticReferenceDensity, which define p' and rho'. These must be on the same mesh as pressure and density, respectively. Also, these fields are meant to be time-independent and set up as vertical profiles (i.e. constant in the horizontal direction).
+         ##
+         ## Note 2: Unlike all the other hydrostatic/geostrophic balance options in Fluidity (i.e. subtract_out_hydrostatic_level under the linear incompressible EoS option, or with HydrostaticPressure or GeostrophicPressure fields), the hydrostatic pressure is not subtracted from the Pressure field itself. In other words, the Pressure field that gets solved for (and output in the .vtu files) is still the combined Pressure (p = p' + p''), and the hydrostatic pressure p' is only subtracted in the momentum equation.
+         element subtract_out_reference_profile {
             empty
-         }  
+         }?
       }   
    )
    

=== modified file 'schemas/equation_of_state.rng'
--- schemas/equation_of_state.rng	2012-05-30 07:26:05 +0000
+++ schemas/equation_of_state.rng	2013-07-10 19:09:26 +0000
@@ -162,6 +162,20 @@
           <empty/>
         </element>
       </choice>
+      <optional>
+        <element name="subtract_out_reference_profile">
+          <a:documentation>This splits up the Density and Pressure fields into a hydrostatic (reference) component (') and a perturbed component (''). 
+The hydrostatic (reference) components, denoted p' and rho', should satisfy the balance:
+grad(p') = rho'*g
+
+Enabling this option will subtract the hydrostatic components, specified here, from the pressure and density used in the pressure gradient and buoyancy terms in the momentum equation. This helps to maintain hydrostatic balance and prevent spurious oscillations in the pressure field when using unbalanced finite element pairs.
+
+Note 1: You will also need to create two prescribed scalar fields, called HydrostaticReferencePressure and HydrostaticReferenceDensity, which define p' and rho'. These must be on the same mesh as pressure and density, respectively. Also, these fields are meant to be time-independent and set up as vertical profiles (i.e. constant in the horizontal direction).
+
+Note 2: Unlike all the other hydrostatic/geostrophic balance options in Fluidity (i.e. subtract_out_hydrostatic_level under the linear incompressible EoS option, or with HydrostaticPressure or GeostrophicPressure fields), the hydrostatic pressure is not subtracted from the Pressure field itself. In other words, the Pressure field that gets solved for (and output in the .vtu files) is still the combined Pressure (p = p' + p''), and the hydrostatic pressure p' is only subtracted in the momentum equation.</a:documentation>
+          <empty/>
+        </element>
+      </optional>
     </element>
   </define>
 </grammar>

=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2013-05-31 07:10:50 +0000
+++ schemas/prognostic_field_options.rnc	2013-07-10 19:09:26 +0000
@@ -1200,6 +1200,7 @@
             }?
          }?
       }?,
+
       ## SurfaceTension
       element tensor_field {
          attribute name { "SurfaceTension" },

=== added directory 'tests/inlet_velocity_bc_compressible'
=== added file 'tests/inlet_velocity_bc_compressible/Makefile'
--- tests/inlet_velocity_bc_compressible/Makefile	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible/Makefile	2013-07-10 19:09:26 +0000
@@ -0,0 +1,17 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=100.0 -- 0.0 10000.0 line
+
+run:
+	@echo **********Running simulation
+	../../bin/fluidity -v2 -l inlet_velocity_bc_compressible_1d.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+

=== added file 'tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible.xml'
--- tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible.xml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible.xml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,63 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>inlet_velocity_bc_compressible</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>  
+    <variable name="gas_velocity_max_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_1d.stat")
+gas_velocity_max_1d = s["Gas"]["Velocity%magnitude"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_velocity_min_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_1d.stat")
+gas_velocity_min_1d = s["Gas"]["Velocity%magnitude"]["min"][-1]
+    </variable>
+
+    <variable name="gas_density_max_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_1d.stat")
+gas_density_max_1d = s["Gas"]["Density"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_density_min_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_1d.stat")
+gas_density_min_1d = s["Gas"]["Density"]["min"][-1]
+    </variable>
+
+    <variable name="solvers_converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Gas::Velocity is at least 300 m/s everywhere in the 1D test case" language="python">
+assert(abs(gas_velocity_max_1d) &lt;= 800.0 and abs(gas_velocity_min_1d) &gt;= 300.0)
+    </test>
+
+    <test name="Gas::Density is between 0.05 and 0.25 kg/m^3 everywhere in the 1D test case" language="python">
+assert(abs(gas_density_max_1d) &lt;= 0.25 and abs(gas_density_min_1d) &gt;= 0.05)
+    </test>
+
+    <test name="Solvers converged" language="python">
+assert(solvers_converged)
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>

=== added file 'tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible_1d.flml'
--- tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible_1d.flml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible/inlet_velocity_bc_compressible_1d.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,537 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">inlet_velocity_bc_compressible_1d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="GeostrophicMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.01</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">600.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.5</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.8</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="1" dim1="dim" rank="1">-1</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+        <subtract_out_reference_profile/>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="DensityMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  rho = p/((gamma-1.0)*e)
+  return rho</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.181985215</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density"/>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1687200</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="HydrostaticReferencePressure" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="HydrostaticReferenceDensity" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  rho = p/((gamma-1.0)*e)
+  return rho</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added directory 'tests/inlet_velocity_bc_compressible_without_gravity'
=== added file 'tests/inlet_velocity_bc_compressible_without_gravity/Makefile'
--- tests/inlet_velocity_bc_compressible_without_gravity/Makefile	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible_without_gravity/Makefile	2013-07-10 19:09:26 +0000
@@ -0,0 +1,21 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=100.0 -- 0.0 10000.0 line
+	@echo **********Creating 2D mesh
+	gmsh -2 -o inlet_velocity_bc_compressible_without_gravity_pseudo1d.msh src/inlet_velocity_bc_compressible_without_gravity_pseudo1d.geo
+	../../bin/gmsh2triangle --2d inlet_velocity_bc_compressible_without_gravity_pseudo1d.msh
+
+run:
+	@echo **********Running simulation
+	../../bin/fluidity inlet_velocity_bc_compressible_without_gravity_1d.flml
+	../../bin/fluidity inlet_velocity_bc_compressible_without_gravity_pseudo1d.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+

=== added file 'tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity.xml'
--- tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity.xml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity.xml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,95 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>inlet_velocity_bc_compressible_without_gravity</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>  
+    <variable name="gas_velocity_max_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_1d.stat")
+gas_velocity_max_1d = s["Gas"]["Velocity%magnitude"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_velocity_min_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_1d.stat")
+gas_velocity_min_1d = s["Gas"]["Velocity%magnitude"]["min"][-1]
+    </variable>
+
+    <variable name="gas_velocity_max_pseudo1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_pseudo1d.stat")
+gas_velocity_max_pseudo1d = s["Gas"]["Velocity%magnitude"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_velocity_min_pseudo1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_pseudo1d.stat")
+gas_velocity_min_pseudo1d = s["Gas"]["Velocity%magnitude"]["min"][-1]
+    </variable>
+
+    <variable name="gas_density_max_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_1d.stat")
+gas_density_max_1d = s["Gas"]["Density"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_density_min_1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_1d.stat")
+gas_density_min_1d = s["Gas"]["Density"]["min"][-1]
+    </variable>
+
+    <variable name="gas_density_max_pseudo1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_pseudo1d.stat")
+gas_density_max_pseudo1d = s["Gas"]["Density"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_density_min_pseudo1d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_compressible_without_gravity_pseudo1d.stat")
+gas_density_min_pseudo1d = s["Gas"]["Density"]["min"][-1]
+    </variable>
+
+    <variable name="solvers_converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Gas::Velocity is 300 m/s everywhere in the 1D test case" language="python">
+assert(abs(gas_velocity_max_1d - 300) &lt; 1 and abs(gas_velocity_min_1d - 300) &lt; 1)
+    </test>
+
+    <test name="Gas::Velocity is 300 m/s everywhere in the pseudo-1D test case" language="python">
+assert(abs(gas_velocity_max_pseudo1d - 300) &lt; 1 and abs(gas_velocity_min_pseudo1d - 300) &lt; 1)
+    </test>
+
+    <test name="Gas::Density is 0.181985215 kg/m^3 everywhere in the 1D test case" language="python">
+assert(abs(gas_density_max_1d - 0.181985215) &lt; 1e-2 and abs(gas_density_min_1d - 0.181985215) &lt; 1e-2)
+    </test>
+
+    <test name="Gas::Density is 0.181985215 kg/m^3 everywhere in the pseudo-1D test case" language="python">
+assert(abs(gas_density_max_pseudo1d - 0.181985215) &lt; 1e-2 and abs(gas_density_min_pseudo1d - 0.181985215) &lt; 1e-2)
+    </test>
+
+    <test name="Solvers converged" language="python">
+assert(solvers_converged)
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>

=== added file 'tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_1d.flml'
--- tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_1d.flml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_1d.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,436 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">inlet_velocity_bc_compressible_without_gravity_1d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="GeostrophicMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">10.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.01</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">600.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.5</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1.0e-6</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">101325</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">101325</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="DensityMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.757877001</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.181985215</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">1.85e-5</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density"/>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1687200</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_pseudo1d.flml'
--- tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_pseudo1d.flml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible_without_gravity/inlet_velocity_bc_compressible_without_gravity_pseudo1d.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,430 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">inlet_velocity_bc_compressible_without_gravity_pseudo1d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="inlet_velocity_bc_compressible_without_gravity_pseudo1d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">10.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">600.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-7</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.5</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1.0e-6</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">101325</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">333</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">101325</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.757877001</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Vent">
+          <surface_ids>
+            <integer_value shape="1" rank="1">999</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.181985215</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix>
+                <use_submesh/>
+              </lump_mass_matrix>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <stress_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Vent">
+          <surface_ids>
+            <integer_value shape="1" rank="1">999</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+              <y_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">111</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">222</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="2 2" dim1="dim" rank="2">1.85e-5 1.85e-5 1.85e-5 1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density"/>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Vent">
+          <surface_ids>
+            <integer_value shape="1" rank="1">999</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1687200</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added directory 'tests/inlet_velocity_bc_compressible_without_gravity/src'
=== added file 'tests/inlet_velocity_bc_compressible_without_gravity/src/inlet_velocity_bc_compressible_without_gravity_pseudo1d.geo'
--- tests/inlet_velocity_bc_compressible_without_gravity/src/inlet_velocity_bc_compressible_without_gravity_pseudo1d.geo	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_compressible_without_gravity/src/inlet_velocity_bc_compressible_without_gravity_pseudo1d.geo	2013-07-10 19:09:26 +0000
@@ -0,0 +1,19 @@
+dx = 100.0;
+Point(1) = {0.0, 0.0, 0.0, dx};
+
+Extrude {dx,0,0} {
+  Point{1}; Layers{dx/dx};
+}
+Extrude {0,10000,0} {
+  Line{1}; Layers{10000/dx};
+}
+
+// Top
+Physical Line(333) = {2};
+// Sides
+Physical Line(111) = {3};
+Physical Line(222) = {4};
+// Bottom
+Physical Line(999) = {1};
+
+Physical Surface(1000) = {5};

=== added directory 'tests/inlet_velocity_bc_incompressible'
=== added file 'tests/inlet_velocity_bc_incompressible/Makefile'
--- tests/inlet_velocity_bc_incompressible/Makefile	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_incompressible/Makefile	2013-07-10 19:09:26 +0000
@@ -0,0 +1,17 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=100.0 -- 0.0 10000.0 line
+
+run:
+	@echo **********Running simulation
+	../../bin/fluidity -v2 -l inlet_velocity_bc_incompressible.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+

=== added file 'tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.flml'
--- tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.flml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,381 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">inlet_velocity_bc_incompressible</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="GeostrophicMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.01</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">20.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.25</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.8</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="1" dim1="dim" rank="1">-1</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">1.23</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = 101325 - 1.23*g*h
+  return p</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = 101325 - 1.23*g*h
+  return p</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <stress_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.xml'
--- tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.xml	1970-01-01 00:00:00 +0000
+++ tests/inlet_velocity_bc_incompressible/inlet_velocity_bc_incompressible.xml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>inlet_velocity_bc_incompressible</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="short" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>  
+    <variable name="gas_velocity_max" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_incompressible.stat")
+gas_velocity_max = s["Gas"]["Velocity%1"]["max"][-1]
+    </variable>
+  
+    <variable name="gas_velocity_min" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("inlet_velocity_bc_incompressible.stat")
+gas_velocity_min = s["Gas"]["Velocity%1"]["min"][-1]
+    </variable>
+
+    <variable name="solvers_converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Gas::Velocity is 300 m/s everywhere" language="python">
+assert(abs(gas_velocity_max - 300) &lt; 1e-7 and abs(gas_velocity_min - 300) &lt; 1e-7)
+    </test>
+     
+    <test name="Solvers converged" language="python">
+assert(solvers_converged)
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>

=== added directory 'tests/mphase_inlet_velocity_bc_compressible'
=== added file 'tests/mphase_inlet_velocity_bc_compressible/Makefile'
--- tests/mphase_inlet_velocity_bc_compressible/Makefile	1970-01-01 00:00:00 +0000
+++ tests/mphase_inlet_velocity_bc_compressible/Makefile	2013-07-10 19:09:26 +0000
@@ -0,0 +1,17 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=1000.0 -- 0.0 10000.0 line
+
+run:
+	@echo **********Running simulation
+	../../bin/fluidity mphase_inlet_velocity_bc_compressible.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+

=== added file 'tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml'
--- tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,863 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">mphase_inlet_velocity_bc_compressible</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">100.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">50</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-4</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.8</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="1" dim1="dim" rank="1">-1.0</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+        <subtract_out_reference_profile/>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">101325</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="DensityMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  rho = p/((gamma-1.0)*e)
+  return rho</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.181985215</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density">
+            <discretisation_options>
+              <spatial_discretisation>
+                <control_volumes>
+                  <face_value name="FirstOrderUpwind"/>
+                </control_volumes>
+              </spatial_discretisation>
+              <temporal_discretisation>
+                <theta>
+                  <real_value rank="0">1.0</real_value>
+                </theta>
+                <control_volumes/>
+              </temporal_discretisation>
+            </discretisation_options>
+          </density>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1687200</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="HydrostaticReferencePressure" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  return p</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="HydrostaticReferenceDensity" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp
+  T0 = 288.15 # K
+  R = 287.058 # Specific gas constant in J/Kg
+  g = 9.8 # m/s**2
+  h = X[0] # m
+  p0 = 101325
+  e = T0*1406.0
+  gamma = 1.33
+  p = p0*exp(-(g*h)/((gamma-1.0)*e))
+  rho = p/((gamma-1.0)*e)
+  return rho</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <multiphase_properties>
+      <effective_conductivity>
+        <real_value rank="0">0.5</real_value>
+      </effective_conductivity>
+      <specific_heat>
+        <real_value rank="0">1406</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <material_phase name="Particles">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2400</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Gas" field_name="Pressure"/>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">300</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.5</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density">
+            <discretisation_options>
+              <spatial_discretisation>
+                <control_volumes>
+                  <face_value name="FirstOrderUpwind"/>
+                </control_volumes>
+              </spatial_discretisation>
+              <temporal_discretisation>
+                <theta>
+                  <real_value rank="0">1.0</real_value>
+                </theta>
+                <control_volumes/>
+              </temporal_discretisation>
+            </discretisation_options>
+          </density>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">274895.1</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1144800</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-7</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1.0e-4</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <multiphase_properties>
+      <particle_diameter>
+        <real_value rank="0">2e-4</real_value>
+      </particle_diameter>
+      <specific_heat>
+        <real_value rank="0">954</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <multiphase_interaction>
+    <fluid_particle_drag>
+      <drag_correlation name="wen_yu"/>
+    </fluid_particle_drag>
+    <heat_transfer>
+      <heat_transfer_coefficient name="gunn"/>
+    </heat_transfer>
+  </multiphase_interaction>
+</fluidity_options>

=== added file 'tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.xml'
--- tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.xml	1970-01-01 00:00:00 +0000
+++ tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.xml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,48 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>mphase_inlet_velocity_bc_compressible</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>  
+    <variable name="gas_velocity_max" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_inlet_velocity_bc_compressible.stat")
+gas_velocity_max = s["Gas"]["Velocity%1"]["max"][-1]
+    </variable>
+    <variable name="gas_velocity_min" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_inlet_velocity_bc_compressible.stat")
+gas_velocity_min = s["Gas"]["Velocity%1"]["min"][-1]
+    </variable>
+    <variable name="solvers_converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Gas::Velocity is at most ~300 m/s" language="python">
+assert(abs(gas_velocity_max) &lt; 305)
+    </test>
+    <test name="Gas::Velocity is at least 75 m/s" language="python">
+# Note: the velocity won't be 300 m/s everywhere because of the effects of changing density and pressure.
+assert(abs(gas_velocity_min) &gt; 75)
+    </test>
+    <test name="Solvers converged" language="python">
+assert(solvers_converged)
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>

=== added directory 'tests/mphase_strong_pressure_bc_compressible'
=== added file 'tests/mphase_strong_pressure_bc_compressible/Makefile'
--- tests/mphase_strong_pressure_bc_compressible/Makefile	1970-01-01 00:00:00 +0000
+++ tests/mphase_strong_pressure_bc_compressible/Makefile	2013-07-10 19:09:26 +0000
@@ -0,0 +1,19 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=100.0 -- 0.0 10000.0 line
+
+run:
+	@echo **********Running P2-P1 simulation
+	../../bin/fluidity mphase_strong_pressure_bc_compressible_p2p1.flml
+	@echo **********Running P0-P1 simulation
+	../../bin/fluidity mphase_strong_pressure_bc_compressible_p0p1.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+

=== added file 'tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible.xml'
--- tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible.xml	1970-01-01 00:00:00 +0000
+++ tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible.xml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,73 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>mphase_strong_pressure_bc_compressible</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>  
+    <variable name="gas_velocity_max_p2p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p2p1.stat")
+gas_velocity_max_p2p1 = s["Gas"]["Velocity%1"]["max"][-1]
+    </variable>
+    <variable name="gas_pressure_max_p2p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p2p1.stat")
+gas_pressure_max_p2p1 = s["Gas"]["Pressure"]["max"][-1]
+    </variable>
+    <variable name="gas_pressure_min_p2p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p2p1.stat")
+gas_pressure_min_p2p1 = s["Gas"]["Pressure"]["min"][-1]
+    </variable>
+    <variable name="gas_velocity_max_p0p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p0p1.stat")
+gas_velocity_max_p0p1 = s["Gas"]["Velocity%1"]["max"][-1]
+    </variable>
+    <variable name="gas_pressure_max_p0p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p0p1.stat")
+gas_pressure_max_p0p1 = s["Gas"]["Pressure"]["max"][-1]
+    </variable>
+    <variable name="gas_pressure_min_p0p1" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_strong_pressure_bc_compressible_p0p1.stat")
+gas_pressure_min_p0p1 = s["Gas"]["Pressure"]["min"][-1]
+    </variable>
+    <variable name="solvers_converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Gas::Velocity is always zero everywhere (P2-P1 simulation)" language="python">
+assert(abs(gas_velocity_max_p2p1) &lt; 1e-7)
+    </test>    
+    <test name="Gas::Pressure is always 101325 Pa everywhere (P2-P1 simulation)" language="python">
+assert(abs(gas_pressure_max_p2p1) - 101325 &lt; 1.0e-7 and abs(gas_pressure_min_p2p1) - 101325 &lt; 1.0e-7)
+    </test>
+    <test name="Gas::Velocity is always zero everywhere (P0-P1 simulation)" language="python">
+assert(abs(gas_velocity_max_p0p1) &lt; 1e-7)
+    </test>    
+    <test name="Gas::Pressure is always 101325 Pa everywhere (P0-P1 simulation)" language="python">
+assert(abs(gas_pressure_max_p0p1) - 101325 &lt; 1.0e-7 and abs(gas_pressure_min_p0p1) - 101325 &lt; 1.0e-7)
+    </test>
+    <test name="Solvers converged" language="python">
+assert(solvers_converged)
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>

=== added file 'tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml'
--- tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,735 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">mphase_strong_pressure_bc_compressible_p0p1</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">0</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="GeostrophicMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">600.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.25</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+  </timestepping>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">101325</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">101325</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="DensityMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.757877001</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.757877001</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density">
+            <discretisation_options>
+              <spatial_discretisation>
+                <control_volumes>
+                  <face_value name="FirstOrderUpwind"/>
+                </control_volumes>
+              </spatial_discretisation>
+              <temporal_discretisation>
+                <theta>
+                  <real_value rank="0">1.0</real_value>
+                </theta>
+                <control_volumes/>
+              </temporal_discretisation>
+            </discretisation_options>
+          </density>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <multiphase_properties>
+      <effective_conductivity>
+        <real_value rank="0">0.5</real_value>
+      </effective_conductivity>
+      <specific_heat>
+        <real_value rank="0">1406</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <material_phase name="Particles">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2400</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Gas" field_name="Pressure"/>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+              <tensor_form/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.5</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density">
+            <discretisation_options>
+              <spatial_discretisation>
+                <control_volumes>
+                  <face_value name="FirstOrderUpwind"/>
+                </control_volumes>
+              </spatial_discretisation>
+              <temporal_discretisation>
+                <theta>
+                  <real_value rank="0">1.0</real_value>
+                </theta>
+                <control_volumes/>
+              </temporal_discretisation>
+            </discretisation_options>
+          </density>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">274895.1</real_value>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">1.0e-7</real_value>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <multiphase_properties>
+      <particle_diameter>
+        <real_value rank="0">2e-4</real_value>
+      </particle_diameter>
+      <specific_heat>
+        <real_value rank="0">954</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml'
--- tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml	2013-07-10 19:09:26 +0000
@@ -0,0 +1,742 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">mphase_strong_pressure_bc_compressible_p2p1</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="GeostrophicMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="PressureMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">600.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">0.25</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="VelocityMesh"/>
+      </courant_number>
+      <minimum_timestep>
+        <real_value rank="0">0.001</real_value>
+      </minimum_timestep>
+      <increase_tolerance>
+        <real_value rank="0">1.1</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+  </timestepping>
+  <material_phase name="Gas">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.33</real_value>
+          </ratio_specific_heats>
+        </stiffened_gas>
+      </compressible>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_compressible_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">101325</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">101325</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="DensityMesh"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.757877001</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Right">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">0.757877001</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind_petrov_galerkin>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind_petrov_galerkin>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <stress_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="1" dim1="dim" rank="1">0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Left">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_asymmetric>
+                <constant>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.85e-5</real_value>
+                </constant>
+              </anisotropic_asymmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density">
+            <discretisation_options>
+              <spatial_discretisation>
+                <control_volumes>
+                  <face_value name="FirstOrderUpwind"/>
+                </control_volumes>
+              </spatial_discretisation>
+              <temporal_discretisation>
+                <theta>
+                  <real_value rank="0">1.0</real_value>
+                </theta>
+                <control_volumes/>
+              </temporal_discretisation>
+            </discretisation_options>
+          </density>
+        </equation>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FirstOrderUpwind"/>
+            <diffusion_scheme name="ElementGradient"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">405138.9</real_value>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <multiphase_properties>
+      <effective_conductivity>
+        <real_value rank="0">0.5</real_value>
+      </effective_conductivity>
+      <specific_heat>
+        <real_value rank="0">1406</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <material_phase name="Particles">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2400</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Gas" field_name="Pressure"/>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind_petrov_galerkin>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind_petrov_galerkin>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>