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

[Christian Jacobs]

Proposing merge of compressible-multiphase branch revisions r3955 to r4050 (inclusive) into trunk.

Summary of changes:
===================

- Added multiphase support for the compressible projection method (currently only for CG Pressure).

- Small changes to the logic in Momentum_Equation.F90. Replaced the use_compressible_projection logical with compressible_eos which is .true. if option_count("/material_phase/equation_of_state/compressible") > 0, and .false. otherwise.

- Added multiphase support for the InternalEnergy equation. Note that the whole state array is now passed to solve_field_equation_cg.

- Implemented the heat transfer term by Gunn (1978).

- Added three P2-P1 multiphase MMS tests for model verification: one without InternalEnergy fields, one with InternalEnergy fields, and one with InternalEnergy fields and heat transfer. The tests will become longtests after merging into the trunk.

- The manual now has documentation on the compressible multiphase model, including the multiphase version of the InternalEnergy equation.

- Added a diagnostic field called CompressibleContinuity which computes the residual of the compressible multiphase continuity equation. This is similar to the SumVelocityDivergence field for incompressible multiphase flow which computes \sum{div(vfrac*u)}.

- Added two gas-solid shock tube tests to help to validate the model: mphase_rogue_shock_tube_dense_bed_glass and mphase_rogue_shock_tube_dense_bed_nylon.

- Added another gas-solid shock tube test (mphase_dusty_gas_shock_tube) which uses the setup by Miura & Glass (1982).

- Implemented two new correlations for fluid-particle drag. The one by Wen & Yu (1966) can support higher Reynolds numbers (compared to the existing Stokes' law correlation). The one by Ergun (1952) can support dense flows where the volume fraction of the particle phase is greater than 0.2, like in the dense bed shock tube tests (and in the pyroclastic flow tests which will hopefully come later).

- Changed the structure of the schema. Users now explicitly select which interaction terms they wish to include, as well as the drag correlation, under "/multiphase_interaction" in Diamond. These options are stored in the new files multiphase_interaction.rnc and multiphase_interaction.rng.

- Added a test for the Wen & Yu (1966) drag correlation. This is similar to mphase_stokes_law.

- Added a test to check if inlet velocity boundary conditions work with the incompressible multiphase flow model.

- Added a test for the flux boundary condition.

[Stephan Kramer (s-kramer)]

Hey there. I've only looked at Advection_Diffusion_CG.F90 yet. Looking great so far, the maths of the changes in the i.e. equation seem to make sense to me. Here it goes:

* Currently the multiphase logical is only set within the
equation_type==FIELD_EQUATION_INTERNALENERGY case. Can we either:
- always set multiphase to be true or false, or
- rename it to mean multiphase_internal_energy_equation
As it is now suggested that multiphase can be used independently of
the value of equation_type.

* Just being picky here: can we get rid of the nullify(vfrac) in line 506. It seems
to suggest that vfrac is used later on, which it's not.

* Advection_Diffusion_CG.F90, line 517:

Why does the behaviour change if we're not multiphase? Do you maybe mean:

if (.not. multiphase .or. have_option(....compressible) then

* Advection_Diffusion_CG.F90, line 794 having an allocate within an element
loop is not ideal. It's not efficient, and it's also not thread-safe. I can
see the problem with following the usual approach (nvfrac might
not be defined), but I think it's better to do some kludge with allocating
a dummy nvfrac or something rather.

* Advection_Diffusion_CG.F90, lines 1492 and 1502: I don't think you intended to
add (istate) there.

* Not too excited about having to pass in an array of states
to advection_diffusion_cg now tbh. One way to keep the damage limited would be
to move the call to add_heat_transfer up one level inside
solve_field_equation_cg. And keep everything else just limited to a
single state, i.e. not passing istate all the way down and changing state ->
state(istate) *everywhere*.

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

4051. By Christian Jacobs on 2012-07-29

Improvements to Advection_Diffusion_CG.F90 based on Stephan's feedback.

4052. By Christian Jacobs on 2012-07-29

Moving the call to add_heat_transfer up one level.

4053. By Christian Jacobs on 2012-07-29

Merging trunk revisions (from r4011 to r4019 inclusive) into branch.

[Christian Jacobs (ctjacobs)]

Thanks for the feedback Stephan - improvements have now been made to Advection_Diffusion_CG.F90. I have decided to always set the 'multiphase' logical rather than rename it to something involving a specific equation type, since it might be used later on when adding multiphase support to the k-epsilon model.

[Stephan Kramer (s-kramer)]

Some more comments trickling through:

* assemble/Compressible_Projection.F90, line 463: dnvfrac_t is not used

* assemble/Compressible_Projection.F90, line 650,660: Why are we setting this
multiphase variable there, it's not used afaics and a bit confusing...

* assemble/Diagnostic_Fields_Matrices.F90, line 372: you query the timestep now
but it doesn't seem to be used? Also in the same routine you've added vfrac
and temp_vfrac that aren't used.

* calculate_compressible_continuity() looks good to me. Is it tested? Would a
better maybe be compressible_continuity_residual? As it's not directly clear
what "compressible_continuity" as a diagnostic quantity is...

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

4054. By Christian Jacobs on 2012-08-03

Some more improvements based on Stephan's feedback.

4055. By Christian Jacobs on 2012-08-03

Added a test for the CompressibleContinuityResidual diagnostic field.

4056. By Christian Jacobs on 2012-08-03

More tests for the CompressibleContinuityResidual diagnostic field.

[Christian Jacobs (ctjacobs)]

Thanks again - improvements made and the (renamed) CompressibleContinuityResidual diagnostic field is now tested within the MMS tests by checking that it converges to the analytical solution at the expected order.

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

4057. By Christian Jacobs on 2012-08-03

Merging trunk revisions (from r4020 to r4022 inclusive) into branch.

[Stephan Kramer (s-kramer)]

Some final, small remarks:

* assemble/Divergence_Matrix.F90: is there a reason you're not using state%option_path? I prefer that as it makes future redesigns of the optionstree easier.

* assemble/Momentum_Equation.F90, line 471: does this comment still apply? if not, how is this fixed?

* assemble/Momentum_Equation.F90, line 1049: why do we all of a sudden extract density from state there? AFAIK Density is not a required field, for instance for Boussinesq. Similar for line 1585: does this not break if we don't have
a Density, but do have a prognostic Velocity and Pressure?

* assemble/Multiphase.F90, add_fluid_particle_drag_element(): doing string comparisons at the element level is not ideal. In particular with OPTION_PATH_LEN (8192 or so?) length string, the compiler might do something
sensible, but I'm not sure and I don't like second-guessing. A solution would be to move the select case up one level, and set an integer with some "enum"-like parameters DRAG_CORRELATION_TYPE_ERGUN=3, etc.

* assemble/Multiphase.F90, line 625 and 632: again here, why not use state%option_path?

* assemble/Multiphase.F90, add_heat_transfer(): looks great and I'm sure is 100% correct - good to see it's well tested though!

Great stuff, well documented, lots of tests. Next time, could we maybe merge this in smaller bits?

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

4058. By Christian Jacobs on 2012-08-08

Using state%option_path instead.

4059. By Christian Jacobs on 2012-08-08

- Don't try to extract the Density field here, as it's not necessary.
- Another change to use state%option_path instead.

4060. By Christian Jacobs on 2012-08-08

Call FLExit if use_theta_pg or use_theta_divergence are true when running a compressible multiphase flow simulation.

4061. By Christian Jacobs on 2012-08-08

Introduced integer parameters to represent the different drag correlations.

4062. By Christian Jacobs on 2012-08-08

Merging trunk revision r4023 into branch.

[Stephan Kramer (s-kramer)]

Awesome!

review: Approve

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

4063. By Christian Jacobs on 2012-08-08

Corrected a typo introduced when changing to state%option_path.

4064. By Christian Jacobs on 2012-08-08

- Corrected a few typos in the XML files.
- Set the mphase_rogue_shock_tube_dense_bed_glass test's length to "medium".

4065. By Christian Jacobs on 2012-08-08

Line continuation fix. The code now compiles on the Intel compiler.

4066. By Christian Jacobs on 2012-08-09

Fixed an assert statement. Only noticed this when debugging was enabled.

4067. By Christian Jacobs on 2012-08-09

Updates from "make makefiles".

[Christian Jacobs (ctjacobs)]

The build queue now has a green light on buildbot. Thanks for the review.

=== modified file 'assemble/Advection_Diffusion_CG.F90'
--- assemble/Advection_Diffusion_CG.F90	2012-04-16 16:41:59 +0000
+++ assemble/Advection_Diffusion_CG.F90	2012-08-09 00:46:20 +0000
@@ -47,6 +47,7 @@
   use upwind_stabilisation
   use sparsity_patterns_meshes
   use porous_media
+  use multiphase_module
   use sparse_tools_petsc
   use colouring
 #ifdef _OPENMP
@@ -117,16 +118,19 @@
   logical :: move_mesh
   ! Include porosity?
   logical :: include_porosity
+  ! Are we running a multiphase flow simulation?
+  logical :: multiphase
 
 contains
 
-  subroutine solve_field_equation_cg(field_name, state, dt, velocity_name, iterations_taken)
+  subroutine solve_field_equation_cg(field_name, state, istate, dt, velocity_name, iterations_taken)
     !!< Construct and solve the advection-diffusion equation for the given
     !!< field using a continuous Galerkin discretisation. Based on
     !!< Advection_Diffusion_DG and Momentum_CG.
     
     character(len = *), intent(in) :: field_name
-    type(state_type), intent(inout) :: state
+    type(state_type), dimension(:), intent(inout) :: state
+    integer, intent(in) :: istate
     real, intent(in) :: dt
     character(len = *), optional, intent(in) :: velocity_name
     integer, intent(out), optional :: iterations_taken
@@ -138,16 +142,22 @@
     ewrite(1, *) "In solve_field_equation_cg"
     
     ewrite(2, *) "Solving advection-diffusion equation for field " // &
-      & trim(field_name) // " in state " // trim(state%name)
+      & trim(field_name) // " in state " // trim(state(istate)%name)
 
-    call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state)
+    call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state(istate))
     
     call profiler_tic(t, "assembly")
-    call assemble_advection_diffusion_cg(t, matrix, rhs, state, dt, velocity_name = velocity_name)    
+    call assemble_advection_diffusion_cg(t, matrix, rhs, state(istate), dt, velocity_name = velocity_name)    
+    
+    ! 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)
+    end if
     call profiler_toc(t, "assembly")
-
+    
     call profiler_tic(t, "solve_total")
-    call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state, &
+    call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state(istate), &
                                       iterations_taken = iterations_taken)
     call profiler_toc(t, "solve_total")
 
@@ -234,6 +244,11 @@
     type(scalar_field), pointer :: density, olddensity
     character(len = FIELD_NAME_LEN) :: density_name
     type(scalar_field), pointer :: pressure
+    
+    ! Volume fraction fields for multiphase flow simulation
+    type(scalar_field), pointer :: vfrac
+    type(scalar_field) :: nvfrac ! Non-linear version
+
     ! Porosity field
     type(scalar_field) :: porosity_theta
         
@@ -373,6 +388,7 @@
        call allocate(porosity_theta, t%mesh, field_type=FIELD_TYPE_CONSTANT)
        call set(porosity_theta, 1.0)
     end if
+
     
     ! Step 2: Pull options out of the options tree
     
@@ -462,6 +478,21 @@
       stabilisation_scheme = STABILISATION_NONE
     end if
     
+    ! PhaseVolumeFraction for multiphase flow simulations
+    if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+       multiphase = .true.
+       vfrac => extract_scalar_field(state, "PhaseVolumeFraction")
+       call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
+       call zero(nvfrac)
+       call get_nonlinear_volume_fraction(state, nvfrac)
+      
+       ewrite_minmax(nvfrac)
+    else
+       multiphase = .false.
+       call allocate(nvfrac, t%mesh, "DummyNonlinearPhaseVolumeFraction", field_type=FIELD_TYPE_CONSTANT)
+       call set(nvfrac, 1.0)
+    end if
+      
     call allocate(dummydensity, t%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
     call set(dummydensity, 1.0)
     ! find out equation type and hence if density is needed or not
@@ -479,19 +510,28 @@
       if(move_mesh) then
         FLExit("Haven't implemented a moving mesh energy equation yet.")
       end if
+      
+      ! Get old and current densities
       call get_option(trim(t%option_path)//'/prognostic/equation[0]/density[0]/name', &
                       density_name)
       density=>extract_scalar_field(state, trim(density_name))
       ewrite_minmax(density)
-      
       olddensity=>extract_scalar_field(state, "Old"//trim(density_name))
       ewrite_minmax(olddensity)
       
-      call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", &
-                      density_theta)
+      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)
+      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
+      end if
                       
       pressure=>extract_scalar_field(state, "Pressure")
       ewrite_minmax(pressure)
+
     case default
       FLExit("Unknown field equation type for cg advection diffusion.")
     end select
@@ -531,7 +571,7 @@
               positions, old_positions, new_positions, &
               velocity, grid_velocity, &
               source, absorption, diffusivity, &
-              density, olddensity, pressure, porosity_theta, &
+              density, olddensity, pressure, porosity_theta, nvfrac, &
               supg_element(thread_num+1))
       end do element_loop
       !$OMP END DO
@@ -545,6 +585,7 @@
     ! which must be included before dirichlet BC's.
     if (add_src_directly_to_rhs) call addto(rhs, source)
     
+    
     ! Step 4: Boundary conditions
     
     if( &
@@ -571,7 +612,7 @@
         call assemble_advection_diffusion_face_cg(i, t_bc_types(i), t, t_bc, t_bc_2, &
                                                   matrix, rhs, &
                                                   positions, velocity, grid_velocity, &
-                                                  density, olddensity)
+                                                  density, olddensity, nvfrac)
       end do
     
       call deallocate(t_bc)
@@ -586,6 +627,7 @@
     ewrite_minmax(rhs)
     
     call deallocate(velocity)
+    call deallocate(nvfrac)
     call deallocate(dummydensity)
     if (stabilisation_scheme == STABILISATION_SUPG) then
        do i = 1, num_threads
@@ -645,7 +687,7 @@
                                       positions, old_positions, new_positions, &
                                       velocity, grid_velocity, &
                                       source, absorption, diffusivity, &
-                                      density, olddensity, pressure, porosity_theta, supg_shape)
+                                      density, olddensity, pressure, porosity_theta, nvfrac, supg_shape)
     integer, intent(in) :: ele
     type(scalar_field), intent(in) :: t
     type(csr_matrix), intent(inout) :: matrix
@@ -661,6 +703,7 @@
     type(scalar_field), intent(in) :: olddensity
     type(scalar_field), intent(in) :: pressure
     type(scalar_field), intent(in) :: porosity_theta
+    type(scalar_field), intent(in) :: nvfrac
     type(element_type), intent(inout) :: supg_shape
 
     integer, dimension(:), pointer :: element_nodes
@@ -669,6 +712,9 @@
     real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)) :: drho_t
     real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t 
     real, dimension(ele_loc(positions, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: dug_t 
+    ! Derivative of shape function for nvfrac field
+    real, dimension(ele_loc(nvfrac, ele), ele_ngi(nvfrac, ele), mesh_dim(nvfrac)) :: dnvfrac_t
+    
     real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)) :: j_mat 
     type(element_type) :: test_function
     type(element_type), pointer :: t_shape
@@ -741,6 +787,16 @@
       end if
     end if
     
+    if(have_advection .and. multiphase .and. (equation_type==FIELD_EQUATION_INTERNALENERGY)) then
+      ! If the field and nvfrac meshes are different, then we need to
+      ! compute the derivatives of the nvfrac shape functions.    
+      if(ele_shape(nvfrac, ele) == t_shape) then
+         dnvfrac_t = dt_t
+      else
+         call transform_to_physical(positions, ele, ele_shape(nvfrac, ele), dshape=dnvfrac_t)
+      end if
+    end if
+
     ! Step 2: Set up test function
     
     select case(stabilisation_scheme)
@@ -763,13 +819,13 @@
     ! Step 3: Assemble contributions
     
     ! Mass
-    if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
+    if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
     
     ! Advection
     if(have_advection) call add_advection_element_cg(ele, test_function, t, &
                                         velocity, grid_velocity, diffusivity, &
-                                        density, olddensity, &
-                                        dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto)
+                                        density, olddensity, nvfrac, &
+                                        dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto)
         
     ! Absorption
     if(have_absorption) call add_absorption_element_cg(ele, test_function, t, absorption, detwei, matrix_addto, rhs_addto)
@@ -784,8 +840,8 @@
     
     ! Pressure
     if(equation_type==FIELD_EQUATION_INTERNALENERGY) call add_pressurediv_element_cg(ele, test_function, t, &
-                                                                                  velocity, pressure, &
-                                                                                  du_t, detwei, rhs_addto)
+                                                                                  velocity, pressure, nvfrac, &
+                                                                                  du_t, dnvfrac_t, detwei, rhs_addto)
     
     ! Step 4: Insertion
             
@@ -795,12 +851,13 @@
 
   end subroutine assemble_advection_diffusion_element_cg
   
-  subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
+  subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
     integer, intent(in) :: ele
     type(element_type), intent(in) :: test_function
     type(scalar_field), intent(in) :: t
     type(scalar_field), intent(in) :: density, olddensity
     type(scalar_field), intent(in) :: porosity_theta    
+    type(scalar_field), intent(in) :: nvfrac    
     real, dimension(ele_ngi(t, ele)), intent(in) :: detwei, detwei_old, detwei_new
     real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto
     real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto
@@ -827,6 +884,8 @@
       else
         if (include_porosity) then
           mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*porosity_theta_at_quad)        
+        else if(multiphase) then
+          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*ele_val_at_quad(nvfrac, ele))
         else
           mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad)
         end if
@@ -877,8 +936,8 @@
   
   subroutine add_advection_element_cg(ele, test_function, t, &
                                 velocity, grid_velocity, diffusivity, &
-                                density, olddensity, &
-                                dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto)
+                                density, olddensity, nvfrac, &
+                                dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto)
     integer, intent(in) :: ele
     type(element_type), intent(in) :: test_function
     type(scalar_field), intent(in) :: t
@@ -886,10 +945,12 @@
     type(vector_field), pointer :: grid_velocity
     type(tensor_field), intent(in) :: diffusivity
     type(scalar_field), intent(in) :: density, olddensity
+    type(scalar_field), intent(in) :: nvfrac
     real, dimension(ele_loc(t, ele), ele_ngi(t, ele), mesh_dim(t)), intent(in) :: dt_t
     real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t
     real, dimension(:, :, :) :: dug_t
     real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)), intent(in) :: drho_t
+    real, dimension(:, :, :), intent(in) :: dnvfrac_t
     real, dimension(ele_ngi(t, ele)), intent(in) :: detwei
     real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)), intent(in) :: j_mat 
     real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto
@@ -904,6 +965,10 @@
     real, dimension(velocity%dim, ele_ngi(density, ele)) :: densitygrad_at_quad
     real, dimension(ele_ngi(density, ele)) :: udotgradrho_at_quad
     
+    real, dimension(ele_ngi(t, ele)) :: nvfrac_at_quad
+    real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad
+    real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad
+    
     assert(have_advection)
     
     t_shape => ele_shape(t, ele)
@@ -922,6 +987,10 @@
       densitygrad_at_quad = density_theta*ele_grad_at_quad(density, ele, drho_t) &
                            +(1.-density_theta)*ele_grad_at_quad(olddensity, ele, drho_t)
       udotgradrho_at_quad = sum(densitygrad_at_quad*velocity_at_quad, 1)
+      
+      if(multiphase) then
+         nvfrac_at_quad = ele_val_at_quad(nvfrac, ele)
+      end if
     end select
                 
     if(integrate_advection_by_parts) then
@@ -931,12 +1000,26 @@
       !    /                                        /
       select case(equation_type)
       case(FIELD_EQUATION_INTERNALENERGY)
-        advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad)
+        if(multiphase) then      
+           advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad*nvfrac_at_quad)
+        else
+           advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad)
+        end if
+           
         if(abs(1.0 - beta) > epsilon(0.0)) then
           velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
-          advection_mat = advection_mat &
-                    - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
-                                                                      +udotgradrho_at_quad)* detwei)
+          
+          if(multiphase) then
+          
+             nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
+             udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1)
+             
+             advection_mat = advection_mat - (1.0-beta) * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) &
+                                                           + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) &
+                                                           + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) )
+          else
+             advection_mat = advection_mat - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad + udotgradrho_at_quad)*detwei)
+          end if
         end if
       case default
         advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei)
@@ -953,12 +1036,31 @@
       !  /                                 /
       select case(equation_type)
       case(FIELD_EQUATION_INTERNALENERGY)
-        advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad)
+        
+        if(multiphase) then     
+           ! vfrac*rho*nu*grad(internalenergy)
+           advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad*nvfrac_at_quad)
+        else
+           advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad)
+        end if
+        
         if(abs(beta) > epsilon(0.0)) then
           velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
-          advection_mat = advection_mat &
-                    + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
-                                                                +udotgradrho_at_quad)*detwei)
+          
+          if(multiphase) then
+             ! advection_mat + internalenergy*div(vfrac*rho*nu)
+             ! Split up div(vfrac*rho*nu) = vfrac*rho*div(nu) + nu*grad(vfrac*rho) = vfrac*rho*div(nu) + nu*(vfrac*grad(rho) + rho*grad(nvfrac))
+             
+             nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
+             udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1)
+             
+             advection_mat = advection_mat + beta * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) &
+                                                     + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) &
+                                                     + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) )
+          else
+             advection_mat = advection_mat + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
+                                                              +udotgradrho_at_quad)*detwei)
+          end if
         end if
       case default
         advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei)
@@ -1060,26 +1162,39 @@
     
   end subroutine add_diffusivity_element_cg
   
-  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, du_t, detwei, rhs_addto)
+  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, dnvfrac_t, detwei, rhs_addto)
   
     integer, intent(in) :: ele
     type(element_type), intent(in) :: test_function
     type(scalar_field), intent(in) :: t
     type(vector_field), intent(in) :: velocity
     type(scalar_field), intent(in) :: pressure
-    real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t
+    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))
     
-    rhs_addto = rhs_addto - &
-                shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei)
+    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) )
+    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
     
   end subroutine add_pressurediv_element_cg
   
-  subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity)
+  subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity, nvfrac)
     integer, intent(in) :: face
     integer, intent(in) :: bc_type
     type(scalar_field), intent(in) :: t
@@ -1092,6 +1207,7 @@
     type(vector_field), pointer :: grid_velocity
     type(scalar_field), intent(in) :: density
     type(scalar_field), intent(in) :: olddensity
+    type(scalar_field), intent(in) :: nvfrac
     
     integer, dimension(face_loc(t, face)) :: face_nodes
     real, dimension(face_ngi(t, face)) :: detwei
@@ -1125,7 +1241,7 @@
     
     ! Advection
     if(have_advection .and. integrate_advection_by_parts) &
-      call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto)
+      call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto)
     
     ! Diffusivity
     if(have_diffusivity) call add_diffusivity_face_cg(face, bc_type, t, t_bc, t_bc_2, detwei, matrix_addto, rhs_addto)
@@ -1138,7 +1254,7 @@
     
   end subroutine assemble_advection_diffusion_face_cg
   
-  subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto)
+  subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto)
     integer, intent(in) :: face
     integer, intent(in) :: bc_type
     type(scalar_field), intent(in) :: t
@@ -1147,6 +1263,7 @@
     type(vector_field), pointer :: grid_velocity
     type(scalar_field), intent(in) :: density
     type(scalar_field), intent(in) :: olddensity
+    type(scalar_field), intent(in) :: nvfrac
     real, dimension(face_ngi(t, face)), intent(in) :: detwei
     real, dimension(mesh_dim(t), face_ngi(t, face)), intent(in) :: normal
     real, dimension(face_loc(t, face), face_loc(t, face)), intent(inout) :: matrix_addto
@@ -1171,8 +1288,12 @@
     case(FIELD_EQUATION_INTERNALENERGY)
       density_at_quad = density_theta*face_val_at_quad(density, face) &
                        +(1.0-density_theta)*face_val_at_quad(olddensity, face)
-
-      advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad)
+                       
+      if(multiphase) then
+         advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad * face_val_at_quad(nvfrac, face))
+      else
+         advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad)
+      end if
     case default
       
       advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1))
 
=== modified file 'assemble/Compressible_Projection.F90'
--- assemble/Compressible_Projection.F90	2012-01-20 15:54:44 +0000
+++ assemble/Compressible_Projection.F90	2012-08-09 00:46:20 +0000
@@ -39,6 +39,7 @@
   use fefields, only: compute_cv_mass
   use state_fields_module
   use upwind_stabilisation
+  use multiphase_module
   implicit none 
 
   ! Buffer for output messages.
@@ -56,7 +57,10 @@
   integer :: nu_bar_scheme
   real :: nu_bar_scale
 
-contains
+  !! Are we running a multiphase flow simulation?
+  logical :: multiphase
+
+ contains
 
   subroutine assemble_compressible_projection_cv(state, cmc, rhs, dt, theta_pg, theta_divergence, cmcget)
 
@@ -367,11 +371,12 @@
 
   end subroutine assemble_mmat_compressible_projection_cv
   
-  subroutine assemble_compressible_projection_cg(state, cmc, rhs, dt, theta_pg, theta_divergence, cmcget)
+  subroutine assemble_compressible_projection_cg(state, istate, cmc, rhs, dt, theta_pg, theta_divergence, cmcget)
 
     ! inputs:
     ! bucket full of fields
     type(state_type), dimension(:), intent(inout) :: state
+    integer, intent(in) :: istate
 
     type(csr_matrix), intent(inout) :: cmc
     type(scalar_field), intent(inout) :: rhs
@@ -380,17 +385,25 @@
     real, intent(in) :: theta_pg, theta_divergence
     logical, intent(in) :: cmcget
 
-    if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
-    
-      call assemble_1mat_compressible_projection_cg(state(1), cmc, rhs, dt, &
-                                                    theta_pg, theta_divergence, cmcget)
-      
+    if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+       multiphase = .true.
+       call assemble_1mat_compressible_projection_cg(state(istate), cmc, rhs, dt, &
+                                                      theta_pg, theta_divergence, cmcget)
     else
-      
-        FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
-      
-    end if
-    
+       multiphase = .false.
+
+       if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
+      
+          call assemble_1mat_compressible_projection_cg(state(1), cmc, rhs, dt, &
+                                                      theta_pg, theta_divergence, cmcget)
+         
+       else
+         
+          FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
+         
+       end if
+
+   end if    
 
   end subroutine assemble_compressible_projection_cg
 
@@ -441,6 +454,11 @@
     logical :: have_absorption, have_source
     integer :: stat
 
+    !! Multiphase variables
+    ! Volume fraction fields
+    type(scalar_field), pointer :: vfrac
+    type(scalar_field) :: nvfrac
+
     ! =============================================================
     ! Subroutine to construct the matrix CT_m (a.k.a. C1/2/3T).
     ! =============================================================
@@ -470,6 +488,15 @@
   
       velocity=>extract_vector_field(state, "Velocity")
       nonlinearvelocity=>extract_vector_field(state, "NonlinearVelocity") ! maybe this should be updated after the velocity solve?
+
+      ! Get the non-linear PhaseVolumeFraction field if multiphase
+      if(multiphase) then
+         vfrac => extract_scalar_field(state, "PhaseVolumeFraction")
+         call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
+         call zero(nvfrac)
+         call get_nonlinear_volume_fraction(state, nvfrac)
+         ewrite_minmax(nvfrac)
+      end if
       
       pressure => extract_scalar_field(state, "Pressure")
   
@@ -552,7 +579,11 @@
         ! Important note: with SUPG the test function derivatives have not been
         ! modified.
   
-        ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*drhodp_at_quad)
+        if(multiphase) then
+           ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*ele_val_at_quad(nvfrac, ele)*drhodp_at_quad)
+        else
+           ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*drhodp_at_quad)
+        end if
         !       /
         ! rhs = |test_shape* &
         !       /
@@ -560,8 +591,14 @@
         ! +(absorption)*(drhodp*theta*(eospressure - (pressure + atmospheric_pressure)) 
         !                - theta*density - (1-theta)*olddensity)
         ! +source)dV
-        ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) &
-                                                       +(olddensity_at_quad - density_at_quad)))
+
+        if(multiphase) then
+            ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*(ele_val_at_quad(nvfrac, ele))*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) &
+                                                       +(olddensity_at_quad - density_at_quad)))
+        else
+            ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) &
+                                                       +(olddensity_at_quad - density_at_quad)))
+        end if
         
         if(have_source) then
           ele_rhs = ele_rhs + shape_rhs(test_shape, detwei*ele_val_at_quad(source, ele))
@@ -587,6 +624,10 @@
   
       call deallocate(drhodp)
       call deallocate(eospressure)
+
+      if(multiphase) then
+         call deallocate(nvfrac)
+      end if
     
     end if
 
@@ -597,17 +638,29 @@
     type(state_type), dimension(:), intent(inout) :: state
     
     type(scalar_field), pointer :: density
-    
-    if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
-    
-      density=>extract_scalar_field(state(1),'Density')
-      
-      if(have_option(trim(density%option_path)//"/prognostic")) then
-        
-        call compressible_eos(state(1), density=density)
-      
-      end if
-    
+
+    integer :: istate
+    
+    if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+       do istate=1,size(state)
+          density=>extract_scalar_field(state(istate),'Density')
+          
+          if(have_option(trim(density%option_path)//"/prognostic")) then
+            call compressible_eos(state(istate), density=density)
+          end if
+       end do
+    else
+       if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
+       
+          density=>extract_scalar_field(state(1),'Density')
+          
+          if(have_option(trim(density%option_path)//"/prognostic")) then
+         
+             call compressible_eos(state(1), density=density)
+         
+          end if
+      
+       end if
     end if
   
   end subroutine update_compressible_density
@@ -620,8 +673,8 @@
     do i=0, option_count("/material_phase")-1
       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")) then
-        FLExit("With a DG pressure you cannot have use_compressible_projection")
+        .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")
       end if
     end do
 
 
=== modified file 'assemble/Diagnostic_Fields_Matrices.F90'
--- assemble/Diagnostic_Fields_Matrices.F90	2012-05-27 16:09:33 +0000
+++ assemble/Diagnostic_Fields_Matrices.F90	2012-08-09 00:46:20 +0000
@@ -40,7 +40,7 @@
   use spud
   use parallel_fields, only: zero_non_owned
   use divergence_matrix_cv, only: assemble_divergence_matrix_cv
-  use divergence_matrix_cg, only: assemble_divergence_matrix_cg
+  use divergence_matrix_cg, only: assemble_divergence_matrix_cg, assemble_compressible_divergence_matrix_cg
   use gradient_matrix_cg, only: assemble_gradient_matrix_cg
   use parallel_tools
   use sparsity_patterns, only: make_sparsity
@@ -58,7 +58,8 @@
   
   public :: calculate_divergence_cv, calculate_divergence_fe, &
             calculate_div_t_cv, calculate_div_t_fe, &
-            calculate_grad_fe, calculate_sum_velocity_divergence
+            calculate_grad_fe, calculate_sum_velocity_divergence, &
+            calculate_compressible_continuity_residual
 
 contains
 
@@ -356,9 +357,9 @@
 
       ! Allocate memory for matrices and sparsity patterns
       call allocate(ctfield, sum_velocity_divergence%mesh, name="CTField")
-      call zero (ctfield)
+      call zero(ctfield)
       call allocate(temp, sum_velocity_divergence%mesh, name="Temp")
-      
+
       ! Require a mass matrix type if not CV tested equation.
       if (.not. test_with_cv_dual) then
          mass_sparsity=make_sparsity(sum_velocity_divergence%mesh, sum_velocity_divergence%mesh, "MassSparsity")
@@ -382,27 +383,29 @@
 
          ! Allocate sparsity patterns, C^T matrix and C^T RHS for current state
          divergence_sparsity=make_sparsity(sum_velocity_divergence%mesh, u%mesh, "DivergenceSparsity")
+
          allocate(ct_m)
          call allocate(ct_m, divergence_sparsity, (/1, u%dim/), name="DivergenceMatrix" )
          call allocate(ct_rhs, sum_velocity_divergence%mesh, name="CTRHS")
-         
+
          ! Reassemble C^T matrix here
          if (test_with_cv_dual) then
             call assemble_divergence_matrix_cv(ct_m, state(i), ct_rhs=ct_rhs, &
                               test_mesh=sum_velocity_divergence%mesh, field=u)
          else
-            if(i==1) then ! Construct the mass matrix (just do this once)                      
+            if(i==1) then ! Construct the mass matrix (just do this once)
                call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, &
-                                 test_mesh=sum_velocity_divergence%mesh, field=u, &
-                                 option_path=sum_velocity_divergence%option_path, div_mass=mass)
+                              test_mesh=sum_velocity_divergence%mesh, field=u, &
+                              option_path=sum_velocity_divergence%option_path, div_mass=mass)
             else
                call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, &
-                                 test_mesh=sum_velocity_divergence%mesh, field=u, &
-                                 option_path=sum_velocity_divergence%option_path)
-            end if            
+                              test_mesh=sum_velocity_divergence%mesh, field=u, &
+                              option_path=sum_velocity_divergence%option_path)
+            end if   
+
          end if
 
-         ! Construct the linear system of equations to solve for \sum{div(vfrac*u)}
+         ! Construct the linear system of equations
          call zero(temp)
          call mult(temp, ct_m, u)
          call addto(temp, ct_rhs, -1.0)
@@ -417,7 +420,8 @@
 
       end do
 
-      ! Solve for sum_velocity_divergence ( = \sum{div(vfrac*u)} )            
+      ! Solve for sum_velocity_divergence
+      ! ( = \sum{div(vfrac*u)} for incompressible multiphase flows )            
       if (test_with_cv_dual) then
          ! get the cv mass matrix
          cv_mass => get_cv_mass(state(1), sum_velocity_divergence%mesh)
@@ -442,4 +446,148 @@
          
   end subroutine calculate_sum_velocity_divergence
 
+  
+  subroutine calculate_compressible_continuity_residual(state, compressible_continuity_residual)
+      !!< Calculates the residual of the continity equation used in compressible multiphase flow simulations:
+      !!< vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) }
+      
+      type(state_type), dimension(:), intent(inout) :: state
+      type(scalar_field), pointer :: compressible_continuity_residual, density, olddensity, vfrac
+      type(scalar_field) :: drhodt
+
+      ! Local variables
+      type(vector_field), pointer :: u, x
+      integer :: i, stat, ele
+      logical :: dg
+
+      type(csr_sparsity) :: divergence_sparsity
+      type(block_csr_matrix), pointer :: ct_m
+
+      type(csr_sparsity) :: mass_sparsity
+      type(csr_matrix) :: mass
+      type(scalar_field) :: ctfield, ct_rhs, temp
+      
+      type(element_type) :: test_function
+      type(element_type), pointer :: compressible_continuity_residual_shape
+      integer, dimension(:), pointer :: compressible_continuity_residual_nodes
+      real, dimension(:), allocatable :: detwei
+      real, dimension(:), allocatable :: drhodt_addto
+      
+      real :: dt
+
+      ewrite(1,*) 'Entering calculate_compressible_continuity_residual'
+
+      ! Allocate memory for matrices and sparsity patterns
+      call allocate(ctfield, compressible_continuity_residual%mesh, name="CTField")
+      call zero(ctfield)
+      call allocate(temp, compressible_continuity_residual%mesh, name="Temp")
+            
+      mass_sparsity=make_sparsity(compressible_continuity_residual%mesh, compressible_continuity_residual%mesh, "MassSparsity")
+      call allocate(mass, mass_sparsity, name="MassMatrix")
+      call zero(mass)
+      
+      call get_option("/timestepping/timestep", dt)
+      
+      ! Sum up over the div's
+      do i = 1, size(state)
+         u => extract_vector_field(state(i), "Velocity", stat)
+
+         ! If there's no velocity then cycle
+         if(stat/=0) cycle
+         ! If this is an aliased velocity then cycle
+         if(aliased(u)) cycle
+         ! If the velocity isn't prognostic then cycle
+         if(.not.have_option(trim(u%option_path)//"/prognostic")) cycle
+
+         ! If velocity field is prognostic, begin calculations below
+         x => extract_vector_field(state(i), "Coordinate")
+
+         ! Allocate sparsity patterns, C^T matrix and C^T RHS for current state
+         divergence_sparsity=make_sparsity(compressible_continuity_residual%mesh, u%mesh, "DivergenceSparsity")
+         allocate(ct_m)
+         call allocate(ct_m, divergence_sparsity, (/1, u%dim/), name="DivergenceMatrix")
+         call allocate(ct_rhs, compressible_continuity_residual%mesh, name="CTRHS")
+
+         ! Reassemble C^T matrix here
+         if(i==1) then ! Construct the mass matrix (just do this once)
+            call assemble_compressible_divergence_matrix_cg(ct_m, state, i, ct_rhs, div_mass=mass)
+         else
+            call assemble_compressible_divergence_matrix_cg(ct_m, state, i, ct_rhs)
+         end if   
+
+         if(have_option("/material_phase::"//trim(state(i)%name)//"/equation_of_state/compressible")) then
+            ! Get the time derivative term for the compressible phase's density, vfrac_c * d(rho_c)/dt
+
+            call allocate(drhodt, compressible_continuity_residual%mesh, name="drhodt")
+                  
+            density => extract_scalar_field(state(i), "Density", stat)
+            olddensity => extract_scalar_field(state(i), "OldDensity", stat)
+            vfrac => extract_scalar_field(state(i), "PhaseVolumeFraction")
+            
+            ! Assumes Density and OldDensity are on the same mesh as Pressure,
+            ! as it should be according to the manual.
+            call zero(drhodt)
+            
+            dg = continuity(compressible_continuity_residual) < 0
+            
+            element_loop: do ele = 1, element_count(compressible_continuity_residual)
+
+               if(.not.dg .or. (dg .and. element_owned(compressible_continuity_residual,ele))) then
+               
+                  allocate(detwei(ele_ngi(compressible_continuity_residual, ele)))
+                  allocate(drhodt_addto(ele_loc(compressible_continuity_residual, ele)))
+               
+                  compressible_continuity_residual_nodes => ele_nodes(compressible_continuity_residual, ele)
+                  compressible_continuity_residual_shape => ele_shape(compressible_continuity_residual, ele)
+                  test_function = compressible_continuity_residual_shape
+                  
+                  call transform_to_physical(x, ele, detwei=detwei)
+                  
+                  drhodt_addto = shape_rhs(test_function, detwei*(ele_val_at_quad(vfrac,ele)*(ele_val_at_quad(density,ele) - ele_val_at_quad(olddensity,ele))/dt))
+                  
+                  call addto(drhodt, compressible_continuity_residual_nodes, drhodt_addto)
+                  
+                  deallocate(detwei)
+                  deallocate(drhodt_addto)
+               end if
+
+            end do element_loop
+            
+            call addto(ctfield, drhodt)    
+            
+            call deallocate(drhodt)
+            
+         end if
+
+         ! Construct the linear system of equations
+         call zero(temp)
+         call mult(temp, ct_m, u)
+         call addto(temp, ct_rhs, -1.0)
+
+         ! Now add it to the sum
+         call addto(ctfield, temp)
+
+         call deallocate(ct_m)
+         deallocate(ct_m)
+         call deallocate(ct_rhs)
+         call deallocate(divergence_sparsity)
+
+      end do
+
+      ! Solve for compressible_continuity_residual      
+      call zero(compressible_continuity_residual)
+      call petsc_solve(compressible_continuity_residual, mass, ctfield)
+      ewrite_minmax(compressible_continuity_residual)
+         
+      ! Deallocate memory
+      call deallocate(ctfield)
+      call deallocate(temp)
+      call deallocate(mass_sparsity)
+      call deallocate(mass)
+
+      ewrite(1,*) 'Exiting calculate_compressible_continuity_residual'
+         
+  end subroutine calculate_compressible_continuity_residual
+  
+  
 end module diagnostic_fields_matrices
 
=== modified file 'assemble/Diagnostic_fields_wrapper.F90'
--- assemble/Diagnostic_fields_wrapper.F90	2012-02-18 15:45:48 +0000
+++ assemble/Diagnostic_fields_wrapper.F90	2012-08-09 00:46:20 +0000
@@ -585,6 +585,21 @@
             FLExit("The SumVelocityDivergence field is only used in multiphase simulations.")
          end if
        end if
+       
+       s_field => extract_scalar_field(state(i), "CompressibleContinuityResidual", stat)
+       if(stat == 0) then
+         ! Check that we are running a compressible multiphase simulation
+         if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1 .and. option_count("/material_phase/equation_of_state/compressible") > 0) then 
+            diagnostic = have_option(trim(s_field%option_path)//"/diagnostic")
+            if(diagnostic .and. .not.(aliased(s_field))) then
+               if(recalculate(trim(s_field%option_path))) then
+                  call calculate_compressible_continuity_residual(state, s_field)
+               end if
+            end if
+         else
+            FLExit("The CompressibleContinuityResidual field is only used in compressible multiphase simulations.")
+         end if
+       end if
 
        ! end of fields that cannot be called through the generic
        ! calculate_diagnostic_variable interface, i.e. - those that need things
 
=== modified file 'assemble/Divergence_Matrix_CG.F90'
--- assemble/Divergence_Matrix_CG.F90	2012-05-27 16:09:33 +0000
+++ assemble/Divergence_Matrix_CG.F90	2012-08-09 00:46:20 +0000
@@ -58,7 +58,10 @@
   integer :: nu_bar_scheme
   real :: nu_bar_scale
 
-contains
+  !! Are we running a multiphase flow simulation?
+  logical :: multiphase
+
+  contains
 
     subroutine assemble_divergence_matrix_cg(CT_m, state, ct_rhs, & 
                                              test_mesh, field, option_path, &
@@ -119,7 +122,6 @@
       logical :: l_get_ct
 
       !! Multiphase variables
-      logical :: multiphase
       ! Volume fraction fields
       type(scalar_field), pointer :: vfrac
       type(scalar_field) :: nvfrac
@@ -355,40 +357,55 @@
 
     end subroutine assemble_divergence_matrix_cg
 
-    subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs)
+    subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass)
       
       ! inputs/outputs
       ! bucket full of fields
       type(state_type), dimension(:), intent(inout) :: state
-
+      integer, intent(in) :: istate
       ! the compressible divergence matrix
       type(block_csr_matrix), pointer :: ctp_m
 
       type(scalar_field), intent(inout), optional :: ct_rhs
 
-      if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
-      
-        call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state(1), ct_rhs)
-        
+      type(csr_matrix), intent(inout), optional :: div_mass
+
+      if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
+         multiphase = .true.
+         call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass)
       else
-      
-        FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
-        
+         multiphase = .false.
+
+         if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
+         
+            call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, 1, ct_rhs, div_mass)
+         
+         else
+
+            FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
+         
+         end if
+
       end if
+
+
     
     end subroutine assemble_compressible_divergence_matrix_cg
 
-    subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs) 
+    subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass) 
 
       ! inputs/outputs
       ! bucket full of fields
-      type(state_type), intent(inout) :: state
+      type(state_type), dimension(:), intent(inout) :: state
+      integer, intent(in) :: istate
 
       ! the compressible divergence matrix
       type(block_csr_matrix), pointer :: ctp_m
 
       type(scalar_field), intent(inout), optional :: ct_rhs
 
+      type(csr_matrix), intent(inout), optional :: div_mass
+
       ! local
       type(mesh_type), pointer :: test_mesh
 
@@ -408,6 +425,8 @@
       real, dimension(:), allocatable :: density_at_quad, olddensity_at_quad
       real, dimension(:,:), allocatable :: density_grad_at_quad, nlvelocity_at_quad
 
+      real, dimension(:,:), allocatable :: div_mass_mat
+
       ! loop integers
       integer :: ele, sele, dim
 
@@ -425,7 +444,16 @@
       integer, dimension(:), allocatable :: density_bc_type
       type(scalar_field) :: density_bc
       
-      integer :: stat
+      integer :: i, stat
+
+      !! Multiphase variables
+      ! Volume fraction fields
+      type(scalar_field), pointer :: vfrac
+      type(scalar_field) :: nvfrac
+      type(element_type), pointer :: nvfrac_shape
+      ! Transformed gradient function for the non-linear PhaseVolumeFraction. 
+      real, dimension(:, :, :), allocatable :: dnvfrac_t
+      logical :: is_compressible_phase ! Is this the (single) compressible phase in a multiphase simulation?
 
       ! =============================================================
       ! Subroutine to construct the matrix ctp_m (a.k.a. C1/2/3TP).
@@ -433,16 +461,42 @@
 
       ewrite(2,*) 'In assemble_1mat_compressible_divergence_matrix_cg'
 
-      coordinate=> extract_vector_field(state, "Coordinate")
-      
-      density => extract_scalar_field(state, "Density")
-      olddensity => extract_scalar_field(state, "OldDensity")
-      
-      pressure => extract_scalar_field(state, "Pressure")
-      
-      velocity=>extract_vector_field(state, "Velocity")
-      nonlinearvelocity=>extract_vector_field(state, "NonlinearVelocity") ! maybe this should be updated after the velocity solve?
-      
+      coordinate=> extract_vector_field(state(istate), "Coordinate")
+      
+      density => extract_scalar_field(state(istate), "Density")
+      olddensity => extract_scalar_field(state(istate), "OldDensity")
+      
+      if(have_option(trim(state(istate)%option_path)//"/equation_of_state/compressible")) then
+         is_compressible_phase = .true.
+      else
+         is_compressible_phase = .false.
+
+         ! Find the compressible phase and extract it's density to form rho_c * div(vfrac_i*u_i), where _c and _i represent
+         ! the compressible and incompressible phases respectively.
+         do i = 1, size(state)
+            if(have_option(trim(state(i)%option_path)//"/equation_of_state/compressible")) then
+               density => extract_scalar_field(state(i), "Density")
+               olddensity => extract_scalar_field(state(i), "OldDensity")
+            end if
+         end do
+      end if
+
+      pressure => extract_scalar_field(state(istate), "Pressure")
+      
+      velocity=>extract_vector_field(state(istate), "Velocity")
+      nonlinearvelocity=>extract_vector_field(state(istate), "NonlinearVelocity") ! maybe this should be updated after the velocity solve?
+      
+      ! Get the non-linear PhaseVolumeFraction field if multiphase
+      if(multiphase) then
+         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)
+      else
+         nullify(vfrac)
+      end if
+
       integrate_by_parts=have_option(trim(complete_field_path(density%option_path, stat=stat))//&
           &"/spatial_discretisation/continuous_galerkin/advection_terms/integrate_advection_by_parts")&
           .or. have_option(trim(complete_field_path(velocity%option_path, stat=stat))//&
@@ -474,7 +528,7 @@
       end if
 
       call get_option(trim(complete_field_path(density%option_path, stat=stat))//&
-          &"/temporal_discretisation/theta", theta)
+                         &"/temporal_discretisation/theta", theta)
       call get_option("/timestepping/timestep", dt)
       
       test_mesh => pressure%mesh
@@ -494,8 +548,14 @@
                olddensity_at_quad(ele_ngi(density, 1)), &
                nlvelocity_at_quad(nonlinearvelocity%dim, ele_ngi(nonlinearvelocity, 1)), &
                density_grad_at_quad(field%dim, ele_ngi(density,1)), &
-               j_mat(field%dim, field%dim, ele_ngi(density, 1)))
+               j_mat(field%dim, field%dim, ele_ngi(density, 1)), &
+               div_mass_mat(ele_loc(test_mesh, 1), ele_loc(test_mesh, 1)))
       
+      if(multiphase) then
+         ! We will need grad(nvfrac) if we are not integrating by parts below
+         allocate(dnvfrac_t(ele_loc(nvfrac,1), ele_ngi(nvfrac,1), field%dim))
+      end if
+
       do ele=1, element_count(test_mesh)
 
         test_nodes=>ele_nodes(test_mesh, ele)
@@ -509,7 +569,7 @@
         olddensity_at_quad = ele_val_at_quad(olddensity, ele)
         
         nlvelocity_at_quad = ele_val_at_quad(nonlinearvelocity, ele)
-        
+
         if(any(stabilisation_scheme == (/STABILISATION_STREAMLINE_UPWIND, STABILISATION_SUPG/))) then
           call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, &
                                      detwei = detwei, j = j_mat)
@@ -517,7 +577,7 @@
           call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, detwei=detwei)
         end if
         
-        if(.not.integrate_by_parts) then
+        if(.not.integrate_by_parts .or. (multiphase .and. integrate_by_parts .and. .not.is_compressible_phase)) then
           ! transform the field (velocity) derivatives into physical space
           call transform_to_physical(coordinate, ele, field_shape, dshape=dfield_t)
           
@@ -545,16 +605,53 @@
 
 
         if(integrate_by_parts) then
+
             ! if SUPG is fixed for P>1 then this dtest_t should be updated
-            ele_mat = -dshape_shape(dtest_t, field_shape, &
-                       detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
+            if(multiphase .and. .not.is_compressible_phase) then
+               density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+&
+                                   (1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t))
+
+               ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) - shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad)
+
+            else if(multiphase .and. is_compressible_phase) then
+               ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
+         
+            else
+         
+               ele_mat = -dshape_shape(dtest_t, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
+            end if
         else
             density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+&
                                    (1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t))
             
-            ele_mat = shape_dshape(test_shape, dfield_t, &
+            if(multiphase) then
+
+               ! 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)
+                  call transform_to_physical(coordinate, ele, nvfrac_shape, dshape=dnvfrac_t)
+               else
+                  dnvfrac_t = dfield_t
+               end if
+
+               ! Split up the divergence term div(rho*vfrac*u) = rho*div(u*vfrac) + vfrac*u*grad(rho)
+               ! = (rho*vfrac*div(u) + rho*u*grad(vfrac)) + u*vfrac*grad(rho)
+
+               ! First assemble rho*div(u*vfrac). This is the incompressible phase's divergence matrix.
+               ele_mat = shape_dshape(test_shape, dfield_t, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)*ele_val_at_quad(nvfrac, ele)) + shape_shape_vector(test_shape, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad), ele_grad_at_quad(nvfrac, ele, dnvfrac_t))
+
+               ! If the phase is compressible, then we now complete the assembly of div(rho*vfrac*u) below.
+               if(is_compressible_phase) then
+                  ele_mat = ele_mat + shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad)
+               end if
+
+            else
+               ele_mat = shape_dshape(test_shape, dfield_t, &
                                    detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) + &
                       shape_shape_vector(test_shape, field_shape, detwei, density_grad_at_quad)
+            end if
+
         end if
         
         ! Stabilisation does not return the right shape for this operator!
@@ -568,6 +665,11 @@
         do dim = 1, field%dim
             call addto(ctp_m, 1, dim, test_nodes, field_nodes, ele_mat(dim,:,:))
         end do
+
+        if(present(div_mass)) then
+           div_mass_mat = shape_shape(test_shape, test_shape, detwei)
+           call addto(div_mass, test_nodes, test_nodes, div_mass_mat)         
+        end if
         
         call deallocate(test_shape)
         
@@ -618,8 +720,14 @@
               &                          detwei_f=detwei_bdy,&
               &                          normal=normal_bdy) 
 
-          ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
-                                           detwei_bdy*density_bdy, normal_bdy)
+
+          if(multiphase) then
+            ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
+                                             detwei_bdy*density_bdy*face_val_at_quad(nvfrac, sele), normal_bdy)
+          else
+            ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
+                                             detwei_bdy*density_bdy, normal_bdy)
+          end if
 
           do dim = 1, field%dim
             if((field_bc_type(dim, sele)==1).and.present(ct_rhs)) then
@@ -641,6 +749,11 @@
         deallocate(test_nodes_bdy, field_nodes_bdy)
 
       end if
+
+      if(multiphase) then
+         deallocate(dnvfrac_t)
+         call deallocate(nvfrac)
+      end if
       
     end subroutine assemble_1mat_compressible_divergence_matrix_cg
 
 
=== modified file 'assemble/Makefile.dependencies'
--- assemble/Makefile.dependencies	2012-08-02 13:03:52 +0000
+++ assemble/Makefile.dependencies	2012-08-09 00:46:20 +0000
@@ -79,12 +79,12 @@
    ../include/boundary_conditions_from_options.mod ../include/colouring.mod \
    ../include/elements.mod ../include/fdebug.h ../include/field_options.mod \
    ../include/fields.mod ../include/fldebug.mod \
-   ../include/global_parameters.mod ../include/petsc_solve_state_module.mod \
-   ../include/porous_media.mod ../include/profiler.mod \
-   ../include/quadrature.mod ../include/sparse_tools.mod \
-   ../include/sparse_tools_petsc.mod ../include/sparsity_patterns_meshes.mod \
-   ../include/spud.mod ../include/state_module.mod \
-   ../include/upwind_stabilisation.mod
+   ../include/global_parameters.mod ../include/multiphase_module.mod \
+   ../include/petsc_solve_state_module.mod ../include/porous_media.mod \
+   ../include/profiler.mod ../include/quadrature.mod \
+   ../include/sparse_tools.mod ../include/sparse_tools_petsc.mod \
+   ../include/sparsity_patterns_meshes.mod ../include/spud.mod \
+   ../include/state_module.mod ../include/upwind_stabilisation.mod
 
 ../include/advection_diffusion_dg.mod: Advection_Diffusion_DG.o
 	@true
@@ -154,10 +154,10 @@
    Compressible_Projection.F90 ../include/equation_of_state.mod \
    ../include/fdebug.h ../include/fefields.mod ../include/field_options.mod \
    ../include/fields.mod ../include/fldebug.mod \
-   ../include/global_parameters.mod ../include/sparse_matrices_fields.mod \
-   ../include/sparse_tools.mod ../include/spud.mod \
-   ../include/state_fields_module.mod ../include/state_module.mod \
-   ../include/upwind_stabilisation.mod
+   ../include/global_parameters.mod ../include/multiphase_module.mod \
+   ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \
+   ../include/spud.mod ../include/state_fields_module.mod \
+   ../include/state_module.mod ../include/upwind_stabilisation.mod
 
 ../include/coriolis_module.mod: Coriolis.o
 	@true
 
=== modified file 'assemble/Momentum_Equation.F90'
--- assemble/Momentum_Equation.F90	2012-04-16 14:55:55 +0000
+++ assemble/Momentum_Equation.F90	2012-08-09 00:46:20 +0000
@@ -92,7 +92,7 @@
       logical :: assemble_schur_auxiliary_matrix
 
       ! Do we want to use the compressible projection method?
-      logical :: use_compressible_projection
+      logical :: compressible_eos
       ! Are we doing a full Schur solve?
       logical :: full_schur
       ! Are we lumping mass or assuming consistent mass?
@@ -282,7 +282,7 @@
             multiphase = .false.
          end if
          ! Do we have fluid-particle drag (for multi-phase simulations)?
-         have_fp_drag = option_count("/material_phase/multiphase_properties/particle_diameter") > 0
+         have_fp_drag = have_option("/multiphase_interaction/fluid_particle_drag")
 
          ! Get the pressure p^{n}, and get the assembly options for the divergence and CMC matrices
          ! find the first non-aliased pressure
@@ -388,7 +388,7 @@
                ! get the CV tested pressure gradient matrix (i.e. the divergence matrix)
                ! if required with a different unique name. Note there is no need
                ! to again decide reassemble_ct_m as ctp_m for this case is assembled when ct_m is.
-               if ((.not. use_compressible_projection) .and. cg_pressure_cv_test_continuity) then
+               if ((.not. compressible_eos) .and. cg_pressure_cv_test_continuity) then
                   ctp_m(istate)%ptr => get_velocity_divergence_matrix(state(istate), ct_m_name = "CVTestedVelocityDivergenceMatrix")
                end if
 
@@ -467,17 +467,16 @@
                end select
             end if
 
-            if (has_boundary_condition(u, "free_surface") .or. use_compressible_projection) then
-               ! this needs fixing for multiphase theta_pg could in principle be chosen
+            if (has_boundary_condition(u, "free_surface") .or. compressible_eos) then
+               ! This needs fixing for multiphase. theta_pg could in principle be chosen
                ! per phase but then we need an array and we'd have to include theta_pg
                ! in cmc_m, i.e. solve for theta_div*dt*dp instead of theta_div*theta_pg*dt*dp
                ! theta_div can only be set once, but where and what is the default?
-               if (multiphase) then
-                 FLExit("Multiphase does not work with a free surface or compressible flow.")
+               if (has_boundary_condition(u, "free_surface") .and. multiphase) then
+                 FLExit("Multiphase does not work with a free surface.")
                end if
 
-
-               ! With free surface or compressible-projection pressures are at integer
+               ! With free surface or compressible-projection, pressures are at integer
                ! time levels and we apply a theta-weighting to the pressure gradient term
                ! Also, obtain theta-weighting to be used in divergence term
                call get_option( trim(u%option_path)//'/prognostic/temporal_discretisation/theta', &
@@ -491,6 +490,15 @@
                ewrite(2,*) "theta_pg: ", theta_pg
                ewrite(2,*) "Velocity divergence is evaluated at n+theta_divergence"
                ewrite(2,*) "theta_divergence: ", theta_divergence
+               
+               ! Note: Compressible multiphase simulations work, but only when use_theta_pg and use_theta_divergence
+               ! are false. This needs improving - see comment above.
+               if(compressible_eos .and. multiphase .and. (use_theta_pg .or. use_theta_divergence)) then
+                  ewrite(-1,*) "Currently, for compressible multiphase flow simulations, the"
+                  ewrite(-1,*) "temporal_discretisation/theta and temporal_discretisation/theta_divergence values"
+                  ewrite(-1,*) "for each Velocity field must be set to 1.0."
+                  FLExit("Multiphase does not work when use_theta_pg or use_theta_divergence are true.")
+               end if
             else
                ! Pressures are, as usual, staggered in time with the velocities
                use_theta_pg=.false.
@@ -514,6 +522,7 @@
             ! Allocation of big_m
             if(dg(istate)) then
                call allocate_big_m_dg(state(istate), big_m(istate), u)
+
                if(subcycle(istate)) then
                   u_sparsity => get_csr_sparsity_firstorder(state, u%mesh, u%mesh)
                   ! subcycle_m currently only contains advection, so diagonal=.true.
@@ -531,8 +540,8 @@
             ! Initialise the big_m, ct_m and ctp_m matrices
             call zero(big_m(istate))
             if(reassemble_ct_m) then
-               call zero(ct_m(istate)%ptr)               
-               if ((.not. use_compressible_projection) .and. cg_pressure_cv_test_continuity) then
+               call zero(ct_m(istate)%ptr)         
+               if ((.not. compressible_eos) .and. cg_pressure_cv_test_continuity) then
                   call zero(ctp_m(istate)%ptr)
                end if
             end if
@@ -580,7 +589,7 @@
                   call subtract_geostrophic_pressure_gradient(mom_rhs(istate), state(istate))
                end if
             else
-               ! This call will form the ct_rhs, which for use_compressible_projection
+               ! This call will form the ct_rhs, which for compressible_eos
                ! or cg_pressure_cv_test_continuity is formed for a second time later below.
                call construct_momentum_cg(u, p, density, x, &
                      big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, &
@@ -635,7 +644,7 @@
 
             call profiler_tic(p, "assembly")
             if(cv_pressure) then
-               ! This call will form the ct_rhs, which for use_compressible_projection
+               ! This call will form the ct_rhs, which for compressible_eos
                ! is formed for a second time later below.
                call assemble_divergence_matrix_cv(ct_m(istate)%ptr, state(istate), ct_rhs=ct_rhs(istate), &
                                              test_mesh=p%mesh, field=u, get_ct=reassemble_ct_m)
@@ -643,7 +652,7 @@
 
             ! Assemble divergence matrix C^T.
             ! At the moment cg does its own ct assembly. We might change this in the future.
-            ! This call will form the ct_rhs, which for use_compressible_projection
+            ! This call will form the ct_rhs, which for compressible_eos
             ! or cg_pressure_cv_test_continuity is formed for a second time later below.
             if(dg(istate) .and. .not. cv_pressure) then
                call assemble_divergence_matrix_cg(ct_m(istate)%ptr, state(istate), ct_rhs=ct_rhs(istate), &
@@ -683,7 +692,11 @@
                
                ! Set up the left C matrix in CMC
                
-               if(use_compressible_projection) then
+               if(compressible_eos) then
+                  ! Note: If we are running a compressible multiphase simulation then the C^T matrix for each phase becomes:
+                  ! rho*div(vfrac*u) for each incompressible phase
+                  ! rho*div(vfrac*u) + vfrac*u*grad(rho) for the single compressible phase.
+
                   allocate(ctp_m(istate)%ptr)
                   call allocate(ctp_m(istate)%ptr, ct_m(istate)%ptr%sparsity, (/1, u%dim/), name="CTP_m")
                   ! NOTE that this is not optimal in that the ct_rhs
@@ -691,7 +704,7 @@
                   if(cv_pressure) then
                      call assemble_compressible_divergence_matrix_cv(ctp_m(istate)%ptr, state, ct_rhs(istate))
                   else
-                     call assemble_compressible_divergence_matrix_cg(ctp_m(istate)%ptr, state, ct_rhs(istate))
+                     call assemble_compressible_divergence_matrix_cg(ctp_m(istate)%ptr, state, istate, ct_rhs(istate))
                   end if               
                else                  
                   ! Incompressible scenario
@@ -708,7 +721,7 @@
                   end if
                end if
                
-               if (use_compressible_projection .or. cg_pressure_cv_test_continuity) then
+               if (compressible_eos .or. cg_pressure_cv_test_continuity) then
                   if (have_rotated_bcs(u)) then
                      if (dg(istate)) then
                        call zero_non_owned(u)
@@ -1041,7 +1054,7 @@
 
                      call profiler_toc(u, "assembly")
 
-                     if(use_compressible_projection) then
+                     if(compressible_eos) then
                         call deallocate(ctp_m(istate)%ptr)
                         deallocate(ctp_m(istate)%ptr)
                      end if
@@ -1252,15 +1265,12 @@
 
          ewrite(1,*) 'Entering get_pressure_options'
 
-
          ! Are we using a compressible projection?
-         use_compressible_projection = have_option(trim(p%option_path)//&
-                                       "/prognostic/scheme&
-                                       &/use_compressible_projection_method")
+         compressible_eos = option_count("/material_phase/equation_of_state/compressible") > 0
 
          reassemble_all_cmc_m = have_option(trim(p%option_path)//&
                      "/prognostic/scheme/update_discretised_equation") .or. &
-                     use_compressible_projection
+                     compressible_eos
 
          reassemble_all_ct_m = have_option(trim(p%option_path)//&
                      "/prognostic/scheme/update_discretised_equation")
@@ -1579,15 +1589,15 @@
          call deallocate(kmk_rhs)
 
          cmc_m => extract_csr_matrix(state(istate), "PressurePoissonMatrix", stat)
-
-         if(use_compressible_projection) then
+         
+         if(compressible_eos .and. have_option(trim(state(istate)%option_path)//'/equation_of_state/compressible')) then
             call allocate(compress_projec_rhs, p%mesh, "CompressibleProjectionRHS")
 
             if(cv_pressure) then
                call assemble_compressible_projection_cv(state, cmc_m, compress_projec_rhs, dt, &
                                                       theta_pg, theta_divergence, reassemble_cmc_m)
             else
-               call assemble_compressible_projection_cg(state, cmc_m, compress_projec_rhs, dt, &
+               call assemble_compressible_projection_cg(state, istate, cmc_m, compress_projec_rhs, dt, &
                                                       theta_pg, theta_divergence, reassemble_cmc_m)
             end if
 
@@ -1720,7 +1730,7 @@
          call addto(p, delta_p, scale=1.0/(theta_pg*dt))
          ewrite_minmax(p)
 
-         if(use_compressible_projection) then
+         if(compressible_eos) then
             call update_compressible_density(state)
          end if
 
 
=== modified file 'assemble/Multiphase.F90'
--- assemble/Multiphase.F90	2011-09-03 00:57:43 +0000
+++ assemble/Multiphase.F90	2012-08-09 00:46:20 +0000
@@ -46,7 +46,7 @@
       private
       public :: get_phase_submaterials, get_nonlinear_volume_fraction, &
                 calculate_diagnostic_phase_volume_fraction, &
-                add_fluid_particle_drag
+                add_fluid_particle_drag, add_heat_transfer
 
    contains
 
@@ -268,7 +268,7 @@
          integer :: ele
          type(element_type) :: test_function
          type(element_type), pointer :: u_shape
-         integer, dimension(:), pointer :: u_ele
+         integer, dimension(:), pointer :: u_nodes
          logical :: dg
              
          type(vector_field), pointer :: velocity_fluid
@@ -281,14 +281,20 @@
          integer :: i, dim
          logical :: not_found ! Error flag. Have we found the fluid phase?
          integer :: istate_fluid, istate_particle
-         
+   
+         ! Types of drag correlation
+         integer, parameter :: DRAG_CORRELATION_TYPE_STOKES = 1, DRAG_CORRELATION_TYPE_WEN_YU = 2, DRAG_CORRELATION_TYPE_ERGUN = 3
          
          ewrite(1, *) "Entering add_fluid_particle_drag"
+         
+         ! Let's check whether we actually have at least one particle phase.
+         if(option_count("/material_phase/multiphase_properties/particle_diameter") == 0) then
+            FLExit("Fluid-particle drag enabled but no particle_diameter has been specified for the particle phase(s).")
+         end if
                
          ! Get the timestepping options
          call get_option("/timestepping/timestep", dt)
-         call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", &
-                        theta)
+         call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", theta)
                                           
          ! For the big_m matrix. Controls whether the off diagonal entries are used    
          block_mask = .false.
@@ -300,15 +306,13 @@
          dg = continuity(u) < 0
 
          ! Is this phase a particle phase?
-         is_particle_phase = have_option("/material_phase["//int2str(istate-1)//&
-                          &"]/multiphase_properties/particle_diameter")
+         is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter")
               
          ! Retrieve the index of the fluid phase in the state array.
          not_found = .true.
          if(is_particle_phase) then    
             do i = 1, size(state)
-               if(.not.have_option("/material_phase["//int2str(i-1)//&
-                        &"]/multiphase_properties/particle_diameter")) then
+               if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then
 
                   velocity_fluid => extract_vector_field(state(i), "Velocity")
                   ! Aliased material_phases will also not have a particle_diameter,
@@ -359,6 +363,8 @@
                type(tensor_field), pointer :: viscosity_fluid
                type(scalar_field) :: nvfrac_fluid, nvfrac_particle
                real :: d ! Particle diameter
+               character(len=OPTION_PATH_LEN) :: drag_correlation_name
+               integer :: drag_correlation
 
                ! Get the necessary fields to calculate the drag force
                velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
@@ -371,8 +377,7 @@
                   density_particle => extract_scalar_field(state(istate_particle), "Density")
                   viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
          
-                  call get_option("/material_phase["//int2str(istate_particle-1)//&
-                           &"]/multiphase_properties/particle_diameter", d)
+                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
 
                   ! Calculate the non-linear approximation to the PhaseVolumeFractions
                   call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
@@ -387,13 +392,25 @@
                   nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity")
                   oldu_fluid => extract_vector_field(state(istate_fluid), "OldVelocity")
                   oldu_particle => extract_vector_field(state(istate_particle), "OldVelocity")
-      
+                  
+                  call get_option("/multiphase_interaction/fluid_particle_drag/drag_correlation/name", drag_correlation_name)
+                  select case(trim(drag_correlation_name))
+                     case("stokes")
+                        drag_correlation = DRAG_CORRELATION_TYPE_STOKES
+                     case("wen_yu")
+                        drag_correlation = DRAG_CORRELATION_TYPE_WEN_YU
+                     case("ergun")
+                        drag_correlation = DRAG_CORRELATION_TYPE_ERGUN
+                     case("default")
+                        FLAbort("Unknown correlation for fluid-particle drag")
+                  end select
+                  
                   ! ----- Volume integrals over elements -------------           
                   call profiler_tic(u, "element_loop")
                   element_loop: do ele = 1, element_count(u)
 
                      if(.not.dg .or. (dg .and. element_owned(u,ele))) then
-                        u_ele=>ele_nodes(u, ele)
+                        u_nodes => ele_nodes(u, ele)
                         u_shape => ele_shape(u, ele)
                         test_function = u_shape                         
 
@@ -403,7 +420,7 @@
                                                             density_fluid, density_particle, &
                                                             nu_fluid, nu_particle, &
                                                             oldu_fluid, oldu_particle, &
-                                                            viscosity_fluid, d)
+                                                            viscosity_fluid, d, drag_correlation)
                      end if
 
                   end do element_loop
@@ -421,7 +438,7 @@
                                                       density_fluid, density_particle, &
                                                       nu_fluid, nu_particle, &
                                                       oldu_fluid, oldu_particle, &
-                                                      viscosity_fluid, d)
+                                                      viscosity_fluid, d, drag_correlation)
                                                          
                integer, intent(in) :: ele
                type(element_type), intent(in) :: test_function
@@ -436,6 +453,7 @@
                type(vector_field), intent(in) :: oldu_fluid, oldu_particle
                type(tensor_field), intent(in) :: viscosity_fluid    
                real, intent(in) :: d ! Particle diameter 
+               integer, intent(in) :: drag_correlation
                
                ! Local variables
                real, dimension(ele_ngi(u,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
@@ -447,14 +465,14 @@
                
                real, dimension(ele_loc(u, ele), ele_ngi(u, ele), x%dim) :: du_t
                real, dimension(ele_ngi(u, ele)) :: detwei
-               real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs_mat
+               real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs
                real, dimension(ele_loc(u, ele), ele_loc(u, ele)) :: interaction_big_m_mat
                real, dimension(u%dim, ele_loc(u,ele)) :: rhs_addto
                real, dimension(u%dim, u%dim, ele_loc(u,ele), ele_loc(u,ele)) :: big_m_tensor_addto
                
-               real, dimension(ele_ngi(u,ele)) :: particle_re ! Particle Reynolds number
-               real, dimension(ele_ngi(u,ele)) :: drag_coefficient
-               real, dimension(ele_ngi(u,ele)) :: magnitude ! |v_f - v_p|
+               real, dimension(ele_ngi(u,ele)) :: particle_re_gi ! Particle Reynolds number
+               real, dimension(ele_ngi(u,ele)) :: drag_coefficient_gi
+               real, dimension(ele_ngi(u,ele)) :: magnitude_gi ! |v_f - v_p|
                
                real, dimension(ele_ngi(u,ele)) :: K
                real, dimension(ele_ngi(u,ele)) :: drag_force_big_m
@@ -476,34 +494,55 @@
          
                ! Compute the magnitude of the relative velocity
                do gi = 1, ele_ngi(u,ele)
-                  magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
+                  magnitude_gi(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
                end do
 
                ! Compute the particle Reynolds number
                ! (Assumes isotropic viscosity for now)
-               particle_re = (vfrac_fluid_gi*density_fluid_gi*magnitude*d) / viscosity_fluid_gi(1,1,:)
+               particle_re_gi = (vfrac_fluid_gi*density_fluid_gi*magnitude_gi*d) / viscosity_fluid_gi(1,1,:)
            
                ! Compute the drag coefficient
-               do gi = 1, ele_ngi(u,ele)
-                  if(particle_re(gi) < 1000) then
-                     drag_coefficient(gi) = (24.0/particle_re(gi))
-                     ! For the Wen & Yu (1966) drag term, this should be:
-                     ! drag_coefficient(gi) = (24.0/particle_re(gi))*(1.0+0.15*particle_re(gi)**0.687)
-                  else
-                     drag_coefficient(gi) = 0.44
-                  end if
-               end do
+               select case(drag_correlation)
+                  case(DRAG_CORRELATION_TYPE_STOKES)
+                     ! Stokes drag correlation
+                     do gi = 1, ele_ngi(u,ele)
+                        if(particle_re_gi(gi) < 1000) then
+                           drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi))
+                        else
+                           drag_coefficient_gi(gi) = 0.44
+                        end if
+                     end do
+                     
+                  case(DRAG_CORRELATION_TYPE_WEN_YU)
+                     ! Wen & Yu (1966) drag correlation
+                     do gi = 1, ele_ngi(u,ele)
+                        if(particle_re_gi(gi) < 1000) then
+                           drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi))*(1.0+0.15*particle_re_gi(gi)**0.687)
+                        else
+                           drag_coefficient_gi(gi) = 0.44
+                        end if
+                     end do
+                     
+                  case(DRAG_CORRELATION_TYPE_ERGUN)
+                     ! No drag coefficient is needed here.                  
+               end select
                       
-               ! Don't let the drag_coefficient be NaN
+               ! Don't let the drag_coefficient_gi be NaN
                do gi = 1, ele_ngi(u,ele)
-                  if(particle_re(gi) == 0) then
-                     drag_coefficient(gi) = 1e16
+                  if(particle_re_gi(gi) == 0) then
+                     drag_coefficient_gi(gi) = 1e16
                   end if
                end do
            
-               ! For the Wen & Yu (1966) drag term, this should be:
-               ! K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d*vfrac_fluid_gi**2.7)
-               K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d)
+               select case(drag_correlation)
+                  case(DRAG_CORRELATION_TYPE_STOKES)
+                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d)
+                  case(DRAG_CORRELATION_TYPE_WEN_YU)
+                     ! Wen & Yu (1966) drag term
+                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d*vfrac_fluid_gi**2.7)
+                  case(DRAG_CORRELATION_TYPE_ERGUN)
+                     K = 150.0*((vfrac_particle_gi**2)*viscosity_fluid_gi(1,1,:))/(vfrac_fluid_gi*(d**2)) + 1.75*(vfrac_particle_gi*density_fluid_gi*magnitude_gi/d)
+               end select               
                
                if(is_particle_phase) then
                   drag_force_big_m = -K
@@ -519,7 +558,7 @@
                
                ! Form the element interaction/drag matrix
                interaction_big_m_mat = shape_shape(test_function, u_shape, detwei*drag_force_big_m)
-               interaction_rhs_mat = shape_vector_rhs(test_function, drag_force_rhs, detwei)
+               interaction_rhs = shape_vector_rhs(test_function, drag_force_rhs, detwei)
               
                ! Add contribution  
                big_m_tensor_addto = 0.0            
@@ -531,16 +570,328 @@
                end if
                do dim = 1, u%dim
                   big_m_tensor_addto(dim, dim, :, :) = big_m_tensor_addto(dim, dim, :, :) - dt*theta*interaction_big_m_mat
-                  rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs_mat(dim,:)
+                  rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs(dim,:)
                end do
                
                ! Add the contribution to mom_rhs
-               call addto(mom_rhs, u_ele, rhs_addto) 
+               call addto(mom_rhs, u_nodes, rhs_addto) 
                ! Add to the big_m matrix
-               call addto(big_m, u_ele, u_ele, big_m_tensor_addto, block_mask=block_mask)
+               call addto(big_m, u_nodes, u_nodes, big_m_tensor_addto, block_mask=block_mask)
 
             end subroutine add_fluid_particle_drag_element
 
       end subroutine add_fluid_particle_drag
 
+      
+      !! Multiphase energy interaction term Q_i
+      !! to be added to the RHS of the internal energy equation.
+      subroutine add_heat_transfer(state, istate, internal_energy, matrix, rhs)
+         !!< This computes the inter-phase heat transfer term.
+         !!< Only between fluid and particle phase pairs.
+         !!< Uses the empirical correlation by Gunn (1978).
+         
+         type(state_type), dimension(:), intent(inout) :: state     
+         integer, intent(in) :: istate
+         type(scalar_field), intent(in) :: internal_energy         
+         type(csr_matrix), intent(inout) :: matrix
+         type(scalar_field), intent(inout) :: rhs
+         
+         ! Local variables              
+         integer :: ele
+         type(element_type) :: test_function
+         type(element_type), pointer :: internal_energy_shape
+         integer, dimension(:), pointer :: internal_energy_nodes
+         logical :: dg
+             
+         type(vector_field), pointer :: x, velocity_fluid
+         
+         real :: dt, theta
+         
+         logical :: is_particle_phase
+         logical :: not_found ! Error flag. Have we found the fluid phase?
+         integer :: i, istate_fluid, istate_particle
+         
+         
+         ewrite(1, *) "Entering add_heat_transfer"
+               
+         ! Get the timestepping options
+         call get_option("/timestepping/timestep", dt)
+         call get_option(trim(internal_energy%option_path)//"/prognostic/temporal_discretisation/theta", &
+                        theta)
+         
+         ! Get the coordinate field from state(istate)
+         x => extract_vector_field(state(istate), "Coordinate")
+         ewrite_minmax(x)
+         assert(x%dim == mesh_dim(internal_energy))
+         assert(ele_count(x) == ele_count(internal_energy))
+
+         ! Are we using a discontinuous Galerkin discretisation?
+         dg = continuity(internal_energy) < 0
+
+         ! Is this phase a particle phase?
+         is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter")
+              
+         ! Retrieve the index of the fluid phase in the state array.
+         not_found = .true.
+         if(is_particle_phase) then    
+            do i = 1, size(state)
+               if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then
+
+                  velocity_fluid => extract_vector_field(state(i), "Velocity")
+                  ! Aliased material_phases will also not have a particle_diameter,
+                  ! so here we make sure that we don't count these as the fluid phase
+                  if(.not.aliased(velocity_fluid)) then
+                     istate_fluid = i
+                     if(.not.not_found) then
+                        FLExit("Heat transfer term does not currently support more than one fluid phase.")
+                     end if
+                     not_found = .false.
+                  end if
+
+               end if
+            end do
+         else
+            istate_fluid = istate
+            not_found = .false.
+         end if
+
+         if(not_found) then
+            FLExit("No fluid phase found for the heat transfer term.")
+         end if
+
+         ! If we have a fluid-particle pair, then assemble the heat transfer term
+         if(is_particle_phase) then
+            call assemble_heat_transfer(istate_fluid, istate)
+         else
+            state_loop: do i = 1, size(state)
+               if(i /= istate_fluid) then
+                  call assemble_heat_transfer(istate_fluid, i)
+               end if
+            end do state_loop
+         end if
+         
+         ewrite(1, *) "Exiting add_heat_transfer"
+         
+         contains
+
+            subroutine assemble_heat_transfer(istate_fluid, istate_particle)
+
+               integer, intent(in) :: istate_fluid, istate_particle
+
+               type(scalar_field), pointer :: vfrac_fluid, vfrac_particle
+               type(scalar_field), pointer :: density_fluid, density_particle
+               type(vector_field), pointer :: velocity_fluid, velocity_particle
+               type(scalar_field), pointer :: internal_energy_fluid, internal_energy_particle
+               type(scalar_field), pointer :: old_internal_energy_fluid, old_internal_energy_particle
+               type(vector_field), pointer :: nu_fluid, nu_particle ! Non-linear approximation to the Velocities
+               type(tensor_field), pointer :: viscosity_fluid
+               type(scalar_field) :: nvfrac_fluid, nvfrac_particle
+               real :: d ! Particle diameter
+               real :: k ! Effective gas conductivity
+               real :: C_fluid, C_particle ! Specific heat of the fluid and particle phases at constant volume
+               real :: gamma ! Ratio of specific heats for compressible phase
+               integer :: kstat, cstat_fluid, cstat_particle, gstat
+
+               ! Get the necessary fields to calculate the heat transfer term
+               velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
+               velocity_particle => extract_vector_field(state(istate_particle), "Velocity")
+               if(.not.aliased(velocity_particle)) then ! Don't count the aliased material_phases
+                  
+                  vfrac_fluid => extract_scalar_field(state(istate_fluid), "PhaseVolumeFraction")
+                  vfrac_particle => extract_scalar_field(state(istate_particle), "PhaseVolumeFraction")
+                  density_fluid => extract_scalar_field(state(istate_fluid), "Density")
+                  density_particle => extract_scalar_field(state(istate_particle), "Density")
+                  viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
+         
+                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
+                           
+                  call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/effective_conductivity", k, kstat)
+                           
+                  call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/specific_heat", C_fluid, cstat_fluid)
+                  
+                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/specific_heat", C_particle, cstat_particle)
+                  
+                  call get_option(trim(state(istate_fluid)%option_path)//"/equation_of_state/compressible/stiffened_gas/ratio_specific_heats", gamma, gstat)
+                           
+                  if(kstat /= 0) then
+                     FLExit("For inter-phase heat transfer, an effective_conductivity needs to be specified for the fluid phase.")
+                  end if                  
+                  if(cstat_fluid /= 0 .or. cstat_particle /= 0) then
+                     FLExit("For inter-phase heat transfer, a specific_heat needs to be specified for each phase.")
+                  end if
+                  if(gstat /= 0) then
+                     FLExit("For inter-phase heat transfer, ratio_specific_heats needs to be specified for the compressible phase.")
+                  end if
+
+                  ! Calculate the non-linear approximation to the PhaseVolumeFractions
+                  call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
+                  call allocate(nvfrac_particle, vfrac_particle%mesh, "NonlinearPhaseVolumeFraction")
+                  call zero(nvfrac_fluid)
+                  call zero(nvfrac_particle)
+                  call get_nonlinear_volume_fraction(state(istate_fluid), nvfrac_fluid)
+                  call get_nonlinear_volume_fraction(state(istate_particle), nvfrac_particle)
+                  
+                  ! Get the non-linear approximation to the Velocities
+                  nu_fluid => extract_vector_field(state(istate_fluid), "NonlinearVelocity")
+                  nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity")
+                  
+                  ! Get the current and old internal energy fields
+                  internal_energy_fluid => extract_scalar_field(state(istate_fluid), "InternalEnergy")
+                  internal_energy_particle => extract_scalar_field(state(istate_particle), "InternalEnergy")
+                  old_internal_energy_fluid => extract_scalar_field(state(istate_fluid), "OldInternalEnergy")
+                  old_internal_energy_particle => extract_scalar_field(state(istate_particle), "OldInternalEnergy")
+      
+                  ! ----- Volume integrals over elements -------------           
+                  call profiler_tic(internal_energy, "element_loop")
+                  element_loop: do ele = 1, element_count(internal_energy)
+
+                     if(.not.dg .or. (dg .and. element_owned(internal_energy,ele))) then
+                        internal_energy_nodes => ele_nodes(internal_energy, ele)
+                        internal_energy_shape => ele_shape(internal_energy, ele)
+                        test_function = internal_energy_shape                         
+
+                        call add_heat_transfer_element(ele, test_function, internal_energy_shape, &
+                                                      x, internal_energy, matrix, rhs, &
+                                                      nvfrac_fluid, nvfrac_particle, &
+                                                      density_fluid, density_particle, &
+                                                      nu_fluid, nu_particle, &
+                                                      internal_energy_fluid, &
+                                                      internal_energy_particle, &
+                                                      old_internal_energy_fluid, &
+                                                      old_internal_energy_particle, &
+                                                      viscosity_fluid, d, k, C_fluid, &
+                                                      C_particle, gamma)
+                     end if
+
+                  end do element_loop
+                  call profiler_toc(internal_energy, "element_loop")
+
+                  call deallocate(nvfrac_fluid)
+                  call deallocate(nvfrac_particle)
+               end if
+
+            end subroutine assemble_heat_transfer
+         
+            subroutine add_heat_transfer_element(ele, test_function, internal_energy_shape, &
+                                                x, internal_energy, matrix, rhs, &
+                                                vfrac_fluid, vfrac_particle, &
+                                                density_fluid, density_particle, &
+                                                nu_fluid, nu_particle, &
+                                                internal_energy_fluid, &
+                                                internal_energy_particle, &
+                                                old_internal_energy_fluid, &
+                                                old_internal_energy_particle, &
+                                                viscosity_fluid, d, k, C_fluid, &
+                                                C_particle, gamma)
+                                                         
+               integer, intent(in) :: ele
+               type(element_type), intent(in) :: test_function
+               type(element_type), intent(in) :: internal_energy_shape
+               type(vector_field), intent(in) :: x
+               type(scalar_field), intent(in) :: internal_energy
+               type(csr_matrix), intent(inout) :: matrix
+               type(scalar_field), intent(inout) :: rhs
+                    
+               type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle
+               type(scalar_field), intent(in) :: density_fluid, density_particle
+               type(vector_field), intent(in) :: nu_fluid, nu_particle
+               type(scalar_field), intent(in) :: internal_energy_fluid, internal_energy_particle
+               type(scalar_field), intent(in) :: old_internal_energy_fluid, old_internal_energy_particle
+               type(tensor_field), intent(in) :: viscosity_fluid    
+               real, intent(in) :: d, k, C_fluid, C_particle, gamma
+               
+               ! Local variables
+               real, dimension(ele_ngi(x,ele)) :: internal_energy_fluid_gi, internal_energy_particle_gi
+               real, dimension(ele_ngi(x,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
+               real, dimension(ele_ngi(x,ele)) :: density_fluid_gi, density_particle_gi
+               real, dimension(x%dim, ele_ngi(x,ele)) :: nu_fluid_gi, nu_particle_gi
+               real, dimension(x%dim, x%dim, ele_ngi(x,ele)) :: viscosity_fluid_gi
+               
+               real, dimension(ele_loc(internal_energy,ele)) :: old_internal_energy_val
+               
+               real, dimension(ele_ngi(x,ele)) :: detwei
+               real, dimension(ele_ngi(x,ele)) :: coefficient_for_matrix, coefficient_for_rhs
+               real, dimension(ele_loc(internal_energy,ele)) :: rhs_addto
+               real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: matrix_addto
+               
+               real, dimension(ele_ngi(x,ele)) :: particle_Re ! Particle Reynolds number
+               real, dimension(ele_ngi(x,ele)) :: Pr ! Prandtl number
+               real, dimension(ele_ngi(x,ele)) :: particle_Nu ! Particle Nusselt number
+               real, dimension(ele_ngi(x,ele)) :: velocity_magnitude ! |v_f - v_p|
+               
+               real, dimension(ele_ngi(x,ele)) :: Q ! heat transfer term = Q*(T_p - T_f)
+               real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: heat_transfer_matrix
+               real, dimension(ele_loc(internal_energy,ele)) :: heat_transfer_rhs
+               
+               integer ::  gi
+               
+               ! Compute detwei
+               call transform_to_physical(x, ele, detwei)
+               
+               ! Get the values of the necessary fields at the Gauss points
+               vfrac_fluid_gi = ele_val_at_quad(vfrac_fluid, ele)
+               vfrac_particle_gi = ele_val_at_quad(vfrac_particle, ele)
+               density_fluid_gi = ele_val_at_quad(density_fluid, ele)
+               density_particle_gi = ele_val_at_quad(density_particle, ele)
+               nu_fluid_gi = ele_val_at_quad(nu_fluid, ele)
+               nu_particle_gi = ele_val_at_quad(nu_particle, ele)
+               viscosity_fluid_gi = ele_val_at_quad(viscosity_fluid, ele)  
+               internal_energy_fluid_gi = ele_val_at_quad(internal_energy_fluid, ele)
+               internal_energy_particle_gi = ele_val_at_quad(internal_energy_particle, ele)
+         
+               ! Compute the magnitude of the relative velocity
+               do gi = 1, ele_ngi(x,ele)
+                  velocity_magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
+               end do
+
+               ! Compute the particle Reynolds number
+               ! (Assumes isotropic viscosity for now)
+               particle_Re = (density_fluid_gi*velocity_magnitude*d) / viscosity_fluid_gi(1,1,:)
+           
+               ! Compute the Prandtl number
+               ! (Assumes isotropic viscosity for now)
+               ! Note: C_fluid (at constant volume) multiplied by gamma = C_fluid at constant pressure
+               Pr = C_fluid*gamma*viscosity_fluid_gi(1,1,:)/k
+               
+               particle_Nu = (7.0 - 10.0*vfrac_fluid_gi + 5.0*vfrac_fluid_gi**2)*(1.0 + 0.7*(particle_Re**0.2)*(Pr**(1.0/3.0))) + &
+                              (1.33 - 2.4*vfrac_fluid_gi + 1.2*vfrac_fluid_gi**2)*(particle_Re**0.7)*(Pr**(1.0/3.0))
+
+               Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d**2)
+               
+               ! Note that the transfer term is defined in terms of temperatures (T_fluid and T_particle)
+               ! Let's convert the temperatures to internal energy (E) using E_i = C_i*T_i,
+               ! where C is the specific heat of phase i at constant volume.
+               if(is_particle_phase) then
+                  coefficient_for_matrix = -Q/C_particle
+                  coefficient_for_rhs = -Q*(-internal_energy_fluid_gi/C_fluid)
+               else
+                  coefficient_for_matrix = -Q/C_fluid
+                  coefficient_for_rhs = Q*(internal_energy_particle_gi/C_particle)
+               end if
+               
+               ! Form the element heat transfer matrix and RHS
+               heat_transfer_matrix = shape_shape(test_function, internal_energy_shape, detwei*coefficient_for_matrix)
+               heat_transfer_rhs = shape_rhs(test_function, coefficient_for_rhs*detwei)
+              
+               ! Add contribution  
+               matrix_addto = 0.0            
+               rhs_addto = 0.0
+               if(is_particle_phase) then
+                  old_internal_energy_val = ele_val(old_internal_energy_particle, ele)
+               else
+                  old_internal_energy_val = ele_val(old_internal_energy_fluid, ele)
+               end if
+
+               matrix_addto = matrix_addto - dt*theta*heat_transfer_matrix
+               rhs_addto = rhs_addto + matmul(heat_transfer_matrix, old_internal_energy_val) + heat_transfer_rhs
+               
+               ! Add to the internal energy equation's RHS
+               call addto(rhs, internal_energy_nodes, rhs_addto) 
+               ! Add to the matrix
+               call addto(matrix, internal_energy_nodes, internal_energy_nodes, matrix_addto)
+
+            end subroutine add_heat_transfer_element
+
+      end subroutine add_heat_transfer
+      
    end module multiphase_module
 
=== modified file 'main/Fluids.F90'
--- main/Fluids.F90	2012-05-29 13:43:29 +0000
+++ main/Fluids.F90	2012-08-09 00:46:20 +0000
@@ -699,7 +699,7 @@
                 else if(have_option(trim(field_optionpath_list(it)) // &
                      & "/prognostic/spatial_discretisation/continuous_galerkin")) then
 
-                   call solve_field_equation_cg(field_name_list(it), state(field_state_list(it)), dt)
+                   call solve_field_equation_cg(field_name_list(it), state, field_state_list(it), dt)
                 else
 
                    ewrite(2, *) "Not solving scalar field " // trim(field_name_list(it)) // " in state " // trim(state(field_state_list(it))%name) //" in an advdif-like subroutine."
 
=== modified file 'main/Shallow_Water.F90'
--- main/Shallow_Water.F90	2011-07-22 12:39:08 +0000
+++ main/Shallow_Water.F90	2012-08-09 00:46:20 +0000
@@ -241,7 +241,7 @@
           call project_cartesian_to_local(state(1), v_field)
        end if 
 
-       call execute_timestep(state(1), dt)
+       call execute_timestep(state, 1, dt)
 
        call set_prescribed_field_values(state, exclude_interpolated=.true., &
             exclude_nonreprescribed=.true., time=current_time + dt)
@@ -528,9 +528,10 @@
 
     end subroutine insert_time_in_state
 
-    subroutine execute_timestep(state, dt)
+    subroutine execute_timestep(state, istate, dt)
       implicit none
-      type(state_type), intent(inout) :: state
+      type(state_type), dimension(:), intent(inout) :: state
+      integer, intent(in) :: istate
       real, intent(in) :: dt
 
       !! Layer thickness
@@ -548,20 +549,20 @@
       real :: energy
       logical :: have_source
 
-      !Pull the fields out of state
-      D=>extract_scalar_field(state, "LayerThickness")
-      U=>extract_vector_field(state, "LocalVelocity")
+      !Pull the fields out of state(istate)
+      D=>extract_scalar_field(state(istate), "LayerThickness")
+      U=>extract_vector_field(state(istate), "LocalVelocity")
       have_source = .false.
-      if (has_vector_field(state, "LocalVelocitySource")) then
+      if (has_vector_field(state(istate), "LocalVelocitySource")) then
         have_source = .true.
-        source => extract_vector_field(state, "LocalVelocitySource")
+        source => extract_vector_field(state(istate), "LocalVelocitySource")
       end if
       dim = U%dim
 
       call execute_timestep_setup(D,U,d_rhs,u_rhs,advecting_u, &
            old_u,old_d,delta_d,delta_u)
 
-      call insert(state,advecting_u,"NonlinearVelocity")
+      call insert(state(istate),advecting_u,"NonlinearVelocity")
 
       call set(advecting_u,u)
       call set(old_u,u)
@@ -577,20 +578,20 @@
             do d1 = 1, dim
                velocity_cpt = extract_scalar_field(u,d1)
                old_velocity_cpt = extract_scalar_field(old_u,d1)
-               call insert(state,velocity_cpt,"VelocityComponent")
-               call insert(state,old_velocity_cpt,"VelocityComponentOld")
-               call solve_advection_diffusion_dg("VelocityComponent", state)
+               call insert(state(istate),velocity_cpt,"VelocityComponent")
+               call insert(state(istate),old_velocity_cpt,"VelocityComponentOld")
+               call solve_advection_diffusion_dg("VelocityComponent", state(istate))
             end do
          end if
          !pressure advection
          if(.not.exclude_pressure_advection) then
-            call solve_field_equation_cg("LayerThickness", state, dt, &
+            call solve_field_equation_cg("LayerThickness", state, 1, dt, &
                  "NonlinearVelocity")
          end if
 
          if(hybridized) then
             call solve_hybridized_helmholtz(&
-                 &state,&
+                 &state(istate),&
                  &U_out=U_rhs,D_out=D_rhs,&
                  &compute_cartesian=.true.,&
                  &check_continuity=.true.,output_dense=.false.)
@@ -623,8 +624,8 @@
             !Construct explicit parts of h rhs in wave equation
             call get_d_rhs(d_rhs,u_rhs,D,U,div_mat,big_mat,D0,dt,theta)
 
-            if (has_scalar_field(state, "LayerThicknessSource")) then
-               d_src => extract_scalar_field(state, "LayerThicknessSource")
+            if (has_scalar_field(state(istate), "LayerThicknessSource")) then
+               d_src => extract_scalar_field(state(istate), "LayerThicknessSource")
                call allocate(md_src, d_src%mesh, "MassMatrixTimesLayerThicknessSource")
                call zero(md_src)
                call mult(md_src, h_mass_mat, d_src)
 
=== modified file 'manual/bibliography.bib'
--- manual/bibliography.bib	2011-11-25 12:32:35 +0000
+++ manual/bibliography.bib	2012-08-09 00:46:20 +0000
@@ -912,6 +912,16 @@
   year		= {2005}
 }
 
+@Article{         ergun1952,
+  author        = {Ergun, S.},
+  journal       = {Chemical Engineering Progress},
+  number        = {2},
+  pages         = {89-94},
+  title         = {Fluid flow through packed columns},
+  volume        = {48},
+  year          = {1952}
+}
+
 @Article{	  erturk2005,
   author	= {E. Erturk and T. C. Corke and C. Gokcol},
   journal	= {International Journal for Numerical Methods in Fluids},
@@ -1298,6 +1308,16 @@
   year		= {1994}
 }
 
+@Article{         gunn1978,
+  author        = {Gunn, D. J.},
+  journal       = {International Journal of Heat and Mass Transfer},
+  number        = {4},
+  pages         = {467-476},
+  title         = {Transfer of heat or mass to particles in fixed and fluidised beds},
+  volume        = {21},
+  year          = {1978}
+}
+
 @Book{		  haidvogelbook,
   author	= {D. B. Haidvogel and A. Beckmann},
   publisher	= {Imperial {C}ollege {P}ress},
@@ -2877,6 +2897,15 @@
   year		= 2008
 }
 
+@Article{ wen_yu_1966,
+  title         = {Mechanics of Fluidization},
+  author        = {C.~Y. Wen and Y.~H. Yu},
+  journal       = {Chem. Eng. Prog. Symp. Ser.},
+  volume        = 62,
+  number        = 100,
+  year          = 1966
+}
+
 @Book{            wilcox1998turbulence,
   author	= {Wilcox, D.C.},
   publisher	= {DCW industries La Canada, CA},
 
=== modified file 'manual/configuring_fluidity.tex'
--- manual/configuring_fluidity.tex	2012-07-26 15:38:42 +0000
+++ manual/configuring_fluidity.tex	2012-08-09 00:46:20 +0000
@@ -719,7 +719,7 @@
 
 When configuring ocean problems (or single material/single phase fluids
 problems), only one \option{/material\_phase} is required.  Multi-material
-and multi-phase problems will require one \option{/material\_phase} for each
+and multiphase problems will require one \option{/material\_phase} for each
 phase or material in the problem.
 
 Note that you must give each of your material phases a name.
@@ -1184,7 +1184,8 @@
 \index{conjugate gradient}
 \index{multigrid}
 \index{PETSc}
-As described in Section \ref{ND_Linear_solvers}, for the solution of large sparse linear systems, the so called iterative methods are usually employed. These methods avoid having to explicitly construct the inverse of the matrix, which is generally dense and therefore costly to compute (both in memory and computer time). FLUIDITY is linked to PETSc: a suite of data structures and routines for the scalable (parallel) solution of scientific applications modelled by partial differential equations.  It employs the MPI standard for parallelism. FLUIDITY therefore supports any iterative method provided by the PETSc library (http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html --- available methods may depend on the PETSc library installed on your system). Examples include Conjugate Gradient (CG), GMRES and FGMRES (Flexible GMRES). Some options are explicitly listed under \option{solver/iterative\_method}, for example CG: \option{solver/iterative\_method::cg}, whereas others can be selected entering the name of the chosen scheme in \option{solver/iterative\_method}.
+As described in Section \ref{ND_Linear_solvers}, for the solution of large sparse linear systems, the so called iterative methods are usually employed. These methods avoid having to explicitly construct the inverse of the matrix, which is generally dense and therefore costly to compute (both in memory and computer time). FLUIDITY is linked to PETSc: a suite of data structures and routines for the scalable (parallel) solution of scientific applications modelled by partial differential equations.  It employs the MPI standard for parallelism. FLUIDITY therefore supports any iterative method provided by the PETSc library (http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html --- available methods may depend on the PETSc library installed on your system). Examples include Conjugate Gradient (CG), GMRES and FGMRES (Flexible GMRES). Some options are explicitly listed under \option{solver/iterative\_method}, for example CG: \option{solver/iterative\_method::cg}, whereas 
+others can be selected entering the name of the chosen scheme in \option{solver/iterative\_method}.
 
 \subsection{Preconditioner}
 \index{linear solvers!preconditioners}
@@ -2313,7 +2314,7 @@
 Models with one prognostic velocity field per \option{material\_phase} are referred to as multiple phase models. The use of these multiple velocity fields permits the inter-penetration and interaction between different phases.
 
 \subsubsection{Simulation requirements}
-In a multi-phase simulation, each \option{material\_phase} requires:
+In a multiphase simulation, each \option{material\_phase} requires:
 \begin{itemize}
  \item an equation of state, and
  \item a PhaseVolumeFraction scalar field.
@@ -2328,13 +2329,15 @@
 \alpha_N = 1 - \sum_{i=1}^{N-1}\alpha_i\textrm{.}
 \end{equation}
 
-\subsubsection{Inter-phase interactions}
-Inter-phase interactions can be included under \option{../multiphase\_properties} in Diamond. Currently, \fluidity\ only supports fluid-particle drag between the continuous phase and dispersed phase(s), given by:
-\begin{equation}\label{eq:multiphase_drag_term}
+\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}.
+
+The drag force using the Stokes drag correlation is given by:
+\begin{equation}\label{eq:stokes_drag_force}
 \mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d},
 \end{equation}
 where $f$ and $p$ denote the fluid (i.e. continuous) and particle (i.e. dispersed) phases respectively, and $d$ is the diameter of a single particle in the dispersed phase. The drag coefficient $C_D$ is defined as:
-\begin{equation}\label{eq:drag_coefficient}
+\begin{equation}\label{eq:stokes_drag_coefficient}
 C_D = \frac{24}{\mathrm{Re}},
 \end{equation}
 with
@@ -2343,17 +2346,52 @@
 \end{equation}
 Note that $\mu_f$ denotes the isotropic viscosity of the fluid (i.e. continuous) phase.
 
-Within a dispersed phase, entering a value for $d$ in the \option{../multiphase\_properties/particle\_diameter} option enables fluid-particle interaction between itself and the continuous phase.
+With the drag correlation by \cite{wen_yu_1966}, $\mathbf{F}_D$ and $C_D$ become:
+\begin{equation}\label{eq:wen_yu_drag_force}
+\mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d\alpha_f^{2.7}},
+\end{equation}
+\begin{equation}\label{eq:wen_yu_drag_coefficient}
+C_D = \frac{24}{\mathrm{Re}}\left(1.0 + 0.15\mathrm{Re}^{0.687}\right).
+\end{equation}
+In contrast to the Stokes drag correlation, the \cite{wen_yu_1966} drag correlation is more suitable for larger particle Reynolds number flows. 
+
+For dense multiphase flows with $\alpha_p > 0.2$, the drag correlation by \cite{ergun1952} is often the most appropriate:
+\begin{equation}\label{eq:ergun_drag_force}
+\mathbf{F}_D = \left(150\frac{\alpha_p^2\mu_f}{\alpha_f d_p^2} + 1.75\frac{\alpha_p\rho_f|\mathbf{u}_f-\mathbf{u}_p|}{d_p}\right)\left(\mathbf{u}_f-\mathbf{u}_p\right).
+\end{equation}
+
+Note that within each dispersed phase, a value for $d$ must be specified in \option{../multiphase\_properties/particle\_diameter} regardless of which drag correlation is used.
+
+\subsubsection{Inter-phase energy interaction}
+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),
+\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},
+\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.
+
+To include this term, switch on the \option{/multiphase\_interaction/heat\_transfer} option in Diamond.
 
 \subsubsection{Current limitations}
 \begin{itemize}
- \item Boussinesq and fully compressible multi-phase simulations are not yet supported.
+ \item Boussinesq multiphase simulations are not yet supported.
  \item The momentum equation for each \option{material\_phase} can only be discretised in non-conservative form.
  \item The stress term must be in tensor form, and only works with an isotropic viscosity.
- \item \option{bassi\_rebay} and \option{compact\_discontinuous\_galerkin} are currently the only DG viscosity schemes available for multi-phase flow simulations.
+ \item \option{bassi\_rebay} and \option{compact\_discontinuous\_galerkin} are currently the only DG viscosity schemes available for multiphase flow simulations.
  \item Discontinuous \option{PhaseVolumeFraction} fields are not yet supported.
- \item The Pressure field only supports \option{continuous\_galerkin} and \option{control\_volume} discretisations.
- \item Prescribed velocity fields cannot yet be used in multi-phase simulations.
+ \item For compressible multiphase flow simulations, the Pressure field only supports a \option{continuous\_galerkin} discretisation.
+ \item Prescribed velocity fields cannot yet be used in multiphase simulations.
  \item Fluid-particle drag can currently only support one continuous (i.e. fluid) phase.
 \end{itemize}
 
@@ -2396,7 +2434,8 @@
 
 This section describes how to configure \fluidity to simulate single phase incompressible Darcy flow.
 
-First, the \option{Porosity} and \option{Permeability} fields found under the option element \option{/porous\_media} require configuring. Both fields can be either \option{prescribed} or \option{diagnostic}. The \option{Permeability} field can be either a \option{scalar\_field} or a \option{vector\_field}. The \option{Porosity} field will only be used in the transport of scalar fields and the metric tensor. The \option{Permeability} field will only be used in forming the absorption coefficient in the Darcy velocity equation. The use of these two fields should be considered when selecting their associated mesh. If \option{Porosity} is to be included in the transport of scalar fields (and the metric tensor) that are discretised using control volumes, due to the way this method is coded, the mesh associated with the \option{Porosity} must be either element wise (meaning discontinuous order zero) or the same order as field to be transported. It is therefore recommended that the \option{Porosity} field always be associated with a mesh that is discontinuous order zero. It is recommended that the \option{Permeability} field also be associated with a mesh that is discontinuous order zero, due to the recommended Darcy velocity - pressure element pair (being p0p1cv or p0p1 with the continuity tested with the p1 control volume dual space)
+First, the \option{Porosity} and \option{Permeability} fields found under the option element \option{/porous\_media} require configuring. Both fields can be either \option{prescribed} or \option{diagnostic}. The \option{Permeability} field can be either a \option{scalar\_field} or a \option{vector\_field}. The \option{Porosity} field will only be used in the transport of scalar fields and the metric tensor. The \option{Permeability} field will only be used in forming the absorption coefficient in the Darcy velocity equation. The use of these two fields should be considered when selecting their associated mesh. If \option{Porosity} is to be included in the transport of scalar fields (and the metric tensor) that are discretised using control volumes, due to the way this method is coded, the mesh associated with the \option{Porosity} must be either element wise (meaning discontinuous order zero) or the same order as field to be transported. It is therefore recommended that the \option{Porosity} field always be 
+associated with a mesh that is discontinuous order zero. It is recommended that the \option{Permeability} field also be associated with a mesh that is discontinuous order zero, due to the recommended Darcy velocity - pressure element pair (being p0p1cv or p0p1 with the continuity tested with the p1 control volume dual space)
 
 Second, the prognostic \option{Velocity} vector field of the (only) phase is now taken to be the Darcy velocity \ref{sec:porous_media_darcy_flow_equations}. There is no option to select Darcy flow but this can be achieved via a careful selection of the \option{Velocity} field options tree given by:
 \begin{itemize}
 
=== modified file 'manual/model_equations.tex'
--- manual/model_equations.tex	2012-07-26 15:38:42 +0000
+++ manual/model_equations.tex	2012-08-09 00:46:20 +0000
@@ -727,7 +727,7 @@
 
 \subsection{Multi-material simulations}
 \index{multi-material flow}
-The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multi-phase approach. The multi-material approach is implemented within \fluidity, and the multi-phase approach (discussed in the next subsection) is currently under development.
+The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multiphase approach. The multi-material approach is implemented within \fluidity, and the multiphase approach (discussed in the next subsection) is currently under development.
 In situations where the model can resolve physical mixing of immiscible materials, or where there is no mixing, only one velocity field (and hence one momentum equation) is required to describe the flow of all materials. The \emph{multi-material} approach, considers all materials to be immiscible materials separated by a sharp interface.
 
 In a multi-material approach, the various forms of the conservation equations can be solved for multiple material flows if the point-wise mass density $\rho$ in the equations is defined as the bulk density at each point.  If the flow comprises $n$ materials and the volume fraction of the $i^{th}$ material is denoted $\phi_i$ then the bulk density is given by:
@@ -749,14 +749,17 @@
 \end{equation}
 where the volume fraction field at time zero must be specified.
 
-\subsection{Multi-phase simulations}
+\subsection{Multiphase simulations}
 \label{sec:multiphase_equations}
-\index{multi-phase flow}
-Multi-phase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body.
+\index{multiphase flow}
+Multiphase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body.
 
 Further to the above definition, each phase is classed as either \textit{continuous} or \textit{dispersed}, where a continuous phase is a connected liquid or gas substance in which dispersed phases (comprising a finite number of solid particles, liquid droplets and/or gas bubbles) may be immersed \citep{croweEtAl1998}.
 
-To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. Thus, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are:
+To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. The form of the governing equations depends on whether we are dealing with incompressible or compressible flow.
+
+\subsubsection{Incompressible flow}
+For an incompressible multiphase flow, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are:
 \begin{equation}
 \sum_{i=1}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) = 0,
 \end{equation}
@@ -769,6 +772,19 @@
 \end{equation}
 where $\mathbf{u}_i$, $\rho_i$, $\mu_i$ and $\alpha_i$ are the velocity, density, isotropic viscosity 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}.
 
+\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:
+\begin{equation}
+\alpha_1\frac{\partial\rho_1}{\partial t} + \nabla\cdot\left(\alpha_1\rho_1\mathbf{u}_1\right) + \rho_1\left(\sum_{i=2}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right)\right) = 0,
+\end{equation}
+where the compressible phase has an index of $i=1$ here.
+
+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
+\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}.
+
 \subsection{Porous Media Darcy Flow}
 \label{sec:porous_media_darcy_flow_equations}
 \index{porous media darcy flow}
 
=== modified file 'manual/visualisation_and_diagnostics.tex'
--- manual/visualisation_and_diagnostics.tex	2012-03-06 10:47:48 +0000
+++ manual/visualisation_and_diagnostics.tex	2012-08-09 00:46:20 +0000
@@ -66,6 +66,7 @@
 Limitations: Requires the diagnostic \option{scalar\_field::IsopycnalCoordinate} and (therefore) is not parallelised.
 \item[BulkMaterialPressure:]Calculates the bulk material pressure based on the MaterialDensity and MaterialVolumeFraction (and MaterialInternalEnergy if appropriate) for the equation of state of all materials.
 \item[CFLNumber:]CFLNumber as defined on the co-ordinate mesh. It is calculated as $\triangle t \vec{u} \mat{J}^{-1}$ where $\triangle t$ is the timestep, $\vec{u}$ the velocity and $\mat{J}$ the Jacobian. 
+\item[CompressibleContinuity:]Computes the residual of the compressible multiphase continuity equation. Used only in compressible multiphase flow simulations. 
 \item[ControlVolumeCFLNumber:]CFL Number as defined on a control volume mesh. It is calculated as $\triangle t \frac{1}{V} \int _{cv_{face}} \vec{u}$ where $\triangle t$ is the timestep and $\vec{u}$ the velocity. The integral is taken over the faces of the control volume and $V$ is the volume of the control volume.
 \item[ControlVolumeDivergence:]Divergence of the velocity field where the divergence operator is defined using the control volume $\mathrm{C}^\mathrm{T}$ matrix. This assumes that the test space is discontinuous control volumes.
 \item[CVMaterialDensityCFLNumber:]Courant Number as defined on a control volume mesh and incorporating the MaterialDensity. Requires a MaterialDensity field!
 
=== modified file 'schemas/fluidity_options.rnc'
--- schemas/fluidity_options.rnc	2012-06-27 14:47:10 +0000
+++ schemas/fluidity_options.rnc	2012-08-09 00:46:20 +0000
@@ -15,6 +15,7 @@
 include "embedded_models.rnc"
 include "flredecomp.rnc"
 include "porous_media.rnc"
+include "multiphase_interaction.rnc"
 include "equation_of_state.rnc"
 
 start =
@@ -1387,6 +1388,17 @@
                comment,
                element particle_diameter { 
                      real 
+               }?,
+               ## If this is the fluid phase, 
+               ## the effective conductivity is required 
+               ## for the inter-phase heat transfer term.
+               element effective_conductivity { 
+                     real 
+               }?,
+               ## The specific heat (at constant volume) is required 
+               ## for the inter-phase heat transfer term.
+               element specific_heat { 
+                     real
                }?
             }?,
             sediment?
@@ -1570,7 +1582,8 @@
          reduced_model?,
          porous_media_model?,
          embedded_models?,
-         flredecomp?
+         flredecomp?,
+         multiphase_interaction?
       }
    )      
 
@@ -4681,7 +4694,27 @@
                   }
                }
             )
-         }              
+         }|
+         
+         ## CompressibleContinuityResidual
+         ##
+         ## Computes the residual of the continuity equation used in compressible multiphase flow simulations
+         ## i.e. vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) }
+         ## where _c and _i denote the compressible and incompressible phases respectively.
+         element scalar_field {
+            attribute rank { "0" },
+            attribute name { "CompressibleContinuityResidual" },
+            (
+               element diagnostic {
+                  velocity_mesh_choice,
+                  internal_algorithm,
+                  diagnostic_scalar_field_no_adapt,
+                  element solver {
+                     linear_solver_options_sym
+                  }
+               }
+            )
+         }    
 
 # Insert new diagnostic scalar fields here using the template:
 #        element scalar_field {
 
=== modified file 'schemas/fluidity_options.rng'
--- schemas/fluidity_options.rng	2012-08-02 13:03:52 +0000
+++ schemas/fluidity_options.rng	2012-08-09 00:46:20 +0000
@@ -16,6 +16,7 @@
   <include href="embedded_models.rng"/>
   <include href="flredecomp.rng"/>
   <include href="porous_media.rng"/>
+  <include href="multiphase_interaction.rng"/>
   <include href="equation_of_state.rng"/>
   <start>
     <element name="fluidity_options">
@@ -1852,6 +1853,21 @@
                   <ref name="real"/>
                 </element>
               </optional>
+              <optional>
+                <element name="effective_conductivity">
+                  <a:documentation>If this is the fluid phase, 
+the effective conductivity is required 
+for the inter-phase heat transfer term.</a:documentation>
+                  <ref name="real"/>
+                </element>
+              </optional>
+              <optional>
+                <element name="specific_heat">
+                  <a:documentation>The specific heat (at constant volume) is required 
+for the inter-phase heat transfer term.</a:documentation>
+                  <ref name="real"/>
+                </element>
+              </optional>
             </element>
           </optional>
           <optional>
@@ -2185,6 +2201,9 @@
       <optional>
         <ref name="flredecomp"/>
       </optional>
+      <optional>
+        <ref name="multiphase_interaction"/>
+      </optional>
     </element>
   </start>
   <define name="sediment">
@@ -5935,6 +5954,27 @@
           </element>
         </element>
       </element>
+      <element name="scalar_field">
+        <a:documentation>CompressibleContinuityResidual
+
+Computes the residual of the continuity equation used in compressible multiphase flow simulations
+i.e. vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) }
+where _c and _i denote the compressible and incompressible phases respectively.</a:documentation>
+        <attribute name="rank">
+          <value>0</value>
+        </attribute>
+        <attribute name="name">
+          <value>CompressibleContinuityResidual</value>
+        </attribute>
+        <element name="diagnostic">
+          <ref name="velocity_mesh_choice"/>
+          <ref name="internal_algorithm"/>
+          <ref name="diagnostic_scalar_field_no_adapt"/>
+          <element name="solver">
+            <ref name="linear_solver_options_sym"/>
+          </element>
+        </element>
+      </element>
     </choice>
     <!--
       Insert new diagnostic scalar fields here using the template:
 
=== added file 'schemas/multiphase_interaction.rnc'
--- schemas/multiphase_interaction.rnc	1970-01-01 00:00:00 +0000
+++ schemas/multiphase_interaction.rnc	2012-08-09 00:46:20 +0000
@@ -0,0 +1,58 @@
+# Phase interaction options for Fluidity's multiphase flow model.
+
+multiphase_interaction = 
+   (
+      ## Phase interaction options for Fluidity's multiphase flow model.
+      element multiphase_interaction {
+         comment,
+         fluid_particle_drag?,
+         heat_transfer?
+      }
+   ) 
+
+fluid_particle_drag =
+   (
+      ## Fluid-particle drag term.
+      element fluid_particle_drag {
+         drag_correlation
+      }
+   )
+
+drag_correlation =
+   (
+      ## Stokes drag correlation: 24/Re_p
+      ## where Re_p is the particle Reynolds number.
+      element drag_correlation {
+         attribute name { "stokes" },
+         comment
+      }|
+      ## Fluid-particle drag term by Wen & Yu (1966).
+      element drag_correlation {
+         attribute name { "wen_yu" },
+         comment
+      }|
+      ## Fluid-particle drag term by Ergun (1952).
+      element drag_correlation {
+         attribute name { "ergun" },
+         comment
+      }
+   )
+   
+heat_transfer =
+   (
+      ## Heat transfer term for the 
+      ## multiphase internal energy equation.
+      ## Note: Only for fluid-particle phase pairs.
+      element heat_transfer {
+         heat_transfer_coefficient
+      }
+   )
+
+heat_transfer_coefficient =
+   (
+      ## Heat transfer coefficient by Gunn (1978).
+      element heat_transfer_coefficient {
+         attribute name { "gunn" },
+         comment
+      }
+   )
 \ No newline at end of file
 
=== added file 'schemas/multiphase_interaction.rng'
--- schemas/multiphase_interaction.rng	1970-01-01 00:00:00 +0000
+++ schemas/multiphase_interaction.rng	2012-08-09 00:46:20 +0000
@@ -0,0 +1,65 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<grammar xmlns:a="http://relaxng.org/ns/compatibility/annotations/1.0" xmlns="http://relaxng.org/ns/structure/1.0">
+  <!-- Phase interaction options for Fluidity's multiphase flow model. -->
+  <define name="multiphase_interaction">
+    <element name="multiphase_interaction">
+      <a:documentation>Phase interaction options for Fluidity's multiphase flow model.</a:documentation>
+      <ref name="comment"/>
+      <optional>
+        <ref name="fluid_particle_drag"/>
+      </optional>
+      <optional>
+        <ref name="heat_transfer"/>
+      </optional>
+    </element>
+  </define>
+  <define name="fluid_particle_drag">
+    <element name="fluid_particle_drag">
+      <a:documentation>Fluid-particle drag term.</a:documentation>
+      <ref name="drag_correlation"/>
+    </element>
+  </define>
+  <define name="drag_correlation">
+    <choice>
+      <element name="drag_correlation">
+        <a:documentation>Stokes drag correlation: 24/Re_p
+where Re_p is the particle Reynolds number.</a:documentation>
+        <attribute name="name">
+          <value>stokes</value>
+        </attribute>
+        <ref name="comment"/>
+      </element>
+      <element name="drag_correlation">
+        <a:documentation>Fluid-particle drag term by Wen &amp; Yu (1966).</a:documentation>
+        <attribute name="name">
+          <value>wen_yu</value>
+        </attribute>
+        <ref name="comment"/>
+      </element>
+      <element name="drag_correlation">
+        <a:documentation>Fluid-particle drag term by Ergun (1952).</a:documentation>
+        <attribute name="name">
+          <value>ergun</value>
+        </attribute>
+        <ref name="comment"/>
+      </element>
+    </choice>
+  </define>
+  <define name="heat_transfer">
+    <element name="heat_transfer">
+      <a:documentation>Heat transfer term for the 
+multiphase internal energy equation.
+Note: Only for fluid-particle phase pairs.</a:documentation>
+      <ref name="heat_transfer_coefficient"/>
+    </element>
+  </define>
+  <define name="heat_transfer_coefficient">
+    <element name="heat_transfer_coefficient">
+      <a:documentation>Heat transfer coefficient by Gunn (1978).</a:documentation>
+      <attribute name="name">
+        <value>gunn</value>
+      </attribute>
+      <ref name="comment"/>
+    </element>
+  </define>
+</grammar>
 
=== added directory 'tests/flux_bc'
=== added file 'tests/flux_bc/Makefile'
--- tests/flux_bc/Makefile	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/Makefile	2012-08-09 00:46:20 +0000
@@ -0,0 +1,22 @@
+preprocess:
+	@echo **********Creating meshes
+	gmsh -2 -o flux_bc_2d.msh src/flux_bc_2d.geo
+	gmsh -3 -o flux_bc_3d.msh src/flux_bc_3d.geo
+	../../bin/gmsh2triangle --2d flux_bc_2d.msh
+	../../bin/gmsh2triangle flux_bc_3d.msh
+
+run:
+	@echo **********Running 2D simulation
+	../../bin/fluidity flux_bc_2d.flml 
+	@echo **********Running 3D simulation
+	../../bin/fluidity flux_bc_3d.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.face
+	rm -f matrixdump* *checkpoint* *.log* *.err*
+
 
=== added file 'tests/flux_bc/flux_bc.xml'
--- tests/flux_bc/flux_bc.xml	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/flux_bc.xml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,44 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<!DOCTYPE testproblem SYSTEM "regressiontest.dtd">
+
+<testproblem>
+   
+  <name>flux_bc</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+    <!-- This checks the integral of the PhaseVolumeFraction field at t = 1 second. -->
+    <!-- The volumetric flux (r) is set to 0.5 m^3/s/m^2, and the flux boundary condition -->
+    <!-- enforces d/dt \int_{volume} vfrac = \int_{surface} r -->
+  </problem_definition>
+  
+  <variables>
+    <variable name="vfrac_integral_2d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("flux_bc_2d.stat")
+vfrac_integral_2d = s["Tephra"]["PhaseVolumeFraction"]["integral"]
+    </variable>
+    
+    <variable name="vfrac_integral_3d" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("flux_bc_3d.stat")
+vfrac_integral_3d = s["Tephra"]["PhaseVolumeFraction"]["integral"]
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Integral of PhaseVolumeFraction after 1 second = 0.15 in 2D simulation." language="python">
+assert abs(vfrac_integral_2d[10] - 0.15) &lt; 1.0e-9
+    </test>
+
+    <test name="Integral of PhaseVolumeFraction after 1 second = 0.045 in 3D simulation." language="python">
+assert abs(vfrac_integral_3d[10] - 0.045) &lt; 1.0e-9
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>
 
=== added file 'tests/flux_bc/flux_bc_2d.flml'
--- tests/flux_bc/flux_bc_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/flux_bc_2d.flml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,396 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">flux_bc_2d</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="flux_bc_2d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1.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">0.1</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-12</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.8</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="2" dim1="dim" rank="1">0.0 -1.0</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Tephra">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2340</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <reference_node>
+          <integer_value rank="0">1</integer_value>
+        </reference_node>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="CoordinateMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <bassi_rebay/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Sides">
+          <surface_ids>
+            <integer_value shape="1" rank="1">666</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Bottom">
+          <surface_ids>
+            <integer_value shape="1" rank="1">444</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">333</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </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.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FiniteElement">
+              <limit_face_value>
+                <limiter name="Sweby"/>
+              </limit_face_value>
+            </face_value>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">333</integer_value>
+          </surface_ids>
+          <type name="flux">
+            <constant>
+              <real_value rank="0">0.5</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Bottom">
+          <surface_ids>
+            <integer_value shape="1" rank="1">444</integer_value>
+          </surface_ids>
+          <type name="zero_flux"/>
+        </boundary_conditions>
+        <boundary_conditions name="Sides">
+          <surface_ids>
+            <integer_value shape="1" rank="1">666</integer_value>
+          </surface_ids>
+          <type name="zero_flux"/>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/flux_bc/flux_bc_3d.flml'
--- tests/flux_bc/flux_bc_3d.flml	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/flux_bc_3d.flml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,405 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">flux_bc_3d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="flux_bc_3d">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1.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">0.1</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-12</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.8</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="3" dim1="dim" rank="1">0.0 -1.0 0.0</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="Tephra">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2340.0</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <reference_node>
+          <integer_value rank="0">1</integer_value>
+        </reference_node>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="CoordinateMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <viscosity_scheme>
+              <bassi_rebay/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0.5</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">444</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Sides_X">
+          <surface_ids>
+            <integer_value shape="1" rank="1">777</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Sides_Z">
+          <surface_ids>
+            <integer_value shape="1" rank="1">666</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <z_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </z_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Bottom">
+          <surface_ids>
+            <integer_value shape="1" rank="1">333</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <apply_weakly/>
+            <align_bc_with_cartesian>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </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.001</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FiniteElement">
+              <limit_face_value>
+                <limiter name="Sweby"/>
+              </limit_face_value>
+            </face_value>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">444</integer_value>
+          </surface_ids>
+          <type name="flux">
+            <constant>
+              <real_value rank="0">0.5</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Others">
+          <surface_ids>
+            <integer_value shape="3" rank="1">333 666 777</integer_value>
+          </surface_ids>
+          <type name="zero_flux"/>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>
 
=== added directory 'tests/flux_bc/src'
=== added file 'tests/flux_bc/src/flux_bc_2d.geo'
--- tests/flux_bc/src/flux_bc_2d.geo	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/src/flux_bc_2d.geo	2012-08-09 00:46:20 +0000
@@ -0,0 +1,19 @@
+Point(1) = {0.0, 0.0, 0.0, 0.05};
+Point(2) = {0.0, 1.0, 0.0, 0.05};
+Point(3) = {0.3, 1.0, 0.0, 0.05};
+Point(4) = {0.3, 0.0, 0.0, 0.05};
+Line(1) = {3,4};
+Line(2) = {4,1};
+Line(3) = {1,2};
+Line(4) = {2,3};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+// Top of the box
+Physical Line(333) = {4};
+// Box sides
+Physical Line(666) = {3,1};
+// Box bottom
+Physical Line(444) = {2};
+// This is just to ensure all the interior
+// elements get written out. 
+Physical Surface(10) = {6};
 
=== added file 'tests/flux_bc/src/flux_bc_3d.geo'
--- tests/flux_bc/src/flux_bc_3d.geo	1970-01-01 00:00:00 +0000
+++ tests/flux_bc/src/flux_bc_3d.geo	2012-08-09 00:46:20 +0000
@@ -0,0 +1,69 @@
+Point(1) = {0.0, 0.0, 0.0, 0.05};
+Point(2) = {0.3, 0.0, 0.0, 0.05};
+Point(3) = {0.3, 0.0, 0.3, 0.05};
+Point(4) = {0.0, 0.0, 0.3, 0.05};
+Point(5) = {0.0, 1.0, 0.0, 0.05};
+Point(6) = {0.3, 1.0, 0.0, 0.05};
+Point(7) = {0.3, 1.0, 0.3, 0.05};
+Point(8) = {0.0, 1.0, 0.3, 0.05};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Line(5) = {5,6};
+Line(6) = {6,7};
+Line(7) = {7,8};
+Line(8) = {8,5};
+
+Line(9) = {5,1};
+Line(10) = {6,2};
+Line(11) = {7,3};
+Line(12) = {8,4};
+
+// Bottom surface
+Line Loop(5) = {1,2,3,4};
+
+// Side surface, nearest to the x axis
+Line Loop(6) = {1,-10,-5,9};
+
+// Side surface, furthest from the x axis
+Line Loop(7) = {-3,-11,7,12};
+
+// Side surface, nearest to the z axis
+Line Loop(8) = {-4,-12,8,9};
+
+// Side surface, furthest from the z axis
+Line Loop(9) = {2,-11,-6,10};
+
+// Top surface
+Line Loop(10) = {5,6,7,8};
+
+Plane Surface(20) = {5};
+Plane Surface(21) = {6};
+Plane Surface(22) = {7};
+Plane Surface(23) = {8};
+Plane Surface(24) = {9};
+Plane Surface(25) = {10};
+
+// Bottom of the box
+Physical Line(333) = {5};
+// Sides of the box perpendicular to the z axis
+Physical Line(666) = {6,7};
+// Sides of the box perpendicular to the x axis
+Physical Line(777) = {8,9};
+// Top of the box
+Physical Line(444) = {10};
+
+// This is just to ensure all the interior
+// elements get written out. 
+Physical Surface(333) = {20};
+Physical Surface(666) = {21,22};
+Physical Surface(777) = {23,24};
+Physical Surface(444) = {25};
+
+Surface Loop(800) = {22, 20, 21, 24, 25, 23};
+Volume(900) = {800};
+Physical Volume(999) = {900};
+
 
=== added directory 'tests/mphase_dusty_gas_shock_tube'
=== added file 'tests/mphase_dusty_gas_shock_tube/Makefile'
--- tests/mphase_dusty_gas_shock_tube/Makefile	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/Makefile	2012-08-09 00:46:20 +0000
@@ -0,0 +1,18 @@
+preprocess:
+	@echo **********Creating 1D mesh
+	../../bin/interval --dx=0.00271002710027 -- -1.35501355014 1.35501355014 line
+
+run:
+	@echo **********Running simulation
+	../../bin/fluidity -v2 -l mphase_dusty_gas_shock_tube.flml
+	../../bin/fluidity -v2 -l single_phase_frozen_flow_test.flml
+
+input: clean preprocess
+
+clean:
+	rm -f *.stat *.steady_state*
+	rm -f *.d.* *.vtu
+	rm -f *.msh
+	rm -f *.ele *.edge *.node *.poly *.bound
+	rm -f matrixdump* *.log* *.err*
+
 
=== added file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml'
--- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,733 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">mphase_dusty_gas_shock_tube</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">multiphase</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">5</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0.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">0.000001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">0.00037772998222</real_value>
+      <comment>This corresponds to a non-dimensional time of tau = 4. Multiplying this by the reference time of l/u_ref gives 0.00037772998222 seconds.
+
+l = (4/3)*(rho_p/rho_g)*d
+u_ref = sqrt(p/rho_g)</comment>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">3</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+  </timestepping>
+  <material_phase name="Air">
+    <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="preonly"/>
+          <preconditioner name="lu"/>
+          <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>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <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>
+        <initial_condition name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  if(X[0] &lt;= 0):
+    return 12.3
+  else:
+    return 1.23</string_value>
+          </python>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <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="preonly"/>
+          <preconditioner name="lu"/>
+          <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="1" dim1="dim" rank="1">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Ends">
+          <surface_ids>
+            <integer_value shape="2" rank="1">1 2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">1.83e-5</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <diagnostic>
+        <mesh name="CoordinateMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="CompressibleContinuityResidual" rank="0">
+      <diagnostic>
+        <mesh name="PressureMesh"/>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <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>
+      </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>
+              <streamline_upwind>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind>
+            </stabilisation>
+            <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="preonly"/>
+          <preconditioner name="lu"/>
+          <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">
+          <constant>
+            <real_value rank="0">205894.308943</real_value>
+            <comment>InternalEnergy from the formula: (3/2)*R*Temperature = (3/2)*287.04*300 K = 203252
+InternalEnergy needed for a Pressure of 1.013e5 = 205894.308943</comment>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <multiphase_properties>
+      <effective_conductivity>
+        <real_value rank="0">0.026</real_value>
+      </effective_conductivity>
+      <specific_heat>
+        <real_value rank="0">758</real_value>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <material_phase name="Dust">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">2500</real_value>
+          </reference_density>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <aliased material_phase_name="Air" 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/>
+            </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="preonly"/>
+          <preconditioner name="lu"/>
+          <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="1" dim1="dim" rank="1">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="Ends">
+          <surface_ids>
+            <integer_value shape="2" rank="1">1 2</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="PhaseVolumeFraction" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <control_volumes>
+            <face_value name="FiniteElement">
+              <limit_face_value>
+                <limiter name="Sweby"/>
+              </limit_face_value>
+            </face_value>
+            <diffusion_scheme name="BassiRebay"/>
+          </control_volumes>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="preonly"/>
+          <preconditioner name="lu"/>
+          <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):
+  if(X[0] &lt;= 0):
+    return 4.92e-8
+  else:
+    return 4.92e-4</string_value>
+            <comment>alpha_p = 0.000492 in dusty gas region.
+
+We know that the ratio of gas density (rho_g) to particle mass concentration (sigma_p) is assumed to be 1.0, so rho_g = sigma_p = 1.23 kg/m^3.
+The normalised mass concentration is 1.0 in the dusty gas region, so we divide 1.23 by the particle density (2500 kg/m^3) get obtain the particle phase volume fraction.</comment>
+          </python>
+        </initial_condition>
+        <boundary_conditions name="All">
+          <surface_ids>
+            <integer_value shape="2" rank="1">1 2</integer_value>
+          </surface_ids>
+          <type name="zero_flux"/>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="ParticleReynoldsNumber" rank="0">
+      <diagnostic>
+        <algorithm name="particle_reynolds_number" material_phase_support="multiple">
+          <depends>
+            <string_value lines="1">Dust::Velocity,Air::Velocity,Air::PhaseVolumeFraction,Air::Density</string_value>
+          </depends>
+          <particle_diameter>
+            <real_value rank="0">10e-6</real_value>
+          </particle_diameter>
+          <continuous_phase_name>Air</continuous_phase_name>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="CFLNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="InternalEnergy" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="InternalEnergy">
+          <density name="Density"/>
+        </equation>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <streamline_upwind>
+                <nu_bar_optimal/>
+                <nu_scale name="0.5">
+                  <real_value shape="1" rank="0">0.5</real_value>
+                </nu_scale>
+              </streamline_upwind>
+            </stabilisation>
+            <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="preonly"/>
+          <preconditioner name="lu"/>
+          <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">
+          <constant>
+            <real_value rank="0">205894.308943</real_value>
+            <comment>InternalEnergy from the formula: (3/2)*R*Temperature = (3/2)*287.04*300 K = 203252
+InternalEnergy needed for a Pressure of 1.013e5 = 205894.308943</comment>
+          </constant>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <multiphase_properties>
+      <particle_diameter>
+        <real_value rank="0">10e-6</real_value>
+      </particle_diameter>
+      <specific_heat>
+        <real_value rank="0">758</real_value>
+        <comment>1406</comment>
+      </specific_heat>
+    </multiphase_properties>
+  </material_phase>
+  <multiphase_interaction>
+    <fluid_particle_drag>
+      <drag_correlation name="wen_yu"/>
+    </fluid_particle_drag>
+    <heat_transfer>
+      <heat_transfer_coefficient name="gunn"/>
+    </heat_transfer>
+  </multiphase_interaction>
+</fluidity_options>
 
=== added file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml'
--- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,156 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+
+<testproblem>
+
+  <!-- A 1D gas-solid shock tube problem. The setup is the same as the one given in the paper by Miura & Glass (1982).  -->
+  <!-- This test uses results at approximately t = 3.77e-4 seconds (normalised time of tau = 4 in the paper by Miura & Glass (1982)). -->
+   
+  <name>mphase_dusty_gas_shock_tube</name>
+  <owner userid="ctj10"/>
+  <tags>flml</tags>
+  
+  <problem_definition length="medium" nprocs="1">
+    <command_line>make run</command_line>
+  </problem_definition>
+  
+  <variables>
+    <variable name="time" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_dusty_gas_shock_tube.stat")
+time=s["ElapsedTime"]["value"]
+    </variable>
+    
+    <variable name="positions" language="python">
+import vtktools
+
+file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu')
+file.GetFieldNames()
+p = file.GetLocations()
+positions = p[:,1]
+    </variable>
+    
+    <variable name="air_velocity" language="python">
+import vtktools
+
+file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu')
+file.GetFieldNames()
+air_velocity = file.GetVectorField('Air::Velocity')
+    </variable>
+    
+    <variable name="dust_vfrac" language="python">
+import vtktools
+
+file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu')
+file.GetFieldNames()
+dust_vfrac = file.GetScalarField('Dust::PhaseVolumeFraction')
+    </variable>
+    
+    <variable name="pressure" language="python">
+import vtktools
+
+file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu')
+file.GetFieldNames()
+pressure = file.GetScalarField('Air::Pressure')
+    </variable>
+    
+    <variable name="max_dust_density" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_dusty_gas_shock_tube.stat")
+max_dust_density=s["Dust"]["Density"]["max"][-1]
+    </variable>
+    
+    <variable name="min_dust_density" language="python">
+from fluidity_tools import stat_parser
+s = stat_parser("mphase_dusty_gas_shock_tube.stat")
+min_dust_density=s["Dust"]["Density"]["min"][-1]
+    </variable>
+  </variables>
+    
+  <pass_tests>
+    <test name="Solvers have converged" language="python">
+import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files
+assert solvers_converged
+    </test>
+    
+    <test name="Finish time is 3.77e-4 seconds" language="python">
+assert (time[-1] &gt;= 3.77e-4)
+assert (time[-1] &lt; (3.77e-4 + 0.000002))
+    </test>    
+    
+    <test name="Normalised Gas::Velocity is 0 for normalised x in [-50, -10]" language="python">
+nnodes = len(air_velocity)
+
+vals = []
+for i in range(nnodes):
+   if(positions[i]/0.0271002710027 &lt; -10.0):
+      assert abs(air_velocity[i,0]/286.980353992) &lt; 1.0e-12
+    </test>
+    
+    <test name="Normalised Gas::Velocity is in [0, 0.9] for normalised x in [-10, 10]" language="python">
+nnodes = len(air_velocity)
+
+vals = []
+for i in range(nnodes):
+   if( (positions[i]/0.0271002710027 &gt; -10.0) and (positions[i]/0.0271002710027 &lt; 10.0) ):
+      assert abs(air_velocity[i,0]/286.980353992) &lt; 0.9
+      <!-- -1.0e-9 is needed to deal with small under-shoots -->
+      assert abs(air_velocity[i,0]/286.980353992) &gt; -1.0e-12
+    </test>
+    
+    <test name="Normalised Gas::Velocity is 0 for normalised x in [10, 50]" language="python">
+nnodes = len(air_velocity)
+
+vals = []
+for i in range(nnodes):
+   if(positions[i]/0.0271002710027 &gt; 10.0):
+      assert abs(air_velocity[i,0]/286.980353992) &lt; 1.0e-12
+    </test>
+    
+    <test name="Normalised Pressure is 10.0 for normalised x in [-50, -10]" language="python">
+nnodes = len(pressure)
+
+vals = []
+for i in range(nnodes):
+   if(positions[i]/0.0271002710027 &lt; -10.0):
+      assert abs(pressure[i]/1.013e5 - 10.0) &lt; 1.0e-6
+    </test>
+    
+    <test name="Normalised Pressure is in [1.0, 10.0] for normalised x in [-10, 10]" language="python">
+nnodes = len(pressure)
+<!-- Note: values of 0.9 and 10.1 have been used instead of 1.0 and 10.0 because of small over-shoots/under-shoots -->
+vals = []
+for i in range(nnodes):
+   if( (positions[i]/0.0271002710027 &gt; -10.0) and (positions[i]/0.0271002710027 &lt; 10.0) ):
+      assert abs(pressure[i]/1.013e5) &lt; 10.1
+      assert abs(pressure[i]/1.013e5) &gt; 0.9
+    </test>
+    
+    <test name="Normalised Pressure is 1.0 for normalised x in [10, 50]" language="python">
+nnodes = len(air_velocity)
+
+vals = []
+for i in range(nnodes):
+   if(positions[i]/0.0271002710027 &gt; 10.0):
+      assert abs(pressure[i]/1.013e5) &gt;= 1.0
+    </test>
+    
+    <test name="Mass concentration of Dust is between 0 and 2.5" language="python">
+nnodes = len(air_velocity)
+
+vals = []
+for i in range(nnodes):
+   assert abs(dust_vfrac[i]*2500.0/1.23) &lt; 2.5
+   assert abs(dust_vfrac[i]*2500.0/1.23) &gt; 0
+    </test>
+    
+    <test name="Dust density stays constant" language="python">
+assert abs(max_dust_density - min_dust_density) &lt; 1.0e-16
+    </test>
+  </pass_tests>
+  
+  <warn_tests>
+  </warn_tests>
+  
+</testproblem>
 
=== added file 'tests/mphase_dusty_gas_shock_tube/plot.py'
--- tests/mphase_dusty_gas_shock_tube/plot.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/plot.py	2012-08-09 00:46:20 +0000
@@ -0,0 +1,93 @@
+import pylab
+import numpy
+import shocktube
+import vtktools
+params = {'text.fontsize': 6,
+         'legend.fontsize': 8,
+         'xtick.labelsize': 8,
+         'ytick.labelsize': 8,
+         'lines.markersize': 6,
+         'axes.titlesize': 'small'}
+pylab.rcParams.update(params)
+
+filename='mphase_dusty_gas_shock_tube_378.vtu' # This corresponds to approximately t = 3.77e-4 seconds (normalised time of tau = 4 in the paper by Miura & Glass (1982))
+vt=vtktools.vtu(filename)
+
+# Time and coordinates
+t=vt.GetScalarField('Air::Time')[0]
+xyz=vt.GetLocations()
+x=xyz[:,0]
+
+# Solution fields
+p=vt.GetScalarField('Air::Pressure')
+uvw_air=vt.GetVectorField('Air::Velocity')
+u_air=uvw_air[:,0]
+uvw_dust=vt.GetVectorField('Dust::Velocity')
+u_dust=uvw_dust[:,0]
+rho_air=vt.GetScalarField('Air::Density')
+rho_dust=vt.GetScalarField('Dust::Density')
+ie_air=vt.GetScalarField('Air::InternalEnergy')
+ie_dust=vt.GetScalarField('Dust::InternalEnergy')
+vfrac=vt.GetScalarField('Dust::PhaseVolumeFraction')
+
+# Plot the numerical results.
+# Note that we have divided the results by the reference velocity/pressure/density/internal energy to normalise them.
+# x is normalised by the reference length.
+pylab.figure(1)
+pylab.subplot(321)
+pylab.plot( x/0.0271002710027, p/1.013e5,'b.', label="Numerical")
+
+pylab.subplot(322)
+print len(x)
+print len(u_air)
+pylab.plot( x[:len(x)]/0.0271002710027, u_air/286.980353992,'b.', label="Numerical (Air)")
+pylab.plot( x[:len(x)]/0.0271002710027, u_dust/286.980353992,'r.', label="Numerical (Particles)")
+
+pylab.subplot(323)
+pylab.plot( x/0.0271002710027, rho_air/1.23,'b.', label="Numerical")
+
+pylab.subplot(324)
+pylab.plot( x/0.0271002710027, rho_dust,'b.', label="Numerical")
+pylab.title('Density of Particles')
+pylab.legend(loc=2)
+
+pylab.subplot(325)
+pylab.plot( x/0.0271002710027, ie_air/82357.7235772,'b.', label="Numerical (Air)")
+pylab.plot( x/0.0271002710027, ie_dust/82357.7235772,'r.', label="Numerical (Dust)")
+
+# Multiply the PhaseVolumeFraction by the Dust Density to get the mass concentration, and then normalise.
+pylab.subplot(326)
+pylab.plot( x/0.0271002710027, vfrac*2500.0/1.23,'b.', label="Numerical")
+pylab.title('Mass Concentration from PhaseVolumeFraction of Particles')
+pylab.legend(loc=2)
+  
+# Now work out the frozen flow ('analytical') solutions
+sol=numpy.array([shocktube.solution(xi,4.0) for xi in x/0.0271002710027])
+p=sol[:,0]
+u=sol[:,1]
+rho=sol[:,2]
+ie=p/rho/(shocktube.gamma-1.0)
+  
+pylab.subplot(321)
+pylab.plot( x/0.0271002710027, p,'-g', label='Frozen flow')
+pylab.title('Normalised Pressure')
+pylab.legend(loc=2)
+  
+pylab.subplot(322)
+pylab.plot( x/0.0271002710027, u,'g-', label='Frozen flow (Air)')
+pylab.title('Normalised Velocity')
+pylab.legend(loc=2)
+
+pylab.subplot(323)
+pylab.plot( x/0.0271002710027, rho,'g-', label='Frozen flow')
+pylab.title('Normalised Density of Air')
+pylab.legend(loc=2)
+
+pylab.subplot(325)
+pylab.plot( x/0.0271002710027, ie,'g-', label='Frozen flow')
+pylab.title('Normalised Internal Energy of Air')
+pylab.legend(loc=2)
+
+pylab.savefig('mphase_dusty_gas_shock_tube.png')
+
+pylab.show()
 
=== added file 'tests/mphase_dusty_gas_shock_tube/shocktube.py'
--- tests/mphase_dusty_gas_shock_tube/shocktube.py	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/shocktube.py	2012-08-09 00:46:20 +0000
@@ -0,0 +1,157 @@
+# bit of python to compute analytical solution for the shock tube problem
+# solution taken from "Mathematical and Computational Methods for Compressible Flow"
+# by Miloslav Feistauer, Jiri Felcman and Ivan Straskraba, 2003
+# eqn numbers below refer to this books
+from math import sqrt
+from scipy.optimize import newton
+
+# ratio of specific heats
+gamma=1.4
+
+# left state:
+ul=0.
+rhol=10.0
+iel=2.5
+
+# right state:
+ur=0.
+rhor=1.0
+ier=2.5
+
+# this eos is assumed in the eqns below (so can't change this):
+def p_eos(ie, rho):
+  return rho*ie*(gamma-1.0)
+  
+def rho_eos(ie, p):
+  return p/ie/(gamma-1.0)
+
+
+# derived quantities, left:
+pl=p_eos(iel,rhol)
+El=rhol*iel
+al=sqrt(gamma*pl/rhol)
+# derived quantities, right:
+pr=p_eos(ier,rhor)
+Er=rhor*ier
+ar=sqrt(gamma*pr/rhor)
+
+
+def F1l(p):
+  """Function F1l defines difference between us and ul:
+     us=ul+F1l(p)    eqn. (3.1.159)"""
+     
+  # eqn. (3.1.160)
+  if p<=pl:
+    return 2.*al/(gamma-1.)*(1.-(p/pl)**((gamma-1.)/(2.*gamma)))
+  else:
+    return -(p-pl)*sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl))
+    
+def F1lprime(p):
+  # derivative of the above
+  if p<=pl:
+    return -al/gamma*(p/pl)**((gamma-1.)/(2.*gamma)-1.)
+  else:
+    return   -1.*sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl)) + \
+      -(p-pl)/2./sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl))
+
+def F3r(p):
+  """Function F3r defines difference between us and ur:
+     us=ur+F3r(p)    eqn. (3.1.161)"""
+     
+  # eqn. (3.1.162)
+  if p<=pr:
+    return -2.*ar/(gamma-1.)*(1.-(p/pr)**((gamma-1.)/(2.*gamma)))
+  else:
+    return (p-pr)*sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr))
+    
+def F3rprime(p):
+  # derivative of the above
+  if p<=pr:
+    return ar/gamma*(p/pr)**((gamma-1.)/(2.*gamma)-1.)
+  else:
+    return      sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr)) + \
+      (p-pr)/2./sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr))
+
+def F(p):
+  # eqn (3.1.165)
+  return F3r(p)-F1l(p)+ur-ul
+  
+def Fprime(p):
+  return F3rprime(p)-F1lprime(p)
+  
+# inital guess:
+p=(pl+pr)/2.
+
+# solve eqn (3.1.164): F(p*)=0, to find p* the pressure between the u-a and u+a waves
+# (is constant over contact discontinuity at the u wave)
+ps=newton(F, p, fprime=Fprime, maxiter=500)
+print "should be zero: F(p*) = ", F(ps)
+print "p* =", ps
+us=ul+F1l(ps)
+print "u* =", us
+print "should be equal to:", ur+F3r(ps)
+
+# the star region is divide by the contact discontinuity which gives a discontinuity in
+# density rhosl/rhosr and internal energy iesl/iesr
+
+# left:
+if ps<pl: # rarefaction
+  rhosl=rhol*(ps/pl)**(1/gamma) # eqn (3.1.172)
+else: # shock
+  rhosl=rhol*( (gamma-1.)/(gamma+1.)*pl/ps +1. )/(pl/ps+(gamma-1.)/(gamma+1.)) # eqn (3.1.176)
+  s1=ul-al*sqrt((gamma+1.)/2./gamma*ps/pl + (gamma-1.)/2./gamma) # shock speed, eqn (3.1.177)
+iesl=ps/rhosl/(gamma-1.)
+asl=sqrt(gamma*ps/rhosl)
+
+# right:
+if ps<pr: # rarefaction
+  rhosr=rhor*(ps/pr)**(1/gamma) # eqn (3.1.178)
+else: # shock
+  rhosr=rhor*( ps/pr + (gamma-1.)/(gamma+1.) )/( (gamma-1.)/(gamma+1.)*ps/pr +1. ) # eqn (3.1.182)
+  s3=ur+ar*sqrt((gamma+1.)/2./gamma*ps/pr + (gamma-1.)/2./gamma) # shock speed, eqn (3.1.183)
+iesr=ps/rhosr/(gamma-1.)
+asr=sqrt(gamma*pr/rhosr)
+
+def solution(x,t):
+  if x/t<us:    
+    # before the contact discontinuity:
+    
+    if ps<pl:
+      # u-a is a rarefaction wave
+      if x/t<ul-al: # left
+        return (pl, ul, rhol)
+      elif x/t<us-asl: # within the rarefaction wave
+        p=pl*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/al*(ul-x/t))**(2*gamma/(gamma-1.)) # eqn (3.1.97)
+        u=2./(gamma+1.)*(al + (gamma-1.)/2.*ul + x/t) # eqn (3.1.97)
+        rho=rhol*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/al*(ul-x/t))**(2./(gamma-1.)) # eqn (3.1.97)
+        return (p, u, rho)
+      else: # after the rarefaction before the contact disc.
+        return (ps,us,rhosl)
+    else:
+      # u-a is a shock wave
+      if x/t<s1: # left
+        return (pl, ul, rhol)
+      else: # between u-a shock and contact disc.
+        return (psl, us, rhosl)
+        
+  else:
+    # after the contact discontinuity:
+        
+    if ps<pr:
+      # u+a is a rarefaction wave
+      if x/t>ur+ar: # right
+        return (pr, ur, rhor)
+      elif x/t>usr+asr: # within the rarefaction wave
+        p=pr*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/ar*(ur-x/t))**(2*gamma/(gamma-1.)) # eqn (3.1.98)
+        u=2./(gamma+1.)*(-ar + (gamma-1.)/2.*ur + x/t) # eqn (3.1.98)
+        rho=rhor*(2./(gamma+1.) - (gamma-1.)/(gamma+1.)/ar*(ur-x/t))**(2./(gamma-1.)) # eqn (3.1.98)
+        return (p, u, rho)
+      else: # after the contact disc. before the rarefaction
+        return (ps,us, rhosr)
+    else:
+      # u+a is a shock wave
+      if x/t>s3: # right
+        return (pr, ur, rhor)
+      else: # between contact disc. and u-a shock
+        return (ps, us, rhosr)
+
 
=== added 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	1970-01-01 00:00:00 +0000
+++ tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml	2012-08-09 00:46:20 +0000
@@ -0,0 +1,379 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">single_phase_frozen_flow_test</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">1</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="line">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="DensityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">5</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0.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">0.000001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">0.00037772998222</real_value>
+      <comment>This corresponds to a non-dimensional time of tau = 4. Multiplying this by the reference time of l/u_ref gives 0.00037772998222 seconds.
+
+l = (4/3)*(rho_p/rho_g)*d
+u_ref = sqrt(p/rho_g)</comment>
+    </finish_time>
+    <nonlinear_iterations>
+      <integer_value rank="0">3</integer_value>
+      <tolerance>
+        <real_value rank="0">1.0e-9</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </nonlinear_iterations>
+  </timestepping>
+  <material_phase name="Air">
+    <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="preonly"/>
+          <preconditioner name="lu"/>
+          <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>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="Density" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>