Merge lp:~asc/fluidity/extrudeup into lp:~asc/fluidity/fluidity_for_extrude

[Adam Candy]

[Adam Candy (asc)]

Requires documenting in the manual.

[lp:~asc/fluidity/extrudeup]

3770. By j-bull08 <email address hidden> on 2011-10-31

small change to 3d BFS example to improve the mesh.
Removing unnecessary mesh files.
Sorry for lateness - usual internet problems up here. I have tested these changes.

3771. By j-bull08 <email address hidden> on 2011-10-31

missing image. Thanks Hannah.

3772. By Hannah Hiester on 2011-10-31

Manual and example Makefile updates
- Makefiles just edited to clean up pdf generated in post processing
- Link between metric advection in adaptivity chapter and examples chapter put it
- Tidying up vtktools table (if someone could have a quick scan that would be helpful)

3773. By Tim Bond on 2011-10-31

Commenting out request for missing PDF file so that mediumtests will run;
please could this be supplied and the line uncommented.

3774. By j-bull08 <email address hidden> on 2011-11-01

Last missing image file. Thanks Tim for bringng this to my attention.

3775. By j-bull08 <email address hidden> on 2011-11-01

undoing Tim's commenting-out of image file; this has been committed.

3776. By Matthew Piggott on 2011-11-01

Fix a figure reference.

[Adam Candy (asc)]

The calculation of top_max in the extrude routines needs to be the max across all partitions - how can I ensure this is negotiated across all nodes?

review: Needs Information

[Stephan Kramer (s-kramer)]

What about:

  use parallel_tools
  call allmax(top_max)

[lp:~asc/fluidity/extrudeup]

3777. By Jon Hill on 2011-11-02

Getting rid of the undefined references in the param section of the manual

[Adam Candy (asc)]

Perfect - thanks.

[Adam Candy (asc)]

The extrude branch has been updated.

- Added back in the removed preprocessing step
- Ensured the align_with_geoid option works when the domain is spread over multiple partitions
- Added more accurate control over the position of the horizontal layers
- Removed some debugging ewrites

This is now complete, except for a description in the manual.
I'll send out an update of developments to everyone.

review: Needs Resubmitting

[Adam Candy (asc)]

That previous comment was obviously supposed to be an 'Approve' from me.

review: Approve

[Adam Candy (asc)]

OK, parallel_tools added and the dependencies regenerated - just a manual entry left now, and that can wait.

review: Approve

[lp:~asc/fluidity/extrudeup]

3778. By Jon Hill on 2011-11-02

Fixing missing citations. This should make medium tests green now

3779. By Alexandros Avdis on 2011-11-02

Adding book title to the sagaut1998 bibliography entry.

3780. By Patrick Farrell on 2011-11-02

Merge the Debian packaging branch. If you want to try it out, go to your source
tree and type

debuild -uc -us

and then the packages will be in the parent directory.

Approved by our man in Amsterdam, S. C. Kramer:
https://code.launchpad.net/~pefarrell/fluidity/debian-packaging/+merge/81048

3781. By Stephan Kramer on 2011-11-02

Some updates on the meshes chapter.

3782. By Alexandros Avdis on 2011-11-02

Added link to the Gmsh tutorial in the ``Meshes for Fluidity'' chapter of the manual. Corrected a few typos.

3783. By Adam Candy on 2011-11-02

Update to the serialtest make directive so that 'make serialtest' on a fresh checkout passes.

3784. By Adam Candy on 2011-11-02

The backward_facing_step_3d example postprocess make directive failed for me.
It assumes the folder numpy_data exists - this change to the makefile should fix it.

3785. By Adam Candy on 2011-11-03

Reinstated version numbering for trunk checkouts.

3786. By Tim Bond on 2011-11-03

Small update adding a clean-debian target to Makefile.in to remove files
generated by the debian 'debuild' process.

3787. By j-bull08 <email address hidden> on 2011-11-03

Missing reference data for BFS 2d

3788. By Tim Bond on 2011-11-03

Thanks to Lawrence for spotting this bug.

Where there is a line such as:

 [ -f foo ] && bar

in the makefile, and foo doesn't exist, [ -f foo ] returns an error
which gets passed through to the end of the command and which make
then interprets as a build error, and exists. Adding in a || true at
the end of the line fixes this.

3789. By Stephan Kramer on 2011-11-03

The final version of the fluidity version saga?

The version information is obtained from bzr as before except when
FLUIDITY_VERSION is passed to configure (used in debian/rules of the debian
package) in which case it is simply set to that value.

3790. By Patrick Farrell on 2011-11-03

python setup.py install now installs fluidity.diagnostics too. Fixes
lp:885642.

3791. By Patrick Farrell on 2011-11-03

Make install should depend on fltools.

3792. By Simon Mouradian on 2011-11-03

Removes warning when 'make clean' with python 2.7.2

3793. By Patrick Farrell on 2011-11-03

Don't pass --install-layout=deb unless we're building the Debian package.

3794. By Patrick Farrell on 2011-11-03

When building the Debian package, only use one processor.

3795. By Tim Bond on 2011-11-03

We want, ultimately, to release a package as '4.1'. Since we suspect that
newness of packages is based on a standard sort, anything appended to 4.1
as a version will probably end up superceding rather than being superceded
by 4.1. Thus, 4.1beta1 wouldn't upgrade as intended when release happens.

Thus, renumbering the changelog to base it on 4.0.999 for the moment, which
can be neatly bumped to 4.1 at release. Note that this also has the 'maverick'
change in, which is required by launchpad so it knows which builder to send the
package to.

3796. By Stephan Kramer on 2011-11-03

Adding fluidity_tools.py to python-fluidity package.

3797. By Patrick Farrell on 2011-11-03

Several changes to the Debian packaging:
 * configure --enable-2d-adaptivity
 * install the fluidity schemas into /usr/share/fluidity/schemas
 * register that schema with diamond in /usr/share/diamond/schemata/flml

This means that Diamond should Just Work for users of the binary package.

3798. By Stephan Kramer on 2011-11-03

Fix of fluidity_tools commit.

3799. By Hannah Hiester on 2011-11-03

Small edits to the adaptivity chapter and including publishers for some of the citations that were making bibtex complain

3800. By Jonathan Bull on 2011-11-03

Switching to strong Dirichlet BC, as the weak log law sometimes blows up.

3801. By Alexandros Avdis on 2011-11-03

Making notations in Parameterisations (LES) chapter consistent.

3802. By Tim Bond on 2011-11-03

The debian package depended on fluidity-dev for both build *and* runtime,
which meant that an install of runtime fluidity required a huge overhead
of dev packages. Since libraries should be handled by the shlib directive
and we *probably* only want to load enough that the standard examples can
be run, this is a first attempt at a (relatively) minimal set of pacakges
for the runtime fluidity to depend on.

I'm taking 'dependencies' a little bit loosely here, namely assuming:

* gmsh and triangle[-bin] are depended on for mesh handling
* diamond is depended on for configuration
* python is depended on for runtime examples, including python-lxml,
    python-scipy, python-numpy, python-vtk, python-matplotlib,
    python-sympy, and python-fluidity. This list was compiled from
    working through all the python calls in the examples directory.

Going with this as an initial attempt at a good set, I'll push the new
package through launchpad and take a 'suck it and see' approach with the
examples on a clean maverick workstation build.

3803. By Jonathan Bull on 2011-11-03

capping the run at 100s to prevent blow-up.

3804. By Patrick Farrell on 2011-11-04

Make the package install shallow_water and burgers_equation too

3805. By Patrick Farrell on 2011-11-04

Apply Lawrence's fix for bug lp:886062. It is utterly daft
that Python changed the behaviour of one of the core methods of
the elementtree module sometime between maverick and natty.

3806. By Tim Bond on 2011-11-04

For the purposes of buildbot, and testing of the make install, we want
testharness to be part of the make install process. Noted, though, that
this is a fairly location-specific program which isn't of particular use
in the debian package, so a neater way of including it for buildbot but
excluding it for debian might be beneficial.

3807. By Tim Bond on 2011-11-04

We now generate a schemas/flml which needs to be cleaned up by make distclean

3808. By Adam Candy on 2011-11-04

Fixes to the backward-facing step 2d reference simulation, as it was still crashing.
With this, preprocess, run and postprocess all run without error.

Still waiting on the other two backward-facing step runs.

3809. By Simon Funke on 2011-11-04

Merge the python interface from the libspud trunk which fixes "set_option".
In addition this commit updates optimality.py to be compatible with these libspud changes.

3810. By David Ham on 2011-11-04

Cause FLLINKER to be used for the hybridized helmholz solver so as to hopefully placate ifort.

3811. By Hannah Hiester on 2011-11-04

Manual and example postprocessing script updates

For the tides_in_the_Mediterranean_Sea example renaming the vtu produced is postprocessing to reflect the sea which it is modelling.

Manual updates:
- Adaptivity chapter are mostly minor or references to other sections. A section on the cost of adaptivity has also been included - would be useful if someone could have a quick scan.

- Configuring fluidity chapter are again mostly minor and also included remaining general adaptivity options.

- Examples chapter are to the timings for flow past a sphere (now takes a little longer) and tides in the Mediterranean Sea (now is much faster). Also updated the RMS differences in tides in the Mediterranean Sea example to the current values (the change is very small but figured it was best to have as much consistency as posible between the numbers that come out of the postprocessing and the manual)

3812. By Alexandros Avdis on 2011-11-04

Improvements to the LES documentation in the fluidity manual, notation is now
consistent, the manual lists and explains the options of the SGS models.

3813. By Christian Jacobs on 2011-11-05

Minor improvements to the diagnostic algorithms section in visualisation_and_diagnostics.tex.

3814. By Alexandros Avdis on 2011-11-05

minor corrections to LES documentation in the Fluidity manual.

3815. By Tim Bond on 2011-11-06

r3810 went half way to fixing the Intel compiler issues with hybridized
helmholz. This is built twice, so needs an additional fix as per this commit.

3816. By Tim Bond on 2011-11-06

Mirroring the previous commit to hybridized helmholz with shallow_water also
needing fixing for the Intel build.

3817. By Tim Bond on 2011-11-06

Removing a large amount of rogue whitespace in an attempt to placate Intel.

3818. By Tim Bond on 2011-11-06

Adding .msh files to the k-epsilon_mms_cg_p1p1 test to remedy failures
resulting from sensitivity to gmsh versions. Tested on oneiric where there were
failures, and appears to do the trick.

3819. By Hannah Hiester on 2011-11-06

Updating to rm -f Sphere_drag.pdf to ensure it doesn't complain if the file doesn't exist

3820. By Tim Bond on 2011-11-06

Bringing in the examples as a tarball to 'make install' so they can be shipped
with the binary package. This also requires a little bit of munging before they
are tarred up, changing relative paths to binaries/scripts into just the names
in the assumption that they will be on the path once installed.

Also adding in vtktools to the install process as it's required by the
examples.

3821. By Tim Bond on 2011-11-07

Adding a one-page-primer for Fluidity. This describes how to add the Fluidity
repository, apt-get Fluidity, uncompress the examples, and run top_hat.

It's presently inserted as a single page following the Overview.

Someone casting a quick eye over this page would be appreciated.

3822. By Tim Bond on 2011-11-07

Adding to commit of r3821, it's helpful if I remember to add the relevant new
LaTeX file having for the one page primer.

3823. By Tim Bond on 2011-11-07

Manual updates to:

* Add -y to apt commands so the user does not have to respond to prompts
* Add in binary package instructions to the 'getting started' chapter
* Remove some of the verbosity from the 'getting started' chapter

3824. By Tim Bond on 2011-11-07

Version updates, reflecting the upgrade of fluidity itself from 4.0 to 4.1 and
an associated overhaul of the external libraries chapter to reflect supporting
software which has been upgraded since our last release. Build information for
supporting software has been brought up to date with the wiki page, and
software versions have been updated to reflect Ubuntu Natty versions and
versions as tested through buildbot or at partner institutions.

3825. By Adam Candy on 2011-11-07

Update to the run times for the backward-facing step exmaples.
(Changes limited to the manual only).

3826. By Adam Candy on 2011-11-07

Fix the spacing in the \fluidity command.

3827. By Axelle Vire on 2011-11-07

Correct typo in getting started: source code should be two words.

3828. By Axelle Vire on 2011-11-07

Fix typo and line cut in Appendix C.

3829. By Hannah Hiester on 2011-11-07

Updates to example xmls to ensure they pass with testharness (donein agreement with example maintainers)

- flow_past_sphere_Re1000: relaxing tolerances by 1%
- lock_exchange: error in the python corrected
- top_hat: relaxing tolerance on DG bounds

3830. By Axelle Vire on 2011-11-07

Fix \fluidity spacing in Appendix D. An mbox was using, hence creating double spacing.

3831. By Tim Bond on 2011-11-07

Updating package changelog in preparation for 4.1-release

3832. By Stephan Kramer on 2011-11-07

Adding explanation for 2+1D configuration options.

3833. By Tim Bond on 2011-11-07

Small update to the tarball information to reflect merging examples into the
main tarball (to make life easier for people).

3834. By Axelle Vire on 2011-11-07

Make chapter 2 (UK) correct.

3835. By Axelle Vire on 2011-11-07

Fix a couple of typos in chapters 2 and 9.

3836. By Florian Rathgeber on 2011-11-07

Fixes to section 1.2 of the manual

3837. By Florian Rathgeber on 2011-11-07

Suggest manual readers to run make -jN rather than make -j

3838. By Stephan Kramer on 2011-11-07

In the exercises for the 1d advection example, suggest a different alternative cg configuration that beats the crap out of dg. Thanks to CianW.

3839. By Axelle Vire on 2011-11-07

Various corrections to chapter 4.

3840. By Axelle Vire on 2011-11-07

Fix typo in chapter 5.

3841. By Axelle Vire on 2011-11-07

Various corrections to chapters 8 and 10.

3842. By Axelle Vire on 2011-11-07

Tiny fixes in chapter 10.

3843. By Tim Bond on 2011-11-07

Clarifying our support policy of aiming to maintain packages and support for
Ubuntu releases for 18 months. This is mandated for non-LTS by the timescale of
Ubuntu support, but given rapid ageing of software is reckoned in Fluidity time
to be reasonable for LTS as well.

3844. By Axelle Vire on 2011-11-07

Fix some capital letters in the bib.

3845. By Axelle Vire on 2011-11-07

One more for the bib.

3846. By Tim Bond on 2011-11-08

Commit purely addressing documentation, orthogonal to functional codebase.

Adding in a short paragraph to the 'Getting Started' chapter noting that unit
testing is broken under Intel as a result of compiler issues, and that we have
no intention of fixing this at our end. Since Intel seem to have no intention
of fixing it at their end either, users are hereby advised not to make unittest
on an Intel build of Fluidity.

3847. By Patrick Farrell on 2011-11-09

This test has been broken since r3671, but failing silently. The test
Makefile had

input: meshes clean

where meshes made the meshes, and then clean deleted them straight after!

3848. By Patrick Farrell on 2011-11-09

And the same again for python_shape_dshape_test_2d.

3849. By Patrick Farrell on 2011-11-09

galerkin -> Galerkin in the text.

3850. By Axelle Vire on 2011-11-10

When constructing the supermesh, the region_ids allocated the intersections need to be explicitely deallocated. This fix a memory leak that was becoming significant when modelling fluid-solid interactions, where the supermesh is constructed at every timestep.

3851. By Axelle Vire on 2011-11-10

This will revert my previous commit, which I meant to do on a branch. Even if this is a small fix, let's do it the Right Way.

3852. By Axelle Vire on 2011-11-11

Merge branch. This is a small change in femtools/Supermesh.F90 that deallocates the region_ids associated to each intersection when constructing the supermesh. This fixes a memory leak that was becoming significant when modelling fluid-solid interactions, where the supermesh is constructed at every timestep. This passes all tests locally.

[Adam Candy (asc)]

Test failures fixed.
Now builds and passes on buildbot:
http://buildbot-ocean.ese.ic.ac.uk:8080/builders/x86_64-extrudeup-gcc4

[lp:~asc/fluidity/extrudeup]

3853. By Paul Woodhams on 2011-11-14

Tidying up the edge-list callback. Removed some redundant code, fixed some MPI data types and added a global calculation of the minimum edge weight.

3854. By Stephan Kramer on 2011-11-17

Bug fix in petsc_readnsolve, for the case where no mesh needs to be read. A dummy mesh is created, but this was done with a bit of a broken hack, now fixed.

3855. By Cian Wilson on 2011-11-17

Improved error handling for remaps.

3856. By Tim Bond on 2011-11-21

Fix for bug lp:892228 making the schemata install honour prefix and datadir.

3857. By Cian Wilson on 2011-11-21

Remove spurious allocate that introduced a memory leak in r3855.

3858. By Patrick Farrell on 2011-11-22

Merge the dont-install-diamond branch.

3859. By Tim Bond on 2011-11-22

A one-line build-system-only fix adding a '-f' to an rm in make install, a
direct bug leading on from r3858 and highlighted by the Intel builder.

3860. By Brendan Tollit on 2011-11-22

Summary:

This merge will correct the sign and number of additions of the weak velocity Dirichlet BC integrals when using a CV tested continuity equation for incompressible flow.

Code Changes:

- assemble/Divergence_Matrix_CV.F90. The variable ct_rhs_local is formed and added (with a minus) for each dimension in the procedure assemble_divergence_matrix_cv.

- schemas/prognostic_field_options.rnc and .rng. A mesh choice is permitted for the Diffusivity of a scalar field. Required for the some of the tests of the fix above and is a useful general addition.

Test Cases Removed:

- tests/darcy_p1dgp2_pressBCinlet_1d

- tests/darcy_p1dgp2_velBCinlet_1d

Test Cases Added:

- tests/darcy_p1dgp2_pressBCinlet. The one removed but now in 1,2 and 3d.

- tests/darcy_p1dgp2_velBCinlet. The one removed but now in 1,2 and 3d.

- tests/darcy_p0p1cv_velBCinlet. Checks the sign of the weak velocity Dirichlet BC aligned with cartesian axis. Single phase.

- tests/darcy_p0p1cv_velBCinlet_1solidphase. Checks the sign of the weak velocity Dirichlet BC aligned with cartesian axis. Two phase with the second representing the solid porosity.

- tests/darcy_p0p1cv_velBCinlet_rotated_domain. Checks the sign of the weak velocity Dirichlet BC NOT aligned with cartesian axis. Single phase.

- tests/mms_ns_p0p1cv_weak_velBC_steady. A mimic of the standard mms_ns test cases using p0p1cv with a weak velocity Dirichlet BC. Here the domain was rotated from a square to a diamond such as to have non Cartesian aligned velocity Dirichlet BC on all surfaces. Checks L2 velocity and pressure convergence is between 0.8 and 1.0. This is a short test.

- tests/mms_ns_p1dgp1cv_weak_velBC_steady. Similar to above but using p1dgp1cv. Checks L2 velocity and pressure convergence is between 1.5 and 2.15. This is a medium test.

3861. By Stephan Kramer on 2011-11-23

Using get_coordinates_remapped_to_surface() to fix https://bugs.launchpad.net/bugs/893893

3862. By Brendan Tollit on 2011-11-23

This merge will add the capability when discretising the pressure
with continuous galerkin to test the continuity equation with the
control volume dual function space (rather than the CG function space).

This is tested on a range of cases checking for capturing balances and
getting the correct order of convergence.

This is documented in the manual (theory and usage).

It is described in the manual and schema that this will only work
for incompressible flow and not work with the free surface,
wetting and drying and implicit solids two way coupling models.
Also the pressure mesh must use the Lagranian element type. This is
also restricted to only work with the available control volume shape
functions.

This will produce a non symmetric pressure correction matrix which
must be considered when selecting the pressure solver options.

3863. By j-bull08 <email address hidden> on 2011-11-24

removing irritating ewrite in smoothing module (used for dynamic LES and Helmholtz smoothing algorithms)

3864. By Tim Bond on 2011-11-25

Commit entirely in documentation targetting bug lp:891733 and bringing in the
revision format from the root Makefile.in to the manual Makefile. Also updating
the front page to refer to 'branch' rather than 'trunk version'.

3865. By Tim Bond on 2011-11-25

Commit purely in documentation, providing fixes for the errors noted in
bug lp:885261. Some duplicate bibliography entries removed, some missing
fields inserted, a pair of missing braces added, and one entry reclassified
from article to book.

3866. By Tim Bond on 2011-11-25

Commit entirely in documentation targetting bug lp:891733 and in addition to
commit for r3864; I had failed to note that LaTeX could complain about the
presence of non-alnum characters. Thus, adding the 'tr' back in to strip all
punctuation. Thanks to Brendan for spotting this one.

3867. By Christian Jacobs on 2011-11-25

Minor updates to the mphase_rising_body test case. Changed the length from 'special' to 'short'. Passes testing on the ctjacobs-multiphase buildbot queue.

3868. By Cian Wilson on 2011-11-26

Remove erroneous dimension from Velocity. Wasn't affecting the test in any way.

3869. By Tim Bond on 2011-11-28

Addressing bug lp:894770 by removing unused tests (either marked 'special'
or which don't have a .xml).

3870. By Jon Hill on 2011-11-29

Fixing the FLML fiels for the GLS KP test such that the wind surface stress returned from python is the correct dimension (1 not 2)

3871. By Patrick Farrell on 2011-11-30

When setting a vector field from Python, check the length of the iterable
returned by the user, and check that it is the length we expect. Fixes
lp:800694.

3872. By Tim Bond on 2011-12-01

Change to documentation only - preparatory manual work for the changes to
Diamond build and install methods.

Note that this does not address (and we have not in the past addressed) the
needs of diamond for supporting software to exist on the system. This probably
needs to be added to one of the appendices at some point; specifically libxml2
and python-lxml, assuming I'm remembering correctly that they are hard
dependencies.

3873. By Brendan Tollit on 2011-12-02

Add an option to include the Source field directly to the right hand side
for prognostic scalar fields. Adjust CG, DG, FV and CV to
act on this as necessary. Include a very simple two field eigenvalue
(using the time step loop as a power iteration loop) 1d test
case for each discretisation checking order of converegence to test
this new functionality.

This also includes a comment in the schema that variations in
the Source field across subcycled time steps is not considered,
a constant value across the whole time step is assumed (whether
added directly or tested with shape_rhs).

In assemble/Advection_Diffusion_CG.F90 a redundant optional
variable passed in called extra_discretised_source is removed, as
well as the code that added it to the rhs.

3874. By Tim Bond on 2011-12-06

This merge supplies:

* An updated version of libspud
* An install target for diamond in the root fluidity makefile
* A target in the root fluidity makefile to install user schemata

3875. By Cian Wilson on 2011-12-15

Loosen tolerances on this test to hopefully make buildbot green again. Still curious why Gerard's changes made a test that's been stable for years change.

3876. By Cian Wilson on 2011-12-17

Cleaning up and making consistent the exclusion of tags in make test and make mediumtest (also deleting make mediumzoltantest as buildbot didn't use it anymore). Most importantly adding an algencan tag and 2dadapt tags where appropriate to medium tests so that they don't fail when fluidity is not configured with algencan (the norm) or 2d adaptivity (unusual).

3877. By Cian Wilson on 2011-12-17

This missed out on a 2dadapt flag in the merge.

3878. By Tim Bond on 2011-12-19

Minor relaxation on tolerance to allow for Precise.

3879. By Tim Bond on 2011-12-19

A slight relaxation of the tolerance on this test to cater for Precise.

3880. By Jon Hill on 2011-12-19

Editing geopressure_conv test Makefile to use fldecomp in the current checkout, not the centrally installed one

3881. By Brendan Tollit on 2011-12-23

This merge will add the capability for the enforce discrete property solenoidal projection when
the lagrange multiplier is CG to test the divergence of the interpolated field with the CV dual.
This is akin to what was added for the velocity-pressure solve previously via this branch.
This adds an option to the schema, commenting that this will make the matrix non symmetric.

There is also a small correction to an addition of the known velocity divergence to the projection
rhs to use ctp_m rather than ct_m (in this solenoidal projection routine).

This is tested via the solenoidal unittest, with an extra case added to check for CV divergence.
Non used code was also removed from here.

3882. By Brendan Tollit on 2011-12-23

This merge will add a procedure to be able to have non zero pressure Dirichlet
boundary conditions when using a CV discretisation for pressure. As for the CG
pressure BCs the integrals are included for both weak and strong Dirichlet.
For multiphase the phase volume fraction is included via a FE spatial
interpolation, akin to their inclusion in interior domain terms.

Code Additions:

- assemble/Pressure_Dirichlet_BCS_CV.F90. Adds the BC integrals via forming
the same local matrix as the weak velocity BCs but uses it transposed,
which is similar to pressure CG version(s).

Code Changes:

- assemble/Momentum_Equation.F90. If cv_pressure calls the above routine
after assembling the main momentum equation terms.

- assemble/Makefile.in. Update for new file/module (using .o rather than .F90
as spotted by Cian)

- assemble/Makefile.dependencies. Update for new file/module.

Test Cases Added:

- tests/darcy_p0p1cv_pressBCinlet. Simple channel flow model using darcy
flow with a pressure gradient using p0p1cv for 1,2 and 3d. Analytic
and very fast.

- tests/darcy_p0p1cv_pressBCinlet_1solidphase. Same as above but with an
extra phase to represent the porous media for darcy flow. Analytic
and very fast.

- tests/darcy_p0p1cv_pressBCinlet_rotated_domain. Similar to the above
but with a rotated domain as the pressure BCs have a normal component
so worth checking non axis aligned integrals. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_pressBCinlet. Similar to above but with
CG pressure and CV continuity. This is added as not many CG pressure
BC test cases appear to exist. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain. Same as
above but rotated. Analytic and very fast.

- tests/darcy_p0p1_test_cty_cv_velBCinlet. Similar to above but with
weak velocity BC instead. Just to test a combination of two previously
merged branches and the test case is very similar to above.

3883. By Rhodri Davies on 2012-01-01

Merging Stokes Options Check branch into trunk. Thanks to skramer for review.

Happy New Year to all!

3884. By Rhodri Davies on 2012-01-01

Some string continuation fixes which will hopefully placate the cx1 build. Also removed unused variable ncars from populate state.

3885. By Brendan Tollit on 2012-01-02

Relax two tolerances (1.0e-09 to 2.0e-09) from a previous merge that seems to affect the i386 builds
for the 3D cases of these tests.

No idea why this is a problem for 3D (and not 1D or 2D) and for i386 and not gcc.

3886. By Jon Hill on 2012-01-03

Adding a try block around the import of scipy in fluidity state_types.py. This allows use of the majority of state_types functionality without scipy being installed, which is useful for HECToR

3887. By Jon Hill on 2012-01-03

Fix for geopressure_convergence test. This works on my 32 bit laptop

3888. By Patrick Farrell on 2012-01-04

Cian pointed out in lp:904816 that running the Burgers equation and
shallow water adjoint tests modifies some files that exist in bzr.
This commit should hopefully fix this. The issue is very minor (utf8
vs. utf-8 in the header of XML files, and minor variations in gmsh
versions).

3889. By Gerard on 2012-01-09

Bugfix for Bug #900774.

This is a work around for a known bug in the advance option when using stream IO for formatted reads:
http://software.intel.com/en-us/forums/showthread.php?t=101333

The exact intel compiler version may have to be tweaked as the Intel fix has not yet been released.

3890. By Axelle Vire on 2012-01-09

Merge from lp:~avire/fluidity/fix-absorption-cg-bug.
This fixes the bug https://bugs.launchpad.net/fluidity/+bug/897358 and tests it.

3891. By Axelle Vire on 2012-01-09

Change the test_source_absorption from short to medium.

3892. By Simon Funke on 2012-01-10

Loosing tolerance for the turbine tests a bit to make buildbot happy.

3893. By Patrick Farrell on 2012-01-11

Remove the time_averaged_vector_squared tensor diagnostic algorithm, as it
was broken for a while and didn't error; it is untested. Resolves lp:788620.

3894. By Simon Mouradian on 2012-01-11

merge salt finger fix

3895. By Jon Hill on 2012-01-11

Relaxing tolerences for two tests for the 32-bit issue.

3896. By Jon Hill on 2012-01-11

Adding back in a bug fix I did ages ago. It was removed when I incorporated changes from a 3rd party into the biology python module. Must be more careful in future and construct a test of this

3897. By Patrick Farrell on 2012-01-13

Rename libadapt's fdebug.h => fdebug_adapt.h. This fixes the build error caused when
a user has -Ifluidity/include in front of libadapt's -I../include. Fixes lp:915907.

3898. By Patrick Farrell on 2012-01-13

Forgot to commit the Makefile.in.

3899. By Tim Bond on 2012-01-18

DOCUMENTATION CHANGE ONLY.

Updating the manual to reflect the new system of versioning, the new system
of tarball generation, and new system of obtaining tarballs. This removes
references to the stable release, replacing them with 4.1, which will be
formally implemented as of 4.1.3.

3900. By Adam Candy on 2012-01-18

Removal of the named solver reordering options.
Notes have been added to both the schema and Solvers.F90 stating the reordering mechanisms have not been tested, and how to use them if you want to try them out.
This addresses bug lp:818015.

3901. By Tim Bond on 2012-01-25

One-line fix in build system only, targetting bug lp:921464 by adding an
extra install line in the 'make install-diamond' target of Makefile.in
to install dxdiff.

3902. By Patrick Farrell on 2012-01-25

Incorporate the libtool config files from Fedora 15, to update libtool
to a more recent version. Our old libtool was too old to correctly
recognise gfortran 4.6. Fixes lp:918235.

3903. By Gerard on 2012-01-25

Should return 0 rather than non-zero when the -V option is specified.

3904. By Adam Candy on 2012-01-25

A help call to Fluidity also does not now complete with a non-zero error.

3905. By Florian Rathgeber on 2012-01-25

Sync libspud with spud trunk r512.

3906. By Tim Greaves on 2012-01-31

Fix in build-system only - targetting bug lp:921608 filed against
spud but also affecting Fluidity, which should fix the install-diamond
target in the Fluidity makefile. This brings libspud in Fluidity in
sync with r513 of lp:spud.

3907. By Tim Greaves on 2012-02-01

Commit purely in build system, reverting the growing untidyness of conflicting
prefixes trying to make a diamond install work from the libspud in Fluidity.
This goes back to a reconfigure with the user's prefix, then installs, then
reconfigures with the prefix that Fluidity expects to be used internally. This
ought to resolve problems with the spud schemas not getting installed.

3908. By Christian Jacobs on 2012-02-06

Merging fix-multiple-material-phase-python-diagnostic branch revisions r3897 to r3902 inclusive into trunk.

Summary of change:
==================
We now pass the whole states array to calculate_diagnostic_variable if a Python diagnostic algorithm is used to compute the VelocityAbsorption, VelocitySource or Viscosity diagnostic fields. Otherwise, we just pass the submaterials array. Brendan has created a test case for this (pydiag_material_phase_support_multiple_velocity_children).

This resolves bug lp:800126.

3909. By Tim Greaves on 2012-02-07

Commit purely in build system. A path error in the Makefile.in spotted by Simon
Mouradian means that the flml schema was being installed in ~/.schemata rather
than into ~/.diamond/schemata. This fix modifies the path to point to the
correct location.

3910. By Brendan Tollit on 2012-02-08

Summary:

This merge will permit the use of P2 control volumes for a tetrahedral mesh.
This extends the CV faces routine to add the p2cv tet shape functions.
The method to form the CV mass matrix is then changed to use a sub CV
analytic approach (for p1 and p2 elements one can calculate by hand the
sub CV fraction of the parent element) rather than lumping a FEM mass
matrix. Routines are added to get the cv mass and compute the cv mass.
Everywhere we could find that required a cv mass matrix now calls these
routines. Diagnostic routines to output the CV mass and the lumped FE mass
have been added. The capability to include the density in the
GridReynoldsNumber diagnostic has also been added. These last two
changes are included in the manual.

11 test cases are added (basic, advection, darcy, mmat, mms). A couple
of test cases where changed to not output cv stats as they didnt need it.

This permits p1dgp2 for multimaterial (although may well be buggy linked to
other features) and is known to not work with conservative bounded mesh
to mesh interpolation.

This was worked on by Cian Wilson and myself and reviewed/checked by Axelle.

3911. By Tim Greaves on 2012-02-08

Fix for lp:928771 , making an install of libspud's spudtools into the root
Fluidity directory a depenency of 'make mediumtest', and adding a rule to the
root Fluidity Makefile.in to this end.

3912. By Jon Hill on 2012-02-08

Trivial fix to the two sinking velocity test so that the user is not asked which algorithm GMSH uses. Simple command line fix is all that is requried.

3913. By Jon Hill on 2012-02-08

Reverting accidental changes to ocean biology

3914. By Brendan Tollit on 2012-02-09

Relax two tolerances for the max GridReynoldsNumber from 1.0e-02 to
2.0e-02 for these two medium test multimaterial cases, to hopefully
make the i386-gcc4 build green.

3915. By Brendan Tollit on 2012-02-09

Relax tolerances on these two test cases slightly more again from
2.0e-02 to 2.5e-02.

3916. By Tim Greaves on 2012-02-09

One-line fix to build system only; correcting 'FLINKER' to 'FLLINKER'. Thanks
to Zhenhua Lin for spotting this one!

3917. By Adam Candy on 2012-02-09

A fix to ensure 'make schema' succeeds on a fresh checkout.

This commit adds a Makefile directive to install spud-preprocess, which is required by the 'schema' directive.
Note that this new directive does not interfere with the build of libspud, which is performed within its own directory.
It is only called when required by 'make schema'.

3918. By Tim Greaves on 2012-02-10

One-line fix in build system only, replacing a plain copy of bin/ in the
install target of the root makefile with a find-based copy to deal with
the appearance of subdirectories which are not wanted in the final install
but which cause non-zero exit states of the install.

The specific case being dealt with here is the tests/ directory made by the
unittest target of the makefile, as spotted by Lin Zhenhua.

3919. By Tim Greaves on 2012-02-10

(Hopefully) minor fixes to the diamond_validation test, removing the dependency
on gtk.gdk which doesn't look to be necessary in the test (and which breaks on
systems with no DISPLAY available), and adding the local copy of the diamond
source to sys.path to allow diamond.debug and diamond.schema to be imported in
situations where the system doesn't have a central diamond install.

3920. By Tim Greaves on 2012-02-10

Bringing Fluidity's version of libspud up to date with r.515 of the
libspud trunk.

3921. By Adam Candy on 2012-02-10

Correction of a mistake in a logging message of the gyre_2d test.

3922. By Adam Candy on 2012-02-10

A new test case: a P1DG-P2 version of the 2d gyre (gyre_2d).

Fangxin asked me how to set this up yesterday for Travers, so I thought it's worth adding it to the tests.
The diagnostic needed adjusting to account for the existance of negative velocities in the y-coordinate direction along the Western boundary (since we apply the Dirichlet no-slip conditions weakly here).

3923. By Adam Candy on 2012-02-11

Contributes to fixing bug lp:920944.
The netCDF reader used when reading a file specified under ../from_netcdf/file options now gives more feedback when the point requested is not present. It reports the coordinates and filename to facilitate easier debugging.

This fixes the problem reported by Yi HanSoo, where the error from the reader was ambiguous as to which points could not be located in which file. We may need to make similar changes to other instances of NetCDF reads in the code, or better, unify them all into one routine.

Merge request made and approved at:
https://code.launchpad.net/~fluidity-core/fluidity/netcdferrorreporting/+merge/92443

3924. By Adam Candy on 2012-02-11

Further contributions to fixing bug lp:920944 (so it can be closed).
Two further instances of NetCDF readers now error more verbosely when the required data is not found. The reports include the coordinates (in both Cartesian and spherical, where applicable) and filename to facilitate easier debugging.

These additional changes are very similar to the change that has already gone in with the approved merge request at:
https://code.launchpad.net/~fluidity-core/fluidity/netcdferrorreporting/+merge/92443
so I have pushed them straight to the trunk.

3925. By Adam Candy on 2012-02-13

Fix to bug lp:922025. This tidies up the manner in which the pressure initial condition is set from a free-surface height field provided in a netCDF file.

Summary of changes:
 - The special case, legacy method has been removed.
 - If the input format is not recognised, the code will abort with an error detailing the problematic field and the schema path of the invalid format.
 - The 'standard' format is now named 'raw' (instead of 'other'). This is the case where no post-processing is performed on the source data, and the field is initialised by the raw data present in the netCDF file.
 - Included logging messages to indicate which field is being populated and the source netCDF file (logging level 1).
 - Logging of the post-processing applied (logging level 3).
 - For the pressure boundary condition set from a free-surface field in particular, logging level 3 will report the min and max values of the field before and after the conversion.
 - A new section in the manual describing how to initialise the pressure field by specifying the free surface height. Removal of the description of the legacy method to achieve this. Other corrections to nearby sections.

There are currently no tests, examples or longtests that are affected by the above change, but possibly user-generated initialisation flml files that do use the legacy implementation of the pressure initial condition set from free-surface netCDF. In this case, the additional logging added with this change should make it easy to identify the change required to adopt the improved version of this import.

In short, the following example highlights the change necessary:
  The option path:
    /material_phase::Ocean/scalar_field::Pressure/prognostic/initial_condition::WholeMesh/from_netcdf/format/Free-surface height

  should become:
    /material_phase::Ocean/scalar_field::Pressure/prognostic/initial_condition::WholeMesh/free_surface/from_netcdf/format/raw

  (where the file location is specified under the 'from_netcdf' node in each case).

3926. By Adam Candy on 2012-02-13

New test case: netcdf_read_free_surface
This creates a valid netCDF file containing a free surface height in the domain x,y in [-1.2,1.2]. This is read in by Fluidity as a free surface and applied as an initial condition on pressure. From this, Fluidity calculates the free surface as a diagnostic, which is provided in the outputted VTU file, alongside a field containing the analytic function used to create the netCDF file in the first place (here described in height.py). The test compares the difference between these two fields, which should be small, due to the interpolation of the netCDF file.

In summary this case tests:
 - NetCDF reader
 - Reading of a free surface and conversion into a pressure initial condition
 - Creation of the free surface diagnostic from the pressure field (including its extrapolation down from the top surface)
 - Creation of prescribed scalar fields from Python

3927. By Adam Candy on 2012-02-13

A new test case: netcdf_read_errors
This checks the errors reported when an attempt to read an invalid NetCDF file is made. Here 'invalid' refers to a file that is itself an invalid, non-conforming NetCDF file, or is of a format not supported by Fluidity.

Five types of NetCDF files are checked:
 - One which does not provide all the data requested,
 - One is missing a dimension (i.e. it provides only a 1d field),
 - One is missing the main variable,
 - One contains an incorrectly-labelled dimension (as 'lat', not 'x'), and
 - One has the variable incorrectly labelled (as 'height', not 'z').

3928. By Brendan Tollit on 2012-02-15

This merge will tidy up, improve, test and document the single phase porous
media Darcy flow model. This will sync the porous_media options model with
all the darcy_... test cases that have appeared over the past couple of months.
The model that used to be initiated by the porous_media option is effectively
removed as it contained code that was never finished, consistent, tested or
documented and was made by people that are not involved in Fluidity any more.
It is now favoured to define the Darcy velocity absorption field via python
diagnostics rather than a hard coded internal algorithm. Sections are added
to the manual to explain the model equations, discretisation and necessary
options.

Old tests cases are removed and new ones added.

3929. By Brendan Tollit on 2012-02-15

Correct a variable name in a ele_ngi check that was in a #ifdef DDEBUG
in assemble/Advection_Diffusion_CG.F90

3930. By Tim Greaves on 2012-02-16

Minor Intel-specific fixes, reducing line lengths which overrun in the
middle of quoted strings.

3931. By Stephan Kramer on 2012-02-20

Bug-fix for SampleNetCDF which is used in tidal boundary conditions. This bug would occur whenever a node was interpolated that is located exactly on one any of the constant longitutde or latitude lines of the input netcdf grid. In this case instead of interpolating between the value of the 2 nearest nodes, the value was extrapolated.

This code clearly needs some more tests. Coming up...

3932. By Axelle Vire on 2012-02-21

Loose the tolerance on this test slightly (with Brendan's consent) since it makes the trunk fails locally on my Maverick machine.

3933. By Axelle Vire on 2012-02-21

Relax a tolerance on this test with Brendan's consent. The test was passing locally but failed on independent buildbot branches.

3934. By Axelle Vire on 2012-02-22

Correct Brendan's option check regarding the use of CV for continuity with two-way coupling. Apologies that it hasn't been done before the merge.

3935. By Axelle Vire on 2012-02-22

Forgot an option check. Sorry for the noise.

3936. By Patrick Farrell on 2012-02-23

Remove the finest simulation (_E) from the Burgers MMS replay test.
This should hopefully speed up the test so that it falls within
buildbot's 20-minute limit ..

3937. By Brendan Tollit on 2012-02-23

Remove the porosity from the DG_CourantNumber and ControlVolumeCFLNumber
diagnostic fields. Add three new diagnostic fields being
InterstitialVelocityCGCourantNumber, InterstitialVelocityCVCourantNumber
and InterstitialVelocityDGCourantNumber which are only for porous media
flow where the Darcy velocity (= porosity * interstitial velocity) is
solved for.

These are magic named diagnostic fields rather than algorithm names due to
the way the courant number fields are already used.

Everywhere in the options where ControlVolumeCFLNumber can be selected
(such as CV face value schemes) permit the selection of
InterstitialVelocityCVCourantNumber.

For CV (or coupled CV) discretised fields there is a restriction that all
the 'porosity' fields used must be the same with the same theta averaging
also, due to the way this is obtained from one call in Field_Equations_CV.
The default options when used satisfy this constraint.

For adaptive time stepping also permit the selection of
InterstitialVelocityCGCourantNumber. If an InterstitialVelocityCG/CV
CourantNumber is used for adaptive time stepping then always force its
calculation in this routine as the options here may be different to a
field that may exist in the state.

For a scalar field discretised with DG add the small new capability to
select which Courant number field to use for the temporal sub cycling
with the available selections being DG_CourantNumber and
InterstitialVelocityDGCourantNumber. To avoid changing all the test cases,
by default if not activated this option will do as before and use
DG_CourantNumber.

Add a copy of a short darcy flow test case that now discretises the tracer
with DG using subcycling and adaptive time stepping related to the new
interstitial velocity courant number fields.

Adjust all relevant darcy test cases to also test this.

Add a routine in assemble/Porous_Media.F90 to find the theta averaged
porosity field, to avoid a code copy in the adv-diff CG/CV/DG routines.

Add to the manual a description of the three new diagnostic fields and
adjust the configuring Porous media flow section as required.

3938. By Brendan Tollit on 2012-02-23

Small correction to last merge. Only use the relevant variables from
Global parameters and correct some repetetive English in a couple of
FLExit's.

3939. By Brendan Tollit on 2012-02-24

Relax two tolerances for the intel build for this recently added test.

5.0e-03 --> 8.0e-03 (integral of tracer)

5.0e-04 --> 6.0e-04 (min of tracer)

3940. By Jon Hill on 2012-02-29

Shortening gls-KP test by remove one of the four tests within this

3941. By Rhodri Davies on 2012-02-29

Change solver and steady state tolerances to speed this test up. Now runs in less than 1 minute on my system.

3942. By Jon Hill on 2012-02-29

And four became two. Hopefully this passes in less than 10 mins on buildbot

3943. By Tim Greaves on 2012-02-29

One-line fix in buildsystem only targetting bug lp:920446.

This checks for an empty DESTDIR and replaces such with / rather than supplying
an empty --root= which was causing problems with distutils.

Thanks to Patrick, Florian, and Stephan for input.

3944. By Jon Hill on 2012-02-29

Shortening the test again. Should be half the time now

3945. By Christian Jacobs on 2012-03-01

Small changes to the mphase_tephra_settling test to reduce the runtime. Passes locally.

3946. By Tim Greaves on 2012-03-02

Merging upstream updates to spud into the Fluidity trunk; this brings
the local copy in Fluidity in sync with r.520 of the spud trunk.

3947. By Tim Greaves on 2012-03-02

Following the update of spud in Fluidity to r.520 of lp:spud, this commit
updates all options files in the Fluidity trunk, resulting from running
tools/update_options.py.

3948. By Tim Greaves on 2012-03-11

r.3947 failed; the system it ran on had an old version of spud/diamond
installed. The Spud 1.1.4 (lp:spud r.520) packages are now out and this is the
result of running tools/update_options.py on a system with those packages
installed.

3949. By Lawrence Mitchell on 2012-03-14

Avoid memory blowup in update_options

When schema files are read in diamond, the tree-structure is memoised
internally to improve performance. But this comes at the expense of
additional memory use. Previously in update_options we would read the
relevant schema file for every file we were checking. This would
never hit the previously memoised values, since object ids were
different. Instead, build a dictionary of pre-parsed schemas and
just grab the relevant one when checking a file.

3950. By Stephan Kramer on 2012-03-14

This deals with bug https://bugs.launchpad.net/fluidity/+bug/947049

"The script gmsh2triangle tries to identify isolated nodes that are not attached to any elements (isolated nodes) and leaves these out of the output triangle files. This requires a renumbering of the nodes. However this renumbering is only applied to the .ele file but not .edge/.face file. Thus this renumbering only works correctly when the isolated nodes have node numbers that are bigger than any of the boundary nodes. Otherwise the node numbering of the surface mesh and the interior mesh no longer correspond leading to very onobvious errors in fluidity."

Fixed by no longer applying any renumbering. This was decided to be the best fix as gmsh meshes are now also be readable directly by fluidity, in which case isolated nodes would be read as well. So if isolated nodes are present in the input gmsh .msh it's up to the user to fix.

3951. By Simon Mouradian on 2012-03-15

A tool to radially scale .vtu files on the sphere.

3952. By Cian Wilson on 2012-03-16

Updating to new base schema.

3953. By Simon Mouradian on 2012-03-16

fix for create_param_sweep

3954. By Brendan Tollit on 2012-03-20

Precise build fix for the cv_mass_matrix_check_p1 and
cv_mass_matrix_check_p2 test cases.

Issue was actually gmsh. These test cases require a one
element mesh to do some very basic checks of the CV mass
matrix. The gmsh with the precise buildbot doesnt seem to
be able to do that anymore with these .geo files. Therefore
commit the very simple .msh files that do work.

Adjust the test Makefile's as required.

3955. By Brendan Tollit on 2012-03-24

Remove options to use the CV Ultimate limiter from metric advection as the
code does not permit it. Add 2 tests for the advection of a scalar field
using the CV Ultimate limiter as there were none. Both use a P1 FE face
value before limiting. One advects a step change, the other is a copy of
an mms case. Also add the mms case for the CV Sweby limiter for comparison.
For the step advection Ultimate is far better than Sweby but not quite as
sharp as HyperC, which is expected. For the mms test case it was found
necessary to apply a strong BC for the Ultimate limiter case so that
it will converge to a steady state - this however reduces the order
of convergence achieved.

No actual code changed, only options schema and test cases.

tests/mms_tracer_cv_adv_steady/MMS_E_cv.flml is removed as it is not used.

Linked to bugs:

https://bugs.launchpad.net/fluidity/+bug/934250
https://bugs.launchpad.net/fluidity/+bug/934259

3956. By Brendan Tollit on 2012-03-24

Remove two test cases that have a memory leak on Precise. This is
linked to bug report:

https://bugs.launchpad.net/fluidity/+bug/963989

and branch:

lp:~btollit/fluidity/3mat_p1dgp2_test_cty_cv_memory_leak

3957. By Jon Hill on 2012-03-28

Minor chnage to configure to get the HECToR build working

3958. By Stephan Kramer on 2012-03-28

Fix for parallel scalar dg subcycling when configured with debugging. It tried to invert the local mass matrix of non-owned elements, which have not been assembled. Fixed by telling it not to do so. Added further asserts so that similar bugs will also be catched in 2d.

3959. By Stephan Kramer on 2012-03-28

Fix for precisison issues on precise. Cian helpfully provided a robust fix by reverting to an earlier version, but since I'm impatient as ever I preferred to keep the non-quadratic algorithm using lexsort. This means Cian may forever blame me if this breaks again.

3960. By Jon Hill on 2012-03-29

Minor fix to the recent dg subcycling merge. mass should not be altered in this routine, so the last if block should be inv_mass. Will fix HECToR and precise builds

3961. By Jon Hill on 2012-03-29

Loosening tolerences on this test to make it pass on Natty and Precise. A note in the XML says what the tolerences should be

3962. By Simon Funke on 2012-03-30

Loosen test tolerance for turbine_flux_dg_2d to make it pass with precise

3963. By Simon Funke on 2012-03-30

Fix the turbine_flux_dg_2plus1 test by loosing the tolerance slightly.

3964. By Simon Funke on 2012-03-30

Fix turbine_flux_penalty_2plus1 test by slightly loosing the tolerances

3965. By Simon Funke on 2012-03-30

Fix the turbine_flux_penalty_2d test by slightly loosing the tolerances.

3966. By Simon Funke on 2012-03-30

Fix shallow_water_adjoint_default_controls_2d by slightly loosing the tolerance.

3967. By Simon Funke on 2012-03-30

Fixing the wetting_and_drying_balzano1_cg test by choosing a slightly smaller timestep.

3968. By Simon Funke on 2012-03-30

Now this is annoying: Apparently I need to reduce the timestep from 600s to 300s to get the test that worked before running on precise!

3969. By Simon Funke on 2012-03-31

Some more tolerance loosening for turbine_flux_penalty_2plus1 and turbine_flux_dg_2plus1

3970. By Tim Greaves on 2012-04-02

Updating the femtools Makefile.dependencies which appears to have been missed
at r3958

3971. By Tim Greaves on 2012-04-02

Minor commit in build system and tests to ensure that 'make distclean' does in
fact return the build tree to a fully clean state.

3972. By Tim Greaves on 2012-04-02

One-line fix targetting bug lp:850727 which had only been partially
implemented. This should complete the fix.

3973. By Jon Hill on 2012-04-03

One of the MPI_AllReduces is only necessary when using additional adapt iterations. This commit fixes the code so that it only does the collective communication when the additional adapt iterations option is being used.

3974. By Jon Hill on 2012-04-06

Move all C/Fortran tools to using the c_iso_binding module. This prevents random errors that have been occurring on some build tests

3975. By Jon Hill on 2012-04-06

Fixing bug in vertical integration - wrong types being passed over from C

3976. By Rhodri Davies on 2012-04-12

Speed this test up further

3977. By Rhodri Davies on 2012-04-12

Merge in branch that remove leading whitespace in diagnistic dependency lists

3978. By Christian Jacobs on 2012-04-15

Merging lp:~fluidity-core/fluidity/bc_from_field branch revisions r3962 to r3975 inclusive into trunk.

Summary of change:
==================
Users can now enforce boundary conditions that use boundary data from another field. This can be done by selecting "from_field" (from the list of "constant", "python", ..., and "from_field" input types for the BC) and specifying the name of the 'parent' field which contains the boundary data.

A short test case called bc_from_field is included.

3979. By dujuan on 2012-04-16

Combining EnKF framework with ReducedModel

3980. By Tim Greaves on 2012-04-16

Reverting r.3979 until it's been through the review and test process; Juan Du,
please could you make a branch and commit these changes to that branch, then
let me know to make you a buildbot queue for it, and put in a merge request.

Thanks!

3981. By David Robinson <email address hidden> on 2012-04-18

A small edit to include velocity in GetFieldNames

3982. By Alexandros Avdis on 2012-04-25

Correction to rotate2ll.

3983. By Tim Greaves on 2012-04-30

Two minor fixes on the move to precise, entirely in test framework:

* netcdf_read_errors had a bash₋specific construct; making generic to
   allow it to be parsed with dash
* sloshing_tank needs a very small relaxing of one check; confirmed
   with the test owner (Simon Funke)

3984. By Lawrence Mitchell on 2012-05-10

Make linear system assembly local (non-communicating)

This merge brings in non-communicating assembly routines for all
equation forms. It also adds OpenMP threaded assembly for DG and CG
Momentum and CG advection diffusion. To enable the latter, configure
--with-openmp.

For non-communicating assembly, we assemble in halo elements for DG as
well as CG such that each process has complete information for its
local dofs. PETSc matrices are then assembled with the option
MAT_IGNORE_OFF_PROC_ENTRIES, which makes the PETSc assembly a no-op.

OpenMP threading makes the assembly loops iterate in colour order.
The correct colouring is stashed on the topology mesh between adapts.
A number of routines are not thread-safe. If you call them in
threaded code they will now abort, fix it!

3985. By Lawrence Mitchell on 2012-05-14

Fix null dereference when checking halo types

Logical evaluations aren't necessarily short-circuited when an
argument evaluates to .false., so

   if (associated(foo) .and. (foo%bar == baz))

may dereference foo even if it is unassociated. So don't do it.

3986. By Stephan Kramer on 2012-05-16

Removing three places with hard-coded values (two different ones) for the radius of the earth. This was all in dead code. Hard-coding this is not a good idea, as the mesh generator may have used another one.

Two more places remain: femtools/projection.cpp and ocean_forcing/climate_reader.cpp, these are used in determining what layer of 3D input data to use. This will be fixed (instead using the outer radius
of the original input mesh coordinates) in the rewrite of the netcdf readers.

3987. By Stephan Kramer on 2012-05-17

Fixing gmsh reader to read meshes on the sphere (i.e. meshes with topological dimension 2, and coordinate dimension 3).

Also adding a parallel spherical_patch test case that tests this reading, writing out to a checkpoint, and flredecomping to a different number of processes.

Some optimisation and cleaning up of Read_GMSH.F90:
- prevent reading of the gmsh file twice. This was because read_gmsh_simple() called identify_gmsh_file() to establish
various dimensions, which it can only do by reading the entire file. Then an appropriate shape functon was allocated
and read_gmsh_to_field was called which read the file again.
- the only version of read_gmsh_files to remain is the _simple version (as this was the only one in use). It
first reads in all the data, then establishes the necessary dimensions and then allocates the appropriate femtools
objects.
- the dimension of the read gmsh mesh is now established from the topology of the mesh, basically looking at the
highest dimension elements present, i.e. if there's tets or hexes it's 3d, otherwise 2d (1d gmsh reading is untested)
Previously this was done by looking at the z coordinate and seeing if it was small using some arbitray, hard-coded definition
of small. This is obviously not a good idea.
- the dimension of the coordinate field that is returned (X%dim) is in most cases chosen to be the same as
the topological dimension. Only with the /geometry/spherical_earth option this is always 3 so that a 2d spherical
horizontal mesh (to be extruded inside fluidity) can be read. This was broken before. The other codes that use
embedded meshes (shallow water) use a specific option that only works for triangle, so that functionality should
remain unchanged.
- remove identify_gmsh_files()
- because of this remove identify_mesh_files() from the mesh_files module
- in all subroutines of the mesh_files module make format a required argument. Before, if the argument was not present
a "guess" would be made. In some tools where no option tree is present (e.g. test_laplacian) this would always default
to "triangle". These therefore now directly call read_triangle_files. In all other cases the format should be worked
out. For output this can be a bit tricky. For checkpointing the same logic is followed as before (look at the option
under the external mesh). For adaptivity debugging meshes it now always uses gmsh.
- identify_mesh_files() was only used in populate_state, in the case of inactive process (flredecomp from less
to more processes), where inactive processed needed to know the dimensions of the mesh. This is now communicated
via a MPI_Scatter()

3988. By Jon Hill on 2012-05-18

Minor change to fldecomp makefile.in to get fldecomp to compile on HECToR. Almost got a fully working build...

3989. By Stephan Kramer on 2012-05-18

Fixing changed dependencies ('make makefiles') since my gmsh commit.

3990. By Stephan Kramer on 2012-05-18

This adds the -fno-realloc-lhs flag to FCFLAGS and FFLAGS during configure, when using a gfortran that accepts it. This is a work-around for memory leaks caused by gfortran >=4.6, see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=53389

3991. By Jon Hill on 2012-05-22

Adding the SVG Fluidity logo to the repository. Same license as the rest of Fluidity

3992. By Brendan Tollit on 2012-05-23

Add back two multimaterial test cases that are set up to test p1dgp2, with checks for volume fractions and velocity divergence.

These were removed because of a memory leak caused by gfortran which now has a fix in the trunk thanks to Sam and Stephan.

3993. By Jon Hill on 2012-05-27

Adding ability to run psuedo-1D columns in parallel by decmping in the vertical. Note that the user must supply the DistanceToTop and Bottom fields of they want them

3994. By Jon Hill on 2012-05-28

Removing warnings from FLML

3995. By Florian Rathgeber on 2012-05-28

Fix the with_zoltan=no configure option, which was without effect

3996. By Cian Wilson on 2012-05-29

Bug fix for intel compilers that now object to converting null pointers into objects in subroutine calls. Specifically in this case the ct_m divergence matrix.

3997. By Lawrence Mitchell on 2012-05-30

Further reduce recompilations in already built trees

include/version.h is now built with a new script bin/version-info.
Makefiles are updated to take advantage of this, and reworked such
that include/version.h is no longer a .PHONY target.

Version is determined in the following ways (in order of preference):

  - make PACKAGE_VERSION=foo
  - ./configure FLUIDITY_VERSION=foo
  - bzr version-info

Additionally, to avoid spurious rebuilds due to
Diagnostic_Fields_New.F90 changing, only overwrite it if the hash of
the existing version differs from the hash of the putative new
version.

3998. By Tim Greaves on 2012-06-06

One-line fix in build system only; the packager was giving unwanted
leading-path entries in the examples tarball, which this commit addresses.

3999. By Alexandros Avdis on 2012-06-08

Correction of errors affected by an earlier commit - Sould have been experimenting on unbound branch.

4000. By Jon Hill on 2012-06-10

Bug fix to period_averaged_scalar field diagnostic to work with adaptivity

4001. By Tim Greaves on 2012-06-13

Change in test tolerance only to relax for 32-bit Precise test platform.

4002. By Stephan Kramer on 2012-06-13

Loosing tolerance: maximum fs may occur one timestep later. Also corrected fs_right_max test.

4003. By Tim Greaves on 2012-06-13

Changes in documentation only: Updating manual to recent versions of the
supporting software documentation, drawn from the AMCG Wiki at:

 http://amcg.ese.ic.ac.uk/index.php?title=Fluidity_Supporting_Software

4004. By Samuel Parkinson on 2012-06-15

merge LES_fix branch into trunk into trunk - fixes bug #1000813

4005. By Tim Greaves on 2012-06-19

Commit purely in build system, adding a fluidity binary package dependency
on libspud-dev

4006. By Brendan Tollit on 2012-06-21

Add the capability to associate any scalar field with the perturbation density for a
fluids linear EoS to complement the Temperature and Salinity components. This is an
added generalisation of the coding and is useful for collapsed multiphase models
involving bubbles, perhaps other applications also.

The options for EoS are added to their own .rnc/g file and new options added.

The manual is slightly changed to include a description of the new functionality.

A test case is added which copies the case buoyancy_adjustment_dg_two_element and
uses a generic scalar field for the EoS rather than Temperature - gets the same
answer as expected.

4007. By Alexandros Avdis on 2012-06-22

Removed some hard coded instances of the Earth radius. The Earth geoid surface radius is now calculated from the mesh.

4008. By Simon Funke on 2012-06-22

Loosening the tolerance for the shallow_water_adjoint_default_controls_2d will hopefully make buildbot happy

4009. By Stephan Kramer on 2012-06-27

Change assert for descending node order with 1d adaptivity, to a check that is also run without debugging and add a helpful error message.

The descending node order for 1d adaptivity is a bit silly, that could easily be fixed if somebody's up for it.

4010. By Stephan Kramer on 2012-06-28

Fixing a few issues with cdisk_2d:
* the test advects a tracer under a prescribed velocity field. The velocity field however does not satisfy no-normal flow at the boundary, so boundary conditions for the tracer are needed.
* only testing for something that is also satisfied by the initial condition is never a good idea. Next to conservation, we now also test the tracer has moved by the right amount.
* adding a 'solvers converged' check pending a more general solution for this.

4011. By Lawrence Mitchell on 2012-07-17

Add checks for solver convergence to most tests

Tests that run fluidity in a meaningful manner should always test for
convergence of the solvers. This merge attempts to do that for most
short and mediumtests.

4012. By Jon Hill on 2012-07-17

Fix bug in libspud's python interface. Same fix has been applied in libspud trunk also

4013. By Lawrence Mitchell on 2012-07-18

Make petsc libs come before system libs in link line

When we discover a set of library flags in configure, we should
prepend, not append to LIBS. The linker looks at library paths in L2R
order on the command line, so the more specific the path we discover,
the closer to the front of the command line it should be.

4014. By Tim Greaves on 2012-07-18

Changes purely in tests reassigning tests from hhiester to new owners.

4015. By Gerard on 2012-07-20

bugfix: FCFLAGS was being clobbered with non-existent OLD_FCLAGS

4016. By Tim Greaves on 2012-07-25

Changes purely in example XMLs: changing length to 'long' in all tests to
  allow them to successfully run as 'longtests' in buildbot.

4017. By Tim Greaves on 2012-07-25

Small change in the makefile of this example to bring the default
 'make input' called by the testing process in line with the rest
 of the example running on 32 cores.

4018. By Tim Greaves on 2012-07-27

Change purely in build infrastructure of an example: for automated running of
  the BFS3D example, a default NPROCS=8 is supplied.

4019. By Brendan Tollit on 2012-07-27

Add the capability for robin BC for control_volume discretisation
(not coupled_cv) only if the diffusion term uses ElementGradient.

Add a simple analytic test case that need only be 1d as the added
code is independent of dimension.

Adjust comments in the schema for this.

Add small sections to the manual.

4020. By Rhodri Davies on 2012-07-31

Bug fix for cases using P2CV and limit theta. Supplied by Cian.

4021. By Jon Hill on 2012-08-02

New GLS updates which include:
 - Ocean parameterisation which parameterises internal wave breaking, adding some surface TKE into deeper water
 - Updates to assembly routines which are now mostly element-wise
 - Updates to tests
 - Manual updates

4022. By Jon Hill on 2012-08-02

Adding tests for PAR and updating manual

4023. By Rhodri Davies on 2012-08-06

Bug fix for CV, when a diffusion term is present. Again supplied by Cian.

4024. By Christian Jacobs on 2012-08-09

Merging compressible-multiphase branch revisions r3955 to r4067 (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 CompressibleContinuityResidual 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.

4025. By Christian Jacobs on 2012-08-09

Moved the compressible multiphase MMS tests to the longtests directory.

4026. By Christian Jacobs on 2012-08-09

- Updated the tephra_settling example's .flml file to use the new multiphase_interaction schema structure.
- Corrected a small typo in the manual.

4027. By Lawrence Mitchell on 2012-08-13

Fix DG momentum assembly corner case

Since r3984 we would assemble elements in DG momentum if
element_neighbour_owned returned true. However, for some corner cases
element_neighbour_owned need not return true if the element is owned
(consider the case of a process with one owned element). In this
situation we would not correctly assemble all elements. To fix this,
change the test for element assembly to directly check for element
ownership.

At the same time, update the documentation for element_neighbour_owned
to note this issue.

4028. By Rhodri Davies on 2012-08-16

Merge relax_halo_verifies branch into trunk. Thanks to dham for review.

4029. By j-bull08 <email address hidden> on 2012-08-17

Removing duplicate file in 3D BFS example.

4030. By Lawrence Mitchell on 2012-08-17

Fix CG momentum assembly on intel fortran

This only affects the have_les and .not. dynamic_les code path.

4031. By Rhodri Davies on 2012-08-17

Merge zoltan configure check for .mod files into trunk. Thanks to steph for help and jon for review.

4032. By Christian Jacobs on 2012-08-17

Added solvers_converged pass_tests to all multiphase tests.

4033. By Christian Jacobs on 2012-08-19

Added a solvers_converged pass_test to the tephra_settling example.

4034. By Samuel Parkinson on 2012-08-19

Updated Momentum_CG.F90 so that upwind stabilisation can be used with partial or full stress form. The issue is that when using partial or full stress form we tend to set all components of viscosity to the same value. This makes the tensor determinant zero and we cannot invert it. This fix simply sets all non-diagonal terms to zero for the purpose of handling stabilisation when using partial or full stress forms. This isn't ideal as it ignores the fact that we may have non-isotropic viscosities, but these aren't handled properly for these stress forms anyway - so this is a bigger issue that may need fixing later.

The same issue arises when we try to calculate the GridReynoldsNumber diagnostic field and the same fix has been applied here.

4035. By Samuel Parkinson on 2012-08-21

bug fix of SUPG stabilisation in Momentum_CG.F90

4036. By Tim Greaves on 2012-08-22

One-line fix, setting the correct flml filename in the test xml.

4037. By Adam Candy on 2012-08-24

Removal of a runtime-generated test file.
jhill - please check!

4038. By Samuel Parkinson on 2012-08-28

Added options for how the bed shear stress is calculated. Yuo can now either:
- use drag formulation as before
- calculate the shear stress based upon the gradient of the velocity field at the bed
- use the shear stress estimate generated by the wall boundary layer boundary condition algorithm

4039. By Samuel Parkinson on 2012-08-31

change to lock_exchange example so that, when run as a test, it runs to 30 seconds so that the tests are checked properly

4040. By Christian Jacobs on 2012-08-31

Merging compressible-multiphase branch revisions r4068 to r4077 (inclusive) into trunk.

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

- Added multiphase support for stress_form, partial_stress_form and isotropic diagonal_viscosity in Momentum_CG.F90. The stress_form will be important in compressible multiphase flow simulations.

- Added an MMS test (mphase_mms_p2p1_stress_form) for the stress_form changes.

4041. By Christian Jacobs on 2012-08-31

Small edit to the manual to generalise the form of the multiphase stress term.

4042. By Xiaohu Guo on 2012-09-01

  1. thread matrix assembly part for Advection_Diffusion_DG.F90
  2. MASSLUMPED_RT0 diffusion scheme is not supported, simply set threads back to 1 if MASSLUMPED_RT0 diffusion scheme is being used, and there will be a warning message.
  3. thread matrix assembly part for Advection_Diffusion_FV.F90

4043. By Christian Jacobs on 2012-09-02

Merging k_epsilon_extend_and_fix branch into trunk.

Changes to k-epsilon model:
===========================
- The code no longer lumps the mass matrix when calculating source terms for the k and epsilon fields, allowing element pairs other than P1-P1 to be used.
- k-epsilon field calculation moved to the beginning of the time-step and calculated with all other diagnostic fields.
- Added source term to momentum equation (-2/3 grad(rho*k)).
- Updated the manual documentation for the k-epsilon model.
- Users must now use partial-stress or full-stress form for the stress term in the momentum equation. This is the correct formulation when using spatially
varying viscosities. The ScalarEddyViscosity from the model is added to all components of the stress tensor to reflect this.
- Fixed calculation of c_eps_3.
- Added lots of debugging options.
- Improved the options_check subroutine.
- Changed how prescribed source terms are implemented for the k and epsilon fields so that multiple non-linear iterations can be carried out
with prescribed source terms.
- Altered calculation of Reynolds Stresses.
- Added a new equation_type (KEpsilon) for the k-epsilon equation. This includes the density in the advection-diffusion equation for the k and epsilon fields. Also made many other changes to the code to take account of density variations. Supports CG and CV discretisations.
- Made a few preliminary steps for adding multiphase support to the k-epsilon model.
- k and epsilon fields are now clipped at Gauss points rather than nodes.
- Moved element assembly to its own subroutine to avoid dynamic memory assignment.
- Changed k and epsilon diffusivity tensor to (nu + nu_t/sigma_*) instead of (nu_t/sigma_*)
- The turbulent buoyancy term is now calculated from the VelocityBuoyancyDensity rather than a summation of scalar fields that have been identified as
being buoyant. This means that any buoyant field will automatically result in a buoyancy term in the k-epsilon model (unless the entire term
is disabled in debugging_options)
- Removed several unused variables.
- Added the buoyancy term to the k-epsilon model.
- Added damping functions for the low_Re k-epsilon model.
- Changed epsilon low_Re boundary condition to d(eps)/dx = 0 as it's simpler than what we had before and still valid (Wilcox 1998).
- Fixes to high-Re boundary conditions.
- Removed length scale limiting.

Changes to tests and examples:
==============================
- The existing k-epsilon MMS test (k-epsilon_mms_cg_p1p1) has been removed and replaced with one that uses a P2-P1 discretisation and runs to steady state.
- Three additional tests have been added for the low_Re boundary condition, and the LinearMomentum Velocity equation_type with CG and CV discretisations.
- Tests include the turbulent buoyancy term and scalar field turbulent diffusion.
- A unit test has been added for the tensor_inner_product function.
- Fixed the backward facing step 2D example based on the changes to the k-epsilon model.

4044. By Christian Jacobs on 2012-09-02

Moved tests/mms_rans_p2p1_keps_linearmomentum and tests/mms_rans_p2p1_keps_linearmomentum_cv to the longtests directory.

4045. By Christian Jacobs on 2012-09-03

Buildbot's makefiles-x86_64 queue has been failing since r4042. These updates from make makefiles will hopefully fix this.

4046. By Adam Candy on 2012-09-03

Wetting and Drying bug fix - the material phase had an assumed name of 'water'.
(It is now just assumed to be the first material phase)

4047. By Christian Jacobs on 2012-09-04

Bug fix following the merge of k_epsilon_extend_and_fix. k-epsilon boundary conditions are now applied to all states that have k-epsilon enabled, not just to state(1).

4048. By Tim Greaves on 2012-09-05

As discussed in today's Fluidity developer meeting, changing this example to
default to using 64 cores to reduce runtimes for automated testing.

4049. By Rhodri Davies on 2012-09-12

Merge python_state_csr_matrices branch into trunk. Thanks to Stephan and Cian for help and to Patrick for review.

4050. By Florian Rathgeber on 2012-09-14

Purge empty directories from the repository and create them as needed.

There are a few empty directories in the repository, which bzr tracks,
other DVCs do not (git, hg). Instead of relying on those directories
being present, these are created by the respective Makefiles.

4051. By Adam Candy on 2012-09-17

Bug fix for the buoyancy adjustment by vertical mixing scheme.
The Galerkin projection to a continuous space is now honoured.

4052. By Adam Candy on 2012-09-17

Compatibility fix for a newer Gmsh version - for the buoyancy adjustment test cases.

4053. By Adam Candy on 2012-09-17

Tagging of the buoyancy adjustment test cases.

4054. By Samuel Parkinson on 2012-09-18

- added new reentrainment boundary condition based on Garcia's paper from 1991 which is the most generally used relationship.
- added hindered settling and sediment concentration dependent viscosity.
- changed how sediment fields are initialised so that they fit more generally with how other fields work in Fluidity.
- manual updates.
- new mms test for all sediment related code.

4055. By Rhodri Davies on 2012-09-18

Merge field weighted partition branch into trunk. As usual, thanks to Stephan for help and advice. Thanks to Jon for review.

4056. By Alexandros Avdis on 2012-09-18

Merge-in of DG_Neumann branch.

4057. By Tim Greaves on 2012-09-19

A commit which should be entirely limited to the build system within the tests
hierarchy and the 'clean' subsystem, attempting to catch up with housekeeping
and ensuring that 'make distclean' from the fluidity base directory reverts any
changes made by the build and testing process.

Major changes here relate to adding/correcting 'make clean' entries, and moving
some tests from copying/changing committed .bml/.oml/.swml files to copying
templates and then modifying the copies.

4058. By Mark Goffin <email address hidden> on 2012-09-19

This commit fixes local edge numbering for quad elements.

4059. By Tim Greaves on 2012-09-19

Having added in a clean-up line for the generated 'version-info' file, it
appears that tools/Makefile requires version-info to exist for the
install-version-info target, but doesn't know how to make it itself.

For the moment, I've backed out the cleaning line and will leave this for a
more comprehensive fix.

4060. By Tim Greaves on 2012-09-19

Thanks to Lawrence for spotting this; my error in r4057 was removing
tools/version-info with 'make clean' whereas it should only be removed
with 'make distclean'. This commit fixes that, adding a 'distclean' target
to the tools makefile to handle removal of version-info.

4061. By Adam Candy on 2012-09-19

Introduces a diagnostic for the buoyancy adjustment by vertical mixing scheme.
Logging for the scheme has been tightened up and generally improved.

4062. By Mark Goffin <email address hidden> on 2012-09-19

This should fix the tests hh_bdfm1 and hybridized_helmholtz that broke after revision 4058.

4063. By Mark Goffin <email address hidden> on 2012-09-20

A minor fix for degree zero quad elements to fix the test - Stokes_mms_q1p0.

4064. By Samuel Parkinson on 2012-09-20

Hopefully these updates will mean that the make makfiles test passes

4065. By Adam Candy on 2012-09-20

New option in gmsh2triangle to interpret a 3d input mesh as a 2d mesh on a spherical shell.
This avoids having to apply sed to the node file.

An improvement would be to have Fluidity detect this independent of the dimension in the node file.
(as it does with the binary gmsh format?)

4066. By Adam Candy on 2012-09-20

Restrict which pylab functions are loaded in the the radial scale tool.

4067. By Adam Candy on 2012-09-20

radial scale tool - no need for pylab dependencies, available in math.

4068. By Simon Mouradian on 2012-09-21

Updates to the manual for changes to gmsh2triangle

4069. By Adam Candy on 2012-09-22

Bug fix to buoyancy adjustment diagnostic.

4070. By Adam Candy on 2012-09-22

Additional documentation and optimisation for the buoyancy adjustment diffusion diagnostic.

4071. By Adam Candy on 2012-09-22

Optimisation in the previous commit requires this set operation.

4072. By Simon Mouradian on 2012-09-25

Fix for rotating channel example

4073. By Simon Mouradian on 2012-09-25

If dev version of libadjoint is installed, the configure script is now smart enough to realise and fluidity is built without adjoint support.

4074. By Matthew Piggott on 2012-09-27

Update the adaptive driven cavity test flml so it is (more) consistent with the fixed mesh driven cavity example.

4075. By Rhodri Davies on 2012-09-28

Stokes example and place holder in manual.

4076. By Samuel Parkinson on 2012-09-28

Tidied up time discretisation of the k-epsilon model and some other bits. See the manual edits for details on this. The main changes to the code are:

- Calculation of keps diagnostic now done twice during each non-linear iteration. First before advection diffusion of scalar fields, then again before the momentum solve. This makes it possible to have a consistent time-discretisation across the equations.
- Non-linear velocities are used in diagnostics and k and epsilon values are averaged between the most recent value and the value at the end of the last time step using a user defined theta value
- changes to how the isotropic Reynolds stresses are calculated. They are now absorbed in to the modified pressure.
- rewritten low-reynolds model which was wrong before

4077. By Matthew Piggott on 2012-10-01

Example failing with a GMYBAD error. Fiddle with adaptivity and zoltan options based on the TEST version of this flml file.

4078. By Matthew Piggott on 2012-10-01

Update the title of the streamfunction diagnostic in the pass chacks.

4079. By Rhodri Davies on 2012-10-02

Merge in changes to scaled pressure mass matrix assembly.

4080. By Christian Jacobs on 2012-10-02

Small bug fix. Removed the requirement to have a VelocitySource field when using the k-epsilon model. The field was previously needed to add in -(2/3)*rho*k, but this is now part of the pressure term.

4081. By Rhodri Davies on 2012-10-03

Merge Full Projection fix into trunk. Thanks to Stephan for review.

4082. By Michael Lange on 2012-10-04

Adding PETSs 3.3 support and dropping support for 3.0.

4083. By Rhodri Davies on 2012-10-05

XML for stokes example. Should be good to go.

4084. By Matthew Piggott on 2012-10-07

Fix my diamond validation error.

4085. By Rhodri Davies on 2012-10-09

Commit mesh files for debugging purposes only.

4086. By Rhodri Davies on 2012-10-10

Change mesh format to triangle. This should fix test.

4087. By Christian Jacobs on 2012-10-10

The use of pylab.grid("True") in plot_data.py causes errors when a user tries to "make postprocess". Changing this to pylab.grid('on') seems to fix the problem. Thanks to Liping for reporting this problem.

4088. By Matthew Piggott on 2012-10-11

Weaken warn tolerances a little.

4089. By Matthew Piggott on 2012-10-11

Weaken one, but strengthen the other pass test tolerances on this - most quantities comfortably converging at second order but going from 64 squared to 128 squared the erturk_v error convergence rate is only 1.4

4090. By Lawrence Mitchell on 2012-10-15

Add --from-file option to testharness

`testharness --from-file=foo' tells testharness to read the list of
tests to run from `foo' (one per line). Useful for debugging failures
in a small number of related tests.

Unmerged revisions

4090. By Lawrence Mitchell on 2012-10-15

Add --from-file option to testharness

`testharness --from-file=foo' tells testharness to read the list of
tests to run from `foo' (one per line). Useful for debugging failures
in a small number of related tests.

4089. By Matthew Piggott on 2012-10-11

Weaken one, but strengthen the other pass test tolerances on this - most quantities comfortably converging at second order but going from 64 squared to 128 squared the erturk_v error convergence rate is only 1.4

4088. By Matthew Piggott on 2012-10-11

Weaken warn tolerances a little.

4087. By Christian Jacobs on 2012-10-10

The use of pylab.grid("True") in plot_data.py causes errors when a user tries to "make postprocess". Changing this to pylab.grid('on') seems to fix the problem. Thanks to Liping for reporting this problem.

4086. By Rhodri Davies on 2012-10-10

Change mesh format to triangle. This should fix test.

4085. By Rhodri Davies on 2012-10-09

Commit mesh files for debugging purposes only.

4084. By Matthew Piggott on 2012-10-07

Fix my diamond validation error.

4083. By Rhodri Davies on 2012-10-05

XML for stokes example. Should be good to go.

4082. By Michael Lange on 2012-10-04

Adding PETSs 3.3 support and dropping support for 3.0.

4081. By Rhodri Davies on 2012-10-03

Merge Full Projection fix into trunk. Thanks to Stephan for review.

=== modified file '.bzrignore'
--- .bzrignore	2011-04-15 16:18:30 +0000
+++ .bzrignore	2012-10-16 13:22:44 +0000
@@ -2,3 +2,4 @@
 backward_facing_step_2d-high.flml
 examples/backward_facing_step_2d/coarse
 coarse
+matrixdump

=== modified file 'AUTHORS'
--- AUTHORS	2011-10-21 17:23:47 +0000
+++ AUTHORS	2012-10-16 13:22:44 +0000
@@ -6,6 +6,7 @@
 ------------
 
 Imperial College London
+University of Edinburgh
 University of Oxford
 University of Reading
 National Oceanography Centre Southampton
@@ -48,6 +49,7 @@
 Frank Milthaler
 Julian Mindel
 Andrew Mitchell
+Lawrence Mitchell <[e-mail]>
 Rhodri Nelson
 Dimitrios Pavlidis
 Ralph Perpeet
@@ -60,6 +62,7 @@
 Brendan Tollit
 Axelle Vire <[e-mail]>
 Hongbin Wang
+Michele Weiland <[e-mail]>
 Martin Wells
 Matthew Whitworth
 Cian Wilson

=== modified file 'Makefile.in'
--- Makefile.in	2011-10-26 13:12:13 +0000
+++ Makefile.in	2012-10-16 13:22:44 +0000
@@ -27,7 +27,6 @@
 SHELL = @SHELL@
 
 PACKAGE_NAME = @PACKAGE_NAME@
-PACKAGE_VERSION = @PACKAGE_VERSION@
 
 FLUIDITY = $(PACKAGE_NAME)
 
@@ -63,9 +62,11 @@
 
 LINKER  = @LINKER@ @LDFLAGS@
 
-LIBS    = -L./lib -l$(FLUIDITY) @BLAS_LIBS@ @LIB_FEMDEM@ ./lib/libadaptivity.a \
+FLLINKER = @FLLINKER@ @LDFLAGS@ $(FCFLAGS)
+
+LIBS    = -L./lib -l$(FLUIDITY) @LIB_FEMDEM@ ./lib/libadaptivity.a \
            ./lib/libvtkfortran.a ./lib/libspatialindex.a ./lib/libspud.a ./lib/libjudy.a \
-           @LIBS@
+           @LIBS@ @BLAS_LIBS@
 
 OBJS = $(shell find . -mindepth 2 -maxdepth 2 -name '*.c' -o -name '*.cpp' -o -name '*.[fF]' -o -name '*.F90' | grep -v ^./lib | xargs)
 MAKE = @MAKE@ -s
@@ -76,8 +77,6 @@
 
 EVAL = tools/silenteval.sh
 
-VERSION = $(shell bzr revision-info)
-
 ifeq (@MBA2D@,yes)
 LIBMBA2D = lib/libmba2d.a
 endif
@@ -94,10 +93,18 @@
 EXCLUDE_TAGS =
 ifneq (@ZOLTAN@,yes)
 	EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e zoltan
+else
+	EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e nozoltan
 endif
 ifneq (@HAVE_ADJOINT@,yes)
 	EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e adjoint
 endif
+ifneq (@ALGENCAN@,yes)
+	EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e algencan
+endif
+ifneq (@MBA2D@,yes)
+	EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e 2dadapt
+endif
 
 .SUFFIXES: .f90 .F90 .c .cpp .o .a
 
@@ -110,7 +117,7 @@
 sub_system: $(LIBMBA2D) $(LIBMBA3D) $(LIBALGENCAN)  @LIB_FEMDEM@ \
 	lib/libadaptivity.a lib/libvtkfortran.a lib/libspatialindex.a \
 	lib/libspud.a lib/libjudy.a include/version.h
- 
+
 lib/libjudy.a:
 	@echo '    MKDIR lib'; mkdir -p lib
 	@cd libjudy; echo '    MAKE libjudy'; $(MAKE) > make.log ; $(MAKE) install >> make.log ; cd ..
@@ -144,7 +151,7 @@
 ifeq (@SPUD_ONLY@,yes)
 	@cd libspud;  echo '    MAKE libspud'; $(MAKE) ; $(MAKE) install-libspud
 else
-	@cd libspud;  echo '    MAKE libspud'; $(MAKE) ; $(MAKE) install; cd ..; cp lib*/python*/site-packages/libspud.so python/
+	@cd libspud;  echo '    MAKE libspud'; $(MAKE) ; $(MAKE) install-libspud; cd python; python setup.py build; cd ../..; cp libspud/python/build/lib*/libspud.so python/
 endif
 
 libfemtools: sub_system
@@ -160,19 +167,17 @@
 	@cd @FLFEMDEM_PATH@;$(MAKE)
 endif
 
-include/libspud.h: lib/libspud.a
-
-main.o: include/version.h include/libspud.h
-
-include/version.h:
-	@echo "DEFINE __FLUIDITY_VERSION__"
-	@grep "$(VERSION)" include/version.h >/dev/null 2>&1 || \
-	printf "#ifndef _FLUIDITY_VERSION_H\n\
-	#define _FLUIDITY_VERSION_H\n\
-	#define __FLUIDITY_VERSION__ \"$(VERSION)\"\n\
-	#endif /* _FLUIDITY_VERSION_H */\n" > include/version.h
-
-.PHONY: include/version.h
+include/spud.h: lib/libspud.a
+
+main.o: include/version.h include/spud.h
+
+include/version.h: version-info
+	@bin/version-info $@
+
+version-info:
+	@$(MAKE) -C tools install-version-info
+
+.PHONY: version-info
 
 bin/$(FLUIDITY): main.o lib/lib$(FLUIDITY).a
 	@echo "BUILD fluidity"
@@ -191,8 +196,8 @@
 	@echo "    MKDIR bin"
 	@mkdir -p bin
 	@echo "    LD shallow_water"
-	@$(EVAL) $(LINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS)
-	@$(EVAL) $(LINKER) -o bin/shallow_water main/Shallow_Water.o $(LIBS)
+	@$(EVAL) $(FLLINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS)
+	@$(EVAL) $(FLLINKER) -o bin/shallow_water main/Shallow_Water.o $(LIBS)
 
 bin/hybridized_helmholtz_solver: fluidity_library main/Hybridized_Helmholtz_Solver.F90
 	@cd main; $(MAKE) Hybridized_Helmholtz_Solver.o
@@ -200,7 +205,7 @@
 	@echo "    MKDIR bin"
 	@mkdir -p bin
 	@echo "    LD hybridized_helmholtz_solver"
-	@$(EVAL) $(LINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS)
+	@$(EVAL) $(FLLINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS)
 
 bin/burgers_equation: fluidity_library main/Burgers_Equation.F90
 	@cd main; $(MAKE) Burgers_Equation.o
@@ -396,6 +401,9 @@
 	@cd horizontal_adaptivity/tests;   $(MAKE) clean
 	@echo "    CLEAN tools"
 	@cd tools;    $(MAKE) clean
+	@echo "    CLEAN adjoint"
+	@cd adjoint;    $(MAKE) clean
+	@cd adjoint/tests;    $(MAKE) clean
 	@echo "    CLEAN bathymetry"
 	@cd bathymetry;  $(MAKE) clean
 	@echo "    CLEAN ocean_forcing"
@@ -415,9 +423,9 @@
 	@cd adjoint; $(MAKE) clean
 	@echo "    CLEAN main"
 	@cd main; $(MAKE) clean
-	@rm -fr bin/* lib
+	@rm -fr bin lib
 
-clean: clean-light clean-test python_clean
+clean: clean-light clean-test python_clean clean-debian
 	@echo "    CLEAN libvtkfortran"
 	@cd libvtkfortran; $(MAKE) clean
 	@echo "    CLEAN libjudy"
@@ -434,17 +442,24 @@
 	@cd libalgencan; $(MAKE) clean
 	@echo "    CLEAN libwm"
 	@cd libwm; $(MAKE) clean
+	@echo "    CLEAN manual"
+	@cd manual; $(MAKE) clean
 	@CLEAN_LIBADAPT@
-	@rm -fr include/*.mod *.cache core *.o config.status
+	@rm -fr lib*/configure.lineno include/*.mod *.cache core *.o config.status
 	@find ./ \( -name work.pc \) -exec rm {} \;
 	@find ./ \( -name work.pcl \) -exec rm {} \;
 
+clean-debian:
+	@echo "    CLEAN debian"
+	@cd debian; rm -rf files tmp fluidity python-fluidity *.substvars *.debhelper* 
+
 clean-test:      
 	@echo "    CLEAN tests"
-	@cd tests; ../tools/testharness.py --clean >/dev/null
+	@cd tests; PYTHONPATH=../python ../tools/testharness.py --clean >/dev/null
+	@cd tests/data; $(MAKE) clean
 
 distclean: clean
-	@echo "   DISTCLEAN"
+	@echo "    DISTCLEAN"
 	@cd spatialindex-1.5 ; make distclean > /dev/null ; rm -rf	\
 	 make.log regressiontest/mvrtree/.deps				\
 	 regressiontest/rtree/.deps regressiontest/tprtree/.deps >	\
@@ -454,6 +469,7 @@
 	@cd libvtkfortran ; make distclean > /dev/null
 	@cd libspud ; make distclean > /dev/null ; rm -rf		\
 	 libspud.so.0 libtool diamond/build python/build > /dev/null
+	@cd tools ; make distclean > /dev/null
 	@rm -fr `find ./ -name config.status` `find ./ -name		\
 	 config.log` share include/Judy.h include/Wm4Command.h		\
 	 include/Wm4FoundationLIB.h include/Wm4FoundationPCH.h		\
@@ -474,7 +490,7 @@
 	 include/spatialindex include/spud include/spud.h		\
 	 include/spud_enums.h include/tinystr.h include/tinyxml.h	\
 	 include/version.h include/vtk.h				\
-	 preprocessor/check_options.F90					\
+	 preprocessor/check_options.F90 schemas/flml			\
 	 preprocessor/register_diagnostics.F90 python/setup.py >	\
 	 /dev/null
 	@for i in `find ./*/ -name Makefile.in`; do rm -f `echo $$i | sed 's/.in$$//'`; done > /dev/null
@@ -483,30 +499,17 @@
 
 test: fltools bin/$(FLUIDITY) bin/shallow_water serialtest bin/burgers_equation
 
-serialtest:
-ifeq (@MBA2D@,yes)
-        ifeq (@ZOLTAN@,yes)
-	        @cd tests; ../bin/testharness -l short -p serial -e nozoltan -n $(THREADS)
-        else
-	        @cd tests; ../bin/testharness -l short -p serial -e zoltan -n $(THREADS)
-        endif
-else
-        ifeq (@ZOLTAN@,yes)
-	        @cd tests; ../bin/testharness -l short -p serial -e "2dadapt nozoltan" -n $(THREADS)
-        else
-	        @cd tests; ../bin/testharness -l short -p serial -e "2dadapt zoltan" -n $(THREADS)
-        endif
-endif
+serialtest: fltools bin/$(FLUIDITY) bin/shallow_water
+	@cd tests; ../bin/testharness -l short -p serial $(EXCLUDE_TAGS) -n $(THREADS)
 
-mediumtest: bin/burgers_equation
+mediumtest: fltools bin/$(FLUIDITY) manual bin/burgers_equation bin/shallow_water spudtools
 	@cd tests; ../bin/testharness -l medium $(EXCLUDE_TAGS) -n $(THREADS)
 
-mediumzoltantest:
-ifeq (@ZOLTAN@,yes)
-	@cd tests; ../bin/testharness -l medium -t zoltan -n $(THREADS)
-else
-	@echo "Warning: build did not have zoltan enabled, make mediumzoltantest will not run" -n $(THREADS)
-endif
+.PHONY: spudtools
+
+spudtools:
+	@cd libspud ; $(MAKE) install-spudtools
+	@echo "    INSTALL spudtools"
 
 setuputs: 
 	@echo "SETUP tests"
@@ -516,7 +519,6 @@
 	@mkdir -p bin/tests
 	@ln -sf ../../tests/data bin/tests/data
 
-
 build_unittest: fluidity_library libfemtools setuputs
 	@echo "    UNITTEST libspud"
 	@cd libspud/src/tests; $(MAKE) test-binaries; find bin/* -maxdepth 1 -type f -exec ln -sf ../../libspud/src/tests/{} ../../../bin/tests \; 
@@ -540,7 +542,11 @@
 	@echo "RUN bin/tests"
 	@bin/unittestharness bin/tests
 
-schema:
+bin/spud-preprocess:
+	@echo "Installing spudtools"
+	@$(MAKE) -C libspud install-spudtools
+
+schema: bin/spud-preprocess
 	@echo "Rebuilding schema fluidity_options.rng"
 	@./bin/spud-preprocess schemas/fluidity_options.rnc
 
@@ -600,7 +606,31 @@
 	$(MAKE) clean-light
 	@echo "    Congratulations, make makefiles succeeded!"
 
-install:
-	@mkdir -p $(DESTDIR)$(bindir)
-	cp bin/* $(DESTDIR)$(bindir)
-	cd python ; python setup.py install --root="${DESTDIR}" --prefix="$(prefix)"
+install: default fltools bin/shallow_water bin/burgers_equation
+	@mkdir -p $(DESTDIR)$(bindir) $(DESTDIR)$(docdir)/fluidity
+	find bin/ -maxdepth 1 -type f -exec cp '{}' $(DESTDIR)$(bindir) \;
+	rm -f $(DESTDIR)$(bindir)/spud-* $(DESTDIR)$(bindir)/diamond $(DESTDIR)$(bindir)/silenteval $(DESTDIR)$(bindir)/runut
+	[ -f manual/fluidity_manual.pdf ] && cp manual/fluidity_manual.pdf $(DESTDIR)$(docdir)/fluidity || true
+	mkdir -p $(DESTDIR)$(datadir)/fluidity
+	cp -r schemas/ $(DESTDIR)$(datadir)/fluidity/
+	mkdir -p $(DESTDIR)$(datadir)/diamond/schemata
+	cp schemas/flml $(DESTDIR)$(datadir)/diamond/schemata
+	cd python ; python setup.py install --root=$(shell echo ${DESTDIR} | sed 's/^$$/\//') --prefix="$(prefix)" $$FLUIDITY_PYTHON_INSTALL_ARGS
+	cp -r examples/ $(DESTDIR)$(docdir)/fluidity
+	find $(DESTDIR)$(docdir)/fluidity/examples -type f -exec sed -i "s/\.\.\/\.\.\/\.\.\/bin\///" '{}' \;
+	find $(DESTDIR)$(docdir)/fluidity/examples -type f -exec sed -i "s/\.\.\/\.\.\/bin\///" '{}' \;
+	cd $(DESTDIR)$(docdir)/fluidity ; tar -czvf $(DESTDIR)$(docdir)/fluidity/examples.tar.gz examples/
+	rm -rf $(DESTDIR)$(docdir)/fluidity/examples/
+
+install-diamond: 
+	cd libspud; ./configure --prefix=$(DESTDIR)@prefix@; cd ../..
+	cd libspud; $(MAKE) clean; cd ../..
+	cd libspud; $(MAKE) install-diamond; cd ../..
+	cd libspud; $(MAKE) install-dxdiff; cd ../..
+	cd libspud; $(MAKE) install-spudtools; cd ../..
+	cd libspud; ./configure --prefix=$(PWD); cd ../..
+
+install-user-schemata:
+	mkdir -p $(HOME)/.diamond/schemata/
+	echo "Fluidity Markup Language" > $(HOME)/.diamond/schemata/flml
+	echo "$(PWD)/schemas/fluidity_options.rng" >> $(HOME)/.diamond/schemata/flml

=== modified file 'aclocal.m4'
--- aclocal.m4	2011-10-03 17:24:46 +0000
+++ aclocal.m4	2012-10-16 13:22:44 +0000
@@ -319,9 +319,6 @@
 	acx_lapack_ok=noblas
 fi
 
-# add BLAS to libs
-LIBS="$BLAS_LIBS $LIBS"
-
 # Are we linking from C?
 case "$ac_ext" in
   f*|F*) dsyev="dsyev" ;;
@@ -415,10 +412,8 @@
   AC_MSG_ERROR( [You need to set PETSC_DIR to point at your PETSc installation... exiting] )
 fi
 
-
-
 PETSC_LINK_LIBS=`make -s -f petsc_makefile getlinklibs`
-LIBS="$LIBS $PETSC_LINK_LIBS"
+LIBS="$PETSC_LINK_LIBS $LIBS"
 
 PETSC_INCLUDE_FLAGS=`make -s -f petsc_makefile getincludedirs`
 CPPFLAGS="$CPPFLAGS $PETSC_INCLUDE_FLAGS"
@@ -470,37 +465,26 @@
 # petsc tutorial - using the headers in the same way as we do in the code
 AC_LINK_IFELSE(
 [AC_LANG_SOURCE([
-program test
+program test_petsc
 #include "petscversion.h"
 #ifdef HAVE_PETSC_MODULES
   use petsc
-#if PETSC_VERSION_MINOR==0
-  use petscvec 
-  use petscmat 
-  use petscksp 
-  use petscpc
-#endif
 #endif
 implicit none
 #ifdef HAVE_PETSC_MODULES
-#include "finclude/petscvecdef.h"
-#include "finclude/petscmatdef.h"
-#include "finclude/petsckspdef.h"
-#include "finclude/petscpcdef.h"
+#include "finclude/petscdef.h"
 #else
 #include "finclude/petsc.h"
-#if PETSC_VERSION_MINOR==0
-#include "finclude/petscvec.h"
-#include "finclude/petscmat.h"
-#include "finclude/petscksp.h"
-#include "finclude/petscpc.h"
-#endif
 #endif
       double precision  norm
       PetscInt  i,j,II,JJ,m,n,its
       PetscInt  Istart,Iend,ione
       PetscErrorCode ierr
+#if PETSC_VERSION_MINOR>=2
+      PetscBool flg
+#else
       PetscTruth  flg
+#endif
       PetscScalar v,one,neg_one
       Vec         x,b,u
       Mat         A 
@@ -573,25 +557,26 @@
       call MatDestroy(A,ierr)
 
       call PetscFinalize(ierr)
-end program test
+end program test_petsc
 ])],
 [
 AC_MSG_NOTICE([PETSc program succesfully compiled and linked.])
 ],
 [
-cp conftest.F90 test.F90
+cp conftest.F90 test_petsc.F90
 AC_MSG_FAILURE([Failed to compile and link PETSc program.])])
 
 AC_LANG_RESTORE
 
 # finally check we have the right petsc version
 AC_COMPUTE_INT(PETSC_VERSION_MAJOR, "PETSC_VERSION_MAJOR", [#include "petscversion.h"], 
-  [AC_MSG_ERROR([Unknown petsc version])])
-if test "x$PETSC_VERSION_MAJOR" == "x3" ; then
-  AC_MSG_NOTICE([Detected PETSc version 3])
-else
-  AC_MSG_NOTICE([Detected PETSc version "$PETSC_VERSION_MAJOR"])
-  AC_MSG_ERROR([Fluidity needs PETSc version 3])
+  [AC_MSG_ERROR([Unknown petsc major version])])
+AC_COMPUTE_INT(PETSC_VERSION_MINOR, "PETSC_VERSION_MINOR", [#include "petscversion.h"], 
+  [AC_MSG_ERROR([Unknown petsc minor version])])
+AC_MSG_NOTICE([Detected PETSc version "$PETSC_VERSION_MAJOR"."$PETSC_VERSION_MINOR"])
+# if major<3 or minor<1
+if test "0$PETSC_VERSION_MAJOR" -lt 3 -o "0$PETSC_VERSION_MINOR" -lt 1; then
+  AC_MSG_ERROR([Fluidity needs PETSc version >=3.1])
 fi
 
 AC_DEFINE(HAVE_PETSC,1,[Define if you have the PETSc library.])
@@ -617,7 +602,7 @@
 tmpCPPFLAGS=$CPPFLAGS
 AC_LANG_SAVE
 AC_LANG_C
-LIBS="$tmpLIBS -L$ParMetis_LIBS_PATH -lparmetis -lmetis -lm"
+LIBS="$tmpLIBS -L$ParMetis_LIBS_PATH -lparmetis -lmetis -lm $ZOLTAN_DEPS"
 AC_CHECK_LIB(
 	[parmetis],
 	[ParMETIS_V3_AdaptiveRepart],
@@ -1005,7 +990,6 @@
 	echo "note: using $zoltan_INCLUDES_PATH/zoltan.mod"
 fi 
 
-
 AC_LANG_SAVE
 AC_LANG_C
 AC_CHECK_LIB(
@@ -1013,7 +997,28 @@
 	[Zoltan_Initialize],
 	[AC_DEFINE(HAVE_ZOLTAN,1,[Define if you have zoltan library.])],
 	[AC_MSG_ERROR( [Could not link in the zoltan library... exiting] )] )
-# Save variables...
+
+# Small test for zoltan .mod files:
+AC_LANG(Fortran)
+ac_ext=F90
+# In fluidity's makefile we explicitly add CPPFLAGS, temporarily add it to
+# FCFLAGS here for this zoltan test:
+tmpFCFLAGS="$FCFLAGS"
+FCFLAGS="$FCFLAGS $CPPFLAGS"
+AC_LINK_IFELSE(
+[AC_LANG_SOURCE([
+program test_zoltan
+ use zoltan
+end program test_zoltan
+])],
+[
+AC_MSG_NOTICE([Great success! Zoltan .mod files exist and are usable])
+],
+[
+cp conftest.F90 test_zoltan.F90
+AC_MSG_FAILURE([Failed to find zoltan.mod files])])
+# And now revert FCFLAGS
+FCFLAGS="$tmpFCFLAGS"
 AC_LANG_RESTORE
 
 ZOLTAN="yes"
@@ -1051,7 +1056,7 @@
 AC_LANG_C
 AC_CHECK_LIB(
 	[adjoint],
-	[adj_register_equation],
+	[adj_register_forward_source_callback],
 	[AC_DEFINE(HAVE_ADJOINT,1,[Define if you have libadjoint.])HAVE_ADJOINT=yes],
 	[AC_MSG_WARN( [Could not link in libadjoint ... ] );HAVE_ADJOINT=no;LIBS=$bakLIBS] )
 # Save variables...

=== modified file 'assemble/Adapt_State.F90'
--- assemble/Adapt_State.F90	2011-09-28 13:38:45 +0000
+++ assemble/Adapt_State.F90	2012-10-16 13:22:44 +0000
@@ -73,13 +73,12 @@
   use diagnostic_variables
   use diagnostic_fields_wrapper_new, only : calculate_diagnostic_variables_new => calculate_diagnostic_variables
   use pickers
+  use write_gmsh
 #ifdef HAVE_ZOLTAN
   use zoltan_integration
   use mpi_interfaces
 #endif
 
-  use mesh_files
-
   implicit none
 
   private
@@ -923,10 +922,10 @@
       output_positions => extract_vector_field(states(1), "Coordinate")
 
       if(isparallel()) then
-        call write_mesh_files(parallel_filename("first_timestep_adapted_mesh"), output_positions)
+        call write_gmsh_file(parallel_filename("first_timestep_adapted_mesh"), output_positions)
         call write_halos("first_timestep_adapted_mesh", output_positions%mesh)
       else
-        call write_mesh_files("first_timestep_adapted_mesh", output_positions)
+        call write_gmsh_file("first_timestep_adapted_mesh", output_positions)
       end if
 
     end if
@@ -1003,13 +1002,9 @@
     call get_option("/mesh_adaptivity/hr_adaptivity/zoltan_options/element_quality_cutoff", &
        & quality_tolerance, default = 0.6)
 
-    if (max_adapt_iteration .ne. 1) then
-       zoltan_min_adapt_iterations = adapt_iterations()
-       zoltan_max_adapt_iterations = zoltan_min_adapt_iterations + zoltan_additional_adapt_iterations
-    else
-       zoltan_min_adapt_iterations = adapt_iterations()
-       zoltan_max_adapt_iterations = zoltan_min_adapt_iterations
-    end if
+    zoltan_min_adapt_iterations = adapt_iterations()
+    zoltan_max_adapt_iterations = zoltan_min_adapt_iterations + zoltan_additional_adapt_iterations
+
 #endif
 
     finished_adapting = .false.
@@ -1238,7 +1233,11 @@
           ! call zoltan now but we need to pass in both the 2d metric (metric) and the 3d full metric (full_metric)
           ! the first is needed to define the element qualities while the second must be interpolated to the newly
           ! decomposed mesh
-          call zoltan_drive(states, final_adapt_iteration, global_min_quality, metric = metric, full_metric = full_metric)
+          if (zoltan_additional_adapt_iterations .gt. 0) then
+             call zoltan_drive(states, final_adapt_iteration, global_min_quality = global_min_quality, metric = metric, full_metric = full_metric)
+          else
+             call zoltan_drive(states, final_adapt_iteration, metric = metric, full_metric = full_metric)
+          end if
           default_stat%zoltan_drive_call=.true.
 
           ! now we can deallocate the horizontal metric and point metric back at the full metric again
@@ -1246,7 +1245,11 @@
           metric = full_metric
         else
 
-          call zoltan_drive(states, final_adapt_iteration, global_min_quality, metric = metric)
+          if (zoltan_additional_adapt_iterations .gt. 0) then
+             call zoltan_drive(states, final_adapt_iteration, global_min_quality = global_min_quality, metric = metric)
+          else
+             call zoltan_drive(states, final_adapt_iteration, metric = metric)
+          end if
           default_stat%zoltan_drive_call=.true.
 
         end if
@@ -1303,10 +1306,12 @@
                ewrite(2,*) "The next iteration will be final adapt iteration else we'll go over the maximum adapt iterations."
             end if
 
-            if (global_min_quality .le. quality_tolerance) then
-               ewrite(-1,*) "Mesh contains elements with quality below element quality tolerance. May need to increase number of adapt iterations to ensure good quality mesh."
-               ewrite(-1,*) "min_quality = ", global_min_quality
-               ewrite(-1,*) "quality_tolerance = ", quality_tolerance
+            if (zoltan_additional_adapt_iterations .gt. 0) then
+               if (global_min_quality .le. quality_tolerance) then
+                  ewrite(-1,*) "Mesh contains elements with quality below element quality tolerance. May need to increase number of adapt iterations to ensure good quality mesh."
+                  ewrite(-1,*) "min_quality = ", global_min_quality
+                  ewrite(-1,*) "quality_tolerance = ", quality_tolerance
+               end if
             end if
             
             final_adapt_iteration = .true.
@@ -1607,7 +1612,7 @@
       file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no)
       call find_mesh_to_adapt(states(1), mesh)
       positions = get_coordinate_field(states(1), mesh)
-      call write_mesh_files(file_name, positions)
+      call write_gmsh_file(file_name, positions)
       if(isparallel()) then
         file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no, add_parallel = .false.)  ! parallel extension is added by write_halos
         call write_halos(file_name, positions%mesh)

=== modified file 'assemble/Adapt_State_Prescribed.F90'
--- assemble/Adapt_State_Prescribed.F90	2011-07-22 15:34:14 +0000
+++ assemble/Adapt_State_Prescribed.F90	2012-10-16 13:22:44 +0000
@@ -109,8 +109,8 @@
     if(have_option(base_path // "/mesh/from_file")) then
       call get_option(base_path // "/mesh/from_file/format/name", format)
       call get_option("/geometry/quadrature/degree", quad_degree)
-          
-      new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree)
+
+      new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree, format = format)
     else
       ewrite(2, *) "Extracting new mesh from state"
     

=== modified file 'assemble/Adaptivity_1D.F90'
--- assemble/Adaptivity_1D.F90	2010-11-09 14:08:19 +0000
+++ assemble/Adaptivity_1D.F90	2012-10-16 13:22:44 +0000
@@ -95,7 +95,12 @@
 
     call deallocate(sizing)
     
-    assert(descending_coordinate_ordered(new_positions))
+    if (.not. descending_coordinate_ordered(new_positions)) then
+      ewrite(-1,*) "To use 1D adaptivity you need an input mesh for which the ordering of the nodes is such that " // &
+        "their coordinates decrease with increasing node number. If you use the 'interval' script to produce a " // &
+        "mesh, you can achieve this by adding the '--reverse' option."
+      FLExit("To be adapted 1D mesh not in descending order.")
+    end if
     
     ! adapt_1d doesn't build a complete mesh. Build the rest of the mesh.
     assert(ele_count(new_positions) == node_count(new_positions) - 1)

=== modified file 'assemble/Advection_Diffusion_CG.F90'
--- assemble/Advection_Diffusion_CG.F90	2011-09-15 13:37:20 +0000
+++ assemble/Advection_Diffusion_CG.F90	2012-10-16 13:22:44 +0000
@@ -39,13 +39,20 @@
   use boundary_conditions_from_options
   use field_options
   use fldebug
-  use global_parameters, only : FIELD_NAME_LEN, OPTION_PATH_LEN
+  use global_parameters, only : FIELD_NAME_LEN, OPTION_PATH_LEN, COLOURING_CG1
   use profiler
   use spud
   use petsc_solve_state_module
   use state_module
   use upwind_stabilisation
   use sparsity_patterns_meshes
+  use porous_media
+  use multiphase_module
+  use sparse_tools_petsc
+  use colouring
+#ifdef _OPENMP
+  use omp_lib
+#endif
   
   implicit none
   
@@ -99,6 +106,8 @@
   logical :: integrate_advection_by_parts
   ! Source?
   logical :: have_source
+  ! Add source directly to the right hand side?
+  logical :: add_src_directly_to_rhs
   ! Absorption?
   logical :: have_absorption
   ! Diffusivity?
@@ -107,20 +116,23 @@
   logical :: isotropic_diffusivity
   ! Is the mesh moving?
   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, &
-                                     extra_discretised_source, 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
-    type(scalar_field), intent(in), optional :: extra_discretised_source
     integer, intent(out), optional :: iterations_taken
     
     type(csr_matrix) :: matrix
@@ -130,17 +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, &
-                                         extra_discretised_source = extra_discretised_source)    
+    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")
 
@@ -206,31 +223,43 @@
     
   end subroutine set_advection_diffusion_cg_initial_guess
   
-  subroutine assemble_advection_diffusion_cg(t, matrix, rhs, state, dt, velocity_name, &
-                                             extra_discretised_source)
+  subroutine assemble_advection_diffusion_cg(t, matrix, rhs, state, dt, velocity_name)
     type(scalar_field), intent(inout) :: t
     type(csr_matrix), intent(inout) :: matrix
     type(scalar_field), intent(inout) :: rhs
-    type(state_type), intent(in) :: state
+    type(state_type), intent(inout) :: state
     real, intent(in) :: dt
     character(len = *), optional, intent(in) :: velocity_name
-    type(scalar_field), intent(in), optional :: extra_discretised_source
 
-    character(len = FIELD_NAME_LEN) :: lvelocity_name
+    character(len = FIELD_NAME_LEN) :: lvelocity_name, velocity_equation_type
     integer :: i, j, stat
     integer, dimension(:), allocatable :: t_bc_types
     type(scalar_field) :: t_bc, t_bc_2
     type(scalar_field), pointer :: absorption, sinking_velocity, source
     type(tensor_field), pointer :: diffusivity
     type(vector_field) :: velocity
-    type(vector_field), pointer :: gravity_direction, velocity_ptr, grid_velocity
+    type(vector_field), pointer :: gravity_direction, temp_velocity_ptr, velocity_ptr, grid_velocity
     type(vector_field), pointer :: positions, old_positions, new_positions
     type(scalar_field), target :: dummydensity
     type(scalar_field), pointer :: density, olddensity
     character(len = FIELD_NAME_LEN) :: density_name
     type(scalar_field), pointer :: pressure
-      
-    type(element_type) :: supg_element
+    
+    ! 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
+        
+    !! Coloring  data structures for OpenMP parallization
+    type(integer_set), dimension(:), pointer :: colours
+    integer :: clr, nnid, len, ele
+    integer :: num_threads, thread_num
+    !! Did we successfully prepopulate the transform_to_physical_cache?
+    logical :: cache_valid
+
+    type(element_type), dimension(:), allocatable :: supg_element
   
     ewrite(1, *) "In assemble_advection_diffusion_cg"
     
@@ -242,7 +271,12 @@
     else
       lvelocity_name = "NonlinearVelocity"
     end if
-    
+
+#ifdef _OPENMP
+    num_threads = omp_get_max_threads()
+#else
+    num_threads = 1
+#endif
     ! Step 1: Pull fields out of state
     
     ! Coordinate
@@ -281,10 +315,19 @@
     if(have_source) then
       assert(mesh_dim(source) == mesh_dim(t))
       assert(ele_count(source) == ele_count(t))
-    
+      
+      add_src_directly_to_rhs = have_option(trim(source%option_path)//'/diagnostic/add_directly_to_rhs')
+      
+      if (add_src_directly_to_rhs) then 
+         ewrite(2, *) "Adding Source field directly to the right hand side"
+         assert(node_count(source) == node_count(t))
+      end if
+      
       ewrite_minmax(source)
     else
       ewrite(2, *) "No source"
+      
+      add_src_directly_to_rhs = .false.
     end if
     
     ! Absorption
@@ -333,6 +376,19 @@
       isotropic_diffusivity = .false.
       ewrite(2, *) "No diffusivity"
     end if
+
+    ! Porosity
+    if (have_option(trim(complete_field_path(t%option_path))//'/porosity')) then
+       include_porosity = .true.
+       
+       ! get the porosity theta averaged field - this will allocate it
+       call form_porosity_theta(porosity_theta, state, option_path = trim(complete_field_path(t%option_path))//'/porosity')       
+    else
+       include_porosity = .false.
+       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
     
@@ -375,6 +431,9 @@
     ! are we moving the mesh?
     move_mesh = (have_option("/mesh_adaptivity/mesh_movement") .and. have_mass)
     if(move_mesh) then
+      if (include_porosity) then
+         FLExit('Cannot include porosity in CG advection diffusion of a field with a moving mesh')
+      end if
       ewrite(2,*) "Moving the mesh"
       old_positions => extract_vector_field(state, "OldCoordinate")
       ewrite_minmax(old_positions)
@@ -392,6 +451,7 @@
       ewrite(2,*) "Not moving the mesh"
     end if
     
+    allocate(supg_element(num_threads))
     if(have_option(trim(t%option_path) // "/prognostic/spatial_discretisation/continuous_galerkin/stabilisation/streamline_upwind")) then
       ewrite(2, *) "Streamline upwind stabilisation"
       stabilisation_scheme = STABILISATION_STREAMLINE_UPWIND
@@ -406,8 +466,10 @@
       call get_upwind_options(trim(t%option_path) // "/prognostic/spatial_discretisation/continuous_galerkin/stabilisation/streamline_upwind_petrov_galerkin", &
           & nu_bar_scheme, nu_bar_scale)
       ! Note this is not mixed mesh safe (but then nothing really is)
-      ! You actually need 1 supg_element per thread.
-      supg_element=make_supg_element(ele_shape(t,1)) 
+      ! You need 1 supg_element per thread.
+      do i = 1, num_threads
+         supg_element(i)=make_supg_element(ele_shape(t,1))
+      end do
       if(move_mesh) then
         FLExit("Haven't thought about how mesh movement works with stabilisation yet.")
       end if
@@ -416,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
@@ -428,24 +505,62 @@
       olddensity => dummydensity
       density_theta = 1.0
       pressure => dummydensity
+
     case(FIELD_EQUATION_INTERNALENERGY)
       ewrite(2,*) "Solving internal energy equation"
       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(FIELD_EQUATION_KEPSILON)
+      ewrite(2,*) "Solving k-epsilon equation"
+      if(move_mesh) then
+        FLExit("Haven't implemented a moving mesh k-epsilon equation yet.")
+      end if
+      
+      ! Depending on the equation type, extract the density or set it to some dummy field allocated above
+      temp_velocity_ptr => extract_vector_field(state, "Velocity")
+      call get_option(trim(temp_velocity_ptr%option_path)//"/prognostic/equation[0]/name", velocity_equation_type)
+      select case(velocity_equation_type)
+         case("LinearMomentum")
+            density=>extract_scalar_field(state, "Density")
+            olddensity => dummydensity
+            density_theta = 1.0
+         case("Boussinesq")
+            density=>dummydensity
+            olddensity => dummydensity
+            density_theta = 1.0
+         case("Drainage")
+            density=>dummydensity
+            olddensity => dummydensity
+            density_theta = 1.0
+         case default
+            ! developer error... out of sync options input and code
+            FLAbort("Unknown equation type for velocity")
+      end select
+      ewrite_minmax(density)
+
     case default
       FLExit("Unknown field equation type for cg advection diffusion.")
     end select
@@ -455,19 +570,51 @@
     call zero(matrix)
     call zero(rhs)
     
-    do i = 1, ele_count(t)
-      call assemble_advection_diffusion_element_cg(i, t, matrix, rhs, &
-                                        positions, old_positions, new_positions, &
-                                        velocity, grid_velocity, &
-                                        source, absorption, diffusivity, &
-                                        density, olddensity, pressure,&
-                                        supg_element)
-    end do
-
-    ! as part of assembly include the already discretised optional source
-    ! needed before applying direchlet boundary conditions
-    call addto_rhs_extra_discretised_source(rhs, extra_discretised_source = extra_discretised_source)
-
+    call profiler_tic(t, "advection_diffusion_loop_overhead")
+
+#ifdef _OPENMP
+    cache_valid = prepopulate_transform_cache(positions)
+    assert(cache_valid)
+#endif
+
+    call get_mesh_colouring(state, t%mesh, COLOURING_CG1, colours)
+    call profiler_toc(t, "advection_diffusion_loop_overhead")
+
+    call profiler_tic(t, "advection_diffusion_loop")
+
+    !$OMP PARALLEL DEFAULT(SHARED) &
+    !$OMP PRIVATE(clr, len, nnid, ele, thread_num)
+
+#ifdef _OPENMP    
+    thread_num = omp_get_thread_num()
+#else
+    thread_num=0
+#endif
+    colour_loop: do clr = 1, size(colours)
+
+      len = key_count(colours(clr))
+      !$OMP DO SCHEDULE(STATIC)
+      element_loop: do nnid = 1, len
+         ele = fetch(colours(clr), nnid)
+         call assemble_advection_diffusion_element_cg(ele, t, matrix, rhs, &
+              positions, old_positions, new_positions, &
+              velocity, grid_velocity, &
+              source, absorption, diffusivity, &
+              density, olddensity, pressure, porosity_theta, nvfrac, &
+              supg_element(thread_num+1))
+      end do element_loop
+      !$OMP END DO
+
+    end do colour_loop
+    !$OMP END PARALLEL
+
+    call profiler_toc(t, "advection_diffusion_loop")
+
+    ! Add the source directly to the rhs if required 
+    ! which must be included before dirichlet BC's.
+    if (add_src_directly_to_rhs) call addto(rhs, source)
+    
+    
     ! Step 4: Boundary conditions
     
     if( &
@@ -494,7 +641,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)
@@ -509,10 +656,17 @@
     ewrite_minmax(rhs)
     
     call deallocate(velocity)
+    call deallocate(nvfrac)
     call deallocate(dummydensity)
-    if (stabilisation_scheme == STABILISATION_SUPG) &
-         call deallocate(supg_element)
+    if (stabilisation_scheme == STABILISATION_SUPG) then
+       do i = 1, num_threads
+          call deallocate(supg_element(i))
+       end do
+    end if
+    deallocate(supg_element)
 
+    call deallocate(porosity_theta)
+    
     ewrite(1, *) "Exiting assemble_advection_diffusion_cg"
     
   end subroutine assemble_advection_diffusion_cg
@@ -562,7 +716,7 @@
                                       positions, old_positions, new_positions, &
                                       velocity, grid_velocity, &
                                       source, absorption, diffusivity, &
-                                      density, olddensity, pressure, 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
@@ -577,6 +731,8 @@
     type(scalar_field), intent(in) :: density
     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
@@ -585,6 +741,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
@@ -611,6 +770,9 @@
       assert(ele_loc(grid_velocity, ele) == ele_loc(positions, ele))
       assert(ele_ngi(grid_velocity, ele) == ele_ngi(velocity, ele))
     end if
+    if (include_porosity) then
+      assert(ele_ngi(porosity_theta, ele) == ele_ngi(t, ele))    
+    end if
 #endif
 
     matrix_addto = 0.0
@@ -630,7 +792,7 @@
         & dshape = dt_t, detwei = detwei)
     end if
     
-    if(have_advection.or.(equation_type==FIELD_EQUATION_INTERNALENERGY)) then
+    if(have_advection.or.(equation_type==FIELD_EQUATION_INTERNALENERGY).or.equation_type==FIELD_EQUATION_KEPSILON) then
       call transform_to_physical(positions, ele, &
            & ele_shape(velocity, ele), dshape = du_t)
     end if
@@ -645,7 +807,7 @@
       call transform_to_physical(new_positions, ele, detwei=detwei_new)
     end if
     
-    if(have_advection.and.(equation_type==FIELD_EQUATION_INTERNALENERGY)) then
+    if(have_advection.and.(equation_type==FIELD_EQUATION_INTERNALENERGY .or. equation_type==FIELD_EQUATION_KEPSILON)) then
       if(ele_shape(density, ele)==t_shape) then
         drho_t = dt_t
       else
@@ -654,6 +816,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)
@@ -661,7 +833,6 @@
         if(have_diffusivity) then
           call supg_test_function(supg_shape, t_shape, dt_t, ele_val_at_quad(velocity, ele), j_mat, diff_q = ele_val_at_quad(diffusivity, ele), &
             & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
-          test_function = supg_shape
         else
           call supg_test_function(supg_shape, t_shape, dt_t, ele_val_at_quad(velocity, ele), j_mat, &
             & nu_bar_scheme = nu_bar_scheme, nu_bar_scale = nu_bar_scale)
@@ -677,13 +848,13 @@
     ! Step 3: Assemble contributions
     
     ! Mass
-    if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, 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)
@@ -692,12 +863,14 @@
     if(have_diffusivity) call add_diffusivity_element_cg(ele, t, diffusivity, dt_t, detwei, matrix_addto, rhs_addto)
     
     ! Source
-    if(have_source) call add_source_element_cg(ele, test_function, t, source, detwei, rhs_addto)
+    if(have_source .and. (.not. add_src_directly_to_rhs)) then 
+       call add_source_element_cg(ele, test_function, t, source, detwei, rhs_addto)
+    end if
     
     ! 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
             
@@ -707,11 +880,13 @@
 
   end subroutine assemble_advection_diffusion_element_cg
   
-  subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, 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
@@ -720,9 +895,12 @@
     real, dimension(ele_loc(t, ele), ele_loc(t, ele)) :: mass_matrix
     
     real, dimension(ele_ngi(density,ele)) :: density_at_quad
+    real, dimension(ele_ngi(porosity_theta,ele)) :: porosity_theta_at_quad
     
     assert(have_mass)
     
+    if (include_porosity) porosity_theta_at_quad = ele_val_at_quad(porosity_theta, ele)
+    
     select case(equation_type)
     case(FIELD_EQUATION_INTERNALENERGY)
       assert(ele_ngi(density, ele)==ele_ngi(olddensity, ele))
@@ -733,15 +911,28 @@
         ! needs to be evaluated at t+dt
         mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei_new*density_at_quad)
       else
-        mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad)
+        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
       end if
+    case(FIELD_EQUATION_KEPSILON)      
+      density_at_quad = ele_val_at_quad(density, ele)
+      mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad)
     case default
     
       if(move_mesh) then
         ! needs to be evaluated at t+dt
         mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei_new)
       else
-        mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei)
+        if (include_porosity) then
+          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*porosity_theta_at_quad)
+        else 
+          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei)
+        end if
       end if
       
     end select
@@ -777,8 +968,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
@@ -786,10 +977,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
@@ -804,6 +997,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)
@@ -822,6 +1019,15 @@
       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
+
+    case(FIELD_EQUATION_KEPSILON)
+      density_at_quad = ele_val_at_quad(density, ele)
+      densitygrad_at_quad = ele_grad_at_quad(density, ele, drho_t)
+      udotgradrho_at_quad = sum(densitygrad_at_quad*velocity_at_quad, 1)
     end select
                 
     if(integrate_advection_by_parts) then
@@ -831,6 +1037,28 @@
       !    /                                        /
       select case(equation_type)
       case(FIELD_EQUATION_INTERNALENERGY)
+        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)
+          
+          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(FIELD_EQUATION_KEPSILON)
         advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad)
         if(abs(1.0 - beta) > epsilon(0.0)) then
           velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
@@ -853,6 +1081,33 @@
       !  /                                 /
       select case(equation_type)
       case(FIELD_EQUATION_INTERNALENERGY)
+        
+        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)
+          
+          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(FIELD_EQUATION_KEPSILON)
         advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad)
         if(abs(beta) > epsilon(0.0)) then
           velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
@@ -960,26 +1215,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
@@ -992,6 +1260,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
@@ -1025,7 +1294,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)
@@ -1038,7 +1307,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
@@ -1047,6 +1316,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
@@ -1071,8 +1341,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))
@@ -1125,22 +1399,6 @@
     end if
 
   end subroutine add_diffusivity_face_cg
-
-  subroutine addto_rhs_extra_discretised_source(rhs, extra_discretised_source)
-     type(scalar_field), intent(inout) :: rhs
-     type(scalar_field), intent(in),optional :: extra_discretised_source
-     
-     ! include the already discretised source into rhs
-     add_extra_source: if (present(extra_discretised_source)) then 
-         
-        ! assert that the rhs and extra_discretised_source have the same mesh
-        assert(trim(rhs%mesh%name) == trim(extra_discretised_source%mesh%name))
-         
-        call addto(rhs, extra_discretised_source)
-      
-     end if add_extra_source
-  
-  end subroutine  addto_rhs_extra_discretised_source
      
   subroutine solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state, iterations_taken)
     type(scalar_field), intent(in) :: t

=== modified file 'assemble/Advection_Diffusion_DG.F90'
--- assemble/Advection_Diffusion_DG.F90	2011-10-19 10:46:35 +0000
+++ assemble/Advection_Diffusion_DG.F90	2012-10-16 13:22:44 +0000
@@ -52,7 +52,14 @@
   use sparse_matrices_fields
   use sparsity_patterns_meshes
   use diagnostic_fields, only: calculate_diagnostic_variable
-  use global_parameters, only : FIELD_NAME_LEN
+  use global_parameters, only: OPTION_PATH_LEN, FIELD_NAME_LEN, COLOURING_DG2, &
+       COLOURING_DG0
+  use porous_media
+  use colouring
+  use Profiler
+#ifdef _OPENMP
+  use omp_lib
+#endif
 
   implicit none
 
@@ -96,12 +103,15 @@
   ! Boundary condition types:
   ! (the numbers should match up with the order in the 
   !  get_entire_boundary_condition call)
-  integer :: BCTYPE_WEAKDIRICHLET=1, BCTYPE_DIRICHLET=2
+  integer :: BCTYPE_WEAKDIRICHLET=1, BCTYPE_DIRICHLET=2, BCTYPE_NEUMANN=3
 
   logical :: include_mass
   ! are we moving the mesh?
   logical :: move_mesh
 
+  ! Include porosity?
+  logical :: include_porosity
+
   ! Stabilisation schemes.
   integer :: stabilisation_scheme
   integer, parameter :: NONE=0
@@ -133,6 +143,7 @@
   logical :: on_sphere
   ! Vertical diffusion by mixing option
   logical :: have_buoyancy_adjustment_by_vertical_diffusion
+  logical :: have_buoyancy_adjustment_diffusivity
 
 contains
 
@@ -483,6 +494,9 @@
 
     character(len=FIELD_NAME_LEN) :: limiter_name
     integer :: i
+    
+    !! Courant number field name used for temporal subcycling
+    character(len=FIELD_NAME_LEN) :: Courant_number_name
 
     T=>extract_scalar_field(state, field_name)
     T_old=>extract_scalar_field(state, "Old"//field_name)
@@ -521,10 +535,12 @@
     call allocate(rhs_diff, T%mesh, trim(field_name)//" Diffusion RHS")
    
     call construct_advection_diffusion_dg(matrix, rhs, field_name, state, &
-         mass, matrix_diff, rhs_diff, semidiscrete=.true., &
+         mass=mass, diffusion_m=matrix_diff, diffusion_rhs=rhs_diff, semidiscrete=.true., &
          velocity_name=velocity_name)
 
-    call get_dg_inverse_mass_matrix(inv_mass, mass)
+    ! mass has only been assembled only for owned elements, so we can only compute
+    ! its inverse for owned elements
+    call get_dg_inverse_mass_matrix(inv_mass, mass, only_owned_elements=.true.)
     
     ! Note that since theta and dt are module global, these lines have to
     ! come after construct_advection_diffusion_dg.
@@ -543,9 +559,17 @@
             &"/prognostic/temporal_discretisation/discontinuous_galerkin"//&
             &"/maximum_courant_number_per_subcycle", Max_Courant_number)
        
-       s_field => extract_scalar_field(state, "DG_CourantNumber")
-       call calculate_diagnostic_variable(state, "DG_CourantNumber", &
-            & s_field)
+       ! Determine the courant field to use to find the max
+       call get_option(trim(T%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin"//&
+            &"/maximum_courant_number_per_subcycle/courant_number/name", &
+            &Courant_number_name, default="DG_CourantNumber")
+       
+       s_field => extract_scalar_field(state, trim(Courant_number_name))
+       call calculate_diagnostic_variable(state, trim(Courant_number_name), &
+            & s_field, option_path=trim(T%option_path)//&
+            &"/prognostic/temporal_discretisation/discontinuous_galerkin"//&
+            &"/courant_number")
        
        subcycles = ceiling( maxval(s_field%val)/Max_Courant_number)
        call allmax(subcycles)
@@ -714,6 +738,22 @@
     integer, dimension(:), allocatable :: bc_type
 
     type(mesh_type), pointer :: mesh_cg
+    
+    ! Porosity field
+    type(scalar_field) :: porosity_theta
+    
+    !! Add the Source directly to the right hand side?
+    logical :: add_src_directly_to_rhs
+
+
+    type(integer_set), dimension(:), pointer :: colours
+    integer :: len, clr, nnid
+    !! Is the transform_to_physical cache we prepopulated valid
+    logical :: cache_valid
+    integer :: num_threads
+
+    !! Diffusivity to add due to the buoyancy adjustment by vertical mixing scheme
+    type(scalar_field) :: buoyancy_adjustment_diffusivity
 
     ewrite(1,*) "Writing advection-diffusion equation for "&
          &//trim(field_name)
@@ -797,9 +837,19 @@
        call allocate(Source, T%mesh, trim(field_name)//"Source",&
             FIELD_TYPE_CONSTANT)
        call zero(Source)
+       
+       add_src_directly_to_rhs = .false.
     else
        ! Grab an extra reference to cause the deallocate below to be safe.
        call incref(Source)
+      
+       add_src_directly_to_rhs = have_option(trim(Source%option_path)//'/diagnostic/add_directly_to_rhs')
+      
+       if (add_src_directly_to_rhs) then 
+          ewrite(2, *) "Adding Source field directly to the right hand side"
+          assert(node_count(Source) == node_count(T))
+       end if
+
     end if
 
     Absorption=extract_scalar_field(state, trim(field_name)//"Absorption"&
@@ -825,6 +875,18 @@
        include_advection=.true.
     end if
 
+    ! Porosity
+    if (have_option(trim(T%option_path)//'/prognostic/porosity')) then
+       include_porosity = .true.
+
+       ! get the porosity theta averaged field - this will allocate it
+       call form_porosity_theta(porosity_theta, state, option_path = trim(T%option_path)//'/prognostic/porosity')       
+    else
+       include_porosity = .false.
+       call allocate(porosity_theta, T%mesh, field_type=FIELD_TYPE_CONSTANT)
+       call set(porosity_theta, 1.0)
+    end if
+
     ! Retrieve scalar options from the options dictionary.
     if (.not.semi_discrete) then
        call get_option(trim(T%option_path)//&
@@ -842,6 +904,9 @@
            
     move_mesh = (have_option("/mesh_adaptivity/mesh_movement").and.include_mass)
     if(move_mesh) then
+      if (include_porosity) then
+         FLExit('Cannot include porosity in DG advection diffusion of a field with a moving mesh')
+      end if
       ewrite(2,*) 'Moving mesh'
       X_old => extract_vector_field(state, "OldCoordinate")
       X_new => extract_vector_field(state, "IteratedCoordinate")
@@ -875,7 +940,9 @@
     ! BCTYPE_WEAKDIRICHLET=1)
     allocate( bc_type(1:surface_element_count(T)) )
     call get_entire_boundary_condition(T, &
-       & (/"weakdirichlet"/), &
+       & (/"weakdirichlet", &
+       &   "dirichlet    ", &
+       &   "neumann      "/), &
        & bc_value, bc_type)
 
     call zero(big_m)
@@ -884,7 +951,8 @@
     if (present(diffusion_m)) call zero(diffusion_m)
     if (present(diffusion_RHS)) call zero(diffusion_RHS)
     if (have_buoyancy_adjustment_by_vertical_diffusion) then
-      if (have_option(trim(T%option_path)//"/prognostic/buoyancy_adjustment"&
+      ewrite(3,*) "Buoyancy adjustment by vertical mixing: enabled"
+      if (have_option(trim(T%option_path)//"/prognostic/buoyancy_adjustment"//&
           &"/by_vertical_diffusion/project_buoyancy_to_continuous_space")) then    
         buoyancy_from_state = extract_scalar_field(state, "VelocityBuoyancyDensity", stat)
         if (stat/=0) FLAbort('Error extracting buoyancy field.')
@@ -895,7 +963,9 @@
         ! Grab an extra reference to cause the deallocate below to be safe.
         ! Check this is OK
         call lumped_mass_galerkin_projection_scalar(state, buoyancy, buoyancy_from_state)
+        ewrite(3,*) "Buoyancy adjustment by vertical mixing: projecting to continuous space"
       else
+        ewrite(3,*) "Buoyancy adjustment by vertical mixing: no projection"
         buoyancy = extract_scalar_field(state, "VelocityBuoyancyDensity", stat)
         if (stat/=0) FLAbort('Error extracting buoyancy field.')
         call incref(buoyancy)
@@ -915,16 +985,64 @@
       end if
       ! Set direction of mixing diffusion, default is in the y- and z-direction for 2- and 3-d spaces respectively
       ! TODO: Align this direction with gravity local to an element
-    end if
-
-    element_loop: do ELE=1,element_count(T)
-       
+      
+      ! Check if the diagnostic associated with the buoyancy adjustment by vertical mixing scheme is required
+      buoyancy_adjustment_diffusivity = extract_scalar_field(state, "BuoyancyAdjustmentDiffusivity", stat)
+      if (stat==0) then
+        have_buoyancy_adjustment_diffusivity = .true.
+        ewrite(3,*) "Buoynacy adjustment by vertical mixing: Updating BuoyancyAdjustmentDiffusivity field."
+      else
+        have_buoyancy_adjustment_diffusivity = .false.
+      end if
+
+    end if
+
+    if (include_diffusion) then
+       call get_mesh_colouring(state, T%mesh, COLOURING_DG2, colours)
+#ifdef _OPENMP
+      if(diffusion_scheme == MASSLUMPED_RT0) then
+         call omp_set_num_threads(1)
+         ewrite(1,*) "WARNING: hybrid assembly can't support The MASSLUMPED_RT0 scheme yet, &
+         set threads back to 1"
+      endif
+#endif
+    else
+       call get_mesh_colouring(state, T%mesh, COLOURING_DG0, colours)
+    end if
+
+#ifdef _OPENMP
+    cache_valid = prepopulate_transform_cache(X)
+    assert(cache_valid)
+#endif
+
+    call profiler_tic(t, "advection_diffusion_dg_loop")
+
+    !$OMP PARALLEL DEFAULT(SHARED) &
+    !$OMP PRIVATE(clr, nnid, ele, len)
+
+    colour_loop: do clr = 1, size(colours) 
+      len = key_count(colours(clr))
+
+      !$OMP DO SCHEDULE(STATIC)
+      element_loop: do nnid = 1, len
+       ele = fetch(colours(clr), nnid)
        call construct_adv_diff_element_dg(ele, big_m, rhs, big_m_diff,&
             & rhs_diff, X, X_old, X_new, T, U_nl, U_mesh, Source, &
-            Absorption, Diffusivity, bc_value, bc_type, q_mesh, mass, &
-            & buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude) 
+            & Absorption, Diffusivity, bc_value, bc_type, q_mesh, mass, &
+            & buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude, &
+            & buoyancy_adjustment_diffusivity, &
+            & add_src_directly_to_rhs, porosity_theta) 
        
-    end do element_loop
+      end do element_loop
+      !$OMP END DO
+
+    end do colour_loop
+    !$OMP END PARALLEL
+
+    call profiler_toc(t, "advection_diffusion_dg_loop")
+    ! Add the source directly to the rhs if required 
+    ! which must be included before dirichlet BC's.
+    if (add_src_directly_to_rhs) call addto(rhs, Source)
     
     ! Drop any extra field references.
     if (have_buoyancy_adjustment_by_vertical_diffusion) call deallocate(buoyancy)
@@ -934,7 +1052,8 @@
     call deallocate(U_nl)
     call deallocate(U_nl_backup)
     call deallocate(bc_value)
-
+    call deallocate(porosity_theta)
+    
   end subroutine construct_advection_diffusion_dg
 
   subroutine lumped_mass_galerkin_projection_scalar(state, field, projected_field)
@@ -1023,14 +1142,13 @@
          
    end subroutine lumped_mass_galerkin_projection_scalar
 
-
-
-
-  subroutine construct_adv_diff_element_dg(ele, big_m, rhs, big_m_diff,&
+   subroutine construct_adv_diff_element_dg(ele, big_m, rhs, big_m_diff,&
        & rhs_diff, &
        & X, X_old, X_new, T, U_nl, U_mesh, Source, Absorption, Diffusivity,&
        & bc_value, bc_type, &
-       & q_mesh, mass, buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude)
+       & q_mesh, mass, buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude, &
+       & buoyancy_adjustment_diffusivity, &
+       & add_src_directly_to_rhs, porosity_theta)
     !!< Construct the advection_diffusion equation for discontinuous elements in
     !!< acceleration form.
     implicit none
@@ -1058,6 +1176,16 @@
     !! Diffusivity
     type(tensor_field), intent(in) :: Diffusivity
 
+    !! Diffusivity to add due to the buoyancy adjustment by vertical mixing scheme
+    type(scalar_field), intent(inout) :: buoyancy_adjustment_diffusivity
+    
+    !! If adding Source directly to rhs then
+    !! do nothing with it here
+    logical, intent(in) :: add_src_directly_to_rhs
+    
+    !! Porosity theta averaged field
+    type(scalar_field), intent(in) :: porosity_theta
+
     !! Flag for a periodic boundary
     logical :: Periodic_neigh 
 
@@ -1203,15 +1331,16 @@
       assert(ele_loc(U_mesh, ele)==ele_loc(X, ele))
       assert(ele_ngi(U_mesh, ele)==ele_ngi(U_nl, ele))
     end if
-
+    
     dg=continuity(T)<0
     primal = .not.dg
     if(diffusion_scheme == CDG) primal = .true.
     if(diffusion_scheme == IP) primal =.true.
 
-    ! In parallel, we only construct the equations on elements we own.
+    ! In parallel, we only construct the equations on elements we own, or
+    ! those in the L1 halo.
     if (dg) then
-       if (.not.element_owned(T,ele)) then
+       if (.not.(element_owned(T, ele).or.element_neighbour_owned(T, ele))) then
           return
        end if
     end if
@@ -1318,12 +1447,18 @@
                  &* gravity_magnitude * dr**2 * gravity_at_node(i) * drho_dz(:)
          end do
        end if
-       
-       !ewrite(3,*) "mixing_density_values", buoyancysample
-       ewrite(3,*) "mixing_grad_rho", minval(grad_rho(:,:)), maxval(grad_rho(:,:))
-       ewrite(3,*) "mixing_drho_dz", minval(drho_dz(:)), maxval(drho_dz(:))
-       ewrite(3,*) "mixing_coeffs amp dt g dr", mixing_diffusion_amplitude, dt, gravity_magnitude, dr**2
-       ewrite(3,*) "mixing_diffusion", minval(mixing_diffusion(2,2,:)), maxval(mixing_diffusion(2,2,:))
+        
+       if(have_buoyancy_adjustment_diffusivity) then
+        call set(buoyancy_adjustment_diffusivity, T_ele, mixing_diffusion_amplitude * dt&
+              &* gravity_magnitude * dr**2 * maxval(drho_dz(:)))
+         ewrite(4,*) "Buoynacy adjustment diffusivity, ele:", ele, "diffusivity:", mixing_diffusion_amplitude * dt * gravity_magnitude * dr**2 * maxval(drho_dz(:))
+       end if 
+        
+       !! Buoyancy adjustment by vertical mixing scheme debugging statements
+       ewrite(4,*) "mixing_grad_rho", minval(grad_rho(:,:)), maxval(grad_rho(:,:))
+       ewrite(4,*) "mixing_drho_dz", minval(drho_dz(:)), maxval(drho_dz(:))
+       ewrite(4,*) "mixing_coeffs amp dt g dr", mixing_diffusion_amplitude, dt, gravity_magnitude, dr**2
+       ewrite(4,*) "mixing_diffusion", minval(mixing_diffusion(2,2,:)), maxval(mixing_diffusion(2,2,:))
 
        mixing_diffusion_rhs=shape_tensor_rhs(T%mesh%shape, mixing_diffusion, detwei_rho)
        t_mass=shape_shape(T%mesh%shape, T%mesh%shape, detwei_rho)
@@ -1352,7 +1487,12 @@
     if(move_mesh) then
       mass_mat = shape_shape(T_shape, T_shape, detwei_new)
     else
-      mass_mat = shape_shape(T_shape, T_shape, detwei)
+      if (include_porosity) then
+        assert(ele_ngi(T, ele)==ele_ngi(porosity_theta, ele))
+        mass_mat = shape_shape(T_shape, T_shape, detwei*ele_val_at_quad(porosity_theta, ele))      
+      else
+        mass_mat = shape_shape(T_shape, T_shape, detwei)
+      end if
     end if
 
     if (include_advection) then
@@ -1513,9 +1653,11 @@
     end if
 
     ! Source term
-    l_T_rhs=l_T_rhs &
-         + shape_rhs(T_shape, detwei*ele_val_at_quad(Source, ele))
-         
+    if (.not. add_src_directly_to_rhs) then
+       l_T_rhs=l_T_rhs &
+              + shape_rhs(T_shape, detwei*ele_val_at_quad(Source, ele))
+    end if
+        
     if(move_mesh) then
       l_T_rhs=l_T_rhs &
          -shape_rhs(T_shape, ele_val_at_quad(t, ele)*(detwei_new-detwei_old)/dt)
@@ -1627,7 +1769,7 @@
        if(primal) then
           call construct_adv_diff_interface_dg(ele, face, face_2, ni,&
                & centre_vec,& 
-               & big_m, rhs, Grad_T_mat, Div_T_mat, X, T, U_nl,&
+               & big_m, rhs, rhs_diff, Grad_T_mat, Div_T_mat, X, T, U_nl,&
                & bc_value, bc_type, &
                & U_mesh, q_mesh, cdg_switch_in, &
                & primal_fluxes_mat, ele2grad_mat,diffusivity, &
@@ -1646,7 +1788,7 @@
        else
           call construct_adv_diff_interface_dg(ele, face, face_2, ni,&
                & centre_vec,&
-               & big_m, rhs, Grad_T_mat, Div_T_mat, X, T, U_nl,&
+               & big_m, rhs, rhs_diff, Grad_T_mat, Div_T_mat, X, T, U_nl,&
                & bc_value, bc_type, &
                & U_mesh, q_mesh)
        end if
@@ -1716,18 +1858,18 @@
                     finish=start+face_loc(T, face)-1
                     
                     if (i==1) then
-                      ! Wipe out boundary condition's coupling to itself.
-                      Diffusivity_mat(start:finish,:)=0.0
+                       ! Wipe out boundary condition's coupling to itself.
+                       Diffusivity_mat(start:finish,:)=0.0
                     else
                       
-                      ! Add BC into RHS
-                      !
-                      call addto(RHS_diff, local_glno, &
-                           & -matmul(Diffusivity_mat(:,start:finish), &
-                           & ele_val( bc_value, face )))
-                      
-                      ! Ensure it is not used again.
-                      Diffusivity_mat(:,start:finish)=0.0
+                       ! Add BC into RHS
+                       !
+                       call addto(RHS_diff, local_glno, &
+                            & -matmul(Diffusivity_mat(:,start:finish), &
+                            & ele_val( bc_value, face )))
+
+                       ! Ensure it is not used again.
+                       Diffusivity_mat(:,start:finish)=0.0
                       
                     end if
                   
@@ -1753,7 +1895,7 @@
     if (include_diffusion.or.have_buoyancy_adjustment_by_vertical_diffusion) then
        call addto(Big_m_diff, local_glno, local_glno,&
             & Diffusivity_mat*theta*dt)
-       
+
        if (.not.semi_discrete) then
           call addto(RHS_diff, local_glno, &
                & -matmul(Diffusivity_mat, node_val(T, local_glno)))
@@ -2211,7 +2353,7 @@
   end subroutine construct_adv_diff_element_dg
   
   subroutine construct_adv_diff_interface_dg(ele, face, face_2, &
-       ni, centre_vec,big_m, rhs, Grad_T_mat, Div_T_mat, &
+       ni, centre_vec,big_m, rhs, rhs_diff, Grad_T_mat, Div_T_mat, &
        & X, T, U_nl,&
        & bc_value, bc_type, &
        & U_mesh, q_mesh, CDG_switch_in, &
@@ -2224,7 +2366,7 @@
 
     integer, intent(in) :: ele, face, face_2, ni
     type(csr_matrix), intent(inout) :: big_m
-    type(scalar_field), intent(inout) :: rhs
+    type(scalar_field), intent(inout) :: rhs, rhs_diff
     real, dimension(:,:,:), intent(inout) :: Grad_T_mat, Div_T_mat
     ! We pass these additional fields to save on state lookups.
     type(vector_field), intent(in) :: X, U_nl
@@ -2275,7 +2417,7 @@
     real, dimension(:,:,:), allocatable :: kappa_gi
 
     integer :: dim, start, finish
-    logical :: boundary, dirichlet
+    logical :: boundary, dirichlet, neumann
 
     logical :: do_primal_fluxes
 
@@ -2318,9 +2460,12 @@
     ! Boundary nodes have both faces the same.
     boundary=(face==face_2)
     dirichlet=.false.
+    neumann=.false.
     if (boundary) then
        if (bc_type(face)==BCTYPE_WEAKDIRICHLET) then
          dirichlet=.true.
+       elseif (bc_type(face)==BCTYPE_NEUMANN) then
+         neumann=.true.
        end if
     end if
 
@@ -2504,13 +2649,13 @@
        end if
 
        ! Insert advection in RHS.
-       
+
        if (.not.dirichlet) then
           ! For interior interfaces this is the upwinding term. For a Neumann
           ! boundary it's necessary to apply downwinding here to maintain the
           ! surface integral. Fortunately, since face_2==face for a boundary
           ! this is automagic.
-          
+
           if (.not.semi_discrete) then
              call addto(RHS, T_face, &
                   ! Outflow boundary integral.
@@ -2525,7 +2670,7 @@
           call addto(RHS, T_face, &
                -matmul(nnAdvection_in,&
                ele_val(bc_value, face)))
-          
+
           if(.not.semi_discrete) then
              ! The interior integral is still interior!
              call addto(RHS, T_face, &
@@ -2533,9 +2678,14 @@
           end if
 
        end if
-       
+      
     end if
 
+   ! Add non-zero contributions from Neumann boundary conditions (if present)
+   if (neumann) then
+      call addto(RHS_diff, T_face, shape_rhs(T_shape, detwei * ele_val_at_quad(bc_value, face)))
+   end if
+
   contains
 
     subroutine arbitrary_upwind_diffusion

=== modified file 'assemble/Advection_Diffusion_FV.F90'
--- assemble/Advection_Diffusion_FV.F90	2011-04-22 17:28:48 +0000
+++ assemble/Advection_Diffusion_FV.F90	2012-10-16 13:22:44 +0000
@@ -45,8 +45,12 @@
   use spud
   use field_options
   use sparsity_patterns_meshes
-  use global_parameters, only : FIELD_NAME_LEN, OPTION_PATH_LEN
+  use global_parameters, only : FIELD_NAME_LEN, OPTION_PATH_LEN, COLOURING_DG1
   use profiler
+  use colouring
+#ifdef _OPENMP
+  use omp_lib
+#endif
   
   implicit none
 
@@ -59,6 +63,8 @@
   logical :: have_advection
   ! Source?
   logical :: have_source
+  ! Add source directly to the right hand side?
+  logical :: add_src_directly_to_rhs
   ! Absorption?
   logical :: have_absorption
   ! Diffusivity?
@@ -138,7 +144,7 @@
     type(scalar_field), intent(inout) :: t
     type(csr_matrix), intent(inout) :: matrix
     type(scalar_field), intent(inout) :: rhs
-    type(state_type), intent(in) :: state
+    type(state_type), intent(inout) :: state
     
     type(vector_field), pointer :: coordinate, &
                                    old_coordinate, new_coordinate, &
@@ -147,7 +153,14 @@
     type(scalar_field), pointer :: source, absorption
     type(tensor_field), pointer :: diffusivity
     
-    integer :: i, j, ele, stat
+    integer :: i, j, stat
+
+    !! Coloring  data structures for OpenMP parallization
+    type(integer_set), dimension(:), pointer :: colours
+    integer :: clr, nnid, len, ele
+    integer :: thread_num
+    !! Did we successfully prepopulate the transform_to_physical_cache?
+    logical :: cache_valid
     
     ewrite(1,*) "In assemble_advection_diffusion_fv"
     
@@ -160,10 +173,19 @@
     if(have_source) then
       assert(mesh_dim(source) == mesh_dim(t))
       assert(ele_count(source) == ele_count(t))
+      
+      add_src_directly_to_rhs = have_option(trim(source%option_path)//'/diagnostic/add_directly_to_rhs')
+      
+      if (add_src_directly_to_rhs) then 
+         ewrite(2, *) "Adding Source field directly to the right hand side"
+         assert(node_count(source) == node_count(t))
+      end if
     
       ewrite_minmax(source)
     else
       ewrite(2,*) 'No source'
+      
+      add_src_directly_to_rhs = .false.
     end if
 
     ! Absorption
@@ -250,11 +272,45 @@
     call zero(matrix)
     call zero(rhs)
     
-    do ele = 1, ele_count(t)
-      call assemble_advection_diffusion_element_fv(ele, t, matrix, rhs, &
+
+#ifdef _OPENMP
+    cache_valid = prepopulate_transform_cache(coordinate)
+    assert(cache_valid)
+#endif
+
+    call get_mesh_colouring(state, T%mesh, COLOURING_DG1, colours)
+
+    call profiler_tic(t, "advection_diffusion_fv_loop")
+
+    !$OMP PARALLEL DEFAULT(SHARED) &
+    !$OMP PRIVATE(clr, len, nnid, ele, thread_num)
+
+#ifdef _OPENMP    
+    thread_num = omp_get_thread_num()
+#else
+    thread_num=0
+#endif
+
+
+    colour_loop: do clr = 1, size(colours)
+      len = key_count(colours(clr))
+      !$OMP DO SCHEDULE(STATIC)
+      element_loop: do nnid = 1, len
+       ele = fetch(colours(clr), nnid)
+       call assemble_advection_diffusion_element_fv(ele, t, matrix, rhs, &
                                                    coordinate, t_coordinate, &
                                                    source, absorption, diffusivity)
-    end do
+      end do element_loop
+       !$OMP END DO
+
+    end do colour_loop
+    !$OMP END PARALLEL
+
+    call profiler_toc(t, "advection_diffusion_fv_loop")
+
+    ! Add the source directly to the rhs if required 
+    ! which must be included before dirichlet BC's.
+    if (add_src_directly_to_rhs) call addto(rhs, source)
 
     ewrite(2, *) "Applying strong Dirichlet boundary conditions"
     call apply_dirichlet_conditions(matrix, rhs, t, dt)
@@ -314,7 +370,9 @@
     if(have_absorption) call add_absorption_element_fv(ele, t_shape, t, absorption, detwei, matrix_addto(:loc,:loc), rhs_addto(:loc))
     
     ! Source
-    if(have_source) call add_source_element_fv(ele, t_shape, t, source, detwei, rhs_addto(:loc))
+    if(have_source .and. (.not. add_src_directly_to_rhs)) then 
+       call add_source_element_fv(ele, t_shape, t, source, detwei, rhs_addto(:loc))
+    end if
     
     if(have_diffusivity.or.have_advection) then
 

=== modified file 'assemble/Assemble_CMC.F90'
--- assemble/Assemble_CMC.F90	2011-04-22 17:28:48 +0000
+++ assemble/Assemble_CMC.F90	2012-10-16 13:22:44 +0000
@@ -172,25 +172,28 @@
 
     subroutine assemble_scaled_pressure_mass_matrix(state,scaled_pressure_mass_matrix)
 
-      ! This routine assembles the scaled_pressure_mass_matrix. It is scaled
-      ! by the inverse of viscosity.
+      ! This routine assembles the scaled_pressure_mass_matrix at the
+      ! quadrature points. It is scaled by the inverse of viscosity.
 
       type(state_type), intent(in) :: state   
-
       ! Scaled pressure mass matrix - already allocated in Momentum_Eq:
       type(csr_matrix), intent(inout) :: scaled_pressure_mass_matrix
-
       ! Pressure field:
       type(scalar_field), pointer :: pressure
       ! Viscosity tensor:
       type(tensor_field), pointer :: viscosity     
       ! Viscosity component:
       type(scalar_field) :: viscosity_component
-      ! Inverse of Viscosity component
-      type(scalar_field) :: inverse_viscosity_component
       ! Positions:
       type(vector_field), pointer :: positions
-      
+
+      ! Relevant declerations for mass matrix calculation:
+      integer :: ele
+      real, dimension(:), allocatable :: detwei
+      type(element_type), pointer :: p_shape
+      real, dimension(:,:), allocatable :: mass_matrix
+      real, dimension(:), allocatable :: mu_gi
+
       ewrite(1,*) 'Entering assemble_scaled_pressure_mass_matrix'    
 
       ! Positions:
@@ -205,20 +208,24 @@
       ! Extract first component of viscosity tensor from full tensor:
       viscosity_component = extract_scalar_field(viscosity,1,1)
 
-      ! Allocate memory for inverse viscosity component scalar field:
-      call allocate(inverse_viscosity_component, viscosity_component%mesh, name="inverse_viscosity_component")
-
-      ! Invert viscosity:
-      call invert(viscosity_component,inverse_viscosity_component)
-
-      ! Compute pressure mass matrix, scaled by inverse viscosity - note that the
-      ! inverse viscosity is supplied as density here:
-      call compute_mass(positions,pressure%mesh,scaled_pressure_mass_matrix,density=inverse_viscosity_component)
+      ! Initialise and assemble scaled pressure mass matrix:
+      allocate(detwei(ele_ngi(pressure, 1)), &
+               mass_matrix(ele_loc(pressure, 1), ele_loc(pressure, 1)), &
+               mu_gi(ele_ngi(viscosity_component, 1)))
+ 
+      call zero(scaled_pressure_mass_matrix)
+      do ele = 1, ele_count(pressure)
+        p_shape => ele_shape(pressure, ele)
+        mu_gi = ele_val_at_quad(viscosity_component, ele)
+        call transform_to_physical(positions, ele, detwei=detwei)
+        mass_matrix = shape_shape(p_shape, p_shape, detwei/mu_gi)
+        call addto(scaled_pressure_mass_matrix, ele_nodes(pressure, ele),&
+             ele_nodes(pressure, ele), mass_matrix)
+      end do
 
       ewrite_minmax(scaled_pressure_mass_matrix)
 
-      ! Deallocate inverse viscosity component scalar field:
-      call deallocate(inverse_viscosity_component)
+      deallocate(detwei, mass_matrix, mu_gi)
 
     end subroutine assemble_scaled_pressure_mass_matrix
       

=== modified file 'assemble/Biology.F90'
--- assemble/Biology.F90	2011-10-19 17:55:20 +0000
+++ assemble/Biology.F90	2012-10-16 13:22:44 +0000
@@ -39,6 +39,9 @@
   use solvers
   use python_state
   use sparsity_patterns_meshes
+  use fefields
+  use field_options
+
   implicit none
 
   private
@@ -66,6 +69,14 @@
     type(state_type), intent(inout) :: state
 
     character(len=OPTION_PATH_LEN) :: prefix, algorithm
+    ! This is the photosynthetic radiation projected onto the 
+    ! same mesh as the biology fields
+    ! It also takes into account the "active" part of the solar radiation
+    type(scalar_field) :: par_bio
+    ! we use the phytoplankton as the "bio" mesh
+    type(scalar_field), pointer :: phytoplankton, PhotosyntheticRadiation
+    type(vector_field) :: coords
+    integer :: stat
 
     call backup_source_terms(state)
 
@@ -91,6 +102,24 @@
     ! Calculate the light field at every point.
     call solve_light_equation(state, prefix)
 
+    par_bio = extract_scalar_field(state, "_PAR", stat)
+    phytoplankton => extract_scalar_field(state, "Phytoplankton")
+    if (stat /= 0) then
+        ! field does not yet exist: create it
+        call allocate(par_bio,phytoplankton%mesh, name="_PAR")
+        call zero(par_bio)
+        call insert(state, par_bio, par_bio%name)
+        call deallocate(par_bio)
+        par_bio = extract_scalar_field(state, "_PAR", stat)
+    end if
+    PhotosyntheticRadiation => extract_scalar_field(state, "PhotosyntheticRadiation")
+    coords = get_coordinate_field(state, par_bio%mesh)
+    ! project the Photosynthetic radaition field onto the _PAR field
+    call project_field(PhotosyntheticRadiation, par_bio, coords)
+    ! scale it to get the active part
+    call scale(par_bio, 0.43)
+    call deallocate(coords)
+
     ! Calculate the sources and sinks at every point.
     call calculate_biology_from_python(state, prefix, algorithm)
 
@@ -210,13 +239,9 @@
     g=extract_vector_field(state, "GravityDirection")
     light=extract_scalar_field(state, "PhotosyntheticRadiation")
     P=extract_scalar_field(state, "Phytoplankton")
-    
-    if(continuity(light)<0) then
-      ! could be cleverer with this and check if diffusion is assembled
-      sparsity=>get_csr_sparsity_secondorder(state, light%mesh, light%mesh)
-    else
-      sparsity=>get_csr_sparsity_firstorder(state, light%mesh, light%mesh)
-    end if
+
+    ! Only need first order sparsity as we have no diffusion
+    sparsity=>get_csr_sparsity_firstorder(state, light%mesh, light%mesh)
 
     call get_option(trim(light%option_path)//&
          "/prognostic/absorption_coefficients/water", k_w)
@@ -369,23 +394,22 @@
           ! External face.
           face_2=face
        end if
-
-       call construct_light_interface(ele, ele_2, face, face_2, ni,&
+       
+       call construct_light_interface(face, face_2,&
             & light_mat, rhs, X, g, light, &
             & bc_value, bc_type)  
-            
     end do neighbourloop
 
   end subroutine construct_light_element
 
-  subroutine construct_light_interface(ele, ele_2, face, face_2, ni,&
+  subroutine construct_light_interface(face, face_2,&
             & light_mat, rhs, X, g, light, bc_value, bc_type)
     !!< Construct the element boundary integrals on the ni-th face of
     !!< element ele. For continuous discretisation, this is only boundary
     !!< faces. For DG it's all of them.
     implicit none
 
-    integer, intent(in) :: ele, ele_2, face, face_2, ni
+    integer, intent(in) :: face, face_2
     type(csr_matrix), intent(inout) :: light_mat
     type(scalar_field), intent(inout) :: rhs
     ! We pass these additional fields to save on state lookups.
@@ -399,7 +423,7 @@
 
     ! Face objects and numberings.
     type(element_type), pointer ::l_shape, l_shape_2
-    integer, dimension(face_loc(light,face)) :: l_face, l_face_l
+    integer, dimension(face_loc(light,face)) :: l_face
     integer, dimension(face_loc(light,face_2)) :: l_face_2
 
     ! Note that both sides of the face can be assumed to have the same
@@ -428,7 +452,7 @@
     ! Boundary nodes have both faces the same.
     boundary=(face==face_2)
     dirichlet=.false.
-    if (boundary) then
+    if (boundary .and. face < size(bc_type)) then
        if (bc_type(face)==BCTYPE_WEAKDIRICHLET) then
           dirichlet=.true.
        end if

=== modified file 'assemble/Compressible_Projection.F90'
--- assemble/Compressible_Projection.F90	2011-07-18 17:11:34 +0000
+++ assemble/Compressible_Projection.F90	2012-10-16 13:22:44 +0000
@@ -36,9 +36,10 @@
   use field_options
   use equation_of_state, only: compressible_eos, compressible_material_eos
   use global_parameters, only: OPTION_PATH_LEN
-  use fefields, only: compute_lumped_mass
+  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)
 
@@ -107,7 +111,7 @@
     type(scalar_field), pointer :: normalisation, &
                                    density, olddensity
     type(scalar_field), pointer :: pressure
-    type(scalar_field), pointer :: p_lumpedmass
+    type(scalar_field), pointer :: p_cvmass
     type(scalar_field) :: lhsfield, invnorm, absrhs
     
     type(scalar_field), pointer :: source, absorption
@@ -126,14 +130,10 @@
       call get_option(trim(pressure%option_path)//'/prognostic/atmospheric_pressure', &
                       atmospheric_pressure, default=0.0)
       
-      if(pressure%mesh%shape%degree>1) then
-        ! try lumping on the submesh
-        p_lumpedmass => get_lumped_mass_on_submesh(state, pressure%mesh)
-      else
-        ! find the lumped mass
-        p_lumpedmass => get_lumped_mass(state, pressure%mesh)
-      end if
-      ewrite_minmax(p_lumpedmass)
+      ! find the cv mass
+      p_cvmass => get_cv_mass(state, pressure%mesh)
+
+      ewrite_minmax(p_cvmass)
       
       call get_option(trim(pressure%option_path)//"/prognostic/scheme/use_compressible_projection_method/normalisation/name", &
                       normalisation_field, stat=norm_stat)
@@ -170,12 +170,12 @@
 
       call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", theta)
 
-      call set(lhsfield, p_lumpedmass)
+      call set(lhsfield, p_cvmass)
       call scale(lhsfield, drhodp)
       call scale(lhsfield, invnorm)
       call addto_diag(cmc, lhsfield, scale=1./(dt*dt*theta_divergence*theta_pg))
       
-!     rhs = invnorm*p_lumpedmass* &
+!     rhs = invnorm*p_cvmass* &
 !      ( (1./dt)*(olddensity - density + drhodp*(eospressure - (pressure + atmospheric_pressure)))
 !       +(absorption)*(drhodp*theta_pg*(eospressure - (pressure + atmospheric_pressure)) - theta_pg*density - (1-theta_pg)*olddensity)
 !       +source)
@@ -215,7 +215,7 @@
         call addto_diag(cmc, lhsfield, scale=(theta/(dt*theta_divergence*theta_pg)))
       end if
       
-      call scale(rhs, p_lumpedmass)
+      call scale(rhs, p_cvmass)
       call scale(rhs, invnorm)
       
       call deallocate(eospressure)
@@ -253,7 +253,7 @@
 
     type(scalar_field), pointer :: pressure
     type(vector_field), pointer :: positions
-    type(scalar_field) :: lumped_mass, tempfield
+    type(scalar_field) :: cv_mass, tempfield
 
     real :: atmospheric_pressure
 
@@ -274,9 +274,9 @@
       if(cmcget) then
 
         positions=>extract_vector_field(state(1), "Coordinate")
-        call allocate(lumped_mass, pressure%mesh, "LumpedMassField")
+        call allocate(cv_mass, pressure%mesh, "CVMassField")
         call allocate(tempfield, pressure%mesh, "TemporaryAssemblyField")
-        call compute_lumped_mass(positions, lumped_mass)
+        call compute_cv_mass(positions, cv_mass)
 
         allocate(dummy_ones)
         call allocate(dummy_ones, pressure%mesh, "DummyOnesField")
@@ -337,16 +337,16 @@
         end do
 
         call zero(tempfield)
-        tempfield%val = (1./(dt*dt))*lumped_mass%val*normdrhodp%val
+        tempfield%val = (1./(dt*dt))*cv_mass%val*normdrhodp%val
 
         call addto_diag(cmc, tempfield)
 
-        rhs%val = (1./dt)*lumped_mass%val* &
+        rhs%val = (1./dt)*cv_mass%val* &
                           ( &
                             normolddensity%val &
                           - normdensity%val &
                           ) &
-               +(1./dt)*lumped_mass%val* &
+               +(1./dt)*cv_mass%val* &
                           ( &
                             normmatdrhodpp%val &
                           - normdrhodp%val*(pressure%val+atmospheric_pressure) &
@@ -360,7 +360,7 @@
         call deallocate(materialpressure)
         call deallocate(materialdrhodp)
 
-        call deallocate(lumped_mass)
+        call deallocate(cv_mass)
         call deallocate(tempfield)
         call deallocate(dummy_ones)
         deallocate(dummy_ones)
@@ -371,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
@@ -384,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
 
@@ -445,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).
     ! =============================================================
@@ -474,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")
   
@@ -556,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* &
         !       /
@@ -564,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))
@@ -591,6 +624,10 @@
   
       call deallocate(drhodp)
       call deallocate(eospressure)
+
+      if(multiphase) then
+         call deallocate(nvfrac)
+      end if
     
     end if
 
@@ -601,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
@@ -624,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/Coriolis.F90'
--- assemble/Coriolis.F90	2011-03-15 22:19:21 +0000
+++ assemble/Coriolis.F90	2012-10-16 13:22:44 +0000
@@ -31,6 +31,7 @@
 
   use embed_python 
   use fldebug
+  use parallel_tools, only: abort_if_in_parallel_region
   use global_parameters, only : current_time, PYTHON_FUNC_LEN
   use spud
   
@@ -145,7 +146,9 @@
     real:: omega
     
     if (coriolis_option/=NOT_INITIALISED) return
-    
+
+    call abort_if_in_parallel_region
+
     ewrite(1, *) "Initialising Coriolis"
     
     call get_option("/geometry/dimension", coriolis_dim)

=== modified file 'assemble/Diagnostic_Fields_Matrices.F90'
--- assemble/Diagnostic_Fields_Matrices.F90	2011-04-28 11:10:02 +0000
+++ assemble/Diagnostic_Fields_Matrices.F90	2012-10-16 13:22:44 +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
 
@@ -73,8 +74,8 @@
 
       character(len=FIELD_NAME_LEN) :: field_name
 
-      type(scalar_field), pointer :: lumpedmass
-      type(scalar_field) :: inverse_lumpedmass
+      type(scalar_field), pointer :: cvmass
+      type(scalar_field) :: inverse_cvmass
       type(scalar_field) :: ctfield, ct_rhs
 
       call get_option(trim(div%option_path)//"/diagnostic/field_name", field_name)
@@ -85,23 +86,19 @@
       call allocate(ctfield, div%mesh, name="CTField")
       call allocate(ct_rhs, div%mesh, name="CTRHS")
 
-      if(element_degree(div, 1)>1) then
-        lumpedmass => get_lumped_mass_on_submesh(state, div%mesh)
-      else
-        lumpedmass => get_lumped_mass(state, div%mesh)
-      end if
+      cvmass => get_cv_mass(state, div%mesh)
 
       CT_m => get_divergence_matrix_cv(state, test_mesh=div%mesh, field=field, div_rhs=ct_rhs)
 
       call mult(ctfield, CT_m, field)
       call addto(ctfield, ct_rhs, -1.0)
 
-      call allocate(inverse_lumpedmass, lumpedmass%mesh, "InverseLumpedMass")
-      call invert(lumpedmass, inverse_lumpedmass)
+      call allocate(inverse_cvmass, cvmass%mesh, "InverseCVMass")
+      call invert(cvmass, inverse_cvmass)
       call set(div, ctfield)
-      call scale(div, inverse_lumpedmass)
+      call scale(div, inverse_cvmass)
       
-      call deallocate(inverse_lumpedmass)
+      call deallocate(inverse_cvmass)
       call deallocate(ctfield)
       call deallocate(ct_rhs)
 
@@ -188,7 +185,7 @@
       type(vector_field), pointer :: field, x
 
       type(csr_sparsity) :: divergence_sparsity
-      type(block_csr_matrix) :: CT_m
+      type(block_csr_matrix), pointer :: CT_m
 
       character(len=FIELD_NAME_LEN) :: field_name
 
@@ -204,6 +201,7 @@
       call allocate(ctfield, div%mesh, name="CTField")
 
       divergence_sparsity=make_sparsity(div%mesh, field%mesh, "DivergenceSparsity")
+      allocate(CT_m)
       call allocate(CT_m, divergence_sparsity, (/1, field%dim/), name="DivergenceMatrix" )
       call allocate(ct_rhs, div%mesh, name="CTRHS")
 
@@ -223,6 +221,7 @@
       call petsc_solve(div, mass, ctfield)
 
       call deallocate(CT_m)
+      deallocate(CT_m)
       call deallocate(ct_rhs)
       call deallocate(ctfield)
       call deallocate(divergence_sparsity)
@@ -239,7 +238,7 @@
       type(scalar_field), pointer :: field
 
       type(csr_sparsity) :: divergence_sparsity
-      type(block_csr_matrix) :: CT_m
+      type(block_csr_matrix), pointer :: CT_m
 
       character(len=FIELD_NAME_LEN) :: field_name
 
@@ -256,6 +255,7 @@
 
       ! Sparsity of C^T - the transpose of the gradient operator.
       divergence_sparsity=make_sparsity(field%mesh, grad%mesh, "DivergenceSparsity")
+      allocate(CT_m)
       call allocate(CT_m, divergence_sparsity, (/1, grad%dim/), name="DivergenceMatrix" )
 
       mass_sparsity=make_sparsity(grad%mesh, grad%mesh, "MassSparsity")
@@ -273,6 +273,7 @@
 
       call deallocate(divergence_sparsity)
       call deallocate(CT_m)
+      deallocate(CT_m)
       call deallocate(mass_sparsity)
       call deallocate(mass)
       call deallocate(cfield)
@@ -340,23 +341,32 @@
       integer :: i, stat
 
       type(csr_sparsity) :: divergence_sparsity
-      type(block_csr_matrix) :: ct_m
+      type(block_csr_matrix), pointer :: ct_m
 
       type(csr_sparsity) :: mass_sparsity
       type(csr_matrix) :: mass
       type(scalar_field) :: ctfield, ct_rhs, temp
+      type(scalar_field), pointer :: cv_mass
+      
+      logical :: test_with_cv_dual
 
       ewrite(1,*) 'Entering calculate_sum_velocity_divergence'
+         
+      ! Are we testing the divergence with the CV dual mesh
+      test_with_cv_dual = have_option(trim(sum_velocity_divergence%option_path)//'/diagnostic/test_with_cv_dual')   
 
       ! 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")
 
-      mass_sparsity=make_sparsity(sum_velocity_divergence%mesh, sum_velocity_divergence%mesh, "MassSparsity")
-      call allocate(mass, mass_sparsity, name="MassMatrix")
-      call zero(mass)
-
+      ! 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")
+         call allocate(mass, mass_sparsity, name="MassMatrix")
+         call zero(mass)      
+      end if 
+      
       ! Sum up over the div's
       do i = 1, size(state)
          u => extract_vector_field(state(i), "Velocity", stat)
@@ -373,21 +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(i==1) then ! Construct the mass matrix (just do this once)          
-            call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, &
+         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)
+               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)
-         else
-            call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, &
+            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   
+
          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)
@@ -396,23 +414,180 @@
          call addto(ctfield, temp)
 
          call deallocate(ct_m)
+         deallocate(ct_m)
          call deallocate(ct_rhs)
          call deallocate(divergence_sparsity)
 
       end do
 
-      ! Solve for sum_velocity_divergence ( = \sum{div(vfrac*u)} )
-      call zero(sum_velocity_divergence)
-      call petsc_solve(sum_velocity_divergence, mass, ctfield)
-
+      ! 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)
+         call set(sum_velocity_divergence, cv_mass)
+         call invert(sum_velocity_divergence)
+         call scale(sum_velocity_divergence, ctfield)
+      else
+         call zero(sum_velocity_divergence)
+         call petsc_solve(sum_velocity_divergence, mass, ctfield)
+      end if 
+      
       ! Deallocate memory
-      call deallocate(mass_sparsity)
-      call deallocate(mass)
       call deallocate(ctfield)
       call deallocate(temp)
+      
+      if (.not. test_with_cv_dual) then
+         call deallocate(mass_sparsity)
+         call deallocate(mass)      
+      end if 
 
       ewrite(1,*) 'Exiting calculate_sum_velocity_divergence'
          
   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	2011-04-28 11:10:02 +0000
+++ assemble/Diagnostic_fields_wrapper.F90	2012-10-16 13:22:44 +0000
@@ -30,7 +30,7 @@
 module diagnostic_fields_wrapper
   !!< A module to link to diagnostic variable calculations.
 
-  use global_parameters, only:FIELD_NAME_LEN 
+  use global_parameters, only: FIELD_NAME_LEN, timestep
   use fields
   use sparse_matrices_fields
   use field_derivatives
@@ -52,7 +52,6 @@
   use diagnostic_fields_matrices
   use equation_of_state
   use momentum_diagnostic_fields
-  use porous_media, only: calculate_porous_media_absorption
   use spontaneous_potentials, only: calculate_formation_conductivity
   use sediment_diagnostics
   use geostrophic_pressure
@@ -135,6 +134,30 @@
          end if
        end if
 
+       s_field => extract_scalar_field(state(i), "InterstitialVelocityCGCourantNumber", stat)
+       if(stat == 0) then
+         if(recalculate(trim(s_field%option_path))) then
+           call calculate_diagnostic_variable(state(i), "InterstitialVelocityCGCourantNumber", &
+             & s_field)
+         end if
+       end if
+
+       s_field => extract_scalar_field(state(i), "InterstitialVelocityDGCourantNumber", stat)
+       if(stat == 0) then
+         if(recalculate(trim(s_field%option_path))) then
+           call calculate_diagnostic_variable(state(i), "InterstitialVelocityDGCourantNumber", &
+             & s_field)
+         end if
+       end if
+
+       s_field => extract_scalar_field(state(i), "InterstitialVelocityCVCourantNumber", stat)
+       if(stat == 0) then
+         if(recalculate(trim(s_field%option_path))) then
+           call calculate_diagnostic_variable(state(i), "InterstitialVelocityCVCourantNumber", &
+             & s_field)
+         end if
+       end if
+
        s_field => extract_scalar_field(state(i), "KineticEnergyDensity", stat)
        if(stat == 0) then
          if(recalculate(trim(s_field%option_path))) then
@@ -519,17 +542,6 @@
        end if
        ! End of vorticity diagnostics
 
-       ! Start of porous media diagnostics
-       if (have_option("/porous_media")) then
-         v_field => extract_vector_field(state(i), "VelocityAbsorption", stat)
-         if(stat == 0) then
-           if(recalculate(trim(v_field%option_path))) then
-             call calculate_porous_media_absorption(state, i, v_field, stat)
-           end if
-         end if
-       end if
-       ! End of porous media diagnostics
-
        ! Start of spontaneous potentials diagnostics
        if(i == 1) then
          s_field => extract_scalar_field(state(i), "ElectricalConductivity", stat)
@@ -544,7 +556,11 @@
 
        ! Start of sediment diagnostics.
        if (have_option("/material_phase[0]/sediment")) then
+          call calculate_sediment_sinking_velocity(state(i))
           call calculate_sediment_flux(state(i))
+          call calculate_sediment_active_layer_d50(state(i))
+          call calculate_sediment_active_layer_sigma(state(i))
+          call calculate_sediment_active_layer_volume_fractions(state(i))
        end if
        ! End of sediment diagnostics.
 
@@ -573,6 +589,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	2011-09-17 01:24:08 +0000
+++ assemble/Divergence_Matrix_CG.F90	2012-10-16 13:22:44 +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, &
@@ -71,7 +74,7 @@
       type(state_type), intent(inout) :: state
 
       ! the velocity divergence gradient matrices
-      type(block_csr_matrix), intent(inout) :: CT_m
+      type(block_csr_matrix), pointer :: CT_m
 
       type(scalar_field), intent(inout), optional :: ct_rhs
 
@@ -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), intent(inout) :: ctp_m
+      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), intent(inout) :: ctp_m
+      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/Divergence_Matrix_CV.F90'
--- assemble/Divergence_Matrix_CV.F90	2011-09-16 21:18:56 +0000
+++ assemble/Divergence_Matrix_CV.F90	2012-10-16 13:22:44 +0000
@@ -107,7 +107,7 @@
       type(vector_field) :: field_bc
 
       real, dimension(:,:,:), allocatable :: ct_mat_local, ct_mat_local_bdy
-      real, dimension(:), allocatable :: ct_rhs_local
+      real, dimension(:,:), allocatable :: ct_rhs_local
 
       logical :: l_get_ct
 
@@ -316,7 +316,7 @@
         allocate(field_nodes_bdy(field%mesh%faces%shape%loc))
         allocate(test_nodes_bdy(test_mesh%faces%shape%loc))
         allocate(ct_mat_local_bdy(x%dim, test_mesh%faces%shape%loc, field%mesh%faces%shape%loc), &
-                ct_rhs_local(test_mesh%faces%shape%loc))
+                 ct_rhs_local(x%dim, test_mesh%faces%shape%loc))
   
         surface_element_loop: do sele = 1, surface_element_count(test_mesh)
   
@@ -366,11 +366,11 @@
                       if((present(ct_rhs)).and.(field_bc_type(dim, sele)==BC_TYPE_WEAKDIRICHLET)) then
   
                         if(multiphase) then
-                           ct_rhs_local(iloc) = ct_rhs_local(iloc) + &
+                           ct_rhs_local(dim, iloc) = ct_rhs_local(dim, iloc) - &
                                  field_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*nvfrac_gi_f(ggi)*&
                                  normal_bdy(dim,ggi)*field_bc_val(dim, jloc)
                         else
-                           ct_rhs_local(iloc) = ct_rhs_local(iloc) + &
+                           ct_rhs_local(dim, iloc) = ct_rhs_local(dim, iloc) - &
                                  field_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*normal_bdy(dim,ggi)*&
                                  field_bc_val(dim, jloc)
                         end if
@@ -403,7 +403,7 @@
   
             if((present(ct_rhs)).and.(field_bc_type(dim, sele)==BC_TYPE_WEAKDIRICHLET)) then
   
-              call addto(ct_rhs, test_nodes_bdy, ct_rhs_local)
+              call addto(ct_rhs, test_nodes_bdy, ct_rhs_local(dim,:))
   
             elseif(l_get_ct) then
   

=== modified file 'assemble/Field_Equations_CV.F90'
--- assemble/Field_Equations_CV.F90	2011-10-10 21:03:46 +0000
+++ assemble/Field_Equations_CV.F90	2012-10-16 13:22:44 +0000
@@ -46,8 +46,8 @@
   use boundary_conditions
   use boundary_conditions_from_options
   use divergence_matrix_cv, only: assemble_divergence_matrix_cv
-  use global_parameters, only: OPTION_PATH_LEN
-  use fefields, only: compute_lumped_mass
+  use global_parameters, only: OPTION_PATH_LEN, FIELD_NAME_LEN
+  use fefields, only: compute_cv_mass
   use petsc_solve_state_module
   use transform_elements, only: transform_cvsurf_to_physical, &
                                 transform_cvsurf_facet_to_physical
@@ -55,6 +55,7 @@
   use halos
   use field_options
   use state_fields_module
+  use porous_media
 
   implicit none
 
@@ -73,8 +74,12 @@
   logical :: move_mesh = .false.
   ! are we including density?
   logical :: include_density = .false.
+  ! is the Density field prognostic?
+  logical :: prognostic_density = .false.
   ! are we including a souce?
   logical :: include_source = .false.
+  ! Add source directly to the right hand side?
+  logical :: add_src_directly_to_rhs = .false.
   ! are we including an absorption?
   logical :: include_absorption = .false.
   ! are we including diffusion?
@@ -127,15 +132,20 @@
       type(csr_sparsity), pointer :: mesh_sparsity_1, mesh_sparsity, &
                                      mesh_sparsity_x, grad_m_t_sparsity
 
-      ! Right hand side vector, lumped mass matrix, 
+      ! Right hand side vector, cv mass matrix, 
       ! locally iterated field (for advection iterations) 
       ! and local old field (for subcycling)
-      type(scalar_field), pointer :: t_lumpedmass, q_lumpedmass
-      type(scalar_field) :: t_lumpedmass_old, t_lumpedmass_new
-      type(scalar_field) :: rhs, lumpedmass, advit_tfield, l_old_tfield
+      type(scalar_field), pointer :: t_cvmass, q_cvmass, t_abs_src_cvmass
+      type(scalar_field) :: t_cvmass_old, t_cvmass_new
+      type(scalar_field) :: rhs, cvmass, advit_tfield, l_old_tfield
       ! Diffusion contribution to rhs
       type(scalar_field) :: diff_rhs
-
+      
+      ! Porosity field
+      type(scalar_field) :: porosity_theta 
+      type(scalar_field), target :: t_cvmass_with_porosity
+      logical :: include_porosity
+              
       ! local copy of option_path for solution field
       character(len=OPTION_PATH_LEN) :: option_path, tdensity_option_path
 
@@ -167,15 +177,15 @@
       type(cv_options_type) :: tfield_options, tdensity_options
 
       ! a dummy density in case we're solving for Advection
-      type(scalar_field), pointer :: dummyscalar
+      type(scalar_field), pointer :: dummyscalar, dummydensity
       type(vector_field), pointer :: dummyvector
       type(tensor_field), pointer :: dummytensor
       ! somewhere to put strings temporarily
       character(len=FIELD_NAME_LEN) :: tmpstring
       ! what equation type are we solving for?
       integer :: equation_type
-      ! success indicator and iterators
-      integer :: stat, i, j
+      ! success indicator
+      integer :: stat
       ! the courant number field
       type(scalar_field) :: cfl_no
       ! nonlinear and grid velocities
@@ -186,6 +196,9 @@
       type(scalar_field), pointer :: sink
       !! Direction of gravity
       type(vector_field), pointer :: gravity
+      !! Temporary pointer to the material_phase's Velocity field
+      type(vector_field), pointer :: temp_velocity_ptr
+      character(len=FIELD_NAME_LEN) :: velocity_equation_type
 
       ! assume explicitness?
       logical :: explicit
@@ -233,6 +246,13 @@
       call allocate(dummytensor, tfield%mesh, name="DummyTensor", field_type=FIELD_TYPE_CONSTANT)
       call zero(dummytensor)
       dummytensor%option_path = " "
+      
+      ! Allocate a dummy density field set to 1.0 in case we are using KEpsilon, 
+      ! but have a Boussinesq or Drainage equation_type for Velocity.
+      allocate(dummydensity)
+      call allocate(dummydensity, tfield%mesh, name="DummyDensity", field_type=FIELD_TYPE_CONSTANT)
+      call set(dummydensity, 1.0)
+      dummydensity%option_path = " "
 
       ! find out equation type and hence if density is needed or not
       equation_type=equation_type_index(trim(option_path))
@@ -244,6 +264,30 @@
         tdensity=>dummyscalar
         oldtdensity=>dummyscalar
 
+      case(FIELD_EQUATION_KEPSILON)
+        ! Depending on the equation type, extract the density or set it to some dummy field allocated above
+        temp_velocity_ptr => extract_vector_field(state(1), "Velocity")
+        call get_option(trim(temp_velocity_ptr%option_path)//"/prognostic/equation[0]/name", velocity_equation_type)
+        select case(velocity_equation_type)
+           case("LinearMomentum")
+              include_density = .true.
+              tdensity=>extract_scalar_field(state(1), "Density")
+              oldtdensity=>extract_scalar_field(state(1), "OldDensity")
+
+              if(have_option(trim(tdensity%option_path)//"/prognostic")) then
+                 prognostic_density = .true.
+              end if
+           case("Boussinesq")
+              tdensity => dummydensity
+              oldtdensity => dummydensity
+           case("Drainage")
+              tdensity => dummydensity
+              oldtdensity => dummydensity
+           case default
+              ! developer error... out of sync options input and code
+              FLAbort("Unknown equation type for velocity")
+        end select
+
       case(FIELD_EQUATION_CONSERVATIONOFMASS, FIELD_EQUATION_REDUCEDCONSERVATIONOFMASS, &
            FIELD_EQUATION_INTERNALENERGY, FIELD_EQUATION_HEATTRANSFER )
         call get_option(trim(option_path)//'/prognostic/equation[0]/density[0]/name', &
@@ -254,17 +298,28 @@
         ! careful with priority ordering
         tdensity=>extract_scalar_field(state(1), trim(tmpstring))
         ewrite_minmax(tdensity)
+
         ! halo exchange? - not currently necessary when suboptimal halo exchange if density
         ! is solved for with this subroutine and the correct priority ordering.
         oldtdensity=>extract_scalar_field(state(1), "Old"//trim(tmpstring))
         ewrite_minmax(oldtdensity)
+           
+        if(have_option(trim(tdensity%option_path)//"/prognostic")) then
+           prognostic_density = .true.
+        end if
       end select
 
       ! get the density option path
-      if(have_option(trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options')) then
-        tdensity_option_path = trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options'
-      else
-        tdensity_option_path = trim(tdensity%option_path)
+      if(include_density) then
+         if(prognostic_density) then
+            if(have_option(trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options')) then
+               tdensity_option_path = trim(option_path)//'/prognostic/equation[0]/density[0]/discretisation_options'
+            else
+               tdensity_option_path = trim(tdensity%option_path)
+            end if
+         else
+            tdensity_option_path = trim(tfield%option_path)
+         end if
       end if
 
       ! now we can get the options for these fields
@@ -333,6 +388,13 @@
       if(.not.include_source) then
         source=>dummyscalar
       else
+        add_src_directly_to_rhs = have_option(trim(source%option_path)//'/diagnostic/add_directly_to_rhs')
+      
+        if (add_src_directly_to_rhs) then 
+          ewrite(2, *) "Adding Source field directly to the right hand side"
+          assert(node_count(source) == node_count(tfield))
+        end if
+
         ewrite_minmax(source)
       end if
       
@@ -344,7 +406,7 @@
       else
         ewrite_minmax(absorption)
       end if
-
+      
       ! create control volume shape functions
       call get_option("/geometry/quadrature/controlvolume_surface_degree", &
                      quaddegree, default=1)
@@ -470,9 +532,9 @@
           FLExit("Can't be explicit and exclude the mass terms.")
         end if
         
-        ! allocate a local lumpedmass field because it will get modified by bc etc.
-        call allocate(lumpedmass, tfield%mesh, name=trim(field_name)//"LumpedMass")
-        call zero(lumpedmass)
+        ! allocate a local cvmass field because it will get modified by bc etc.
+        call allocate(cvmass, tfield%mesh, name=trim(field_name)//"CVMass")
+        call zero(cvmass)
       end if
 
       if(include_diffusion) then
@@ -485,35 +547,54 @@
 
       ! allocate the rhs of the equation
       call allocate(rhs, tfield%mesh, name=trim(field_name)//"RHS")
-
-      if(tfield%mesh%shape%degree>1) then
-        ! try lumping on the submesh
-        t_lumpedmass => get_lumped_mass_on_submesh(state, tfield%mesh)
-      else
-        ! then find the lumped mass
-        ! (this replaces hart2/3d etc. in old code!!)
-        t_lumpedmass => get_lumped_mass(state, tfield%mesh)
-      end if
-      ewrite_minmax(t_lumpedmass)
+      
+      ! are we including a porosity coefficient on the time term?
+      if (have_option(trim(complete_field_path(tfield%option_path))//'/porosity')) then
+         include_porosity = .true.
+
+         ! get the porosity theta averaged field - this will allocate it
+         call form_porosity_theta(porosity_theta, state(1), &
+             &option_path = trim(complete_field_path(tfield%option_path))//'/porosity')       
+                  
+         call allocate(t_cvmass_with_porosity, tfield%mesh, name="CVMassWithPorosity")
+         call compute_cv_mass(x, t_cvmass_with_porosity, porosity_theta)
+         ewrite_minmax(t_cvmass_with_porosity)
+         
+         call deallocate(porosity_theta)
+      else
+         include_porosity = .false.
+      end if
+
+      ! find the cv mass that is used for the absorption and source terms
+      t_abs_src_cvmass => get_cv_mass(state, tfield%mesh)
+      ewrite_minmax(t_abs_src_cvmass)
+      
+      ! find the cv mass that is used for the time term derivative
+      if (include_porosity) then
+         t_cvmass => t_cvmass_with_porosity
+      else
+         t_cvmass => t_abs_src_cvmass      
+      end if
+      ewrite_minmax(t_cvmass)
 
       move_mesh = have_option("/mesh_adaptivity/mesh_movement")
       if(move_mesh) then
         if(.not.include_advection) then
           FLExit("Moving the mesh but not including advection is not possible yet.")
         end if
+        if (include_porosity) then
+           FLExit("Moving mesh not set up to work when including porosity")
+        end if
         ewrite(2,*) "Moving mesh."
         x_old=>extract_vector_field(state(1), "OldCoordinate")
         x_new=>extract_vector_field(state(1), "IteratedCoordinate")
-        call allocate(t_lumpedmass_old, tfield%mesh, name=trim(field_name)//"OldLumpedMass")
-        call allocate(t_lumpedmass_new, tfield%mesh, name=trim(field_name)//"NewLumpedMass")
-        if(tfield%mesh%shape%degree>1) then
-          FLExit("Lumping on submesh while moving the mesh not set up.")
-        else
-          call compute_lumped_mass(x_old, t_lumpedmass_old)
-          call compute_lumped_mass(x_new, t_lumpedmass_new)
-        end if
-        ewrite_minmax(t_lumpedmass_old)
-        ewrite_minmax(t_lumpedmass_new)
+        call allocate(t_cvmass_old, tfield%mesh, name=trim(field_name)//"OldCVMass")
+        call allocate(t_cvmass_new, tfield%mesh, name=trim(field_name)//"NewCVMass")
+
+        call compute_cv_mass(x_old, t_cvmass_old)
+        call compute_cv_mass(x_new, t_cvmass_new)
+        ewrite_minmax(t_cvmass_old)
+        ewrite_minmax(t_cvmass_new)
         
         ug=>extract_vector_field(state(1), "GridVelocity")
         ewrite_minmax(ug)
@@ -531,12 +612,12 @@
 
       if(include_diffusion.and.(tfield_options%diffusionscheme==CV_DIFFUSION_BASSIREBAY)) then
         if(.not.(tfield%mesh==diffusivity%mesh)) then
-          q_lumpedmass => get_lumped_mass(state, diffusivity%mesh)
+          q_cvmass => get_cv_mass(state, diffusivity%mesh)
         else
-          q_lumpedmass => t_lumpedmass
+          q_cvmass => t_cvmass
         end if
       else
-        q_lumpedmass => t_lumpedmass
+        q_cvmass => t_cvmass
       end if
 
       ! allocate a field to store the locally iterated values in
@@ -576,6 +657,9 @@
           ! yes, we're subcycling
           ! we should have already calculated the courant number (or aborted in the attempt)
           no_subcycles=ceiling(max_cfl/max_sub_cfl)
+          if(include_diffusion.or.include_source.or.include_absorption) then
+            no_subcycles = max(no_subcycles, 1)
+          end if
           if(no_subcycles>1) then
             sub_dt=dt/real(no_subcycles)
             call scale(cfl_no, 1.0/real(no_subcycles))
@@ -626,25 +710,25 @@
                                         t_cvshape_full, t_cvbdyshape_full, &
                                         diff_cvshape_full, diff_cvbdyshape_full, &
                                         state, advu, ug, x, x_tfield, cfl_no, sub_dt, &
-                                        diffusivity, q_lumpedmass, &
+                                        diffusivity, q_cvmass, &
                                         mesh_sparsity_x, grad_m_t_sparsity)
           end if
 
           ! assemble it all into a coherent equation
-          call assemble_field_eqn_cv(M, A_m, lumpedmass, rhs, &
+          call assemble_field_eqn_cv(M, A_m, cvmass, rhs, &
                                     tfield, l_old_tfield, &
                                     tdensity, oldtdensity, tdensity_options, &
                                     source, absorption, tfield_options%theta, &
                                     state, advu, sub_dt, explicit, &
-                                    t_lumpedmass, t_lumpedmass_old, t_lumpedmass_new, & 
+                                    t_cvmass, t_abs_src_cvmass, t_cvmass_old, t_cvmass_new, & 
                                     D_m, diff_rhs)
 
 
           ! Solve for the change in tfield.
           if(explicit) then
-            call apply_dirichlet_conditions(lumpedmass, rhs, tfield, sub_dt)
+            call apply_dirichlet_conditions(cvmass, rhs, tfield, sub_dt)
 
-            delta_tfield%val = rhs%val/lumpedmass%val
+            delta_tfield%val = rhs%val/cvmass%val
           else
             ! apply strong dirichlet boundary conditions (if any)
             ! note that weak conditions (known as control volume boundary conditions)
@@ -662,7 +746,7 @@
 
           call halo_update(tfield)  ! exchange the extended halos
 
-          call test_and_write_advection_convergence(tfield, advit_tfield, x, &
+          call test_and_write_advection_convergence(tfield, advit_tfield, x, t_cvmass, &
                                     filename=trim(state(1)%name)//"__"//trim(tfield%name), &
                                     time=time+sub_dt, dt=sub_dt, it=global_it, adv_it=adv_it, &
                                     subcyc=sub, error=error)
@@ -682,7 +766,7 @@
       call deallocate(rhs)
       if(.not.explicit) call deallocate(A_m)
       if(.not.explicit) call deallocate(M)
-      if(explicit) call deallocate(lumpedmass)
+      if(explicit) call deallocate(cvmass)
       call deallocate(cfl_no)
       call deallocate(x_cvbdyshape)
       call deallocate(x_cvbdyshape_full)
@@ -706,29 +790,34 @@
       deallocate(dummyvector)
       call deallocate(dummytensor)
       deallocate(dummytensor)
+      call deallocate(dummydensity)
+      deallocate(dummydensity)
       if (include_diffusion) then
         call deallocate(D_m)
         call deallocate(diff_rhs)
       end if
       call deallocate(x_tfield)
       if(move_mesh) then
-        call deallocate(t_lumpedmass_new)
-        call deallocate(t_lumpedmass_old)
+        call deallocate(t_cvmass_new)
+        call deallocate(t_cvmass_old)
       end if
       call deallocate(ug_cvshape)
       call deallocate(ug_cvbdyshape)
+      if (include_porosity) then
+        call deallocate(t_cvmass_with_porosity)
+      end if
 
     end subroutine solve_field_eqn_cv
     ! end of solution wrapping subroutines
     !************************************************************************
     !************************************************************************
     ! equation wrapping subroutines
-    subroutine assemble_field_eqn_cv(M, A_m, m_lumpedmass, rhs, &
+    subroutine assemble_field_eqn_cv(M, A_m, m_cvmass, rhs, &
                                     tfield, oldtfield, &
                                     tdensity, oldtdensity, tdensity_options, &
                                     source, absorption, theta, &
                                     state, advu, dt, explicit, &
-                                    lumpedmass, lumpedmass_old, lumpedmass_new, &
+                                    cvmass, abs_src_cvmass, cvmass_old, cvmass_new, &
                                     D_m, diff_rhs)
 
       ! This subroutine assembles the equation
@@ -744,8 +833,8 @@
       ! matrix containing advective terms - to be incorporated
       ! into M during this subroutine
       type(csr_matrix), intent(inout) :: A_m
-      ! rhs of equation
-      type(scalar_field), intent(inout) :: m_lumpedmass, rhs
+      ! explicit lhs and rhs of equation
+      type(scalar_field), intent(inout) :: m_cvmass, rhs
       ! the field we are solving for
       type(scalar_field), intent(inout) :: tfield
       type(scalar_field), intent(inout) :: oldtfield, tdensity, oldtdensity
@@ -762,10 +851,14 @@
       real, intent(in) :: dt
       ! are we assuming this is a fully explicit equation?
       logical, intent(in) :: explicit
+      ! cv mass to use for time derivative term
+      type(scalar_field), intent(in) :: cvmass
       ! moving mesh stuff
-      type(scalar_field), intent(in) :: lumpedmass
-      type(scalar_field), intent(in) :: lumpedmass_old
-      type(scalar_field), intent(in) :: lumpedmass_new
+      type(scalar_field), intent(in) :: cvmass_old
+      type(scalar_field), intent(in) :: cvmass_new
+      ! cv mass to use for absorption and source
+      type(scalar_field), intent(in) :: abs_src_cvmass      
+      
       ! diffusion:
       type(csr_matrix), intent(inout), optional :: D_m
       type(scalar_field), intent(inout), optional :: diff_rhs
@@ -809,7 +902,7 @@
 
       ! zero the "matrices" being assembled
       if(explicit) then
-        call zero(m_lumpedmass)
+        call zero(m_cvmass)
       else
         call zero(M)
       end if
@@ -818,36 +911,36 @@
       if(include_mass) then
         if(move_mesh) then
           if(explicit) then
-            call set(m_lumpedmass, lumpedmass_new)
+            call set(m_cvmass, cvmass_new)
           else
-            call addto_diag(M, lumpedmass_new)
+            call addto_diag(M, cvmass_new)
           end if
         else
           if(explicit) then
-            call set(m_lumpedmass, lumpedmass)
+            call set(m_cvmass, cvmass)
           else
-            call addto_diag(M, lumpedmass)
+            call addto_diag(M, cvmass)
           end if
         end if
       end if
 
       ! allocate some memory for assembly
       call allocate(MT_old, rhs%mesh, name="MT_oldProduct" )
-      if(include_source) then
+      if(include_source .and. (.not. add_src_directly_to_rhs)) then
         call allocate(masssource, rhs%mesh, name="MassSourceProduct" )
-        call set(masssource, lumpedmass)
+        call set(masssource, abs_src_cvmass)
         call scale(masssource, source)
       end if
       if(include_absorption) then
         call allocate(massabsorption, rhs%mesh, name="MassAbsorptionProduct" )
-        call set(massabsorption, lumpedmass)
+        call set(massabsorption, abs_src_cvmass)
         call scale(massabsorption, absorption)
       end if
       
       if(move_mesh) then
         call allocate(massconservation, rhs%mesh, name="MovingMeshMassConservation")
-        call set(massconservation, lumpedmass_old)
-        call addto(massconservation, lumpedmass_new, scale=-1.0)
+        call set(massconservation, cvmass_old)
+        call addto(massconservation, cvmass_new, scale=-1.0)
         call scale(massconservation, 1./dt)
         call scale(massconservation, oldtfield)
       end if
@@ -875,7 +968,7 @@
         end if
 
 
-        if(include_source) call addto(rhs, masssource)
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
 
         if(include_absorption) then
           ! massabsorption has already been added to the matrix so it can now be scaled
@@ -909,7 +1002,7 @@
           call addto(consterm, oldtdensity, -1.0)
           call scale(consterm, oldtfield)
           call scale(consterm, 1./dt)
-          call scale(consterm, lumpedmass)
+          call scale(consterm, cvmass)
 
           call addto(rhs, consterm, -1.0)
 
@@ -920,7 +1013,7 @@
         ! multiply the diagonal of M by the up to date density
         if(explicit) then
           if(include_mass) then
-            call scale(m_lumpedmass, tdensity)
+            call scale(m_cvmass, tdensity)
           end if
         else
           if(include_mass) then
@@ -936,7 +1029,7 @@
           end if
         end if
         
-        if(include_source) call addto(rhs, masssource)
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
 
         if(include_absorption) then
           ! massabsorption has already been added to the matrix so it can now be scaled
@@ -968,7 +1061,7 @@
         ! multiply the diagonal by the previous timesteps density
         if(explicit) then
           if(include_mass) then
-            m_lumpedmass%val = m_lumpedmass%val*(tdensity_theta*tdensity%val+(1.0-tdensity_theta)*oldtdensity%val)
+            m_cvmass%val = m_cvmass%val*(tdensity_theta*tdensity%val+(1.0-tdensity_theta)*oldtdensity%val)
           end if
         else
           if(include_mass) then 
@@ -984,7 +1077,7 @@
           end if
         end if
 
-        if(include_source) call addto(rhs, masssource)
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
 
         if(include_absorption) then
           ! massabsorption has already been added to the matrix so it can now be scaled
@@ -1016,7 +1109,7 @@
         ! multiply the diagonal by the previous timesteps density
         if(explicit) then
           if(include_mass) then
-            m_lumpedmass%val = m_lumpedmass%val*(tdensity_theta*tdensity%val+(1.0-tdensity_theta)*oldtdensity%val)
+            m_cvmass%val = m_cvmass%val*(tdensity_theta*tdensity%val+(1.0-tdensity_theta)*oldtdensity%val)
           end if
         else
           if(include_mass) then 
@@ -1032,7 +1125,7 @@
           end if
         end if
 
-        if(include_source) call addto(rhs, masssource)
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
 
         if(include_absorption) then
           ! massabsorption has already been added to the matrix so it can now be scaled
@@ -1087,22 +1180,67 @@
         ! construct M
         if(explicit) then
           if(include_mass) then
-            call scale(m_lumpedmass, tdensity)
-          end if
-        else
-          if(include_mass) then
-            call mult_diag(M, tdensity)
-          end if
-          if(include_advection) call addto(M, A_m, dt)
-          if(include_absorption) call addto_diag(M, massabsorption, theta*dt)
-        
-          if(include_advection) then
-            call mult(MT_old, A_m, oldtfield)
-            call addto(rhs, MT_old, -1.0)
-          end if
-        end if
-
-        if(include_source) call addto(rhs, masssource)
+            call scale(m_cvmass, tdensity)
+          end if
+        else
+          if(include_mass) then
+            call mult_diag(M, tdensity)
+          end if
+          if(include_advection) call addto(M, A_m, dt)
+          if(include_absorption) call addto_diag(M, massabsorption, theta*dt)
+        
+          if(include_advection) then
+            call mult(MT_old, A_m, oldtfield)
+            call addto(rhs, MT_old, -1.0)
+          end if
+        end if
+
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
+
+        if(include_absorption) then
+          ! massabsorption has already been added to the matrix so it can now be scaled
+          ! by the old field value to add it to the rhs
+          call scale(massabsorption, oldtfield)
+          call addto(rhs, massabsorption, -1.0)
+        end if
+
+        if(include_diffusion) then
+          call mult(MT_old, D_m, oldtfield)
+          call addto(rhs, MT_old, -1.0)
+          call addto(rhs, diff_rhs, -1.0)
+
+          if(.not.explicit) then
+            call addto(M, D_m, theta*dt)
+          end if
+        end if
+        
+        if(move_mesh) then
+          FLExit("Moving mesh with this equation type not yet supported.")
+        end if
+
+      case (FIELD_EQUATION_KEPSILON)
+
+        ! [\rho^{n+1}M + dt*A_m + dt*theta*D_m](T^{n+1}-T^{n})/dt = rhs - [A_m + D_m]*T^{n} - diff_rhs
+
+        ! construct M
+        if(explicit) then
+          if(include_mass) then
+            call scale(m_cvmass, tdensity)
+          end if
+        else
+          if(include_mass) then
+            call mult_diag(M, tdensity)
+          end if
+          if(include_advection) call addto(M, A_m, dt)
+          if(include_absorption) call addto_diag(M, massabsorption, theta*dt)
+        
+          if(include_advection) then
+            call mult(MT_old, A_m, oldtfield)
+            call addto(rhs, MT_old, -1.0)
+          end if
+        end if
+
+        if(include_source .and. (.not. add_src_directly_to_rhs)) call addto(rhs, masssource)
 
         if(include_absorption) then
           ! massabsorption has already been added to the matrix so it can now be scaled
@@ -1127,7 +1265,11 @@
 
       end select
 
-      if(include_source) call deallocate(masssource)
+      ! Add the source directly to the rhs if required 
+      ! which must be included before dirichlet BC's.
+      if (add_src_directly_to_rhs) call addto(rhs, source)
+
+      if(include_source .and. (.not. add_src_directly_to_rhs)) call deallocate(masssource)
       if(include_absorption) call deallocate(massabsorption)
       call deallocate(MT_old)
       if(move_mesh) call deallocate(massconservation)
@@ -1149,7 +1291,7 @@
                                        t_cvshape_full, t_cvbdyshape_full, &
                                        diff_cvshape_full, diff_cvbdyshape_full, &
                                        state, advu, ug, x, x_tfield, cfl_no, dt, &
-                                       diffusivity, q_lumpedmass, &
+                                       diffusivity, q_cvmass, &
                                        mesh_sparsity, grad_m_t_sparsity)
 
       ! This subroutine assembles the advection and diffusion matrices and rhs(s) for
@@ -1210,8 +1352,8 @@
 
       ! the diffusivity tensor
       type(tensor_field), intent(in) :: diffusivity
-      ! the lumped mass = the mass matrix for the auxilliary diffusion equation
-      type(scalar_field), intent(in) :: q_lumpedmass
+      ! the cv mass = the mass matrix for the auxilliary diffusion equation
+      type(scalar_field), intent(in) :: q_cvmass
       ! sparsity pattern for the gradient transposed operator
       type(csr_sparsity), intent(inout) :: grad_m_t_sparsity
 
@@ -1219,6 +1361,7 @@
       ! memory for coordinates, velocity, normals, determinants, nodes
       ! and the cfl number at the gauss pts and nodes
       real, dimension(x%dim,ele_loc(x,1)) :: x_ele
+      real, dimension(x%dim,ele_loc(tfield,1)) :: xt_ele
       real, dimension(x%dim,face_loc(x,1)) :: x_ele_bdy
       real, dimension(x%dim,x_cvshape%ngi) :: x_f
       real, dimension(advu%dim,u_cvshape%ngi) :: u_f
@@ -1268,7 +1411,7 @@
       ! the type of the bc if integrating over domain boundaries
       integer, dimension(:), allocatable :: tfield_bc_type, tdensity_bc_type
       ! fields for the bcs over the entire surface mesh
-      type(scalar_field) :: tfield_bc, tdensity_bc
+      type(scalar_field) :: tfield_bc, tfield_bc2, tdensity_bc
 
       ! local element matrices - allow the assembly of an entire face without multiple calls to csr_pos
       real, dimension(mesh_dim(tfield), ele_loc(tfield,1), ele_loc(tfield,1)) :: grad_mat_local
@@ -1277,6 +1420,11 @@
       real, dimension(mesh_dim(tfield), face_loc(tfield,1)) :: grad_mat_local_bdy, grad_rhs_local_bdy
       real, dimension(face_loc(tfield,1)) :: mat_local_bdy, rhs_local_bdy, div_rhs_local_bdy
       real, dimension(ele_loc(tfield,1)) :: rhs_local
+      
+      ! Data associated with Robin BC's.
+      real, dimension(face_loc(tfield,1)) :: robin_diff_mat_local_bdy
+      real, dimension(face_loc(tfield,1)) :: robin_diff_rhs_local_bdy
+      real, dimension(face_loc(tfield,1)) :: robin_bc_val1_ele_bdy, robin_bc_val2_ele_bdy
 
       ! the auxilliary gradient matrix (assembled as a divergence confusingly)
       type(block_csr_matrix) :: div_m
@@ -1293,7 +1441,7 @@
 
       ! Boundary condition types
       integer, parameter :: BC_TYPE_WEAKDIRICHLET = 1, BC_TYPE_NEUMANN = 2, BC_TYPE_INTERNAL = 3, &
-                            BC_TYPE_ZEROFLUX = 4, BC_TYPE_FLUX = 5
+                            BC_TYPE_ZEROFLUX = 4, BC_TYPE_FLUX = 5, BC_TYPE_ROBIN = 6
 
       ewrite(1, *) "In assemble_advectiondiffusion_m_cv"
 
@@ -1365,6 +1513,7 @@
       ! loop over elements
       element_loop: do ele=1, element_count(tfield)
         x_ele=ele_val(x, ele)
+        xt_ele=ele_val(x_tfield, ele)
         x_f=ele_val_at_quad(x, ele, x_cvshape)
         nodes=>ele_nodes(tfield, ele)
         ! the nodes in this element from the coordinate mesh projected
@@ -1477,7 +1626,7 @@
                                         tfield_face_val, &
                                         oldtfield_face_val, &
                                         tfield_options%theta, dt, udotn, &
-                                        x_ele, tfield_options%limit_theta, &
+                                        xt_ele, tfield_options%limit_theta, &
                                         tfield_ele, oldtfield_ele, &
                                         ftheta=ftheta)
 
@@ -1505,7 +1654,7 @@
                                           tdensity_face_val, &
                                           oldtdensity_face_val, &
                                           tdensity_options%theta, dt, udotn, &
-                                          x_ele, tdensity_options%limit_theta, &
+                                          xt_ele, tdensity_options%limit_theta, &
                                           tdensity_ele, oldtdensity_ele)
 
                       if(assemble_advection_matrix) then
@@ -1637,7 +1786,8 @@
         "neumann      ", &
         "internal     ", &
         "zero_flux    ", &
-        "flux         "/), tfield_bc, tfield_bc_type)
+        "flux         ", &
+        "robin        "/), tfield_bc, tfield_bc_type, boundary_second_value = tfield_bc2)
       if(include_density) then
         allocate(tdensity_bc_type(surface_element_count(tdensity)))
         call get_entire_boundary_condition(tdensity, (/"weakdirichlet"/), tdensity_bc, tdensity_bc_type)
@@ -1719,7 +1869,17 @@
         end if
 
         if(assemble_diffusion) then
-          ghost_gradtfield_ele_bdy = ele_val(tfield_bc, sele)
+          if (tfield_bc_type(sele)==BC_TYPE_ROBIN) then
+             if (.not. tfield_options%diffusionscheme == CV_DIFFUSION_ELEMENTGRADIENT) then
+                FLExit('Can only use CV robin BC with ElementGradient diffusion scheme')
+             end if
+             robin_bc_val1_ele_bdy    = ele_val(tfield_bc, sele)
+             robin_bc_val2_ele_bdy    = ele_val(tfield_bc2, sele)
+             robin_diff_mat_local_bdy = 0.0
+             robin_diff_rhs_local_bdy = 0.0
+          else
+             ghost_gradtfield_ele_bdy = ele_val(tfield_bc, sele)
+           end if
         end if
 
         ! zero small matrices for assembly
@@ -1869,6 +2029,18 @@
                       div_rhs_local_bdy(iloc) = div_rhs_local_bdy(iloc) &
                                   -detwei_bdy(ggi)*ghost_gradtfield_ele_bdy(iloc)
 
+                    else if (tfield_bc_type(sele)==BC_TYPE_ROBIN) then
+                        
+                        ! Add first coeff term to rhs
+                        robin_diff_rhs_local_bdy(iloc) = robin_diff_rhs_local_bdy(iloc) - &
+                                                         detwei_bdy(ggi) * &
+                                                         robin_bc_val1_ele_bdy(iloc)
+                                                
+                        ! Add implicit second coeff * tfield term to matrix
+                        robin_diff_mat_local_bdy(iloc) = robin_diff_mat_local_bdy(iloc) + &
+                                                         detwei_bdy(ggi) * &
+                                                         robin_bc_val2_ele_bdy(iloc)
+                    
                     else
 
                       ! because transform to physical doesn't give the full gradient at a face
@@ -1927,6 +2099,12 @@
 
               call addto(diff_rhs, nodes_bdy, div_rhs_local_bdy)
 
+            else if (tfield_bc_type(sele)==BC_TYPE_ROBIN) then
+              
+               call addto(diff_rhs, nodes_bdy, robin_diff_rhs_local_bdy)
+               
+               call addto_diag(D_m, nodes_bdy, robin_diff_mat_local_bdy)
+              
             else
 
             !    assume zero neumann for the moment
@@ -1948,21 +2126,23 @@
 
       if(assemble_diffusion.and.(tfield_options%diffusionscheme==CV_DIFFUSION_BASSIREBAY)) then
 
-        ! assemble div_m and q_lumpedmass into the final D_m
+        ! assemble div_m and q_cvmass into the final D_m
         call assemble_bassirebay_diffusion_m_cv(D_m, diff_rhs, &
                                      div_m, grad_rhs, &
-                                     diffusivity, q_lumpedmass)
+                                     diffusivity, q_cvmass)
         ! ElementGradient assembles D_m directly so no need for a step like this
 
       end if
 
       deallocate(tfield_bc_type)
       call deallocate(tfield_bc)
+      call deallocate(tfield_bc2)
       call deallocate(tfield_upwind)
       call deallocate(oldtfield_upwind)
       if(include_density) then
         deallocate(tdensity_bc_type)
         call deallocate(tdensity_bc)
+
         call deallocate(tdensity_upwind)
         call deallocate(oldtdensity_upwind)
       end if
@@ -1978,7 +2158,7 @@
 
     subroutine assemble_bassirebay_diffusion_m_cv(D_m, diff_rhs, &
                                        div_m, grad_rhs, &
-                                       diffusivity, q_lumpedmass)
+                                       diffusivity, q_cvmass)
 
       type(csr_matrix), intent(inout) :: D_m
       type(scalar_field), intent(inout) :: diff_rhs
@@ -1987,7 +2167,7 @@
       type(vector_field), intent(inout) :: grad_rhs
 
       type(tensor_field), intent(in) :: diffusivity
-      type(scalar_field), intent(in) :: q_lumpedmass
+      type(scalar_field), intent(in) :: q_cvmass
 
       logical :: isotropic
 
@@ -1996,10 +2176,10 @@
       ! an optimisation that reduces the number of matrix multiplies if we're isotropic
       isotropic=isotropic_field(diffusivity)
 
-      call mult_div_tensorinvscalar_div_T(D_m, div_m, diffusivity, q_lumpedmass, div_m, &
+      call mult_div_tensorinvscalar_div_T(D_m, div_m, diffusivity, q_cvmass, div_m, &
                                           isotropic)
 
-      call mult_div_tensorinvscalar_vector(diff_rhs, div_m, diffusivity, q_lumpedmass, grad_rhs, &
+      call mult_div_tensorinvscalar_vector(diff_rhs, div_m, diffusivity, q_cvmass, grad_rhs, &
                                            isotropic)
 
     end subroutine assemble_bassirebay_diffusion_m_cv
@@ -2115,14 +2295,20 @@
       ! sparsity structure to construct the matrices with
       type(csr_sparsity), pointer :: mesh_sparsity, mesh_sparsity_x
 
-      ! Right hand side vector, lumped mass matrix, 
+      ! Right hand side vector, cv mass matrix, 
       ! locally iterated field (for advection iterations) 
       ! and local old field (for subcycling)
       type(scalar_field), dimension(nfields) :: rhs, advit_tfield
       type(scalar_field_pointer), dimension(nfields) :: l_old_tfield
-      type(scalar_field), dimension(nfields) :: lumpedmass
-      type(scalar_field), pointer :: t_lumpedmass
-      type(scalar_field) :: t_lumpedmass_old, t_lumpedmass_new
+      type(scalar_field), dimension(nfields) :: cvmass
+      type(scalar_field), pointer :: t_abs_src_cvmass
+      type(scalar_field_pointer), dimension(nfields) :: t_cvmass
+      type(scalar_field) :: t_cvmass_old, t_cvmass_new
+      
+      ! Porosity field
+      type(scalar_field) :: porosity_theta 
+      type(scalar_field), dimension(nfields), target :: t_cvmass_with_porosity
+      logical, dimension(nfields) :: include_porosity
 
       ! local copy of option_path for solution field
       character(len=OPTION_PATH_LEN), dimension(nfields) :: option_path
@@ -2281,6 +2467,9 @@
         source(f)%ptr=>extract_scalar_field(state(state_indices(f)), trim(field_name)//"Source", stat=stat)
         if(stat==0) then
           FLExit("Coupled CV broken with Sources")
+          ! If Coupled CV is ever fixed to work with Sources then the 
+          ! option source(f)%option_path//'/diagnostic/add_directly_to_rhs'
+          ! should be accounted for.
         end if
         if(stat/=0) source(f)%ptr=>dummyscalar
         absorption(f)%ptr=>extract_scalar_field(state(state_indices(f)), trim(field_name)//"Absorption", stat=stat)
@@ -2358,23 +2547,40 @@
           call allocate(A_m(f), mesh_sparsity, name=trim(field_name)//int2str(f)//"AdvectionMatrix")
           call zero(A_m(f))
         else
-          call allocate(lumpedmass(f), tfield(1)%ptr%mesh, name=trim(field_name)//"LocalLumpedMass")
-          call zero(lumpedmass(f)) 
+          call allocate(cvmass(f), tfield(1)%ptr%mesh, name=trim(field_name)//"LocalCVMass")
+          call zero(cvmass(f)) 
         end if
 
         ! allocate the rhs of the equation
         call allocate(rhs(f), tfield(f)%ptr%mesh, name=trim(field_name)//int2str(f)//"RHS")
       end do
 
-      if(element_degree(tfield(1)%ptr, 1)>1) then
-        ! try lumping on the submesh
-        t_lumpedmass => get_lumped_mass_on_submesh(state, tfield(1)%ptr%mesh)
-      else
-        ! then find the lumped mass
-        ! (this replaces hart2/3d etc. in old code!!)
-        t_lumpedmass => get_lumped_mass(state, tfield(1)%ptr%mesh)
-      end if
-      ewrite_minmax(t_lumpedmass)
+      ! find the cv mass that is used for the absorption and source terms
+      t_abs_src_cvmass => get_cv_mass(state, tfield(1)%ptr%mesh)
+      ewrite_minmax(t_abs_src_cvmass)
+      
+      ! find the cv mass that is used for the time term derivative - which may include the porosity
+      do f = 1,nfields                  
+         if (have_option(trim(complete_field_path(tfield(f)%ptr%option_path))//'/porosity')) then            
+            include_porosity(f) = .true.
+
+            ! get the porosity theta averaged field - this will allocate it
+            call form_porosity_theta(porosity_theta, state(state_indices(f)), &
+                &option_path = trim(complete_field_path(tfield(f)%ptr%option_path))//'/porosity')       
+                     
+            call allocate(t_cvmass_with_porosity(f), tfield(f)%ptr%mesh, name="CVMassWithPorosity")
+            call compute_cv_mass(x, t_cvmass_with_porosity(f), porosity_theta)
+            
+            call deallocate(porosity_theta)
+            
+            t_cvmass(f)%ptr => t_cvmass_with_porosity(f)                    
+         else 
+            include_porosity(f) = .false.
+
+            t_cvmass(f)%ptr => t_abs_src_cvmass
+         end if
+         ewrite_minmax(t_cvmass(f)%ptr)
+      end do
 
       move_mesh = have_option("/mesh_adaptivity/mesh_movement")
       if(move_mesh) then
@@ -2382,19 +2588,19 @@
         if(.not.include_advection) then
           FLExit("Moving the mesh but not including advection is not possible yet.")
         end if
+        if (any(include_porosity)) then
+           FLExit("Moving mesh not set up to work when including porosity")
+        end if
         ewrite(2,*) "Moving mesh."
         x_old=>extract_vector_field(state(1), "OldCoordinate")
         x_new=>extract_vector_field(state(1), "IteratedCoordinate")
-        call allocate(t_lumpedmass_old, tfield(1)%ptr%mesh, name=trim(field_name)//"OldLumpedMass")
-        call allocate(t_lumpedmass_new, tfield(1)%ptr%mesh, name=trim(field_name)//"NewLumpedMass")
-        if(tfield(1)%ptr%mesh%shape%degree>1) then
-          FLExit("Lumping on submesh while moving the mesh not set up.")
-        else
-          call compute_lumped_mass(x_old, t_lumpedmass_old)
-          call compute_lumped_mass(x_new, t_lumpedmass_new)
-        end if
-        ewrite_minmax(t_lumpedmass_old)
-        ewrite_minmax(t_lumpedmass_new)
+        call allocate(t_cvmass_old, tfield(1)%ptr%mesh, name=trim(field_name)//"OldCVMass")
+        call allocate(t_cvmass_new, tfield(1)%ptr%mesh, name=trim(field_name)//"NewCVMass")
+        
+        call compute_cv_mass(x_old, t_cvmass_old)
+        call compute_cv_mass(x_new, t_cvmass_new)
+        ewrite_minmax(t_cvmass_old)
+        ewrite_minmax(t_cvmass_new)
         
         ug=>extract_vector_field(state(1), "GridVelocity")
         ewrite_minmax(ug)
@@ -2517,18 +2723,18 @@
           do f = 1, nfields
 
             ! assemble it all into a coherent equation
-            call assemble_field_eqn_cv(M(f), A_m(f), lumpedmass(f), rhs(f), &
+            call assemble_field_eqn_cv(M(f), A_m(f), cvmass(f), rhs(f), &
                                       tfield(f)%ptr, l_old_tfield(f)%ptr, &
                                       tdensity(f)%ptr, oldtdensity(f)%ptr, tdensity_options(f), &
                                       source(f)%ptr, absorption(f)%ptr, tfield_options(f)%theta, &
                                       state(state_indices(f):state_indices(f)), advu, sub_dt, explicit(f), &
-                                      t_lumpedmass, t_lumpedmass_old, t_lumpedmass_new)
+                                      t_cvmass(f)%ptr, t_abs_src_cvmass, t_cvmass_old, t_cvmass_new)
 
             ! Solve for the change in tfield.
             if(explicit(f)) then
-              call apply_dirichlet_conditions(lumpedmass(f), rhs(f), tfield(f)%ptr, sub_dt)
+              call apply_dirichlet_conditions(cvmass(f), rhs(f), tfield(f)%ptr, sub_dt)
 
-              delta_tfield(f)%val = rhs(f)%val/lumpedmass(f)%val
+              delta_tfield(f)%val = rhs(f)%val/cvmass(f)%val
             else
               ! apply strong dirichlet boundary conditions (if any)
               ! note that weak conditions (known as control volume boundary conditions)
@@ -2548,7 +2754,7 @@
 
             call halo_update(tfield(f)%ptr)  ! exchange the extended halos
 
-            call test_and_write_advection_convergence(tfield(f)%ptr, advit_tfield(f), x, &
+            call test_and_write_advection_convergence(tfield(f)%ptr, advit_tfield(f), x, t_cvmass(f)%ptr, &
                                       filename=trim(state(state_indices(f))%name)//"__"//trim(tfield(f)%ptr%name), &
                                       time=time+sub_dt, dt=sub_dt, it=global_it, adv_it=adv_it, &
                                       subcyc=sub, error=error(f))
@@ -2574,7 +2780,7 @@
         deallocate(l_old_tfield(f)%ptr)
         call deallocate(rhs(f))
         if(.not.explicit(f)) call deallocate(A_m(f))
-        if(explicit(f)) call deallocate(lumpedmass(f))
+        if(explicit(f)) call deallocate(cvmass(f))
         if(.not.explicit(f)) call deallocate(M(f))
       end do
       call deallocate(cfl_no)
@@ -2591,12 +2797,15 @@
       deallocate(dummyscalar)
       call deallocate(x_tfield)
       if(move_mesh) then
-        call deallocate(t_lumpedmass_new)
-        call deallocate(t_lumpedmass_old)
+        call deallocate(t_cvmass_new)
+        call deallocate(t_cvmass_old)
       end if
       call deallocate(ug_cvshape)
       call deallocate(ug_cvbdyshape)
-
+      do f = 1,nfields
+         if (include_porosity(f)) call deallocate(t_cvmass_with_porosity(f))
+      end do
+      
     end subroutine solve_coupled_cv
 
     ! coupled assembly:
@@ -3264,12 +3473,13 @@
 
     end subroutine initialise_advection_convergence
 
-    subroutine test_and_write_advection_convergence(field, nlfield, coordinates, filename, &
+    subroutine test_and_write_advection_convergence(field, nlfield, coordinates, cv_mass, filename, &
                                                     time, dt, it, subcyc, adv_it, &
                                                     error)
 
        type(scalar_field), intent(inout) :: field, nlfield
        type(vector_field), intent(in) :: coordinates
+       type(scalar_field), intent(in) :: cv_mass
        character(len=*), intent(in) :: filename
        real, intent(in) :: time, dt
        integer, intent(in) :: it, subcyc, adv_it
@@ -3287,7 +3497,7 @@
 
        error = 0.0
        call field_con_stats(field, nlfield, error, &
-                            convergence_norm, coordinates)
+                            convergence_norm, coordinates, cv_mass)
 
        format='(e15.6e3)'
        iformat='(i4)'

=== modified file 'assemble/Full_Projection.F90'
--- assemble/Full_Projection.F90	2011-04-06 15:37:23 +0000
+++ assemble/Full_Projection.F90	2012-10-16 13:22:44 +0000
@@ -44,37 +44,18 @@
     use petsc_solve_state_module
     use boundary_conditions_from_options
 
-#ifdef HAVE_PETSC_MODULES
 #include "petscversion.h"
+#ifdef HAVE_PETSC_MODULES
     use petsc
-#if PETSC_VERSION_MINOR==0
-    use petscvec
-    use petscmat
-    use petscksp
-    use petscpc
-#endif
 #endif
 
     implicit none
     ! Module to provide solvers, preconditioners etc... for full_projection Solver.
     ! Not this is currently tested for Full CMC solves and Stokes flow:
 #ifdef HAVE_PETSC_MODULES
-#if PETSC_VERSION_MINOR==0
-#include "finclude/petscvecdef.h"
-#include "finclude/petscmatdef.h"
-#include "finclude/petsckspdef.h"
-#include "finclude/petscpcdef.h"
-#else
 #include "finclude/petscdef.h"
-#endif
 #else
 #include "finclude/petsc.h"
-#if PETSC_VERSION_MINOR==0
-#include "finclude/petscvec.h"
-#include "finclude/petscmat.h"
-#include "finclude/petscksp.h"
-#include "finclude/petscpc.h"
-#endif
 #endif
     
     private
@@ -85,8 +66,7 @@
     
 !--------------------------------------------------------------------------------------------------------------------
     subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat,&
-      state, inner_mesh, &
-      auxiliary_matrix)
+      state, inner_mesh, auxiliary_matrix)
 !--------------------------------------------------------------------------------------------------------------------
 
       ! Solve Schur complement problem the nice way (using petsc) !!
@@ -161,10 +141,8 @@
 
 !--------------------------------------------------------------------------------------------------------
     subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,solver_option_path, &
-         lstartfromzero,inner_m,div_matrix_comp, &
-         div_matrix_incomp,option_path,preconditioner_matrix,rhs, &
-         state, inner_mesh, &
-         auxiliary_matrix)
+         lstartfromzero,inner_m,div_matrix_comp, div_matrix_incomp,option_path,preconditioner_matrix,rhs, &
+         state, inner_mesh, auxiliary_matrix)
          
 !--------------------------------------------------------------------------------------------------------
 
@@ -230,7 +208,7 @@
       character(len=OPTION_PATH_LEN) :: inner_option_path, inner_solver_option_path
       
       integer reference_node, stat, i
-      logical parallel, have_auxiliary_matrix
+      logical parallel, have_auxiliary_matrix, have_preconditioner_matrix
 
       logical :: apply_reference_node, apply_reference_node_from_coordinates, reference_node_owned
 
@@ -239,7 +217,7 @@
       inner_option_path= trim(option_path)//&
               "/prognostic/scheme/use_projection_method&
               &/full_schur_complement/inner_matrix[0]"
-      
+
       if (have_option(trim(option_path)//'/name')) then
          call get_option(trim(option_path)//'/name', name)
          ewrite(1,*) 'Inside petsc_solve_(block_)csr, solving for: ', trim(name)
@@ -347,18 +325,10 @@
       ! Build Schur complement:
       ewrite(2,*) 'Building Schur complement'                
       if(have_auxiliary_matrix) then
-#if PETSC_VERSION_MINOR==0
-         call MatCreateSchurComplement(inner_M%M,G,G_t_comp,S,A,ierr)
-#else
          call MatCreateSchurComplement(inner_M%M,inner_M%M,G,G_t_comp,S,A,ierr)
-#endif
       else
          myPETSC_NULL_OBJECT=PETSC_NULL_OBJECT
-#if PETSC_VERSION_MINOR==0
-         call MatCreateSchurComplement(inner_M%M,G,G_t_comp,PETSC_NULL_OBJECT,A,ierr)
-#else
          call MatCreateSchurComplement(inner_M%M,inner_M%M,G,G_t_comp,PETSC_NULL_OBJECT,A,ierr)
-#endif
          if (myPETSC_NULL_OBJECT/=PETSC_NULL_OBJECT) then
            FLAbort("PETSC_NULL_OBJECT has changed please report to skramer")
          end if
@@ -386,8 +356,14 @@
       
       ! leaving out petsc_numbering and mesh, so "iteration_vtus" monitor won't work!
 
-      ! Assemble preconditioner matrix in petsc format:
-      pmat=csr2petsc(preconditioner_matrix, petsc_numbering_p, petsc_numbering_p)
+      ! Assemble preconditioner matrix in petsc format (if required):
+      have_preconditioner_matrix=.not.(have_option(trim(option_path)//&
+              "/prognostic/scheme/use_projection_method&
+              &/full_schur_complement/preconditioner_matrix::NoPreconditionerMatrix"))
+
+      if(have_preconditioner_matrix) then
+         pmat=csr2petsc(preconditioner_matrix, petsc_numbering_p, petsc_numbering_p)
+      end if
 
       ! Set up RHS and Solution vectors (note these are loaded later):
       b = PetscNumberingCreateVec(petsc_numbering_p)
@@ -398,8 +374,12 @@
 
       parallel=IsParallel()
 
-      call SetupKSP(ksp,A,pmat,solver_option_path,parallel,&
-        petsc_numbering_p, lstartfromzero)
+      if(have_preconditioner_matrix) then
+         call SetupKSP(ksp,A,pmat,solver_option_path,parallel,petsc_numbering_p, lstartfromzero)
+      else
+         ! If preconditioner matrix is not required, send in A instead:
+         call SetupKSP(ksp,A,A,solver_option_path,parallel,petsc_numbering_p, lstartfromzero)
+      end if
       
       ! Destroy the matrices setup for the schur complement computation. While
       ! these matrices are destroyed here, they are still required for the inner solve,
@@ -410,7 +390,7 @@
       call MatDestroy(G_t_comp,ierr) ! Destroy Compressible Divergence Operator.
       call MatDestroy(G_t_incomp, ierr) ! Destroy Incompressible Divergence Operator.
       call MatDestroy(G, ierr) ! Destroy Gradient Operator (i.e. transpose of incompressible div).
-      call MatDestroy(pmat,ierr) ! Destroy preconditioning matrix if allocated.
+      if(have_preconditioner_matrix) call MatDestroy(pmat,ierr) ! Destroy preconditioning matrix.
       if(have_auxiliary_matrix) call MatDestroy(S,ierr) ! Destroy stabilization matrix
       
       call deallocate( petsc_numbering_u )

=== modified file 'assemble/Geostrophic_Pressure.F90'
--- assemble/Geostrophic_Pressure.F90	2011-07-29 12:30:06 +0000
+++ assemble/Geostrophic_Pressure.F90	2012-10-16 13:22:44 +0000
@@ -110,7 +110,7 @@
     !! Mass option path
     character(len = OPTION_PATH_LEN) :: mass_option_path
     !! Divergence matrix
-    type(block_csr_matrix) :: ct_m
+    type(block_csr_matrix), pointer :: ct_m
     !! RHS terms from integrating the divergence operator by parts
     type(scalar_field) :: ct_rhs
     !! Mass matrix. Only used when not lumping mass for continuous u_mesh.
@@ -785,6 +785,7 @@
     ewrite(2, *) "Boundary conditions field: ", trim(lbcfield%name)
     
     ct_sparsity => get_csr_sparsity_firstorder(state, matrices%p_mesh, matrices%u_mesh)
+    allocate(matrices%ct_m)
     call allocate(matrices%ct_m, ct_sparsity, blocks = (/1, dim/), name = "CT")
     call allocate(matrices%ct_rhs, matrices%p_mesh, "CTRHS")
     
@@ -1197,6 +1198,7 @@
       end select
     end if
     call deallocate(matrices%ct_m)
+    deallocate(matrices%ct_m)
     call deallocate(matrices%ct_rhs)
     
     if(matrices%have_cmc_m) then
@@ -1434,7 +1436,7 @@
     !!< Compute the divergence of a field
   
     type(vector_field), intent(in) :: field
-    type(block_csr_matrix), intent(in) :: ct_m
+    type(block_csr_matrix), pointer :: ct_m
     type(csr_matrix), intent(in) :: mass
     type(scalar_field), intent(inout) :: div
     

=== modified file 'assemble/LES_viscosity.F90'
--- assemble/LES_viscosity.F90	2011-05-10 10:39:15 +0000
+++ assemble/LES_viscosity.F90	2012-10-16 13:22:44 +0000
@@ -54,17 +54,17 @@
     ewrite(2,*) "Initialising optional dynamic LES diagnostic fields"
     ! Filter width
     if(have_filter_width) then
-      tensorfield => extract_tensor_field(state, "DynamicFilterWidth")
+      tensorfield => extract_tensor_field(state, "FilterWidth")
       call zero(tensorfield)
     end if
     ! Strain rate field S1
     if(have_strain) then
-      tensorfield => extract_tensor_field(state, "DynamicStrainRate")
+      tensorfield => extract_tensor_field(state, "StrainRate")
       call zero(tensorfield)
     end if
     ! Filtered strain rate field S2
     if(have_filtered_strain) then
-      tensorfield => extract_tensor_field(state, "DynamicFilteredStrainRate")
+      tensorfield => extract_tensor_field(state, "FilteredStrainRate")
       call zero(tensorfield)
     end if
     ! Eddy viscosity field m_ij
@@ -98,31 +98,31 @@
 
     ! Strain rate field S1
     if(have_strain) then
-      tensorfield => extract_tensor_field(state, "DynamicStrainRate")
+      tensorfield => extract_tensor_field(state, "StrainRate")
       tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), strain_gi, detwei)
       call addto(tensorfield, ele_nodes(nu, ele), tensor_loc)
     end if
 
     ! Filtered strain rate field S2
     if(have_filtered_strain) then
-      tensorfield => extract_tensor_field(state, "DynamicFilteredStrainRate")
+      tensorfield => extract_tensor_field(state, "FilteredStrainRate")
       tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), t_strain_gi, detwei)
       call addto(tensorfield, ele_nodes(nu, ele), tensor_loc)
     end if
 
     ! Filter width
     if(have_filter_width) then
-      tensorfield => extract_tensor_field(state, "DynamicFilterWidth")
+      tensorfield => extract_tensor_field(state, "FilterWidth")
       tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), mesh_size_gi, detwei)
       call addto(tensorfield, ele_nodes(nu, ele), tensor_loc)
     end if
 
   end subroutine les_set_diagnostic_tensor_fields
 
-  subroutine leonard_tensor(mnu, positions, tnu, leonard, alpha, path)
+  subroutine leonard_tensor(nu, positions, tnu, leonard, alpha, path)
 
     ! Unfiltered velocity
-    type(vector_field), pointer               :: mnu
+    type(vector_field), pointer               :: nu
     type(vector_field), intent(in)            :: positions
     ! Filtered velocity
     type(vector_field), pointer               :: tnu
@@ -140,31 +140,30 @@
 
     ! Path is to level above solver options
     lpath = (trim(path)//"/dynamic_les")
-    ewrite(2,*) "path: ", trim(lpath)
-    ewrite(2,*) "alpha: ", alpha
+    ewrite(2,*) "filter factor alpha: ", alpha
 
-    ewrite_minmax(mnu)
+    ewrite_minmax(nu)
     ewrite_minmax(tnu)
 
-    call anisotropic_smooth_vector(mnu, positions, tnu, alpha, lpath)
+    call anisotropic_smooth_vector(nu, positions, tnu, alpha, lpath)
 
-    ewrite_minmax(mnu)
+    ewrite_minmax(nu)
     ewrite_minmax(tnu)
 
     ! Velocity products (ui*uj)
     allocate(ui_uj); allocate(tui_tuj)
-    call allocate(ui_uj, mnu%mesh, "NonlinearVelocityProduct")
-    call allocate(tui_tuj, mnu%mesh, "TestNonlinearVelocityProduct")
+    call allocate(ui_uj, nu%mesh, "NonlinearVelocityProduct")
+    call allocate(tui_tuj, nu%mesh, "TestNonlinearVelocityProduct")
     call zero(ui_uj); call zero(tui_tuj)
 
     ! Other local variables
-    allocate(u_loc(mnu%dim)); allocate(t_loc(mnu%dim, mnu%dim))
+    allocate(u_loc(nu%dim)); allocate(t_loc(nu%dim, nu%dim))
     u_loc=0.0; t_loc=0.0
 
     ! Get cross products of velocities
-    do i=1, node_count(mnu)
+    do i=1, node_count(nu)
       ! Mesh filter ^r
-      u_loc = node_val(mnu,i)
+      u_loc = node_val(nu,i)
       t_loc = outer_product(u_loc, u_loc)
       call set( ui_uj, i, t_loc )
       ! Test filter ^t

=== modified file 'assemble/Makefile.dependencies'
--- assemble/Makefile.dependencies	2011-10-26 08:48:33 +0000
+++ assemble/Makefile.dependencies	2012-10-16 13:22:44 +0000
@@ -31,14 +31,14 @@
    ../include/interpolation_module.mod \
    ../include/intersection_finder_module.mod ../include/mba2d_integration.mod \
    ../include/mba3d_integration.mod ../include/mba_adapt_module.mod \
-   ../include/mesh_files.mod ../include/metric_assemble.mod \
-   ../include/mpi_interfaces.mod ../include/parallel_tools.mod \
-   ../include/pickers.mod ../include/populate_state_module.mod \
+   ../include/metric_assemble.mod ../include/mpi_interfaces.mod \
+   ../include/parallel_tools.mod ../include/pickers.mod \
+   ../include/populate_state_module.mod \
    ../include/project_metric_to_surface_module.mod ../include/quadrature.mod \
    ../include/reserve_state_module.mod ../include/sam_integration.mod \
    ../include/sparse_tools.mod ../include/state_module.mod \
    ../include/tictoc.mod ../include/timeloop_utilities.mod \
-   ../include/zoltan_integration.mod
+   ../include/write_gmsh.mod ../include/zoltan_integration.mod
 
 ../include/adapt_state_prescribed_module.mod: Adapt_State_Prescribed.o
 	@true
@@ -76,11 +76,13 @@
 
 Advection_Diffusion_CG.o ../include/advection_diffusion_cg.mod: \
    Advection_Diffusion_CG.F90 ../include/boundary_conditions.mod \
-   ../include/boundary_conditions_from_options.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/profiler.mod \
-   ../include/quadrature.mod ../include/sparse_tools.mod \
+   ../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/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
 
@@ -89,11 +91,13 @@
 
 Advection_Diffusion_DG.o ../include/advection_diffusion_dg.mod: \
    Advection_Diffusion_DG.F90 ../include/boundary_conditions.mod \
-   ../include/boundary_conditions_from_options.mod ../include/coordinates.mod \
-   ../include/dgtools.mod ../include/diagnostic_fields.mod \
-   ../include/elements.mod ../include/fdebug.h ../include/fefields.mod \
-   ../include/fetools.mod ../include/fields.mod ../include/fldebug.mod \
+   ../include/boundary_conditions_from_options.mod ../include/colouring.mod \
+   ../include/coordinates.mod ../include/dgtools.mod \
+   ../include/diagnostic_fields.mod ../include/elements.mod \
+   ../include/fdebug.h ../include/fefields.mod ../include/fetools.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/shape_functions.mod ../include/slope_limiters_dg.mod \
    ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \
    ../include/sparsity_patterns.mod ../include/sparsity_patterns_meshes.mod \
@@ -106,9 +110,9 @@
 
 Advection_Diffusion_FV.o ../include/advection_diffusion_fv.mod: \
    Advection_Diffusion_FV.F90 ../include/boundary_conditions.mod \
-   ../include/boundary_conditions_from_options.mod ../include/elements.mod \
-   ../include/fdebug.h ../include/fetools.mod ../include/field_options.mod \
-   ../include/fields.mod ../include/fldebug.mod \
+   ../include/boundary_conditions_from_options.mod ../include/colouring.mod \
+   ../include/elements.mod ../include/fdebug.h ../include/fetools.mod \
+   ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \
    ../include/global_parameters.mod ../include/petsc_solve_state_module.mod \
    ../include/profiler.mod ../include/quadrature.mod \
    ../include/shape_functions.mod ../include/sparse_tools.mod \
@@ -131,11 +135,11 @@
 
 Biology.o ../include/biology.mod: Biology.F90 \
    ../include/boundary_conditions.mod ../include/fdebug.h \
-   ../include/fetools.mod ../include/fields.mod \
-   ../include/global_parameters.mod ../include/python_state.mod \
-   ../include/solvers.mod ../include/sparse_tools.mod \
-   ../include/sparsity_patterns_meshes.mod ../include/spud.mod \
-   ../include/state_module.mod
+   ../include/fefields.mod ../include/fetools.mod ../include/field_options.mod \
+   ../include/fields.mod ../include/global_parameters.mod \
+   ../include/python_state.mod ../include/solvers.mod \
+   ../include/sparse_tools.mod ../include/sparsity_patterns_meshes.mod \
+   ../include/spud.mod ../include/state_module.mod
 
 ../include/burgers_assembly.mod: Burgers_Assembly.o
 	@true
@@ -151,17 +155,18 @@
    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
 
 Coriolis.o ../include/coriolis_module.mod: Coriolis.F90 \
    ../include/embed_python.mod ../include/fdebug.h ../include/fldebug.mod \
-   ../include/global_parameters.mod ../include/spud.mod
+   ../include/global_parameters.mod ../include/parallel_tools.mod \
+   ../include/spud.mod
 
 ../include/diagnostic_children.mod: Diagnostic_Children.o
 	@true
@@ -200,10 +205,9 @@
    ../include/geostrophic_pressure.mod ../include/global_parameters.mod \
    ../include/momentum_diagnostic_fields.mod \
    ../include/multimaterial_module.mod ../include/multiphase_module.mod \
-   ../include/parallel_tools.mod ../include/porous_media.mod \
-   ../include/sediment_diagnostics.mod ../include/sparse_matrices_fields.mod \
-   ../include/spontaneous_potentials.mod ../include/spud.mod \
-   ../include/state_module.mod ../include/tidal_module.mod \
+   ../include/parallel_tools.mod ../include/sediment_diagnostics.mod \
+   ../include/sparse_matrices_fields.mod ../include/spontaneous_potentials.mod \
+   ../include/spud.mod ../include/state_module.mod ../include/tidal_module.mod \
    ../include/vorticity_diagnostics.mod
 
 ../include/discrete_properties_module.mod: Discrete_Properties.o
@@ -278,7 +282,7 @@
    ../include/fdebug.h ../include/fefields.mod ../include/field_options.mod \
    ../include/fields.mod ../include/global_parameters.mod ../include/halos.mod \
    ../include/parallel_tools.mod ../include/petsc_solve_state_module.mod \
-   ../include/sparse_matrices_fields.mod \
+   ../include/porous_media.mod ../include/sparse_matrices_fields.mod \
    ../include/sparsity_patterns_meshes.mod ../include/spud.mod \
    ../include/state_fields_module.mod ../include/state_module.mod \
    ../include/transform_elements.mod
@@ -512,21 +516,23 @@
 
 Momentum_CG.o ../include/momentum_cg.mod: Momentum_CG.F90 \
    ../include/boundary_conditions.mod \
-   ../include/boundary_conditions_from_options.mod ../include/coordinates.mod \
-   ../include/coriolis_module.mod ../include/edge_length_module.mod \
-   ../include/elements.mod ../include/fdebug.h ../include/fefields.mod \
-   ../include/fetools.mod ../include/field_derivatives.mod \
-   ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \
+   ../include/boundary_conditions_from_options.mod ../include/colouring.mod \
+   ../include/coordinates.mod ../include/coriolis_module.mod \
+   ../include/edge_length_module.mod ../include/elements.mod \
+   ../include/fdebug.h ../include/fefields.mod ../include/fetools.mod \
+   ../include/field_derivatives.mod ../include/field_options.mod \
+   ../include/fields.mod ../include/fldebug.mod \
    ../include/global_parameters.mod ../include/halos.mod \
    ../include/les_viscosity_module.mod ../include/metric_tools.mod \
    ../include/multiphase_module.mod ../include/petsc_solve_state_module.mod \
-   ../include/rotated_boundary_conditions.mod ../include/smoothing_module.mod \
-   ../include/solvers.mod ../include/sparse_matrices_fields.mod \
-   ../include/sparse_tools.mod ../include/sparse_tools_petsc.mod \
-   ../include/sparsity_patterns_meshes.mod ../include/spud.mod \
-   ../include/state_fields_module.mod ../include/state_matrices_module.mod \
-   ../include/state_module.mod ../include/transform_elements.mod \
-   ../include/upwind_stabilisation.mod ../include/vector_tools.mod
+   ../include/profiler.mod ../include/rotated_boundary_conditions.mod \
+   ../include/smoothing_module.mod ../include/solvers.mod \
+   ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \
+   ../include/sparse_tools_petsc.mod ../include/sparsity_patterns_meshes.mod \
+   ../include/spud.mod ../include/state_fields_module.mod \
+   ../include/state_matrices_module.mod ../include/state_module.mod \
+   ../include/transform_elements.mod ../include/upwind_stabilisation.mod \
+   ../include/vector_tools.mod
 
 ../include/momentum_dg.mod: Momentum_DG.o
 	@true
@@ -558,6 +564,7 @@
    ../include/equation_of_state.mod ../include/fdebug.h \
    ../include/field_priority_lists.mod ../include/fields.mod \
    ../include/fldebug.mod ../include/global_parameters.mod \
+   ../include/initialise_fields_module.mod ../include/k_epsilon.mod \
    ../include/multimaterial_module.mod ../include/multiphase_module.mod \
    ../include/spud.mod ../include/state_module.mod
 
@@ -579,7 +586,8 @@
    ../include/momentum_cg.mod ../include/momentum_dg.mod \
    ../include/momentum_diagnostic_fields.mod ../include/multiphase_module.mod \
    ../include/oceansurfaceforcing.mod ../include/parallel_tools.mod \
-   ../include/petsc_solve_state_module.mod ../include/profiler.mod \
+   ../include/petsc_solve_state_module.mod \
+   ../include/pressure_dirichlet_bcs_cv.mod ../include/profiler.mod \
    ../include/reduced_model_runtime.mod \
    ../include/rotated_boundary_conditions.mod ../include/slope_limiters_dg.mod \
    ../include/solvers.mod ../include/sparse_matrices_fields.mod \
@@ -642,7 +650,23 @@
 	@true
 
 Porous_Media.o ../include/porous_media.mod: Porous_Media.F90 \
-   ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \
+   ../include/fdebug.h ../include/fields.mod ../include/futils.mod \
+   ../include/global_parameters.mod ../include/spud.mod \
+   ../include/state_module.mod
+
+../include/pressure_dirichlet_bcs_cv.mod: Pressure_Dirichlet_BCS_CV.o
+	@true
+
+Pressure_Dirichlet_BCS_CV.o ../include/pressure_dirichlet_bcs_cv.mod: \
+   Pressure_Dirichlet_BCS_CV.F90 ../include/boundary_conditions.mod \
+   ../include/cv_face_values.mod ../include/cv_faces.mod \
+   ../include/cv_fields.mod ../include/cv_options.mod \
+   ../include/cv_shape_functions.mod ../include/cv_upwind_values.mod \
+   ../include/cvtools.mod ../include/diagnostic_fields.mod ../include/fdebug.h \
+   ../include/fetools.mod ../include/field_derivatives.mod \
+   ../include/field_options.mod ../include/fields.mod ../include/futils.mod \
+   ../include/global_parameters.mod ../include/multiphase_module.mod \