Merge lp:~ctjacobs-multiphase/fluidity/compressible-multiphase into lp:fluidity

[Christian Jacobs]

Proposing merge of compressible-multiphase branch revisions r4078 to r4108 (inclusive) into trunk.

Summary of changes:
===================
- Added CV support for the multiphase InternalEnergy equation. This is tested with the new MMS test mphase_mms_p2p1_compressible_ie_p1cv.

- Changed -p*div(vfrac*u) to -p*vfrac*div(u) in the InternalEnergy equation.

- Added the option to include or exclude the -p*vfrac*div(u) term in the schema. This term should only be included in the fluid phase.

- Added an options check to make sure the prognostic PhaseVolumeFraction field is in the particle phase(s) only.

- Included the PhaseVolumeFraction in the Diffusivity term if we are solving the multiphase InternalEnergy equation. This allows the user to add in the heat flux term diff( k/Cv*vfrac*grad(ie) ) by specifying an isotropic diffusivity k/Cv.

- The MMS test mphase_mms_p2p1_compressible_ie_heat_flux will not be merged into the trunk, but the .flml will replace the one in longtests/mphase_mms_p2p1_compressible_ie to test the heat flux term.

[lp:~ctjacobs-multiphase/fluidity/compressible-multiphase]

4109. By Christian Jacobs on 2012-10-23

Improvements to the mphase_rogue_shock_tube_dense_bed_nylon test:
- Switch to stress_form in the Gas phase.
- Include a heat flux term for the Gas phase's InternalEnergy equation.
- Set zero_flux BCs for the prognostic PhaseVolumeFraction so no particles escape.

4110. By Christian Jacobs on 2012-10-25

Added the .msh file.

4111. By Christian Jacobs on 2012-10-26

Improvements to the mphase_rogue_shock_tube_dense_bed_glass test:
- Switch to stress_form in the Gas phase.
- Include a heat flux term for the Gas phase's InternalEnergy equation.
- Set zero_flux BCs for the prognostic PhaseVolumeFraction so no particles escape.

Also added the .msh file.

4112. By Christian Jacobs on 2012-10-26

Corrected the value of k/Cv for the heat flux term.

4113. By Christian Jacobs on 2012-10-26

Switched to stress_form in the Gas phase.

4114. By Christian Jacobs on 2012-10-27

Added documentation on the heat flux term and the include_pressure_term option for the InternalEnergy field.

4115. By Christian Jacobs on 2012-10-27

Merging trunk revisions (from r4095 to r4108 inclusive) into branch.

4116. By Christian Jacobs on 2012-11-06

Merging trunk revisions (from r4109 to r4125 inclusive) into branch.

4117. By Christian Jacobs on 2012-11-28

Merging trunk revisions (from r4126 to r4134 inclusive) into branch.

4118. By Christian Jacobs on 2012-12-15

Corrected a typo in 'make clean' definition.

4119. By Christian Jacobs on 2012-12-15

Added tephra settling paper reference to the section on the tephra_settling example, and also to the section on Fluidity's incompressible multiphase flow model.

4120. By Christian Jacobs on 2012-12-15

Merging trunk revisions (from r4135 to r4144 inclusive) into branch.

4121. By Christian Jacobs on 2012-12-15

Added another compressible multiphase MMS test which, unlike the other MMS tests, uses a prognostic PhaseVolumeFraction field instead of a prescribed one.

4122. By Christian Jacobs on 2012-12-16

Added an MMS test to check the implementation of the heat transfer term by Gunn (1978). This will be moved to the longtests directory once the compressible-multiphase branch is merged with the trunk.

4123. By Christian Jacobs on 2013-01-04

Merging trunk revisions (from r4145 to r4147 inclusive) into branch.

4124. By Christian Jacobs on 2013-01-07

Added a Python script to track the position of the lower and upper front of the 2 cm particle bed in the mphase_rogue_shock_tube_dense_bed_nylon simulation.

The upper front is defined as the position where the particle PhaseVolumeFraction reaches 0.3 (approximately the same value used in the article by Rogue et al. (1998)). For the lower front, a PhaseVolumeFraction value of 0.01 is used since this area is more dilute in comparison to the upper front.

Fluidity's results do not match well with the experimental data, probably because the mesh has been kept quite coarse. Furthermore, Rogue et al. (1998) point out that "the important dilution of the cloud in the lower front area makes the comparison with the experimental position quite ambiguous", so the numerical results should be viewed with caution: http://amcg.ese.ic.ac.uk/~ctj10/images/fronts.png

[Stephan Kramer (s-kramer)]

About assemble/Advection_Diffusion_CG.F90:

I'm a bit confused what you've done with the diffusion term. In the comments you claim to implement div( (k/Cv) * vfrac * grad(ie) ). Now k an d Cv are passed down indeed (missing intent in line 531 btw) but not actually used anywhere.

About the new 'compressible' logical what is this used for? At the moment it only seems to decide whether to include the p div u term. But that's now also controlled by the include_pressure_term logical that is set by a new option. Does this mean existing single phase compressible now need to set this option as well? What happens if we're in the (non-compressible) particle phase and still set this option: at the moment it'll just ignore it, do we want an error? Afaics there are two options to determine whether to include the term: 1) we base it purely on the compressible logical and don't have an option 2) we base it purely on the new option. Hm 2) is probably too error prone and if we combine 2) with an error if the user gets it wrong then we might as well do 1). So maybe having an option wasn't such a good idea after all...thoughts?

[lp:~ctjacobs-multiphase/fluidity/compressible-multiphase]

4125. By Christian Jacobs on 2013-01-26

Merging trunk revisions (from r4148 to r4159 inclusive) into branch.

4126. By Christian Jacobs on 2013-01-26

Updated the bibliography.

4127. By Christian Jacobs on 2013-01-27

The code now includes the p*div(u) term based on the compressible logical only. The include_pressure_term option has been removed.

[Christian Jacobs (ctjacobs)]

The code now includes the p*div(u) term based on the compressible logical only. The include_pressure_term option has been removed.

Regarding the diffusion term, the idea was for the user to give k/Cv in the .flml file. I've now removed the unused k and Cv variables, since k/Cv will be supplied in the 'diffusivity' tensor.

[lp:~ctjacobs-multiphase/fluidity/compressible-multiphase]

4128. By Christian Jacobs on 2013-01-28

Bug fix: use the OldPhaseVolumeFraction field here instead.

4129. By Christian Jacobs on 2013-01-28

Bug fix: there should be brackets surrounding the material_phase index in this options check.

4130. By Christian Jacobs on 2013-02-08

Merging trunk revisions (from r4160 to r4172 inclusive) into branch.

[Stephan Kramer (s-kramer)]

All looks very good. Ship it!

review: Approve

[lp:~ctjacobs-multiphase/fluidity/compressible-multiphase]

4131. By Christian Jacobs on 2013-02-20

Merging trunk revisions (from r4173 to r4177 inclusive) into branch.

=== modified file 'assemble/Advection_Diffusion_CG.F90'
--- assemble/Advection_Diffusion_CG.F90	2012-10-04 10:51:59 +0000
+++ assemble/Advection_Diffusion_CG.F90	2013-02-08 22:48:22 +0000
@@ -118,6 +118,8 @@
   logical :: move_mesh
   ! Include porosity?
   logical :: include_porosity
+  ! Is this material_phase compressible?
+  logical :: compressible = .false.
   ! Are we running a multiphase flow simulation?
   logical :: multiphase
 
@@ -151,8 +153,9 @@
     
     ! Note: the assembly of the heat transfer term is done here to avoid 
     ! passing in the whole state array to assemble_advection_diffusion_cg.
-    if(have_option("/multiphase_interaction/heat_transfer")) then
-       call add_heat_transfer(state, istate, t, matrix, rhs)
+    if(have_option("/multiphase_interaction/heat_transfer") .and. &
+       equation_type_index(trim(t%option_path)) == FIELD_EQUATION_INTERNALENERGY) then
+      call add_heat_transfer(state, istate, t, matrix, rhs)
     end if
     call profiler_toc(t, "assembly")
     
@@ -520,18 +523,22 @@
       olddensity=>extract_scalar_field(state, "Old"//trim(density_name))
       ewrite_minmax(olddensity)
       
-      if(.not.multiphase .or. have_option('/material_phase::'//trim(state%name)//&
-         &'/equation_of_state/compressible')) then         
-         call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", &
-                        density_theta)
+      if(have_option(trim(state%option_path)//'/equation_of_state/compressible')) then         
+         call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", density_theta)
+         compressible = .true.
+         
+         ! We always include the p*div(u) term if this is the compressible phase.
+         pressure=>extract_scalar_field(state, "Pressure")
+         ewrite_minmax(pressure)
       else
          ! Since the particle phase is always incompressible then its Density
          ! will not be prognostic. Just use a fixed theta value of 1.0.
          density_theta = 1.0
+         compressible = .false.
+         
+         ! Don't include the p*div(u) term if this is the incompressible particle phase.
+         pressure => dummydensity
       end if
-                      
-      pressure=>extract_scalar_field(state, "Pressure")
-      ewrite_minmax(pressure)
 
     case(FIELD_EQUATION_KEPSILON)
       ewrite(2,*) "Solving k-epsilon equation"
@@ -860,7 +867,7 @@
     if(have_absorption) call add_absorption_element_cg(ele, test_function, t, absorption, detwei, matrix_addto, rhs_addto)
     
     ! Diffusivity
-    if(have_diffusivity) call add_diffusivity_element_cg(ele, t, diffusivity, dt_t, detwei, matrix_addto, rhs_addto)
+    if(have_diffusivity) call add_diffusivity_element_cg(ele, t, diffusivity, dt_t, nvfrac, detwei, matrix_addto, rhs_addto)
     
     ! Source
     if(have_source .and. (.not. add_src_directly_to_rhs)) then 
@@ -868,9 +875,10 @@
     end if
     
     ! Pressure
-    if(equation_type==FIELD_EQUATION_INTERNALENERGY) call add_pressurediv_element_cg(ele, test_function, t, &
-                                                                                  velocity, pressure, nvfrac, &
-                                                                                  du_t, dnvfrac_t, detwei, rhs_addto)
+    if(equation_type==FIELD_EQUATION_INTERNALENERGY .and. compressible) then
+       call add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, detwei, rhs_addto)
+    end if
+                                                                                  
     
     ! Step 4: Insertion
             
@@ -1187,9 +1195,9 @@
     
   end subroutine add_absorption_element_cg
   
-  subroutine add_diffusivity_element_cg(ele, t, diffusivity, dt_t, detwei, matrix_addto, rhs_addto)
+  subroutine add_diffusivity_element_cg(ele, t, diffusivity, dt_t, nvfrac, detwei, matrix_addto, rhs_addto)
     integer, intent(in) :: ele
-    type(scalar_field), intent(in) :: t
+    type(scalar_field), intent(in) :: t, nvfrac
     type(tensor_field), intent(in) :: diffusivity
     real, dimension(ele_loc(t, ele), ele_ngi(t, ele), mesh_dim(t)), intent(in) :: dt_t
     real, dimension(ele_ngi(t, ele)), intent(in) :: detwei
@@ -1202,9 +1210,22 @@
     assert(have_diffusivity)
     
     diffusivity_gi = ele_val_at_quad(diffusivity, ele)
+
     if(isotropic_diffusivity) then
-      assert(size(diffusivity_gi, 1) > 0)
-      diffusivity_mat = dshape_dot_dshape(dt_t, dt_t, detwei * diffusivity_gi(1, 1, :))
+      assert(size(diffusivity_gi, 1) > 0)    
+          
+      if(multiphase .and. equation_type==FIELD_EQUATION_INTERNALENERGY) then
+         ! This allows us to use the Diffusivity term as the heat flux term
+         ! in the multiphase InternalEnergy equation: div( (k/Cv) * vfrac * grad(ie) ).
+         ! The user needs to input k/Cv for the prescribed diffusivity,
+         ! where k is the effective conductivity and Cv is the specific heat
+         ! at constant volume. We've assumed this will always be isotropic here.
+         ! The division by Cv is needed because the heat flux
+         ! is defined in terms of temperature T = ie/Cv.
+         diffusivity_mat = dshape_dot_dshape(dt_t, dt_t, detwei * diffusivity_gi(1, 1, :) * ele_val_at_quad(nvfrac, ele))
+      else
+         diffusivity_mat = dshape_dot_dshape(dt_t, dt_t, detwei * diffusivity_gi(1, 1, :))
+      end if
     else
       diffusivity_mat = dshape_tensor_dshape(dt_t, diffusivity_gi, dt_t, detwei)
     end if
@@ -1215,7 +1236,7 @@
     
   end subroutine add_diffusivity_element_cg
   
-  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, dnvfrac_t, detwei, rhs_addto)
+  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, detwei, rhs_addto)
   
     integer, intent(in) :: ele
     type(element_type), intent(in) :: test_function
@@ -1224,23 +1245,15 @@
     type(scalar_field), intent(in) :: pressure
     type(scalar_field), intent(in) :: nvfrac
     real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)), intent(in) :: du_t
-    real, dimension(:, :, :), intent(in) :: dnvfrac_t
     real, dimension(ele_ngi(t, ele)), intent(in) :: detwei
     real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto
     
-    real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad
-    real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad
-    
     assert(equation_type==FIELD_EQUATION_INTERNALENERGY)
     assert(ele_ngi(pressure, ele)==ele_ngi(t, ele))
     
     if(multiphase) then
-       ! -p * (div(vfrac*nu)) 
-       nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
-       udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*ele_val_at_quad(velocity, ele), 1)
-             
-       rhs_addto = rhs_addto - ( shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(nvfrac, ele) * ele_val_at_quad(pressure, ele) * detwei) &
-                                 + shape_rhs(test_function, udotgradnvfrac_at_quad * ele_val_at_quad(pressure, ele) * detwei) )
+       ! -p * vfrac * div(nu)
+       rhs_addto = rhs_addto - shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei * ele_val_at_quad(nvfrac, ele))
     else
        rhs_addto = rhs_addto - shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei)
     end if

=== modified file 'assemble/Compressible_Projection.F90'
--- assemble/Compressible_Projection.F90	2012-08-03 11:18:57 +0000
+++ assemble/Compressible_Projection.F90	2013-02-08 22:48:22 +0000
@@ -671,7 +671,7 @@
     integer:: i
 
     do i=0, option_count("/material_phase")-1
-      prognostic_pressure_path="/material_phase"//int2str(i)//"/scalar_field::Pressure/prognostic"
+      prognostic_pressure_path="/material_phase["//int2str(i)//"]/scalar_field::Pressure/prognostic"
       if (have_option(trim(prognostic_pressure_path)//"/spatial_discretisation/discontinuous_galerkin") &
         .and. have_option(trim(prognostic_pressure_path)//"/scheme/use_compressible_projection_method")) then
         FLExit("With a DG pressure you cannot have use_compressible_projection_method")

=== modified file 'assemble/Divergence_Matrix_CG.F90'
--- assemble/Divergence_Matrix_CG.F90	2012-08-08 19:51:41 +0000
+++ assemble/Divergence_Matrix_CG.F90	2013-02-08 22:48:22 +0000
@@ -629,7 +629,7 @@
                ! If the field and nvfrac meshes are different, then we need to
                ! compute the derivatives of the nvfrac shape functions.
                if(.not.(nvfrac%mesh == field%mesh)) then
-                  nvfrac_shape => ele_shape(nvfrac%mesh, ele)
+                  nvfrac_shape => ele_shape(nvfrac, ele)
                   call transform_to_physical(coordinate, ele, nvfrac_shape, dshape=dnvfrac_t)
                else
                   dnvfrac_t = dfield_t

=== modified file 'assemble/Divergence_Matrix_CV.F90'
--- assemble/Divergence_Matrix_CV.F90	2011-11-25 13:09:24 +0000
+++ assemble/Divergence_Matrix_CV.F90	2013-02-08 22:48:22 +0000
@@ -57,7 +57,7 @@
     !************************************************************************
     subroutine assemble_divergence_matrix_cv(CT_m, state, ct_rhs, & 
                                              test_mesh, field, &
-                                             get_ct, exclude_boundaries)
+                                             get_ct, exclude_boundaries, include_vfrac)
 
       ! inputs/outputs
       ! bucket full of fields
@@ -74,6 +74,7 @@
       
       logical, intent(in), optional :: get_ct
       logical, intent(in), optional :: exclude_boundaries
+      logical, intent(in), optional :: include_vfrac
 
       ! local
       ! degree of quadrature over cv faces
@@ -109,7 +110,7 @@
       real, dimension(:,:,:), allocatable :: ct_mat_local, ct_mat_local_bdy
       real, dimension(:,:), allocatable :: ct_rhs_local
 
-      logical :: l_get_ct
+      logical :: l_get_ct, l_include_vfrac
 
       !! Multiphase variables
       logical :: multiphase
@@ -138,6 +139,17 @@
       else
         l_get_ct = .true.
       end if
+      
+      ! In some cases we might not want to include the PhaseVolumeFraction
+      ! field in the divergence matrix, even though we are running a multiphase
+      ! flow simulation.
+      ! For example, in the InternalEnergy equation we just want
+      ! div(u), not div(vfrac*u).
+      if(present(include_vfrac)) then
+        l_include_vfrac = include_vfrac
+      else
+        l_include_vfrac = .true.
+      end if
 
       x=>extract_vector_field(state, "Coordinate")
       allocate(x_ele(x%dim,x%mesh%shape%loc))
@@ -154,7 +166,7 @@
                             quaddegree=quaddegree)
 
       ! Check if we need to multiply through by the non-linear volume fraction
-      if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+      if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1 .and. l_include_vfrac) then
          multiphase = .true.
 
          vfrac => extract_scalar_field(state, "PhaseVolumeFraction")

=== modified file 'assemble/Field_Equations_CV.F90'
--- assemble/Field_Equations_CV.F90	2012-10-04 10:51:59 +0000
+++ assemble/Field_Equations_CV.F90	2013-02-08 22:48:22 +0000
@@ -56,6 +56,7 @@
   use field_options
   use state_fields_module
   use porous_media
+  use multiphase_module
 
   implicit none
 
@@ -93,11 +94,13 @@
   logical :: assemble_advection_matrix = .true.
   ! diffusion?
   logical :: assemble_diffusion = .false.
+  ! Are we running a multiphase flow simulation?
+  logical :: multiphase = .false.
 
 contains
     !************************************************************************
     ! solution wrapping subroutines
-    subroutine solve_field_eqn_cv(field_name, state, global_it)
+    subroutine solve_field_eqn_cv(field_name, state, istate, global_it)
       !!< Construct and solve the advection-diffusion equation for the given
       !!< field using control volumes.
 
@@ -105,6 +108,7 @@
       character(len=*), intent(in) :: field_name
       !! Collection of fields defining system state.
       type(state_type), dimension(:), intent(inout) :: state
+      integer, intent(in) :: istate ! The material_phase index in which the field called field_name is in.
       ! global iteration number - passed in so it can be output to file
       integer, intent(in) :: global_it
 
@@ -199,6 +203,9 @@
       !! Temporary pointer to the material_phase's Velocity field
       type(vector_field), pointer :: temp_velocity_ptr
       character(len=FIELD_NAME_LEN) :: velocity_equation_type
+      ! Volume fraction fields for multiphase flow simulation
+      type(scalar_field), pointer :: vfrac, oldvfrac
+      type(scalar_field) :: nvfrac ! Non-linear version
 
       ! assume explicitness?
       logical :: explicit
@@ -210,20 +217,20 @@
       character(len=OPTION_PATH_LEN), dimension(1) :: tdensity_option_path_array
 
       ewrite(2,*) 'in solve_field_eqn_cv'
-      ewrite(2,*) 'solving for '//field_name//' in state '//trim(state(1)%name)
+      ewrite(2,*) 'solving for '//field_name//' in state '//trim(state(istate)%name)
 
       ! extract lots of fields:
       ! the actual thing we're trying to solve for
-      tfield=>extract_scalar_field(state(1), trim(field_name))
+      tfield=>extract_scalar_field(state(istate), trim(field_name))
       ewrite_minmax(tfield)
       option_path=tfield%option_path
       ! its previous timelevel - if we're doing more than one iteration per timestep!
-      oldtfield=>extract_scalar_field(state(1), "Old"//trim(field_name))
+      oldtfield=>extract_scalar_field(state(istate), "Old"//trim(field_name))
       ! because fluidity resets tfield to oldtfield at the start of every
       ! global iteration we need to undo this so that the control volume faces
       ! are discretised using the most up to date values
       ! therefore extract the iterated values:
-      it_tfield=>extract_scalar_field(state(1), "Iterated"//trim(field_name))
+      it_tfield=>extract_scalar_field(state(istate), "Iterated"//trim(field_name))
       ! and set tfield to them:
       call set(tfield, it_tfield)
       ewrite_minmax(tfield)
@@ -253,6 +260,21 @@
       call allocate(dummydensity, tfield%mesh, name="DummyDensity", field_type=FIELD_TYPE_CONSTANT)
       call set(dummydensity, 1.0)
       dummydensity%option_path = " "
+      
+      ! PhaseVolumeFraction for multiphase flow simulations
+      if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+         multiphase = .true.
+         vfrac => extract_scalar_field(state(istate), "PhaseVolumeFraction")
+         call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
+         call zero(nvfrac)
+         call get_nonlinear_volume_fraction(state(istate), nvfrac)
+         
+         ewrite_minmax(nvfrac)
+         
+         oldvfrac => extract_scalar_field(state(istate), "OldPhaseVolumeFraction")
+      else
+         multiphase = .false.
+      end if
 
       ! find out equation type and hence if density is needed or not
       equation_type=equation_type_index(trim(option_path))
@@ -266,16 +288,18 @@
 
       case(FIELD_EQUATION_KEPSILON)
         ! Depending on the equation type, extract the density or set it to some dummy field allocated above
-        temp_velocity_ptr => extract_vector_field(state(1), "Velocity")
+        temp_velocity_ptr => extract_vector_field(state(istate), "Velocity")
         call get_option(trim(temp_velocity_ptr%option_path)//"/prognostic/equation[0]/name", velocity_equation_type)
         select case(velocity_equation_type)
            case("LinearMomentum")
               include_density = .true.
-              tdensity=>extract_scalar_field(state(1), "Density")
-              oldtdensity=>extract_scalar_field(state(1), "OldDensity")
+              tdensity=>extract_scalar_field(state(istate), "Density")
+              oldtdensity=>extract_scalar_field(state(istate), "OldDensity")
 
               if(have_option(trim(tdensity%option_path)//"/prognostic")) then
                  prognostic_density = .true.
+              else
+                 prognostic_density = .false.
               end if
            case("Boussinesq")
               tdensity => dummydensity
@@ -296,16 +320,18 @@
         ! density needed so extract the type specified in the input
         ! ?? are there circumstances where this should be "Iterated"... need to be
         ! careful with priority ordering
-        tdensity=>extract_scalar_field(state(1), trim(tmpstring))
+        tdensity=>extract_scalar_field(state(istate), trim(tmpstring))
         ewrite_minmax(tdensity)
 
         ! halo exchange? - not currently necessary when suboptimal halo exchange if density
         ! is solved for with this subroutine and the correct priority ordering.
-        oldtdensity=>extract_scalar_field(state(1), "Old"//trim(tmpstring))
+        oldtdensity=>extract_scalar_field(state(istate), "Old"//trim(tmpstring))
         ewrite_minmax(oldtdensity)
            
         if(have_option(trim(tdensity%option_path)//"/prognostic")) then
            prognostic_density = .true.
+        else
+           prognostic_density = .false.
         end if
       end select
 
@@ -326,14 +352,18 @@
       ! handily wrapped in a new type...
       tfield_options=get_cv_options(tfield%option_path, tfield%mesh%shape%numbering%family, mesh_dim(tfield))
       if(include_density) then
-        tdensity_options=get_cv_options(tdensity_option_path, tdensity%mesh%shape%numbering%family, mesh_dim(tdensity),  coefficient_field=.true.)
+        if(prognostic_density) then
+          tdensity_options=get_cv_options(tdensity_option_path, tdensity%mesh%shape%numbering%family, mesh_dim(tdensity),  coefficient_field=.true.)
+        else
+          tdensity_options=tfield_options
+        end if
       end if
 
       ! extract fields from state
       ! the base Coordinate field
-      x=>extract_vector_field(state(1), "Coordinate")
+      x=>extract_vector_field(state(istate), "Coordinate")
       ! the Coordinate field on the same mesh as tfield
-      x_tfield=get_coordinate_field(state(1), tfield%mesh)
+      x_tfield=get_coordinate_field(state(istate), tfield%mesh)
 
       ! are we including the mass (generally yes)?
       include_mass = .not.have_option(trim(tfield%option_path)//"/prognostic/spatial_discretisation/control_volumes/mass_terms/exclude_mass_terms")
@@ -341,7 +371,7 @@
       ! are we inclulding advection (generally yes)?
       include_advection = .not.(tfield_options%facevalue==CV_FACEVALUE_NONE)
       if(include_advection) then
-        nu=>extract_vector_field(state(1), "NonlinearVelocity")
+        nu=>extract_vector_field(state(istate), "NonlinearVelocity")
         ewrite_minmax(nu)
         ! find relative velocity
         allocate(advu)
@@ -354,10 +384,10 @@
             ewrite(2,*) "Removing velocity from ", trim(field_name)
         end if
         ! add in sinking velocity
-        sink=>extract_scalar_field(state(1), trim(field_name)//"SinkingVelocity"&
+        sink=>extract_scalar_field(state(istate), trim(field_name)//"SinkingVelocity"&
             &, stat=stat)
         if(stat==0) then
-          gravity=>extract_vector_field(state(1), "GravityDirection")
+          gravity=>extract_vector_field(state(istate), "GravityDirection")
           ! this may perform a "remap" internally from CoordinateMesh to VelocityMesh
           call addto(advu, gravity, scale=sink)
         end if
@@ -365,7 +395,7 @@
       else
         ewrite(2,*) 'Excluding advection'
         advu => dummyvector
-        if(has_scalar_field(state(1), trim(field_name)//"SinkingVelocity")) then
+        if(has_scalar_field(state(istate), trim(field_name)//"SinkingVelocity")) then
             ewrite(-1,*) "No advection in "//trim(field_name)
             FLExit("But you have a sinking velocity on. Can't have that")
         end if
@@ -374,7 +404,7 @@
       ! do we have a diffusivity - this will control whether we construct an auxiliary
       ! eqn or not (if BassiRebay is selected) or whether we construct gradients (if ElementGradient
       ! is selected)
-      diffusivity=>extract_tensor_field(state(1), trim(field_name)//"Diffusivity", stat=stat)
+      diffusivity=>extract_tensor_field(state(istate), trim(field_name)//"Diffusivity", stat=stat)
       include_diffusion = (stat==0).and.(tfield_options%diffusionscheme/=CV_DIFFUSION_NONE)
       if(.not.include_diffusion) then
         diffusivity => dummytensor
@@ -383,7 +413,7 @@
       end if
       
       ! do we have a source?
-      source=>extract_scalar_field(state(1), trim(field_name)//"Source", stat=stat)
+      source=>extract_scalar_field(state(istate), trim(field_name)//"Source", stat=stat)
       include_source = (stat==0)
       if(.not.include_source) then
         source=>dummyscalar
@@ -399,7 +429,7 @@
       end if
       
       ! do we have an absorption?
-      absorption=>extract_scalar_field(state(1), trim(field_name)//"Absorption", stat=stat)
+      absorption=>extract_scalar_field(state(istate), trim(field_name)//"Absorption", stat=stat)
       include_absorption = (stat==0)
       if(.not.include_absorption) then
         absorption=>dummyscalar
@@ -471,7 +501,7 @@
       option_path_array(1) = trim(option_path)                  ! temporary hack for compiler failure
       tdensity_option_path_array(1) = tdensity_option_path
       call cv_disc_get_cfl_no(option_path_array, &
-                      state(1), tfield%mesh, cfl_no, &
+                      state(istate), tfield%mesh, cfl_no, &
                       tdensity_option_path_array)
 
       ! get the mesh sparsity for the matrices
@@ -553,7 +583,7 @@
          include_porosity = .true.
 
          ! get the porosity theta averaged field - this will allocate it
-         call form_porosity_theta(porosity_theta, state(1), &
+         call form_porosity_theta(porosity_theta, state(istate), &
              &option_path = trim(complete_field_path(tfield%option_path))//'/porosity')       
                   
          call allocate(t_cvmass_with_porosity, tfield%mesh, name="CVMassWithPorosity")
@@ -586,8 +616,8 @@
            FLExit("Moving mesh not set up to work when including porosity")
         end if
         ewrite(2,*) "Moving mesh."
-        x_old=>extract_vector_field(state(1), "OldCoordinate")
-        x_new=>extract_vector_field(state(1), "IteratedCoordinate")
+        x_old=>extract_vector_field(state(istate), "OldCoordinate")
+        x_new=>extract_vector_field(state(istate), "IteratedCoordinate")
         call allocate(t_cvmass_old, tfield%mesh, name=trim(field_name)//"OldCVMass")
         call allocate(t_cvmass_new, tfield%mesh, name=trim(field_name)//"NewCVMass")
 
@@ -596,7 +626,7 @@
         ewrite_minmax(t_cvmass_old)
         ewrite_minmax(t_cvmass_new)
         
-        ug=>extract_vector_field(state(1), "GridVelocity")
+        ug=>extract_vector_field(state(istate), "GridVelocity")
         ewrite_minmax(ug)
 
         ug_cvshape=make_cv_element_shape(cvfaces, ug%mesh%shape)
@@ -702,14 +732,14 @@
             ! diffusion matrix then assemble A_m, D_m and rhs here
             call assemble_advectiondiffusion_m_cv(A_m, rhs, D_m, diff_rhs, &
                                         tfield, l_old_tfield, tfield_options, &
-                                        tdensity, oldtdensity, tdensity_options, &
-                                        cvfaces, x_cvshape, x_cvbdyshape, &
+                                        tdensity, oldtdensity, nvfrac, oldvfrac, &
+                                        tdensity_options, cvfaces, x_cvshape, x_cvbdyshape, &
                                         u_cvshape, u_cvbdyshape, t_cvshape, &
                                         ug_cvshape, ug_cvbdyshape, &
                                         x_cvshape_full, x_cvbdyshape_full, &
                                         t_cvshape_full, t_cvbdyshape_full, &
                                         diff_cvshape_full, diff_cvbdyshape_full, &
-                                        state, advu, ug, x, x_tfield, cfl_no, sub_dt, &
+                                        state(istate:istate), advu, ug, x, x_tfield, cfl_no, sub_dt, &
                                         diffusivity, q_cvmass, &
                                         mesh_sparsity_x, grad_m_t_sparsity)
           end if
@@ -719,10 +749,14 @@
                                     tfield, l_old_tfield, &
                                     tdensity, oldtdensity, tdensity_options, &
                                     source, absorption, tfield_options%theta, &
-                                    state, advu, sub_dt, explicit, &
+                                    state(istate:istate), advu, sub_dt, explicit, &
                                     t_cvmass, t_abs_src_cvmass, t_cvmass_old, t_cvmass_new, & 
-                                    D_m, diff_rhs)
+                                    D_m, diff_rhs, nvfrac, oldvfrac)
 
+          if(have_option("/multiphase_interaction/heat_transfer") .and. &
+             equation_type_index(trim(tfield%option_path)) == FIELD_EQUATION_INTERNALENERGY) then
+            call add_heat_transfer(state, istate, tfield, M, rhs)
+          end if
 
           ! Solve for the change in tfield.
           if(explicit) then
@@ -736,7 +770,7 @@
             call apply_dirichlet_conditions(M, rhs, tfield, sub_dt)
 
             call zero(delta_tfield)
-            call petsc_solve(delta_tfield, M, rhs, state(1))
+            call petsc_solve(delta_tfield, M, rhs, state(istate))
           end if
 
           ! reset tfield to l_old_tfield before applying change
@@ -747,7 +781,7 @@
           call halo_update(tfield)  ! exchange the extended halos
 
           call test_and_write_advection_convergence(tfield, advit_tfield, x, t_cvmass, &
-                                    filename=trim(state(1)%name)//"__"//trim(tfield%name), &
+                                    filename=trim(state(istate)%name)//"__"//trim(tfield%name), &
                                     time=time+sub_dt, dt=sub_dt, it=global_it, adv_it=adv_it, &
                                     subcyc=sub, error=error)
 
@@ -806,6 +840,9 @@
       if (include_porosity) then
         call deallocate(t_cvmass_with_porosity)
       end if
+      if(multiphase) then
+         call deallocate(nvfrac)
+      end if
 
     end subroutine solve_field_eqn_cv
     ! end of solution wrapping subroutines
@@ -818,7 +855,7 @@
                                     source, absorption, theta, &
                                     state, advu, dt, explicit, &
                                     cvmass, abs_src_cvmass, cvmass_old, cvmass_new, &
-                                    D_m, diff_rhs)
+                                    D_m, diff_rhs, nvfrac, oldvfrac)
 
       ! This subroutine assembles the equation
       ! M(T^{n+1}-T^{n})/\Delta T = rhs
@@ -862,6 +899,8 @@
       ! diffusion:
       type(csr_matrix), intent(inout), optional :: D_m
       type(scalar_field), intent(inout), optional :: diff_rhs
+      
+      type(scalar_field), intent(inout), optional :: nvfrac, oldvfrac
 
       ! local memory:
       ! for all equation types:
@@ -1153,38 +1192,50 @@
         ! [\rho^{n+1}M + dt*A_m + dt*theta*D_m](T^{n+1}-T^{n})/dt = rhs - [A_m + D_m]*T^{n} - diff_rhs - (p+atm_p)*CT_m*u
 
         ! construct rhs
-        p=>extract_scalar_field(state(1), "Pressure")
-        ewrite_minmax(p)
-        assert(p%mesh==tfield%mesh)
-        ! halo exchange not necessary as it is done straight after solve
-        call get_option(trim(p%option_path)//'/prognostic/atmospheric_pressure', &
-                              atmospheric_pressure, default=0.0)
-        gradient_sparsity => get_csr_sparsity_firstorder(state, p%mesh, advu%mesh)
-
-        call allocate(CT_m, gradient_sparsity, (/1, advu%dim/), name="DivergenceMatrix" )
-        call assemble_divergence_matrix_cv(CT_m, state(1), &
-                                           test_mesh=p%mesh, field=advu)
-
-        call allocate(pterm, p%mesh, "PressureTerm")
-
-        ! construct the pressure term
-        call mult(pterm, CT_m, advu) 
-                                ! should this really be the advection velocity or just the relative or the nonlinear?
-        pterm%val = pterm%val*(p%val+atmospheric_pressure)
-
-        call addto(rhs, pterm, -1.0)
-
-        call deallocate(CT_m)
-        call deallocate(pterm)
-
+        if(have_option(trim(state(1)%option_path)//'/equation_of_state/compressible')) then
+           p=>extract_scalar_field(state(1), "Pressure")
+           ewrite_minmax(p)
+           assert(p%mesh==tfield%mesh)
+           ! halo exchange not necessary as it is done straight after solve
+           call get_option(trim(p%option_path)//'/prognostic/atmospheric_pressure', &
+                                   atmospheric_pressure, default=0.0)
+           gradient_sparsity => get_csr_sparsity_firstorder(state, p%mesh, advu%mesh)
+
+           call allocate(CT_m, gradient_sparsity, (/1, advu%dim/), name="DivergenceMatrix" )
+           call assemble_divergence_matrix_cv(CT_m, state(1), &
+                                                test_mesh=p%mesh, field=advu, include_vfrac=.false.)
+                                                
+           call allocate(pterm, p%mesh, "PressureTerm")
+            
+           ! construct the pressure term
+           call mult(pterm, CT_m, advu) 
+           if(multiphase) then
+              call scale(pterm, nvfrac) ! We need vfrac*p*div(u) for the multiphase InternalEnergy equation
+           end if
+            
+                                   ! should this really be the advection velocity or just the relative or the nonlinear?
+           pterm%val = pterm%val*(p%val+atmospheric_pressure)
+
+           call addto(rhs, pterm, -1.0)
+
+           call deallocate(CT_m)
+           call deallocate(pterm)
+        end if
+         
         ! construct M
         if(explicit) then
           if(include_mass) then
             call scale(m_cvmass, tdensity)
+            if(multiphase) then
+               call scale(m_cvmass, nvfrac) 
+            end if
           end if
         else
           if(include_mass) then
             call mult_diag(M, tdensity)
+            if(multiphase) then
+               call mult_diag(M, nvfrac) 
+            end if
           end if
           if(include_advection) call addto(M, A_m, dt)
           if(include_absorption) call addto_diag(M, massabsorption, theta*dt)
@@ -1283,8 +1334,8 @@
     ! assembly subroutines 
     subroutine assemble_advectiondiffusion_m_cv(A_m, rhs, D_m, diff_rhs, &
                                        tfield, oldtfield, tfield_options, &
-                                       tdensity, oldtdensity, tdensity_options, &
-                                       cvfaces, x_cvshape, x_cvbdyshape, &
+                                       tdensity, oldtdensity, nvfrac, oldvfrac, &
+                                       tdensity_options, cvfaces, x_cvshape, x_cvbdyshape, &
                                        u_cvshape, u_cvbdyshape, t_cvshape, &
                                        ug_cvshape, ug_cvbdyshape, &
                                        x_cvshape_full, x_cvbdyshape_full, &
@@ -1346,6 +1397,9 @@
       type(scalar_field), intent(in) :: cfl_no
       ! timestep
       real, intent(in) :: dt
+      
+      type(scalar_field), intent(in) :: nvfrac, oldvfrac
+      real, dimension(x_cvshape%ngi) :: nvfrac_gi
 
       ! mesh sparsity for upwind value matrices
       type(csr_sparsity), intent(in) :: mesh_sparsity
@@ -1397,6 +1451,9 @@
 
       ! loop integers
       integer :: ele, sele, iloc, oloc, dloc, face, gi, ggi, dim
+      
+      ! what equation type are we solving for?
+      integer :: equation_type
 
       ! upwind value matrices for the fields and densities
       type(csr_matrix)  :: tfield_upwind, &
@@ -1551,6 +1608,12 @@
           call transform_to_physical(X, ele, x_shape=x_cvshape_full, &
                                     shape=t_cvshape_full, dshape=dt_t)
           diffusivity_gi = ele_val_at_quad(diffusivity, ele, diff_cvshape_full)
+          
+          equation_type = equation_type_index(trim(tfield%option_path))
+          if(multiphase .and. equation_type==FIELD_EQUATION_INTERNALENERGY) then
+            nvfrac_gi = ele_val_at_quad(nvfrac, ele, diff_cvshape_full)
+          end if
+          
         end if
 
         cfl_ele = ele_val(cfl_no, ele)
@@ -1561,6 +1624,11 @@
         if(include_density) then
           tdensity_ele = ele_val(tdensity, ele)
           oldtdensity_ele = ele_val(oldtdensity, ele)
+          
+          if(multiphase) then
+            tdensity_ele = tdensity_ele*ele_val(nvfrac,ele)
+            oldtdensity_ele = oldtdensity_ele*ele_val(oldvfrac,ele)
+          end if
         end if
 
         notvisited=.true.
@@ -1634,28 +1702,25 @@
                       ! do the same for the density but save some effort if it's just a dummy
                       select case (tdensity%field_type)
                       case(FIELD_TYPE_CONSTANT)
-
                           tdensity_face_val = tdensity_ele(iloc)
                           oldtdensity_face_val = oldtdensity_ele(iloc)
 
                       case default
-
                           call evaluate_face_val(tdensity_face_val, oldtdensity_face_val, &
-                                                iloc, oloc, ggi, upwind_nodes, &
-                                                t_cvshape,&
-                                                tdensity_ele, oldtdensity_ele, &
-                                                tdensity_upwind, oldtdensity_upwind, &
-                                                inflow, cfl_ele, &
-                                                tdensity_options)
+                                                   iloc, oloc, ggi, upwind_nodes, &
+                                                   t_cvshape,&
+                                                   tdensity_ele, oldtdensity_ele, &
+                                                   tdensity_upwind, oldtdensity_upwind, &
+                                                   inflow, cfl_ele, &
+                                                   tdensity_options)
 
                       end select
-
                       tdensity_theta_val=theta_val(iloc, oloc, &
-                                          tdensity_face_val, &
-                                          oldtdensity_face_val, &
-                                          tdensity_options%theta, dt, udotn, &
-                                          xt_ele, tdensity_options%limit_theta, &
-                                          tdensity_ele, oldtdensity_ele)
+                                             tdensity_face_val, &
+                                             oldtdensity_face_val, &
+                                             tdensity_options%theta, dt, udotn, &
+                                             xt_ele, tdensity_options%limit_theta, &
+                                             tdensity_ele, oldtdensity_ele)
 
                       if(assemble_advection_matrix) then
                         mat_local(iloc, oloc) = mat_local(iloc, oloc) &
@@ -1731,16 +1796,35 @@
                       end do dimension_loop1
 
                     case(CV_DIFFUSION_ELEMENTGRADIENT)
-
-                      do dloc=1,size(dt_t,1)
-                        ! n_i K_{ij} dT/dx_j
-                        diff_mat_local(iloc,dloc) = diff_mat_local(iloc,dloc) - &
-                          sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*normgi, 1)*detwei(ggi)
-
-                        ! notvisited
-                        diff_mat_local(oloc, dloc) = diff_mat_local(oloc,dloc) - &
-                          sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*(-normgi), 1)*detwei(ggi)
-                      end do
+                    
+                     if(multiphase .and. equation_type==FIELD_EQUATION_INTERNALENERGY) then
+                        ! This allows us to use the Diffusivity term as the heat flux term
+                        ! in the multiphase InternalEnergy equation: div( (k/Cv) * vfrac * grad(ie) ).
+                        ! The user needs to input k/Cv for the prescribed diffusivity,
+                        ! where k is the effective conductivity and Cv is the specific heat
+                        ! at constant volume. The division by Cv is needed because the heat flux
+                        ! is defined in terms of temperature T = ie/Cv.
+                     
+                        do dloc=1,size(dt_t,1)
+                           ! n_i K_{ij} dT/dx_j
+                           diff_mat_local(iloc,dloc) = diff_mat_local(iloc,dloc) - &
+                           sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*normgi, 1)*detwei(ggi)*nvfrac_gi(ggi)
+
+                           ! notvisited
+                           diff_mat_local(oloc, dloc) = diff_mat_local(oloc,dloc) - &
+                           sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*(-normgi), 1)*detwei(ggi)*nvfrac_gi(ggi)
+                        end do
+                     else       
+                        do dloc=1,size(dt_t,1)
+                           ! n_i K_{ij} dT/dx_j
+                           diff_mat_local(iloc,dloc) = diff_mat_local(iloc,dloc) - &
+                           sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*normgi, 1)*detwei(ggi)
+
+                           ! notvisited
+                           diff_mat_local(oloc, dloc) = diff_mat_local(oloc,dloc) - &
+                           sum(matmul(diffusivity_gi(:,:,ggi), dt_t(dloc, ggi, :))*(-normgi), 1)*detwei(ggi)
+                        end do
+                     end if
 
                     end select
                   end if
@@ -1854,17 +1938,30 @@
             if(tdensity_bc_type(sele)==BC_TYPE_WEAKDIRICHLET) then
               ghost_tdensity_ele_bdy=ele_val(tdensity_bc, sele)
             else
-              ghost_tdensity_ele_bdy=face_val(tdensity, sele)
+              if(multiphase) then
+                 ghost_tdensity_ele_bdy=face_val(tdensity, sele)*face_val(nvfrac, sele)
+              else
+                 ghost_tdensity_ele_bdy=face_val(tdensity, sele)
+              end if
             end if
 
             if(tdensity_bc_type(sele)==BC_TYPE_WEAKDIRICHLET) then
               ghost_oldtdensity_ele_bdy=ele_val(tdensity_bc, sele) ! not considering time varying bcs yet
             else
-              ghost_oldtdensity_ele_bdy=face_val(oldtdensity, sele)
+              if(multiphase) then
+                 ghost_oldtdensity_ele_bdy=face_val(oldtdensity, sele)*face_val(oldvfrac, sele)
+              else
+                 ghost_oldtdensity_ele_bdy=face_val(oldtdensity, sele)
+              end if
             end if
 
             tdensity_ele_bdy=face_val(tdensity, sele)
             oldtdensity_ele_bdy=face_val(oldtdensity, sele)
+            
+            if(multiphase) then
+               tdensity_ele_bdy=tdensity_ele_bdy*face_val(nvfrac, sele)
+               oldtdensity_ele_bdy=oldtdensity_ele_bdy*face_val(oldvfrac, sele)
+            end if
           end if
         end if
 

=== modified file 'assemble/Makefile.dependencies'
--- assemble/Makefile.dependencies	2013-01-26 18:26:38 +0000
+++ assemble/Makefile.dependencies	2013-02-08 22:48:22 +0000
@@ -281,8 +281,9 @@
    ../include/diagnostic_variables.mod ../include/divergence_matrix_cv.mod \
    ../include/fdebug.h ../include/fefields.mod ../include/field_options.mod \
    ../include/fields.mod ../include/global_parameters.mod ../include/halos.mod \
-   ../include/parallel_tools.mod ../include/petsc_solve_state_module.mod \
-   ../include/porous_media.mod ../include/sparse_matrices_fields.mod \
+   ../include/multiphase_module.mod ../include/parallel_tools.mod \
+   ../include/petsc_solve_state_module.mod ../include/porous_media.mod \
+   ../include/sparse_matrices_fields.mod \
    ../include/sparsity_patterns_meshes.mod ../include/spud.mod \
    ../include/state_fields_module.mod ../include/state_module.mod \
    ../include/transform_elements.mod

=== modified file 'assemble/Momentum_Equation.F90'
--- assemble/Momentum_Equation.F90	2012-12-12 21:42:17 +0000
+++ assemble/Momentum_Equation.F90	2013-02-08 22:48:22 +0000
@@ -183,7 +183,7 @@
          ! Scaled pressure mass matrix - used for preconditioning full projection solve:
          type(csr_matrix), target :: scaled_pressure_mass_matrix
          type(csr_sparsity), pointer :: scaled_pressure_mass_matrix_sparsity
-         ! Left hand matrix of CMC. For incompressibe flow this points to ct_m as they are identical, 
+         ! Left hand matrix of CMC. For incompressible flow this points to ct_m as they are identical, 
          ! unless for CG pressure with CV tested continuity case when this matrix will be the 
          ! CV divergence tested matrix and ct_m the CG divergence tested matrix (right hand matrix of CMC).
          ! For compressible flow this differs to ct_m in that it will contain the variable density.
@@ -2196,6 +2196,17 @@
                
             end if
             
+            ! Check that each particle phase has a PROGNOSTIC PhaseVolumeFraction field.
+            ! The fluid phase cannot have a prognostic PhaseVolumeFraction as this is not always valid.
+            ! For example, since we do not include the Density in the advection-diffusion equation for the PhaseVolumeFraction,
+            ! solving this equation for the compressible fluid phase would not be correct. The particle phases on the other hand
+            ! are always incompressible where the density is constant.
+            if(have_option("/material_phase["//int2str(i)//"]/multiphase_properties/particle_diameter") .and. &
+               .not.(have_option("/material_phase["//int2str(i)//"]/scalar_field::PhaseVolumeFraction/prognostic") .or. &
+               have_option("/material_phase["//int2str(i)//"]/scalar_field::PhaseVolumeFraction/prescribed"))) then
+               FLExit("All particle phases must have a prognostic/prescribed PhaseVolumeFraction field. The diagnostic PhaseVolumeFraction field should always be in the continuous/fluid phase.")
+            end if
+            
          end do
 
          ewrite(1,*) 'Finished checking momentum discretisation options'

=== modified file 'main/Fluids.F90'
--- main/Fluids.F90	2012-10-02 11:35:46 +0000
+++ main/Fluids.F90	2013-02-08 22:48:22 +0000
@@ -697,8 +697,7 @@
 
                    ! Solve the pure control volume form of the equations
                    call solve_field_eqn_cv(field_name=trim(field_name_list(it)), &
-                        state=state(field_state_list(it):field_state_list(it)), &
-                        global_it=its)
+                        state=state, istate=field_state_list(it), global_it=its)
 
                 else if(have_option(trim(field_optionpath_list(it)) // &
                      & "/prognostic/spatial_discretisation/continuous_galerkin")) then

=== modified file 'manual/bibliography.bib'
--- manual/bibliography.bib	2012-12-04 17:19:40 +0000
+++ manual/bibliography.bib	2013-02-08 22:48:22 +0000
@@ -1672,6 +1672,17 @@
   year		= {2003}
 }
 
+@Article{         Jacobs_etal_2013,
+   author       = {Jacobs, C. T. and Collins, G. S. and Piggott, M. D. and Kramer, S. C. and Wilson, C. R. G.},
+   title        = {Multiphase flow modelling of volcanic ash particle settling in water using adaptive unstructured meshes},
+   journal      = {Geophysical Journal International},
+   year         = {2013},
+   volume       = {192},
+   number       = {2},
+   pages        = {647--665},
+   doi          = {10.1093/gji/ggs059}
+}
+
 @Article{	  jarrin_06,
   author	= {N. Jarrin and S. Benhamadouche and D. Laurence and R. Prosser},
   journal	= {International Journal of Heat and Fluid Flow},

=== modified file 'manual/configuring_fluidity.tex'
--- manual/configuring_fluidity.tex	2012-12-04 17:19:40 +0000
+++ manual/configuring_fluidity.tex	2013-02-08 22:48:22 +0000
@@ -2369,6 +2369,7 @@
 \begin{equation}
 \alpha_N = 1 - \sum_{i=1}^{N-1}\alpha_i\textrm{.}
 \end{equation}
+For compressible multiphase simulations, the diagnostic PhaseVolumeFraction field should always be in the compressible (fluid) phase. This is because we do not include the Density in the advection-diffusion equation for prognostic PhaseVolumeFraction fields, so solving this equation in the compressible phase would not be correct. The particle phases on the other hand are always incompressible where the density is constant.
 
 \subsubsection{Inter-phase momentum interaction}
 Currently, \fluidity\ only supports fluid-particle drag between the continuous phase and dispersed phase(s). This can be enabled by switching on the \option{/multiphase\_interaction/fluid\_particle\_drag} option in Diamond and specifying the drag correlation: Stokes, \cite{wen_yu_1966} or \cite{ergun1952}.
@@ -2407,20 +2408,20 @@
 This subsection considers compressible multiphase flow simulations where each phase has its own InternalEnergy field. Users can choose to include the inter-phase heat transfer term by \cite{gunn1978}, denoted $Q$ in Section \ref{sec:multiphase_equations}, on the RHS of each internal energy equation:
 
 \begin{equation}\label{eq:gunn_heat_transfer_term}
-Q = \frac{6 k \alpha_p \mathrm{Nu}_p}{d_p^2}\left(\frac{e_f}{C_f} - \frac{e_p}{C_p}\right),
+Q = \frac{6 k \alpha_p \mathrm{Nu}_p}{d_p^2}\left(\frac{e_f}{C_{v,f}} - \frac{e_p}{C_{v,p}}\right),
 \end{equation}
 where 
 \begin{equation}
 \mathrm{Nu}_p = \left(7 - 10\alpha_f + 5\alpha_f^2\right)\left(1 + 0.7\mathrm{Re}_p^{0.2}\mathrm{Pr}^{\frac{1}{3}}\right) + \left(1.33 - 2.4\alpha_f + 1.2\alpha_f^2\right)\mathrm{Re}_p^{0.7}\mathrm{Pr}^{\frac{1}{3}},
 \end{equation}
 \begin{equation}
- \mathrm{Pr} = \frac{C_f \gamma \mu_f}{k},
+ \mathrm{Pr} = \frac{C_{v,f} \gamma \mu_f}{k},
 \end{equation}
 and
 \begin{equation}
  \mathrm{Re} = \frac{\rho_f d_p |\mathbf{u}_f-\mathbf{u}_p|}{\mu_f},
 \end{equation}
-are the Nusselt, Prandtl and Reynolds number respectively. $C_i$ denotes the specific heat of phase i at constant volume, and $k$ denotes the effective conductivity of the fluid phase; these must be specified in \option{../multiphase\_properties/effective\_conductivity} and \option{../multiphase\_properties/specific\_heat}. Note that we have written the above in terms of internal energy (rather than temperature) by dividing $e_i$ by the specific heat at constant volume.
+are the Nusselt, Prandtl and Reynolds number respectively. $C_{v,i}$ denotes the specific heat of phase i at constant volume, and $k$ denotes the effective conductivity of the fluid phase; these must be specified in \option{../multiphase\_properties/effective\_conductivity} and \option{../multiphase\_properties/specific\_heat}. Note that we have written the above in terms of internal energy (rather than temperature) by dividing $e_i$ by the specific heat at constant volume.
 
 To include this term, switch on the \option{/multiphase\_interaction/heat\_transfer} option in Diamond.
 
@@ -2463,6 +2464,11 @@
 \item in the presence of viscosity \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/stress\_form} is required and all components of the anisotropic symmetric Viscosity tensor should be filled out.  This functionality is only available with continuous Galerkin velocities.
 \end{itemize}
 
+\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.
+\end{itemize}
+
 \subsection{Restrictions: discretisation options and element pairs}
 Either continuous Galerkin or control volumes can be used as discretisation options for pressure and density (both fields need to have the same option). When using control volumes pressure and density have to be on the same order of parent mesh.
 

=== modified file 'manual/examples.tex'
--- manual/examples.tex	2012-10-29 14:48:22 +0000
+++ manual/examples.tex	2013-02-08 22:48:22 +0000
@@ -307,7 +307,8 @@
 \begin{equation}
 E_b = \int \rho(z_*) g z_* \, dz_* \, ,
 \end{equation}
-where $g$ is gravity and $\rho$ the density. Most crucially, for a closed system, the background potential energy can only be altered by diapycnal mixing and increases in the background potential energy correspond to mixing in the system. The reference state is calculated using the method of \citet{tseng2001} which uses a probability density function to calculate the value of $z_*$ associated with a given $\rho$ (or here temperature). Here the probability density function is obtained from a set of `mixing bins'. Each mixing bin contains the fraction of fluid in the domain that has temperature within a given range or `class'. A set of mixing bins, the reference state and background potential energy are presented in figure \ref{fig:lock_exchange_mixing}. As the simulation progresses, the mixing and the amount of mixed fluid increases, more rapidly at first as Kelvin--Helmholtz billows form and just after the fronts hit the end walls. Later, as the dynamics become less turbulent and the fluid oscillates back and forth in the tank the mixing is less vigorous. The increase in mixed fluid is compensated for by a decrease in non--mixed fluid.
+where $g$ is gravity and $\rho$ the density. Most crucially, for a closed system, the background potential energy can only be altered by diapycnal mixing and increases in the background potential energy correspond to mixing in the system. The reference state is calculated using the method of \citet{tseng2001} which uses a probability density function to calculate the value of $z_*$ associated with a given $\rho$ (or here temperature). Here the probability density function is obtained from a set of `mixing bins'. Each mixing bin contains the fraction of fluid in the domain that has temperature within a given range or `class'. A set of mixing bins, the reference state and background potential energy are presented in figure \ref{fig:lock_exchange_mixing}. As the simulation progresses, the mixing and the amount of mixed fluid increases, more rapidly at first as Kelvin--Helmholtz billows form and just after the fronts hit the end walls. Later, as the dynamics become less turbulent and the fluid oscillates back 
+and forth in the tank the mixing is less vigorous. The increase in mixed fluid is compensated for by a decrease in non--mixed fluid.
 
 
 \begin{figure}[ht]
@@ -1000,7 +1001,8 @@
 \label{sec:water_collapse}
 
 \subsection{Overview}
-A commonly used validation experiment for multi-material models is that of a collapsing column of liquid, normally water, within an atmosphere or vacuum \citep{lakehal_interface_2002}, also known as the dam break problem.  In the experimental set-up a reservoir of water is held behind an impermeable barrier separating it from the rest of the tank.  The barrier is then quickly removed, allowing the water column to collapse and flood the remaining sections of the tank.  In the numerical analogue the initial condition is generally taken as the trapped water column, still behind the dam.  At the start of the simulation the barrier is imagined to have been removed instantaneously and switching on gravity, $|g| = 9.81$, causes the column to collapse.   Several experimental set-ups have been published and used as comparison and validation tools for numerical models \citep{martin_part_1952, greaves_simulation_2006}.  Those with water depth gauges distributed throughout the tank are particularly useful, allowing the direct comparison of data.  Furthermore, pressure gauges located on the tank walls or on any obstacles within the tank provide another useful validation tool.
+A commonly used validation experiment for multi-material models is that of a collapsing column of liquid, normally water, within an atmosphere or vacuum \citep{lakehal_interface_2002}, also known as the dam break problem.  In the experimental set-up a reservoir of water is held behind an impermeable barrier separating it from the rest of the tank.  The barrier is then quickly removed, allowing the water column to collapse and flood the remaining sections of the tank.  In the numerical analogue the initial condition is generally taken as the trapped water column, still behind the dam.  At the start of the simulation the barrier is imagined to have been removed instantaneously and switching on gravity, $|g| = 9.81$, causes the column to collapse.   Several experimental set-ups have been published and used as comparison and validation tools for numerical models \citep{martin_part_1952, greaves_simulation_2006}.  Those with water depth gauges distributed throughout the tank are particularly useful, allowing the 
+direct comparison of data.  Furthermore, pressure gauges located on the tank walls or on any obstacles within the tank provide another useful validation tool.
 
 In this example, \fluidity\ is used to simulate a simple dam break experiment \citep{zhou_nonlinear_1999}.  The example demonstrates the following functionality:
 
@@ -1018,7 +1020,8 @@
 \subsection{Problem specification}
 The experiment on which this example is based was a simple dam break problem in a $3.22\times2\times1m$ (length $\times$ height $\times$ depth) tank \citep{zhou_nonlinear_1999}.  A reservoir of water $1.2\times0.6\times1m$ (length $\times$ height $\times$ depth) was held behind a barrier at one end of the tank.  Water depth gauges were placed at two points, marked H1 and H2 in Figure \ref{fig:zhouinitial}(a), at $x_1 = 2.725m$ and $2.228m$ respectively.  Additionally, a pressure gauge was located at the point marked P2 in Figure \ref{fig:zhouinitial}(a), at $x_2=0.16m$ on the wall facing the initial water column.
 
-As no variations were introduced in the third dimension, the experiment is reproduced here numerically in two dimensions within the domain $\Omega$: $x_1 \in [0,3.22]$, $x_2 \in [0,2]$ \citep{lee_numerical_2002, colagrossi_numerical_2003, park_volume-of-fluid_2009}.  The two materials (water and air) are distinguished by scalar fields $\alpha_k$ representing their volume fraction, where the volume fraction of air $\alpha_2 = 1-\alpha_1$ and $\alpha_1$ is the volume fraction of water. The initial condition of the water volume fraction is shown in Figure \ref{fig:zhouinitial}(a).  The presence of water is indicated as a blue region and the interface to air is delineated by contours at volume fraction values of $0.025$, $0.5$ and $0.975$.  The densities of the water and air are taken as $ 1,000kg\,m^{-2}$ and $1kg\,m^{-2}$ respectively.  Both fluids are treated inviscidly. As the simulation is inviscid, free slip boundary conditions are imposed on the tank bottom, $x_2=0$, and sides, $x_1=0,\,3.22$.  The top of the tank, $x_2=2$, is left open.
+As no variations were introduced in the third dimension, the experiment is reproduced here numerically in two dimensions within the domain $\Omega$: $x_1 \in [0,3.22]$, $x_2 \in [0,2]$ \citep{lee_numerical_2002, colagrossi_numerical_2003, park_volume-of-fluid_2009}.  The two materials (water and air) are distinguished by scalar fields $\alpha_k$ representing their volume fraction, where the volume fraction of air $\alpha_2 = 1-\alpha_1$ and $\alpha_1$ is the volume fraction of water. The initial condition of the water volume fraction is shown in Figure \ref{fig:zhouinitial}(a).  The presence of water is indicated as a blue region and the interface to air is delineated by contours at volume fraction values of $0.025$, $0.5$ and $0.975$.  The densities of the water and air are taken as $ 1,000kg\,m^{-2}$ and $1kg\,m^{-2}$ respectively.  Both fluids are treated inviscidly. As the simulation is inviscid, free slip boundary conditions are imposed on the tank bottom, $x_2=0$, and sides, $x_1=0,\,3.22$.  The top of 
+the tank, $x_2=2$, is left open.
 
 \begin{figure}[tbp]
 \hspace{1cm}(a)
@@ -1097,7 +1100,8 @@
 \end{center}
 \end{figure}
 
-For validation purposes, the post-processing script provided for this example extracts information from the detector file (pressure) and from the water volume fraction field stored in the vtu files (water depth) and plots these data in comparison with the experimental results. The water depth gauge data is displayed in Figure \ref{fig:zhoudepth} alongside the experimental data.  The numerical results show the total thickness of water at the points H1 and H2, discounting any air bubbles that cross the gauges.  The simulation results show a close similarity to the experimental results with the exception of a small lip of water when the initial water head passes the gauge.  It is unclear what causes this structure, though it may be related to the initial withdrawal of the barrier in the experiment or drag effects from the bottom of the tank.  All previous published attempts to model the experiment also fail to reproduce this initial lip \citep{zhou_nonlinear_1999, lee_numerical_2002, colagrossi_numerical_2003, park_volume-of-fluid_2009}.
+For validation purposes, the post-processing script provided for this example extracts information from the detector file (pressure) and from the water volume fraction field stored in the vtu files (water depth) and plots these data in comparison with the experimental results. The water depth gauge data is displayed in Figure \ref{fig:zhoudepth} alongside the experimental data.  The numerical results show the total thickness of water at the points H1 and H2, discounting any air bubbles that cross the gauges.  The simulation results show a close similarity to the experimental results with the exception of a small lip of water when the initial water head passes the gauge.  It is unclear what causes this structure, though it may be related to the initial withdrawal of the barrier in the experiment or drag effects from the bottom of the tank.  All previous published attempts to model the experiment also fail to reproduce this initial lip \citep{zhou_nonlinear_1999, lee_numerical_2002, colagrossi_numerical_2003, 
+park_volume-of-fluid_2009}.
 
 After $t=1.5s$ the overturning wave starts to pass the water gauges and the match between the experimental results and the numerical simulation deteriorates. As would be expected from such complex behaviour, all previous published attempts have also failed to reproduce the experimental depth gauge data after this point.  However, the broad pattern and average depth observed in the simulation after $t=1.5s$ can be seen in Figure \ref{fig:zhoudepth} to match the experiment reasonably well.
 
@@ -1384,19 +1388,22 @@
 \label{sec:tephra_settling}
 
 \subsection{Overview}
-In this example, \fluidity\ is used to replicate a laboratory experiment of tephra (fine volcanic ash particles) settling through a tank of water \citep{carey1997}. The \cite{carey1997} experiments introduced tephra particles into a $0.3\times 0.3\times 0.7$ m tank, filled with water, from above using a delivery system and a particle disperser. The particles settled through the air in the tank at an approximately constant rate until they landed in the water and began to settle through the water at a much reduced velocity. While the particles in the water were sufficiently dispersed their settling velocity was that predicted by Stokes' flow (a few mm/s).  However, the build-up of particles caused by the air-water interface eventually created a layer of particles and water with a bulk density so great, relative to the density of the particle-poor water beneath, as to be gravitationally unstable and promote the formation of a vertical gravity current (plumes). The settling velocity of these plumes of particles and water was observed to be an order of magnitude greater than the Stokes settling velocity of the individual particles.
-
-In the \fluidity\ simulations, rather than represent each tephra particle individually using a multimaterial approach (which would be prohibitively expensive), a ``dispersed multiphase'' approach is adopted (see \ref{sec:multiphase_equations}), in which one phase represents the water in the tank and the other phase represents all the particles dispersed within the water. Similar to the water column collapse example (\ref{sec:water_collapse}), a volume fraction field is used to distinguish between the continuous phase (water) and the dispersed phase (particles). In an element the volume fraction of particles represents the fraction of the element volume that is occupied by particles, which is typically very small. In contrast to the multimaterial approach, however, the multiphase approach assigns each phase a separate velocity, allowing the dispersed phase (particles) to move through the continuous phase (water). In this example, the dispersed phase is more dense than the continuous phase and so sinks under the influence of gravity.
-
-The simulation replicates one of the \cite{carey1997} experiments (with a characteristic particle size of 48 $\mu$m) by considering a constant influx of the particle phase into a 2D tank of water. The simulation results can be compared to the experiment in terms of the conditions under which plumes of tephra particles form as well as the settling velocities of the individual particles and plumes. The example demonstrates the following functionality:
+In this example, \fluidity\ is used to replicate a laboratory experiment of tephra (fine volcanic ash particles) settling through a tank of water \citep{carey1997}. The experiments by \cite{carey1997} introduced tephra particles into a $0.3\times 0.3\times 0.7$ m tank, filled with water, from above using a delivery system and a particle disperser. The particles settled through the air in the tank at an approximately constant rate until they landed in the water and began to settle through the water at a much reduced velocity. While the particles in the water were sufficiently dispersed their settling velocity was that predicted by Stokes' flow (a few mm/s).  However, the build-up of particles caused by the air-water interface eventually created a layer of particles and water with a bulk density so great, relative to the density of the particle-poor water beneath, as to be gravitationally unstable and promote the formation of a vertical gravity current (plumes). The settling velocity of these plumes of particles 
+and water was observed to be an order of magnitude greater than the Stokes settling velocity of the individual particles.
+
+In the \fluidity\ simulations, rather than represent each tephra particle individually using a multimaterial approach (which would be prohibitively expensive), a ``dispersed multiphase'' approach is adopted (see \ref{sec:multiphase_equations}), in which one phase represents the water in the tank and the other phase represents all the particles dispersed within the water. Similar to the water column collapse example (\ref{sec:water_collapse}), a volume fraction field is used to distinguish between the continuous phase (water) and the dispersed phase (particles). In an element the volume fraction of particles represents the fraction of the element volume that is occupied by particles, which is typically very small. In contrast to the multimaterial approach, however, the multiphase approach assigns each phase a separate velocity, allowing the dispersed phase (particles) to move through the continuous phase (water). In this example, the dispersed phase is more dense than the continuous phase and so sinks under 
+the influence of gravity.
+
+The simulation replicates one of the experiments by \cite{carey1997} (with a characteristic particle size of 48 $\mu$m) by considering a constant influx of the particle phase into a 2D tank of water. The simulation results can be compared to the experiment in terms of the conditions under which plumes of tephra particles form as well as the settling velocities of the individual particles and plumes. The example demonstrates the following functionality:
 \begin{itemize}
 \item 2D flow
 \item Multi-phase flow
 \item Incompressible flow
 \item Adaptive timestepping
 \end{itemize}
+The simulation illustrated here took $\sim$1 hour to run in serial on an Intel Xeon E5430 2.66 GHz processor. 
 
-The simulation illustrated here took $\sim$1 hour to run in serial on an Intel Xeon E5430 2.66 GHz processor.
+For more information about the simulation of volcanic ash particles settling in water using \fluidity, see the paper by \cite{Jacobs_etal_2013}.
 
 \subsection{Problem specification}
 The simulation uses a 0.3 x 0.7 metre domain, replicating the cross-section of the water tank used in the experiments, with an initial characteristic element size of 0.01 metres. No normal flow boundary conditions are weakly imposed along with a zero velocity initial condition. The parameters for density, viscosity, particle diameter and gravity are $\rho_p = 2340\ \mathrm{kgm^{-3}}$, $\rho_l = 1000\ \mathrm{kgm^{-3}}$, $\mu_l = 0.001\ \mathrm{Pa\cdot s}$, $d = 48\times 10^{-6}\ \mathrm{m}$, and $\mathbf{g} = 9.8\ \mathrm{ms^{-2}}$ respectively (with the subscripts $p$ and $l$ denoting the properties of the particle/tephra and liquid phase).

=== modified file 'manual/model_equations.tex'
--- manual/model_equations.tex	2012-12-04 17:19:40 +0000
+++ manual/model_equations.tex	2013-02-08 22:48:22 +0000
@@ -831,7 +831,7 @@
 \nabla\cdot\left(\alpha_i\tautens_i\right) +
 \mathbf{F}_i,
 \end{equation}
-where $\mathbf{u}_i$, $\rho_i$, $\tautens_i$ and $\alpha_i$ are the velocity, density, viscous stress tensor and volume fraction of phase $i$ respectively, and $\mathbf{F}_i$ represents the forces imposed on phase $i$ by the other $N-1$ phases. Details of momentum transfer terms are given in Chapter \ref{chap:configuration}.
+where $\mathbf{u}_i$, $\rho_i$, $\tautens_i$ and $\alpha_i$ are the velocity, density, viscous stress tensor and volume fraction of phase $i$ respectively, and $\mathbf{F}_i$ represents the forces imposed on phase $i$ by the other $N-1$ phases. Details of momentum transfer terms are given in Chapter \ref{chap:configuration}. For more information about the incompressible multiphase flow model implemented in Fluidity, see the paper by \cite{Jacobs_etal_2013}.
 
 \subsubsection{Compressible flow}
 Fluidity currently only supports the case where the continuous phase is compressible. All other phases (the dispersed phases) are incompressible. The momentum equation is the same as the one given above, but the continuity equation becomes:
@@ -842,9 +842,9 @@
 
 An internal energy equation may also need to be solved depending on the equation of state used for the compressible phase. This is given by:
 \begin{equation}
-\alpha_i\rho_i\frac{\partial e_i}{\partial t} + \alpha_i\rho_i\mathbf{u}_i\cdot\nabla e_i = -p\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) + Q_i
+\alpha_i\rho_i\frac{\partial e_i}{\partial t} + \alpha_i\rho_i\mathbf{u}_i\cdot\nabla e_i = -\alpha_i p\nabla\cdot\mathbf{u}_i + \nabla\cdot\left(\alpha_i\frac{k_i}{C_{v,i}}\nabla e_i\right) + Q_i
 \end{equation}
-where $e_i$ is the internal energy and $Q_i$ is an inter-phase energy transfer term described in Chapter \ref{chap:configuration}.
+where $e_i$ is the internal energy, $k_i$ is the effective conductivity, $C_{v,i}$ is the specific heat (at constant volume), and $Q_i$ is an inter-phase energy transfer term described in Chapter \ref{chap:configuration}. Note that, in Fluidity, the pressure term is always neglected in the (incompressible) particle phase's internal energy equation.
 
 \subsection{Porous Media Darcy Flow}
 \label{sec:porous_media_darcy_flow_equations}

=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2013-02-06 18:43:09 +0000
+++ schemas/prognostic_field_options.rnc	2013-02-08 22:48:22 +0000
@@ -1947,14 +1947,12 @@
             attribute name { "ReducedConservationOfMass" },
             equation_coefficients
          }|
-         ## ***UNDER TESTING***
-         ##
          ## Select the equation used to solve for this field.
          ## Internal Energy equation - requires the selection of a Density field.
-         ## ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS
+         ## ONLY WORKS FOR CG and CV DISCRETISATIONS
          ## Solve the internal energy equation for this field.
          ## Requires pressure and velocity fields to be present.
-         ## Uses a nonconservative time discretisation.
+         ## Uses a non-conservative time discretisation.
          element equation {
             attribute name { "InternalEnergy" },
             equation_coefficients
@@ -2002,7 +2000,7 @@
         coefficient_discretisation_options?
       }
    )
-
+   
 coefficient_discretisation_options =
   (
     ## Provide discretisation options for the coefficient density field.

=== modified file 'schemas/prognostic_field_options.rng'
--- schemas/prognostic_field_options.rng	2013-02-06 18:43:09 +0000
+++ schemas/prognostic_field_options.rng	2013-02-08 22:48:22 +0000
@@ -2400,14 +2400,12 @@
         <ref name="equation_coefficients"/>
       </element>
       <element name="equation">
-        <a:documentation>***UNDER TESTING***
-
-Select the equation used to solve for this field.
+        <a:documentation>Select the equation used to solve for this field.
 Internal Energy equation - requires the selection of a Density field.
-ONLY WORKS FOR CONTROL VOLUME DISCRETISATIONS
+ONLY WORKS FOR CG and CV DISCRETISATIONS
 Solve the internal energy equation for this field.
 Requires pressure and velocity fields to be present.
-Uses a nonconservative time discretisation.</a:documentation>
+Uses a non-conservative time discretisation.</a:documentation>
         <attribute name="name">
           <value>InternalEnergy</value>
         </attribute>

=== modified file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml'
--- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	2012-10-17 15:59:53 +0000
+++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	2013-02-08 22:48:22 +0000
@@ -165,7 +165,7 @@
             <mass_terms/>
           </continuous_galerkin>
           <conservative_advection>
-            <real_value rank="0">0.0</real_value>
+            <real_value rank="0">1.0</real_value>
           </conservative_advection>
         </spatial_discretisation>
         <temporal_discretisation>
@@ -215,7 +215,7 @@
             </mass_terms>
             <advection_terms/>
             <stress_terms>
-              <tensor_form/>
+              <stress_form/>
             </stress_terms>
           </continuous_galerkin>
           <conservative_advection>
@@ -268,11 +268,11 @@
         <tensor_field name="Viscosity" rank="2">
           <prescribed>
             <value name="WholeMesh">
-              <isotropic>
+              <anisotropic_asymmetric>
                 <constant>
-                  <real_value rank="0">1.83e-5</real_value>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.83e-5</real_value>
                 </constant>
-              </isotropic>
+              </anisotropic_asymmetric>
             </value>
             <output/>
           </prescribed>

=== modified file 'tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml'
--- tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml	2012-10-17 15:59:53 +0000
+++ tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml	2013-02-08 22:48:22 +0000
@@ -165,7 +165,7 @@
             <mass_terms/>
           </continuous_galerkin>
           <conservative_advection>
-            <real_value rank="0">0.0</real_value>
+            <real_value rank="0">1.0</real_value>
           </conservative_advection>
         </spatial_discretisation>
         <temporal_discretisation>
@@ -215,7 +215,7 @@
             </mass_terms>
             <advection_terms/>
             <stress_terms>
-              <tensor_form/>
+              <stress_form/>
             </stress_terms>
           </continuous_galerkin>
           <conservative_advection>
@@ -268,11 +268,11 @@
         <tensor_field name="Viscosity" rank="2">
           <prescribed>
             <value name="WholeMesh">
-              <isotropic>
+              <anisotropic_asymmetric>
                 <constant>
-                  <real_value rank="0">1.83e-5</real_value>
+                  <real_value symmetric="false" dim2="dim" shape="1 1" dim1="dim" rank="2">1.83e-5</real_value>
                 </constant>
-              </isotropic>
+              </anisotropic_asymmetric>
             </value>
             <output/>
           </prescribed>

=== added directory 'tests/mphase_mms_p2p1_compressible_ie_heat_flux'
=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/MMS_X.flml'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/MMS_X.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/MMS_X.flml	2013-02-08 22:48:22 +0000
@@ -0,0 +1,1296 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">MMS_X</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="src/MMS_X">
+        <format name="gmsh"/>
+        <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="ErrorMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">4</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>continuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">8</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1000.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">999.9</real_value>
+      <comment>gives a max cfl number of approximately 0.1</comment>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1000.0</real_value>
+      <comment>10.0</comment>
+    </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>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1.E-7</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">1.0</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="2" dim1="dim" rank="1">0.707106781 0.707106781</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Phase1">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.4</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/>
+          </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">10000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <boundary_conditions name="Incoming">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.p(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <output>
+          <include_previous_time_step/>
+        </output>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+        <consistent_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">0.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):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.rho1(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.rho1(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.forcing_rho1(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix>
+                <use_submesh/>
+              </lump_mass_matrix>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</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):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.velocity1(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="4" rank="1">7 8 9 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.u1(X)</string_value>
+                </python>
+              </x_component>
+              <y_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.v1(X)</string_value>
+                </python>
+              </y_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.7</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <vector_field name="Source" rank="1">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   density = tools.rho1(X)
+   forcing_velocity1 = tools.forcing_velocity1(X)
+   # Remember to divide through by density here
+   # because we multiply through by Density in
+   # Momentum_CG.F90.
+   return (forcing_velocity1[0]/density, forcing_velocity1[1]/density)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </vector_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>
+          <exclude_from_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>
+          <exclude_from_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>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <advection_terms/>
+            <mass_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>
+        </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">10000</integer_value>
+          </max_iterations>
+          <ignore_all_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie1(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie1(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Diffusivity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.forcing_ie1(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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="InternalEnergyAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="InternalEnergyError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="InternalEnergyAnalytical" source_field_2_name="InternalEnergyProjection" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="InternalEnergyProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="InternalEnergy">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="PressureAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.p(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="PressureProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="Pressure">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="PressureError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="PressureProjection" source_field_2_name="PressureAnalytical" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="DensityAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.rho1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="DensityError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="DensityAnalytical" source_field_2_name="DensityProjection" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="DensityProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="Density">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="VelocityAnalytical" rank="1">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.velocity1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <vector_field name="VelocityProjection" rank="1">
+      <diagnostic>
+        <algorithm source_field_type="vector" material_phase_support="single" name="vector_galerkin_projection" source_field_name="Velocity">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <vector_field name="VelocityError" rank="1">
+      <diagnostic>
+        <algorithm source_field_2_type="vector" name="vector_difference" source_field_1_name="VelocityProjection" source_field_2_name="VelocityAnalytical" material_phase_support="single" source_field_1_type="vector">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+  <material_phase name="Phase2">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">3.0</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Phase1" 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>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix>
+                <use_submesh/>
+              </lump_mass_matrix>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</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):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.velocity2(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="4" rank="1">7 8 9 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.u2(X)</string_value>
+                </python>
+              </x_component>
+              <y_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.v2(X)</string_value>
+                </python>
+              </y_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.7</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <vector_field name="Source" rank="1">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   density = tools.rho2(X)
+   forcing_velocity2 = tools.forcing_velocity2(X)
+   # Remember to divide through by density here
+   # because we multiply through by Density in
+   # Momentum_CG.F90.
+   return (forcing_velocity2[0]/density, forcing_velocity2[1]/density)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </vector_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>
+          <exclude_from_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>
+          <exclude_from_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>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <advection_terms/>
+            <mass_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>
+        </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">10000</integer_value>
+          </max_iterations>
+          <ignore_all_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie2(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie2(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Diffusivity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.forcing_ie2(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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="InternalEnergyAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="InternalEnergyError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="InternalEnergyAnalytical" source_field_2_name="InternalEnergyProjection" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="InternalEnergyProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="InternalEnergy">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="PhaseVolumeFraction" rank="0">
+      <prescribed>
+        <mesh name="CoordinateMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.vfrac2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="InternalEnergyError2" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="InternalEnergy" source_field_2_name="InternalEnergyAnalytical2" material_phase_support="single" source_field_1_type="scalar"/>
+        <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>
+    <scalar_field name="InternalEnergyAnalytical2" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.ie2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <vector_field name="VelocityAnalytical" rank="1">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_tools as tools
+   return tools.velocity2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <vector_field name="VelocityProjection" rank="1">
+      <diagnostic>
+        <algorithm source_field_type="vector" material_phase_support="single" name="vector_galerkin_projection" source_field_name="Velocity">
+          <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">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <vector_field name="VelocityError" rank="1">
+      <diagnostic>
+        <algorithm source_field_2_type="vector" name="vector_difference" source_field_1_name="VelocityProjection" source_field_2_name="VelocityAnalytical" material_phase_support="single" source_field_1_type="vector">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/Makefile'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/Makefile	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/Makefile	2013-02-08 22:48:22 +0000
@@ -0,0 +1,8 @@
+input:	clean
+	./copy_script
+	gmsh -2 -bin src/MMS_A.geo
+	gmsh -2 -bin src/MMS_B.geo
+	gmsh -2 -bin src/MMS_C.geo
+
+clean:
+	rm -f *.vtu *.stat *.s *.log *_A* *_B* *_C* matrix* *~ mms_rans_p2p1_keps.py *.pyc fluidity.*

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/README'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/README	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/README	2013-02-08 22:48:22 +0000
@@ -0,0 +1,25 @@
+CREATING mms_rans_p2p1_keps_tools.py:
+
+edit mms_rans_p2p1_keps.sage as required
+in the terminal, run:
+sage mms_rans_p2p1_keps.sage > mms_rans_p2p1_keps_tools.py
+
+EDITING MMS_X.flml and MMS_X.geo (copyscript):
+
+These are copied to MMS_A, MMS_B, MMS_C, MMS_D cases by copyscript
+DO NOT CHANGE - timestep or references to MMS_X in MMS_X.flml
+DO NOT CHANGE - n = XX in MMS_X.geo
+
+DISABLING PARTS OF THE MODEL:
+
+alternative flml files are included in the alternative_flml folder for debugging failures - these include flml files that have prescribed velocity or k-epsilon fields so that the model can be debugged in stages.
+
+additionally, individual source and absorbtion terms can be disabled within the flml.
+
+FUNCTION PRINTER
+
+run using:
+python function_printer.py AA BB CC DD .. n_rows
+where:
+	AA, BB, CC, DD are names of functions in mms_rans_p2p1_keps_tools.py (any number can be entered)
+	n_rows is the number of rows to display the functions on
\ No newline at end of file

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/copy_script'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/copy_script	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/copy_script	2013-02-08 22:48:22 +0000
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+n=4
+for F in A B C
+do
+    echo creating set $F with n_ele=$n	
+
+    cp MMS_X.flml MMS_$F.flml 
+    sed s/MMS_X/MMS_$F/ MMS_$F.flml > MMS_$F.flml.cp
+
+    #sed s/XX/$n/ src/MMS_X.geo > src/MMS_$F.geo	
+
+    n=$((n*2))
+done
+
+echo setting flml timestep values
+
+sed s/999.9/0.005/ MMS_A.flml.cp > MMS_A.flml
+rm MMS_A.flml.cp
+sed s/999.9/0.0025/ MMS_B.flml.cp > MMS_B.flml
+rm MMS_B.flml.cp
+sed s/999.9/0.00125/ MMS_C.flml.cp > MMS_C.flml
+rm MMS_C.flml.cp

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/function_printer.py'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/function_printer.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/function_printer.py	2013-02-08 22:48:22 +0000
@@ -0,0 +1,40 @@
+from mphase_mms_p2p1_compressible_ie_tools import *
+from numpy import *
+import matplotlib
+import matplotlib.pyplot as plt
+import sys
+
+'''
+run using:
+python function_printer.py AA BB CC DD .. n_rows
+where:
+	AA, BB, CC, DD are names of functions in mms_rans_p2p1_keps_tools.py (any number can be entered)
+	n_rows is the number of rows to display the functions on
+'''
+
+functions = []
+for arg in sys.argv[1:-1]:
+    functions.append(arg)
+n_rows = int(sys.argv[-1])
+
+matplotlib.rcParams['xtick.direction'] = 'out'
+matplotlib.rcParams['ytick.direction'] = 'out'
+plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
+                    wspace=0.2, hspace=0.2)
+   
+res = 50
+X = linspace(0.0, pi, res)
+Y = linspace(0.0, pi, res)
+x = [0,0]
+
+data = empty([len(functions), res, res])
+for z, function in enumerate(functions):
+    for j, x[0] in enumerate(X):
+        for i, x[1] in enumerate(Y):
+            data[z,i,j] = eval(function + '(x)')
+    plt.subplot(n_rows, len(functions)/n_rows + 1, z+1)
+    CS = plt.contour(X, Y, data[z])
+    plt.clabel(CS, inline=1, fontsize=10)
+    plt.title(functions[z])
+    
+plt.show()

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.py'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.py	2013-02-08 22:48:22 +0000
@@ -0,0 +1,141 @@
+# This file was *autogenerated* from the file mphase_mms_p2p1_compressible_ie.sage.
+from sage.all_cmdline import *   # import sage library
+_sage_const_2 = Integer(2); _sage_const_1p0 = RealNumber('1.0'); _sage_const_500p0 = RealNumber('500.0'); _sage_const_1p4 = RealNumber('1.4'); _sage_const_1p5 = RealNumber('1.5'); _sage_const_100p0 = RealNumber('100.0'); _sage_const_0p1 = RealNumber('0.1'); _sage_const_0p5 = RealNumber('0.5'); _sage_const_0p707106781 = RealNumber('0.707106781'); _sage_const_0p7 = RealNumber('0.7'); _sage_const_0p8 = RealNumber('0.8'); _sage_const_3p0 = RealNumber('3.0')
+y = var('y')
+
+def function(phi_0, phi_x, phi_y, phi_xy, 
+             f_sin_x, f_cos_x, f_sin_y, f_cos_y, f_sin_xy, f_cos_xy, 
+             alpha_x, alpha_y, alpha_xy):
+    
+    f_0 = phi_0 
+    f_x = phi_x*(f_sin_x*sin(alpha_x*x) + f_cos_x*cos(alpha_x*x)) 
+    f_y = phi_y*(f_sin_y*sin(alpha_y*y) + f_cos_y*cos(alpha_y*y)) 
+    f_xy = phi_xy*(f_sin_xy*sin(alpha_xy*x*y/pi) + f_cos_xy*cos(alpha_xy*x*y/pi)) 
+    f = f_0 + f_x + f_y + f_xy
+    return f
+             
+rho1 = _sage_const_0p5 *(sin(x*x + y*y) + _sage_const_1p5 )
+rho2 = _sage_const_3p0 
+
+ie1 = _sage_const_0p5 *(cos(x + y) + _sage_const_1p5 )
+ie2 = _sage_const_3p0 *y
+
+u1 = _sage_const_1p0 *(sin(x**_sage_const_2 +y**_sage_const_2 )+_sage_const_0p5 )
+v1 = _sage_const_0p1 *(cos(x**_sage_const_2 +y**_sage_const_2 )+_sage_const_0p5 )
+
+u2 = _sage_const_1p0 *(sin(x**_sage_const_2 +y**_sage_const_2 )+_sage_const_0p5 )
+v2 = _sage_const_0p1 *(cos(x**_sage_const_2 +y**_sage_const_2 )+_sage_const_0p5 )
+ 
+gamma = _sage_const_1p4 
+rho0 = _sage_const_0p5 
+csq = _sage_const_100p0 
+#p = csq*(rho1 - rho0)
+p = (gamma-_sage_const_1p0 )*ie1*rho1
+             
+vfrac1 = _sage_const_0p8 
+vfrac2 = _sage_const_1p0  - vfrac1
+             
+nu1 = _sage_const_0p7 
+nu2 = _sage_const_0p7 
+
+g_x = _sage_const_0p707106781 
+g_y = _sage_const_0p707106781 
+
+tau_xx1 = nu1*diff(u1,x)            
+tau_yy1 = nu1*diff(v1,y)
+tau_xy1 = nu1*diff(u1,y)
+tau_yx1 = nu1*diff(v1,x)  
+
+tau_xx2 = nu2*diff(u2,x)            
+tau_yy2 = nu2*diff(v2,y)
+tau_xy2 = nu2*diff(u2,y)
+tau_yx2 = nu2*diff(v2,x)  
+
+Su1 = vfrac1*rho1*u1*diff(u1,x) + vfrac1*rho1*v1*diff(u1,y) - diff(vfrac1*tau_xx1, x) - diff(vfrac1*tau_xy1, y) - vfrac1*g_x*rho1 + vfrac1*diff(p,x)  
+Sv1 = vfrac1*rho1*u1*diff(v1,x) + vfrac1*rho1*v1*diff(v1,y) - diff(vfrac1*tau_yy1, y) - diff(vfrac1*tau_yx1, x) - vfrac1*g_y*rho1 + vfrac1*diff(p,y)  
+
+Su2 = vfrac2*rho2*u2*diff(u2,x) + vfrac2*rho2*v2*diff(u2,y) - diff(vfrac2*tau_xx2, x) - diff(vfrac2*tau_xy2, y) - vfrac2*g_x*rho2 + vfrac2*diff(p,x)  
+Sv2 = vfrac2*rho2*u2*diff(v2,x) + vfrac2*rho2*v2*diff(v2,y) - diff(vfrac2*tau_yy2, y) - diff(vfrac2*tau_yx2, x) - vfrac2*g_y*rho2 + vfrac2*diff(p,y) 
+
+Srho1 = diff(vfrac1*u1*rho1,x) + diff(vfrac1*v1*rho1,y) + rho1*diff(u2*vfrac2, x) + rho1*diff(v2*vfrac2, y)
+
+k = _sage_const_0p5 
+Cv1 = _sage_const_500p0 
+Cv2 = _sage_const_500p0 
+
+heat_flux_x1 = (k/Cv1)*vfrac1*diff(ie1,x)
+heat_flux_y1 = (k/Cv1)*vfrac1*diff(ie1,y)
+heat_flux_x2 = (k/Cv2)*vfrac2*diff(ie2,x)
+heat_flux_y2 = (k/Cv2)*vfrac2*diff(ie2,y)
+
+Sie1 = vfrac1*rho1*u1*diff(ie1, x) + vfrac1*rho1*v1*diff(ie1, y) + vfrac1*p*diff(u1, x) + vfrac1*p*diff(v1, y) - diff(heat_flux_x1, x) - diff(heat_flux_y1, y)
+Sie2 = vfrac2*rho2*u2*diff(ie2, x) + vfrac2*rho2*v2*diff(ie2, y) - diff(heat_flux_x2, x) - diff(heat_flux_y2, y)
+
+print 'from math import sin, cos, tanh, pi'
+print ''
+print 'def u1(X):'
+print '    return', str(u1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v1(X):'
+print '    return', str(v1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def p(X):'
+print '    return', str(p).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho1(X):'
+print '    return', str(rho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie1(X):'
+print '    return', str(ie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac1(X):'
+print '    return', str(vfrac1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u1(X):'
+print '    return', str(Su1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v1(X):'
+print '    return', str(Sv1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_rho1(X):'
+print '    return', str(Srho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie1(X):'
+print '    return', str(Sie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity1(X):'
+print '   return [u1(X), v1(X)]'
+print ''
+print 'def forcing_velocity1(X):'
+print '   return [forcing_u1(X), forcing_v1(X)]'
+print ''
+print 'def u2(X):'
+print '    return', str(u2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v2(X):'
+print '    return', str(v2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho2(X):'
+print '    return', str(rho2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie2(X):'
+print '    return', str(ie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac2(X):'
+print '    return', str(vfrac2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u2(X):'
+print '    return', str(Su2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v2(X):'
+print '    return', str(Sv2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie2(X):'
+print '    return', str(Sie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity2(X):'
+print '   return [u2(X), v2(X)]'
+print ''
+print 'def forcing_velocity2(X):'
+print '   return [forcing_u2(X), forcing_v2(X)]'
+print ''

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.sage'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.sage	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie.sage	2013-02-08 22:48:22 +0000
@@ -0,0 +1,138 @@
+y = var('y')
+
+def function(phi_0, phi_x, phi_y, phi_xy, 
+             f_sin_x, f_cos_x, f_sin_y, f_cos_y, f_sin_xy, f_cos_xy, 
+             alpha_x, alpha_y, alpha_xy):
+    
+    f_0 = phi_0 
+    f_x = phi_x*(f_sin_x*sin(alpha_x*x) + f_cos_x*cos(alpha_x*x)) 
+    f_y = phi_y*(f_sin_y*sin(alpha_y*y) + f_cos_y*cos(alpha_y*y)) 
+    f_xy = phi_xy*(f_sin_xy*sin(alpha_xy*x*y/pi) + f_cos_xy*cos(alpha_xy*x*y/pi)) 
+    f = f_0 + f_x + f_y + f_xy
+    return f
+             
+rho1 = 0.5*(sin(x*x + y*y) + 1.5)
+rho2 = 3.0
+
+ie1 = 0.5*(cos(x + y) + 1.5)
+ie2 = 3.0*y
+
+u1 = 1.0*(sin(x**2+y**2)+0.5)
+v1 = 0.1*(cos(x**2+y**2)+0.5)
+
+u2 = 1.0*(sin(x**2+y**2)+0.5)
+v2 = 0.1*(cos(x**2+y**2)+0.5)
+ 
+gamma = 1.4
+rho0 = 0.5
+csq = 100.0
+#p = csq*(rho1 - rho0)
+p = (gamma-1.0)*ie1*rho1
+             
+vfrac1 = 0.8
+vfrac2 = 1.0 - vfrac1
+             
+nu1 = 0.7
+nu2 = 0.7
+
+g_x = 0.707106781
+g_y = 0.707106781
+
+tau_xx1 = nu1*diff(u1,x)            
+tau_yy1 = nu1*diff(v1,y)
+tau_xy1 = nu1*diff(u1,y)
+tau_yx1 = nu1*diff(v1,x)  
+
+tau_xx2 = nu2*diff(u2,x)            
+tau_yy2 = nu2*diff(v2,y)
+tau_xy2 = nu2*diff(u2,y)
+tau_yx2 = nu2*diff(v2,x)  
+
+Su1 = vfrac1*rho1*u1*diff(u1,x) + vfrac1*rho1*v1*diff(u1,y) - diff(vfrac1*tau_xx1, x) - diff(vfrac1*tau_xy1, y) - vfrac1*g_x*rho1 + vfrac1*diff(p,x)  
+Sv1 = vfrac1*rho1*u1*diff(v1,x) + vfrac1*rho1*v1*diff(v1,y) - diff(vfrac1*tau_yy1, y) - diff(vfrac1*tau_yx1, x) - vfrac1*g_y*rho1 + vfrac1*diff(p,y)  
+
+Su2 = vfrac2*rho2*u2*diff(u2,x) + vfrac2*rho2*v2*diff(u2,y) - diff(vfrac2*tau_xx2, x) - diff(vfrac2*tau_xy2, y) - vfrac2*g_x*rho2 + vfrac2*diff(p,x)  
+Sv2 = vfrac2*rho2*u2*diff(v2,x) + vfrac2*rho2*v2*diff(v2,y) - diff(vfrac2*tau_yy2, y) - diff(vfrac2*tau_yx2, x) - vfrac2*g_y*rho2 + vfrac2*diff(p,y) 
+
+Srho1 = diff(vfrac1*u1*rho1,x) + diff(vfrac1*v1*rho1,y) + rho1*diff(u2*vfrac2, x) + rho1*diff(v2*vfrac2, y)
+
+k = 0.5
+Cv1 = 500.0
+Cv2 = 500.0
+
+heat_flux_x1 = (k/Cv1)*vfrac1*diff(ie1,x)
+heat_flux_y1 = (k/Cv1)*vfrac1*diff(ie1,y)
+heat_flux_x2 = (k/Cv2)*vfrac2*diff(ie2,x)
+heat_flux_y2 = (k/Cv2)*vfrac2*diff(ie2,y)
+
+Sie1 = vfrac1*rho1*u1*diff(ie1, x) + vfrac1*rho1*v1*diff(ie1, y) + vfrac1*p*diff(u1, x) + vfrac1*p*diff(v1, y) - diff(heat_flux_x1, x) - diff(heat_flux_y1, y)
+Sie2 = vfrac2*rho2*u2*diff(ie2, x) + vfrac2*rho2*v2*diff(ie2, y) - diff(heat_flux_x2, x) - diff(heat_flux_y2, y)
+
+print 'from math import sin, cos, tanh, pi'
+print ''
+print 'def u1(X):'
+print '    return', str(u1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v1(X):'
+print '    return', str(v1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def p(X):'
+print '    return', str(p).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho1(X):'
+print '    return', str(rho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie1(X):'
+print '    return', str(ie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac1(X):'
+print '    return', str(vfrac1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u1(X):'
+print '    return', str(Su1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v1(X):'
+print '    return', str(Sv1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_rho1(X):'
+print '    return', str(Srho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie1(X):'
+print '    return', str(Sie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity1(X):'
+print '   return [u1(X), v1(X)]'
+print ''
+print 'def forcing_velocity1(X):'
+print '   return [forcing_u1(X), forcing_v1(X)]'
+print ''
+print 'def u2(X):'
+print '    return', str(u2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v2(X):'
+print '    return', str(v2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho2(X):'
+print '    return', str(rho2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie2(X):'
+print '    return', str(ie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac2(X):'
+print '    return', str(vfrac2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u2(X):'
+print '    return', str(Su2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v2(X):'
+print '    return', str(Sv2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie2(X):'
+print '    return', str(Sie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity2(X):'
+print '   return [u2(X), v2(X)]'
+print ''
+print 'def forcing_velocity2(X):'
+print '   return [forcing_u2(X), forcing_v2(X)]'
+print ''

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie_tools.py'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie_tools.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/mphase_mms_p2p1_compressible_ie_tools.py	2013-02-08 22:48:22 +0000
@@ -0,0 +1,68 @@
+from math import sin, cos, tanh, pi
+
+def u1(X):
+    return sin(X[0]**2 + X[1]**2) + 0.500
+
+def v1(X):
+    return 0.100*cos(X[0]**2 + X[1]**2) + 0.0500
+
+def p(X):
+    return (0.200*cos(X[0] + X[1]) + 0.300)*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)
+
+def rho1(X):
+    return 0.500*sin(X[0]**2 + X[1]**2) + 0.750
+
+def ie1(X):
+    return 0.500*cos(X[0] + X[1]) + 0.750
+
+def vfrac1(X):
+    return 0.800
+
+def forcing_u1(X):
+    return 2*(0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*X[1]*cos(X[0]**2 + X[1]**2) + 2*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*X[0]*cos(X[0]**2 + X[1]**2) + 0.800*(0.200*cos(X[0] + X[1]) + 0.300)*X[0]*cos(X[0]**2 + X[1]**2) + 2.24*X[0]**2*sin(X[0]**2 + X[1]**2) + 2.24*X[1]**2*sin(X[0]**2 + X[1]**2) - 0.160*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*sin(X[0] + X[1]) - 0.282842712400000*sin(X[0]**2 + X[1]**2) - 2.24*cos(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_v1(X):
+    return -0.200*(0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*X[1]*sin(X[0]**2 + X[1]**2) - 0.200*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*X[0]*sin(X[0]**2 + X[1]**2) + 0.800*(0.200*cos(X[0] + X[1]) + 0.300)*X[1]*cos(X[0]**2 + X[1]**2) + 0.224*X[0]**2*cos(X[0]**2 + X[1]**2) + 0.224*X[1]**2*cos(X[0]**2 + X[1]**2) - 0.160*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*sin(X[0] + X[1]) - 0.0588427124000001*sin(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_rho1(X):
+    return (0.0800*cos(X[0]**2 + X[1]**2) + 0.0400)*X[1]*cos(X[0]**2 + X[1]**2) + 2.00*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[0]*cos(X[0]**2 + X[1]**2) - 0.200*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[1]*sin(X[0]**2 + X[1]**2) + (0.800*sin(X[0]**2 + X[1]**2) + 0.400)*X[0]*cos(X[0]**2 + X[1]**2)
+
+def forcing_ie1(X):
+    return 1.60*(0.200*cos(X[0] + X[1]) + 0.300)*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[0]*cos(X[0]**2 + X[1]**2) - 0.160*(0.200*cos(X[0] + X[1]) + 0.300)*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[1]*sin(X[0]**2 + X[1]**2) - 0.500*(0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*sin(X[0] + X[1]) - 0.500*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*sin(X[0] + X[1]) + 0.000800*cos(X[0] + X[1])
+
+def velocity1(X):
+   return [u1(X), v1(X)]
+
+def forcing_velocity1(X):
+   return [forcing_u1(X), forcing_v1(X)]
+
+def u2(X):
+    return sin(X[0]**2 + X[1]**2) + 0.500
+
+def v2(X):
+    return 0.100*cos(X[0]**2 + X[1]**2) + 0.0500
+
+def rho2(X):
+    return 3.00
+
+def ie2(X):
+    return 3.00*X[1]
+
+def vfrac2(X):
+    return 0.200
+
+def forcing_u2(X):
+    return 2*(0.0600*cos(X[0]**2 + X[1]**2) + 0.0300)*X[1]*cos(X[0]**2 + X[1]**2) + 0.200*(0.200*cos(X[0] + X[1]) + 0.300)*X[0]*cos(X[0]**2 + X[1]**2) + 2*(0.600*sin(X[0]**2 + X[1]**2) + 0.300)*X[0]*cos(X[0]**2 + X[1]**2) + 0.560*X[0]**2*sin(X[0]**2 + X[1]**2) + 0.560*X[1]**2*sin(X[0]**2 + X[1]**2) - 0.0400*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*sin(X[0] + X[1]) - 0.560*cos(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_v2(X):
+    return -0.200*(0.0600*cos(X[0]**2 + X[1]**2) + 0.0300)*X[1]*sin(X[0]**2 + X[1]**2) + 0.200*(0.200*cos(X[0] + X[1]) + 0.300)*X[1]*cos(X[0]**2 + X[1]**2) - 0.200*(0.600*sin(X[0]**2 + X[1]**2) + 0.300)*X[0]*sin(X[0]**2 + X[1]**2) + 0.0560*X[0]**2*cos(X[0]**2 + X[1]**2) + 0.0560*X[1]**2*cos(X[0]**2 + X[1]**2) - 0.0400*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*sin(X[0] + X[1]) + 0.0560*sin(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_ie2(X):
+    return 0.180*cos(X[0]**2 + X[1]**2) + 0.0900
+
+def velocity2(X):
+   return [u2(X), v2(X)]
+
+def forcing_velocity2(X):
+   return [forcing_u2(X), forcing_v2(X)]
+

=== added directory 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src'
=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.geo'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.1};
+Point(2) = {0.7,0.2,0,0.1};
+Point(3) = {0.7,0.8,0,0.1};
+Point(4) = {-0.1,0.8,0,0.1};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_A.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.geo'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.05};
+Point(2) = {0.7,0.2,0,0.05};
+Point(3) = {0.7,0.8,0,0.05};
+Point(4) = {-0.1,0.8,0,0.05};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_B.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.geo'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.025};
+Point(2) = {0.7,0.2,0,0.025};
+Point(3) = {0.7,0.8,0,0.025};
+Point(4) = {-0.1,0.8,0,0.025};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_C.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_X.geo'
--- tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_X.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_heat_flux/src/MMS_X.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,18 @@
+pi = 3.141592653589793;
+n = XX;
+
+Point(1) = {0.0, 0.0, 0.0, 1};
+
+Extrude {pi, 0.0, 0.0} {
+  Point{1}; Layers{n};
+}
+
+Extrude {0.0, pi, 0.0} {
+  Line{1}; Layers{n};
+}
+
+Physical Line(7) = {1};
+Physical Line(8) = {2};
+Physical Line(9) = {3};
+Physical Line(10) = {4}; 
+Physical Surface(6) = {5};

=== added directory 'tests/mphase_mms_p2p1_compressible_ie_p1cv'
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/MMS_X.flml'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/MMS_X.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/MMS_X.flml	2013-02-08 22:48:22 +0000
@@ -0,0 +1,1235 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">MMS_X</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="src/MMS_X">
+        <format name="gmsh"/>
+        <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="ErrorMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">4</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>continuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <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>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">8</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1000.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">999.9</real_value>
+      <comment>gives a max cfl number of approximately 0.1</comment>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1000.0</real_value>
+      <comment>10.0</comment>
+    </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>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1.E-7</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">1.0</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="2" dim1="dim" rank="1">0.707106781 0.707106781</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Phase1">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.4</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/>
+          </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">10000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <output>
+          <include_previous_time_step/>
+        </output>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+        <consistent_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">0.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):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.rho1(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.rho1(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.forcing_rho1(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix>
+                <use_submesh/>
+              </lump_mass_matrix>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</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="All">
+          <surface_ids>
+            <integer_value shape="4" rank="1">7 8 9 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.u1(X)</string_value>
+                </python>
+              </x_component>
+              <y_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.v1(X)</string_value>
+                </python>
+              </y_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.7</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <vector_field name="Source" rank="1">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   density = tools.rho1(X)
+   forcing_velocity1 = tools.forcing_velocity1(X)
+   # Remember to divide through by density here
+   # because we multiply through by Density in
+   # Momentum_CG.F90.
+   return (forcing_velocity1[0]/density, forcing_velocity1[1]/density)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </vector_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>
+          <exclude_from_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>
+          <exclude_from_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="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">10000</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):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie1(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie1(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Diffusivity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.forcing_ie1(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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="InternalEnergyAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="InternalEnergyError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="InternalEnergyAnalytical" source_field_2_name="InternalEnergyProjection" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="InternalEnergyProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="InternalEnergy">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="PressureAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.p(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="PressureProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="Pressure">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="PressureError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="PressureProjection" source_field_2_name="PressureAnalytical" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="DensityAnalytical" rank="0">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.rho1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="DensityError" rank="0">
+      <diagnostic>
+        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="DensityAnalytical" source_field_2_name="DensityProjection" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="DensityProjection" rank="0">
+      <diagnostic>
+        <algorithm source_field_type="scalar" material_phase_support="single" name="scalar_galerkin_projection" source_field_name="Density">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <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="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="VelocityAnalytical" rank="1">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.velocity1(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <vector_field name="VelocityProjection" rank="1">
+      <diagnostic>
+        <algorithm source_field_type="vector" material_phase_support="single" name="vector_galerkin_projection" source_field_name="Velocity">
+          <solver>
+            <iterative_method name="cg"/>
+            <preconditioner name="sor"/>
+            <relative_error>
+              <real_value rank="0">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <vector_field name="VelocityError" rank="1">
+      <diagnostic>
+        <algorithm source_field_2_type="vector" name="vector_difference" source_field_1_name="VelocityProjection" source_field_2_name="VelocityAnalytical" material_phase_support="single" source_field_1_type="vector">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <multiphase_properties>
+      <effective_conductivity>
+        <real_value rank="0">0.5</real_value>
+      </effective_conductivity>
+      <specific_heat>
+        <real_value rank="0">700.0</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <material_phase name="Phase2">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">3.0</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Phase1" field_name="Pressure"/>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="DensityMesh"/>
+        <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>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix>
+                <use_submesh/>
+              </lump_mass_matrix>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_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">1e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</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="All">
+          <surface_ids>
+            <integer_value shape="4" rank="1">7 8 9 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.u2(X)</string_value>
+                </python>
+              </x_component>
+              <y_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.v2(X)</string_value>
+                </python>
+              </y_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.7</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <vector_field name="Source" rank="1">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   density = tools.rho2(X)
+   forcing_velocity2 = tools.forcing_velocity2(X)
+   # Remember to divide through by density here
+   # because we multiply through by Density in
+   # Momentum_CG.F90.
+   return (forcing_velocity2[0]/density, forcing_velocity2[1]/density)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat>
+              <include_in_stat/>
+            </stat>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </vector_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>
+          <exclude_from_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>
+          <exclude_from_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="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">10000</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):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie2(X)</string_value>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">7 10</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <python>
+              <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie2(X)</string_value>
+            </python>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Diffusivity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">0.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <scalar_field name="Source" rank="0">
+          <prescribed>
+            <value name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.forcing_ie2(X)</string_value>
+              </python>
+            </value>
+            <output/>
+            <stat/>
+            <detectors>
+              <exclude_from_detectors/>
+            </detectors>
+          </prescribed>
+        </scalar_field>
+        <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="InternalEnergyAnalytical" rank="0">
+      <prescribed>
+        <mesh name="PressureMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.ie2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prescribed>
+        <mesh name="CoordinateMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.vfrac2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ScalarAbsoluteDifference" rank="0">
+      <diagnostic field_name_b="InternalEnergy" field_name_a="InternalEnergyAnalytical">
+        <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="VelocityAnalytical" rank="1">
+      <prescribed>
+        <mesh name="ErrorMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   import mphase_mms_p2p1_compressible_ie_p1cv_tools as tools
+   return tools.velocity2(X)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <vector_field name="VelocityProjection" rank="1">
+      <diagnostic>
+        <algorithm source_field_type="vector" material_phase_support="single" name="vector_galerkin_projection" source_field_name="Velocity">
+          <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">1e-10</real_value>
+            </relative_error>
+            <max_iterations>
+              <integer_value rank="0">1000</integer_value>
+            </max_iterations>
+            <never_ignore_solver_failures/>
+            <diagnostics>
+              <monitors/>
+            </diagnostics>
+          </solver>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <vector_field name="VelocityError" rank="1">
+      <diagnostic>
+        <algorithm source_field_2_type="vector" name="vector_difference" source_field_1_name="VelocityProjection" source_field_2_name="VelocityAnalytical" material_phase_support="single" source_field_1_type="vector">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="ErrorMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+    <multiphase_properties>
+      <particle_diameter>
+        <real_value rank="0">0.1</real_value>
+      </particle_diameter>
+      <specific_heat>
+        <real_value rank="0">700.0</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/Makefile'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/Makefile	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/Makefile	2013-02-08 22:48:22 +0000
@@ -0,0 +1,8 @@
+input:	clean
+	./copy_script
+	gmsh -2 -bin src/MMS_A.geo
+	gmsh -2 -bin src/MMS_B.geo
+	gmsh -2 -bin src/MMS_C.geo
+
+clean:
+	rm -f *.vtu *.stat *.s *.log *_A* *_B* *_C* matrix* *~ mphase_mms_p2p1_compressible_ie_p1cv.py *.pyc fluidity.*

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/copy_script'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/copy_script	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/copy_script	2013-02-08 22:48:22 +0000
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+n=4
+for F in A B C
+do
+    echo creating set $F with n_ele=$n	
+
+    cp MMS_X.flml MMS_$F.flml 
+    sed s/MMS_X/MMS_$F/ MMS_$F.flml > MMS_$F.flml.cp
+
+    #sed s/XX/$n/ src/MMS_X.geo > src/MMS_$F.geo	
+
+    n=$((n*2))
+done
+
+echo setting flml timestep values
+
+sed s/999.9/0.005/ MMS_A.flml.cp > MMS_A.flml
+rm MMS_A.flml.cp
+sed s/999.9/0.0025/ MMS_B.flml.cp > MMS_B.flml
+rm MMS_B.flml.cp
+sed s/999.9/0.00125/ MMS_C.flml.cp > MMS_C.flml
+rm MMS_C.flml.cp

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.sage'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.sage	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.sage	2013-02-08 22:48:22 +0000
@@ -0,0 +1,156 @@
+y = var('y')
+
+def function(phi_0, phi_x, phi_y, phi_xy, 
+             f_sin_x, f_cos_x, f_sin_y, f_cos_y, f_sin_xy, f_cos_xy, 
+             alpha_x, alpha_y, alpha_xy):
+    
+    f_0 = phi_0 
+    f_x = phi_x*(f_sin_x*sin(alpha_x*x) + f_cos_x*cos(alpha_x*x)) 
+    f_y = phi_y*(f_sin_y*sin(alpha_y*y) + f_cos_y*cos(alpha_y*y)) 
+    f_xy = phi_xy*(f_sin_xy*sin(alpha_xy*x*y/pi) + f_cos_xy*cos(alpha_xy*x*y/pi)) 
+    f = f_0 + f_x + f_y + f_xy
+    return f
+
+ie1 = 0.5*(cos(x + y) + 1.5)
+ie2 = 2.5*(sin(x*x) + 1.5)
+#ie2 = 3.0*y
+
+#u1 = 0.25*(sin(2*x+(x**2)*(y**2))+0.5)
+#v1 = 0.1*(cos(x**2+y**2)+0.5)
+
+#u2 = 1.0*(sin(x**2+y**2)+0.25)
+#v2 = 1.15*(cos(x*y)+0.75)
+
+u1 = 1.0*(sin(x**2+y**2)+0.5)
+v1 = 0.1*(cos(x**2+y**2)+0.5)
+u2 = cos(x**2+y**2)+2.5*x
+v2 = 0.5*x*y 
+
+rho1 = 0.5*(sin(x*x + y*y) + 1.5)
+rho2 = 3.0
+
+ie1 = 1.25*x*y + cos(x + y)
+ie2 = sin(x*y)
+
+
+gamma = 1.4
+rho0 = 0.5
+csq = 100.0
+#p = csq*(rho1 - rho0)
+p = (gamma-1.0)*ie1*rho1
+             
+#vfrac1 = 0.5*cos(x)*cos(x) + 0.1
+vfrac1 = 0.8
+vfrac2 = 1.0 - vfrac1 
+            
+nu1 = 0.7
+nu2 = 0.7
+
+g_x = 0.707106781
+g_y = 0.707106781
+
+tau_xx1 = nu1*diff(u1,x)            
+tau_yy1 = nu1*diff(v1,y)
+tau_xy1 = nu1*diff(u1,y)
+tau_yx1 = nu1*diff(v1,x)  
+
+tau_xx2 = nu2*diff(u2,x)            
+tau_yy2 = nu2*diff(v2,y)
+tau_xy2 = nu2*diff(u2,y)
+tau_yx2 = nu2*diff(v2,x)  
+
+Su1 = vfrac1*rho1*u1*diff(u1,x) + vfrac1*rho1*v1*diff(u1,y) - diff(vfrac1*tau_xx1, x) - diff(vfrac1*tau_xy1, y) - vfrac1*g_x*rho1 + vfrac1*diff(p,x)  
+Sv1 = vfrac1*rho1*u1*diff(v1,x) + vfrac1*rho1*v1*diff(v1,y) - diff(vfrac1*tau_yy1, y) - diff(vfrac1*tau_yx1, x) - vfrac1*g_y*rho1 + vfrac1*diff(p,y)  
+
+Su2 = vfrac2*rho2*u2*diff(u2,x) + vfrac2*rho2*v2*diff(u2,y) - diff(vfrac2*tau_xx2, x) - diff(vfrac2*tau_xy2, y) - vfrac2*g_x*rho2 + vfrac2*diff(p,x)  
+Sv2 = vfrac2*rho2*u2*diff(v2,x) + vfrac2*rho2*v2*diff(v2,y) - diff(vfrac2*tau_yy2, y) - diff(vfrac2*tau_yx2, x) - vfrac2*g_y*rho2 + vfrac2*diff(p,y) 
+
+Srho1 = diff(vfrac1*u1*rho1,x) + diff(vfrac1*v1*rho1,y) + rho1*diff(u2*vfrac2, x) + rho1*diff(v2*vfrac2, y)
+
+# Heat transfer term for the internal energy equations
+k = 0.5
+Cv1 = 500.0
+Cv2 = 500.0
+d = 0.1
+
+Pr = Cv1*gamma*nu1/k
+Re = rho1*d*(sqrt((u1 - u2)**2 + (v1 - v2)**2))/nu1
+Nu = (7.0 - 10.0*vfrac1 + 5.0*vfrac1**2)*(1.0 + 0.7*(Re**0.2)*(Pr**(1.0/3.0))) + (1.33 - 2.4*vfrac1 + 1.2*vfrac1**2)*(Re**0.7)*(Pr**(1.0/3.0))
+Q = (6.0*k*vfrac2*Nu/(d**2))*(ie2/Cv2 - ie1/Cv1)
+
+heat_flux_x1 = (k/Cv1)*vfrac1*diff(ie1,x)
+heat_flux_y1 = (k/Cv1)*vfrac1*diff(ie1,y)
+heat_flux_x2 = (k/Cv2)*vfrac2*diff(ie2,x)
+heat_flux_y2 = (k/Cv2)*vfrac2*diff(ie2,y)
+
+Sie1 = vfrac1*rho1*u1*diff(ie1, x) + vfrac1*rho1*v1*diff(ie1, y) - diff(heat_flux_x1, x) - diff(heat_flux_y1, y) + vfrac1*p*diff(u1, x) + vfrac1*p*diff(v1, y) #- Q
+Sie2 = vfrac2*rho2*u2*diff(ie2, x) + vfrac2*rho2*v2*diff(ie2, y) - diff(heat_flux_x2, x) - diff(heat_flux_y2, y) #+ Q
+
+print 'from math import sin, cos, tanh, pi, sqrt'
+print ''
+print 'def u1(X):'
+print '    return', str(u1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v1(X):'
+print '    return', str(v1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def p(X):'
+print '    return', str(p).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho1(X):'
+print '    return', str(rho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie1(X):'
+print '    return', str(ie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac1(X):'
+print '    return', str(vfrac1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u1(X):'
+print '    return', str(Su1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v1(X):'
+print '    return', str(Sv1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_rho1(X):'
+print '    return', str(Srho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie1(X):'
+print '    return', str(Sie1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity1(X):'
+print '   return [u1(X), v1(X)]'
+print ''
+print 'def forcing_velocity1(X):'
+print '   return [forcing_u1(X), forcing_v1(X)]'
+print ''
+print 'def u2(X):'
+print '    return', str(u2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def v2(X):'
+print '    return', str(v2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def rho2(X):'
+print '    return', str(rho2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def ie2(X):'
+print '    return', str(ie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def vfrac2(X):'
+print '    return', str(vfrac2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''  
+print 'def forcing_u2(X):'
+print '    return', str(Su2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_v2(X):'
+print '    return', str(Sv2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def forcing_ie2(X):'
+print '    return', str(Sie2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
+print ''
+print 'def velocity2(X):'
+print '   return [u2(X), v2(X)]'
+print ''
+print 'def forcing_velocity2(X):'
+print '   return [forcing_u2(X), forcing_v2(X)]'
+print ''

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.xml'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.xml	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv.xml	2013-02-08 22:48:22 +0000
@@ -0,0 +1,142 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+<testproblem>
+  <name>mphase_mms_p2p1_compressible_ie_p1cv</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  <problem_definition length="special" nprocs="1">
+    <command_line>fluidity -v3 MMS_A.flml &gt; MMS_A.log; fluidity -v3 MMS_B.flml &gt; MMS_B.log; fluidity -v3 MMS_C.flml &gt; MMS_C.log; fluidity -v3 MMS_D.flml &gt; MMS_D.log</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="convergence" language="python">
+from fluidity_tools import stat_parser as stat
+from math import log
+import numpy as np
+
+meshes = [['A','B']]#, ['B','C']]
+
+convergence = np.ones(5) * 1e10
+
+print ''
+print 'ORDER OF CONVERGENCE'
+print '-------------------------------------------'
+
+print 'VelocityError:'
+print '-------------------------------------------'
+
+for i, mesh in enumerate(meshes):
+
+    a_error_x = stat("MMS_"+str(mesh[0])+".stat")["Phase1"]["VelocityError%1"]["l2norm"][-1]
+    b_error_x = stat("MMS_"+str(mesh[1])+".stat")["Phase1"]["VelocityError%1"]["l2norm"][-1]
+    a_error_y = stat("MMS_"+str(mesh[0])+".stat")["Phase1"]["VelocityError%2"]["l2norm"][-1]
+    b_error_y = stat("MMS_"+str(mesh[1])+".stat")["Phase1"]["VelocityError%2"]["l2norm"][-1]
+
+    ratio_x = a_error_x / b_error_x
+    ratio_y = a_error_y / b_error_y
+
+    print mesh[0] + '->' + mesh[1] + ': ', [log(ratio_x, 2), log(ratio_y, 2)]
+
+    convergence[0] = min(log(ratio_x, 2), log(ratio_y, 2), convergence[0])
+
+print '-------------------------------------------'
+
+fields = [
+          'PressureError',
+          'InternalEnergyError',
+          'DensityError',
+          ]
+
+for i, field in enumerate(fields):
+    print field
+    print '-------------------------------------------'
+
+    for j, mesh in enumerate(meshes):
+        
+        a_error = stat("MMS_"+str(mesh[0])+".stat")["Phase1"][field]["l2norm"][-1]
+        b_error = stat("MMS_"+str(mesh[1])+".stat")["Phase1"][field]["l2norm"][-1]
+        
+        ratio = a_error / b_error
+        
+        print mesh[0] + '->' + mesh[1] + ': ', log(ratio, 2)
+
+        convergence[i+1] = min(log(ratio, 2), convergence[i+1])
+    
+    print '-------------------------------------------'
+
+print ''
+
+
+fields = [
+          'ScalarAbsoluteDifference',
+          ]
+
+for i, field in enumerate(fields):
+    print field
+    print '-------------------------------------------'
+
+    for j, mesh in enumerate(meshes):
+        
+        a_error = stat("MMS_"+str(mesh[0])+".stat")["Phase2"][field]["l2norm"][-1]
+        b_error = stat("MMS_"+str(mesh[1])+".stat")["Phase2"][field]["l2norm"][-1]
+        
+        ratio = a_error / b_error
+        
+        print mesh[0] + '->' + mesh[1] + ': ', log(ratio, 2)
+
+        convergence[i+4] = min(log(ratio, 2), convergence[i+1])
+    
+    print '-------------------------------------------'
+
+print ''
+
+    </variable>
+    <variable name="a_finish_time" language="python">
+from fluidity_tools import stat_parser as stat
+a_finish_time = stat("MMS_A.stat")["ElapsedTime"]["value"][-1]
+    </variable>
+    <variable name="b_finish_time" language="python">
+from fluidity_tools import stat_parser as stat
+b_finish_time = stat("MMS_B.stat")["ElapsedTime"]["value"][-1]
+    </variable>
+    <variable name="c_finish_time" language="python">
+from fluidity_tools import stat_parser as stat
+c_finish_time = stat("MMS_C.stat")["ElapsedTime"]["value"][-1]
+    </variable>
+    <variable name="d_finish_time" language="python">
+from fluidity_tools import stat_parser as stat
+d_finish_time = stat("MMS_D.stat")["ElapsedTime"]["value"][-1]
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="velocity convergence order > 2.8" language="python">
+assert(convergence[0] &gt; 2.8)
+    </test>
+    <test name="pressure convergence order > 1.9" language="python">
+assert(convergence[1] &gt; 1.9)
+    </test>
+    <test name="ie (phase 1) convergence order > 1.9" language="python">
+assert(convergence[2] &gt; 1.9)
+    </test>
+    <test name="density (phase 1) convergence order > 1.9" language="python">
+assert(convergence[3] &gt; 1.9)
+    </test>
+    <test name="ie (phase 2) convergence order > 1.9" language="python">
+assert(convergence[4] &gt; 1.9)
+    </test>
+    <test name="checking A finished in less than 100 seconds" language="python">
+assert(a_finish_time-100.0 &lt; 1.E-10)
+    </test>
+    <test name="checking B finished in less than 100 seconds" language="python">
+assert(b_finish_time-100.0 &lt; 1.E-10)
+    </test>
+    <test name="checking C finished in less than 100 seconds" language="python">
+assert(c_finish_time-100.0 &lt; 1.E-10)
+    </test>
+    <test name="checking D finished in less than 100 seconds" language="python">
+assert(d_finish_time-100.0 &lt; 1.E-10)
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
+

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv_tools.py'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv_tools.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/mphase_mms_p2p1_compressible_ie_p1cv_tools.py	2013-02-08 22:48:22 +0000
@@ -0,0 +1,68 @@
+from math import sin, cos, tanh, pi, sqrt
+
+def u1(X):
+    return sin(X[0]**2 + X[1]**2) + 0.500
+
+def v1(X):
+    return 0.100*cos(X[0]**2 + X[1]**2) + 0.0500
+
+def p(X):
+    return (0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))
+
+def rho1(X):
+    return 0.500*sin(X[0]**2 + X[1]**2) + 0.750
+
+def ie1(X):
+    return 1.25*X[0]*X[1] + cos(X[0] + X[1])
+
+def vfrac1(X):
+    return 0.800
+
+def forcing_u1(X):
+    return 2*(0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*X[1]*cos(X[0]**2 + X[1]**2) + 2*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*X[0]*cos(X[0]**2 + X[1]**2) + 0.800*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[0]*cos(X[0]**2 + X[1]**2) + 2.24*X[0]**2*sin(X[0]**2 + X[1]**2) + 2.24*X[1]**2*sin(X[0]**2 + X[1]**2) + 0.800*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[1] - 0.400*sin(X[0] + X[1])) - 0.282842712400000*sin(X[0]**2 + X[1]**2) - 2.24*cos(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_v1(X):
+    return -0.200*(0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*X[1]*sin(X[0]**2 + X[1]**2) - 0.200*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*X[0]*sin(X[0]**2 + X[1]**2) + 0.800*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[1]*cos(X[0]**2 + X[1]**2) + 0.224*X[0]**2*cos(X[0]**2 + X[1]**2) + 0.224*X[1]**2*cos(X[0]**2 + X[1]**2) + 0.800*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[0] - 0.400*sin(X[0] + X[1])) - 0.0588427124000001*sin(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_rho1(X):
+    return (0.0800*cos(X[0]**2 + X[1]**2) + 0.0400)*X[1]*cos(X[0]**2 + X[1]**2) + 1.60*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[0]*cos(X[0]**2 + X[1]**2) - 0.160*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[1]*sin(X[0]**2 + X[1]**2) + (0.800*sin(X[0]**2 + X[1]**2) + 0.400)*X[0]*cos(X[0]**2 + X[1]**2) + (0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(-0.400*X[0]*sin(X[0]**2 + X[1]**2) + 0.500) + 0.100*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*X[0]
+
+def forcing_ie1(X):
+    return 1.60*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[0]*cos(X[0]**2 + X[1]**2) - 0.160*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[1]*sin(X[0]**2 + X[1]**2) + (0.100*cos(X[0]**2 + X[1]**2) + 0.0500)*(0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(1.25*X[0] - sin(X[0] + X[1])) + (0.400*sin(X[0]**2 + X[1]**2) + 0.600)*(sin(X[0]**2 + X[1]**2) + 0.500)*(1.25*X[1] - sin(X[0] + X[1])) + 0.00160*cos(X[0] + X[1])
+
+def velocity1(X):
+   return [u1(X), v1(X)]
+
+def forcing_velocity1(X):
+   return [forcing_u1(X), forcing_v1(X)]
+
+def u2(X):
+    return 2.50*X[0] + cos(X[0]**2 + X[1]**2)
+
+def v2(X):
+    return 0.500*X[0]*X[1]
+
+def rho2(X):
+    return 3.00
+
+def ie2(X):
+    return sin(X[0]*X[1])
+
+def vfrac2(X):
+    return 0.200
+
+def forcing_u2(X):
+    return -0.600*X[0]*X[1]**2*sin(X[0]**2 + X[1]**2) + 0.200*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[0]*cos(X[0]**2 + X[1]**2) + 0.560*X[0]**2*cos(X[0]**2 + X[1]**2) + 0.560*X[1]**2*cos(X[0]**2 + X[1]**2) + 0.200*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[1] - 0.400*sin(X[0] + X[1])) - (2.00*X[0]*sin(X[0]**2 + X[1]**2) - 2.50)*(1.50*X[0] + 0.600*cos(X[0]**2 + X[1]**2)) + 0.560*sin(X[0]**2 + X[1]**2) - 0.424264068600000
+
+def forcing_v2(X):
+    return 0.200*(0.500*X[0]*X[1] + 0.400*cos(X[0] + X[1]))*X[1]*cos(X[0]**2 + X[1]**2) + 0.150*X[0]**2*X[1] + 0.200*(0.500*sin(X[0]**2 + X[1]**2) + 0.750)*(0.500*X[0] - 0.400*sin(X[0] + X[1])) + 0.500*(1.50*X[0] + 0.600*cos(X[0]**2 + X[1]**2))*X[1] - 0.424264068600000
+
+def forcing_ie2(X):
+    return 0.300*X[0]**2*X[1]*cos(X[0]*X[1]) + (1.50*X[0] + 0.600*cos(X[0]**2 + X[1]**2))*X[1]*cos(X[0]*X[1]) + 0.000200*X[0]**2*sin(X[0]*X[1]) + 0.000200*X[1]**2*sin(X[0]*X[1])
+
+def velocity2(X):
+   return [u2(X), v2(X)]
+
+def forcing_velocity2(X):
+   return [forcing_u2(X), forcing_v2(X)]
+

=== added directory 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src'
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.geo'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.1};
+Point(2) = {0.7,0.2,0,0.1};
+Point(3) = {0.7,0.8,0,0.1};
+Point(4) = {-0.1,0.8,0,0.1};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_A.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.geo'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.05};
+Point(2) = {0.7,0.2,0,0.05};
+Point(3) = {0.7,0.8,0,0.05};
+Point(4) = {-0.1,0.8,0,0.05};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_B.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.geo'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,15 @@
+Point(1) = {-0.1,0.2,0,0.025};
+Point(2) = {0.7,0.2,0,0.025};
+Point(3) = {0.7,0.8,0,0.025};
+Point(4) = {-0.1,0.8,0,0.025};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Physical Line(7) = {3};
+Physical Line(8) = {2};
+Physical Line(9) = {1};
+Physical Line(10) = {4};
+Physical Surface(12) = {6};

=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.msh'
Binary files tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.msh	1970-01-01 00:00:00 +0000 and tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_C.msh	2013-02-08 22:48:22 +0000 differ
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_X.geo'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_X.geo	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv/src/MMS_X.geo	2013-02-08 22:48:22 +0000
@@ -0,0 +1,18 @@
+pi = 3.141592653589793;
+n = XX;
+
+Point(1) = {0.0, 0.0, 0.0, 1};
+
+Extrude {pi, 0.0, 0.0} {
+  Point{1}; Layers{n};
+}
+
+Extrude {0.0, pi, 0.0} {
+  Line{1}; Layers{n};
+}
+
+Physical Line(7) = {1};
+Physical Line(8) = {2};
+Physical Line(9) = {3};
+Physical Line(10) = {4}; 
+Physical Surface(6) = {5};

=== added directory 'tests/mphase_mms_p2p1_compressible_ie_p1cv_heat_transfer'
=== added file 'tests/mphase_mms_p2p1_compressible_ie_p1cv_heat_transfer/MMS_X.flml'
--- tests/mphase_mms_p2p1_compressible_ie_p1cv_heat_transfer/MMS_X.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_mms_p2p1_compressible_ie_p1cv_heat_transfer/MMS_X.flml	2013-02-08 22:48:22 +0000
@@ -0,0 +1,1216 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">MMS_X</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="src/MMS_X">
+        <format name="gmsh"/>
+        <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="ErrorMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">4</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>continuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <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>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">8</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1000.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">999.9</real_value>
+      <comment>gives a max cfl number of approximately 0.1</comment>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1000.0</real_value>
+      <comment>10.0</comment>
+    </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>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1.E-7</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">1.0</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="2" dim1="dim" rank="1">0.707106781 0.707106781</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Fluid">
+    <equation_of_state>
+      <compressible>
+        <stiffened_gas>
+          <ratio_specific_heats>
+            <real_value rank="0">1.4</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/>
+          </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">10000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <output>
+          <include_previous_time_step/>
+        </output>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>