Merge lp:~g-bhutani12/fluidity/PBM into lp:fluidity

[Gaurav Bhutani]

A fully functional population balance modelling capability has been added to Fluidity, It is based on the famous Direct Quadrature Method of Moments (DQMOM) by Marchisio and Fox (2005). Currently it can solve population balance equations with one internal variable. Support for multiple internal variables will be added at a later stage (1-2 years time). The branch has been tested for adaptivity and parallel runs.

[lp:~g-bhutani12/fluidity/PBM]

4210. By Gaurav Bhutani on 2013-10-21

Removed some old test cases

4211. By Gaurav Bhutani on 2013-10-22

Merge with trunk

4212. By Gaurav Bhutani on 2013-10-24

Removed a few old test cases

4213. By Gaurav Bhutani on 2013-10-27

Changes made to diamond tree in Multiphase options to accept a scalar field for partcle diameter (only accepted a constant earlier). Multiphase.F90 modified to accept a scalar field as particle diameter and use it in drag force term calculations for momentum equations; the file has also been modified for new diamond structure. Multiphase test cases modified to suit modified diamond structure. All tests passed with this new modifications.

4214. By Gaurav Bhutani on 2013-10-28

Included the scalar field -- SauterMeanDia -- in the list of statistics. Modified the comparison conditions for other statistics to exact string comparison. Modified the rnc file to include the scalar field called SauterMeanDia in the list of statistics. SauterMeanDia has been added to the statistics so that it can be calculated after every non-linear iteration in Fluids.F90 and the new sauter mean diameter be used in the momentum calculations.

4215. By Gaurav Bhutani on 2013-10-28

Changes to Fluids.F90 to call dqmom_calculate_moments() and dqmom_calculate_statistics inside the nonlinear loop. This updates the SauterMeanDia after every nonlinear iteration to be used in the drag force term in the momentum equation in next nonlinear iteration.

4216. By Gaurav Bhutani on 2013-10-31

Corrected the PD algorithm implementation in DQMOM.F90

[Tim Greaves (tim-greaves)]

Hi Gaurav -- given this was proposed more than a couple of months ago and looks to be in active development at the moment I'm setting it back to 'in progress'; if you get to the point you want to move for a review and merge, please do set back to 'needs review'.

[Gaurav Bhutani (g-bhutani12)]

Hi Tim,

Yes, I have been working on the branch recently to couple it with Multiphase.F90.
Apologies, I completely forgot to change the status of it to WIP.
Thanks for that.
I will propose it for a merge very soon with all the documentation and test cases in it.

Regards,
Gaurav
________________________________________
From: <email address hidden> [<email address hidden>] on behalf of Tim Greaves [<email address hidden>]
Sent: 01 November 2013 17:41
To: Bhutani, Gaurav
Subject: Re: [Merge] lp:~g-bhutani12/fluidity/PBM into lp:fluidity

Hi Gaurav -- given this was proposed more than a couple of months ago and looks to be in active development at the moment I'm setting it back to 'in progress'; if you get to the point you want to move for a review and merge, please do set back to 'needs review'.
--
https://code.launchpad.net/~g-bhutani12/fluidity/PBM/+merge/183296
You are the owner of lp:~g-bhutani12/fluidity/PBM.

[lp:~g-bhutani12/fluidity/PBM]

4217. By Gaurav Bhutani on 2013-11-05

Included a breakage kernel as used by McCoy and Madras (2003) for the daughter distribution function. Also corrected the aggregation birth term contribution into source to handle negative abscissas and their cube roots.

4218. By Gaurav Bhutani on 2013-11-06

Added a test case for homogeneous population balance with breakage and aggregation. The analytical solution of McCoy and Madras (2003) was used to compare against the numerical results.

4219. By Gaurav Bhutani on 2013-11-15

Added saling factor to scale Sauter Mean Diameter in case of non-dimensional implementation of PBE.

4220. By Gaurav Bhutani on 2013-11-20

Added modified Schiller & Naumann (1933) interphase drag coefficient to Multiphase.F90 as implemented in ANSYS CFX. Changes made to diamond tree for the option for this new drag coefficient. Slight change in Multiphase.F90 to change defualt to default in a switch-case statement.

4221. By Gaurav Bhutani on 2013-11-25

Slight change to Schiller-Naumann (1933) drag implementation - multiplied by the volume fraction of continuous phase in drag force term as done by Abbasi et. al. (2013).

4222. By Gaurav Bhutani on 2013-12-06

PBM: Added a sum_aggregation kernel as implemented by Abbasi et. al. (2013). Also added the population balance statistic called MeanDia10 (d_10). DQMOM.F90 and fluidity_oprions.rnc modified for the two changes.

4223. By Gaurav Bhutani on 2013-12-12

Modified Multiphase_Diagnostics.F90 to allow particle_reynolds_number to take a scalar field for particle diameter. schemas/diagnostic_algorithms.rnc was also modified for the same reason. flml in 4 test cases in multiphase were modified as per the schema modifications.

Unmerged revisions

4223. By Gaurav Bhutani on 2013-12-12

Modified Multiphase_Diagnostics.F90 to allow particle_reynolds_number to take a scalar field for particle diameter. schemas/diagnostic_algorithms.rnc was also modified for the same reason. flml in 4 test cases in multiphase were modified as per the schema modifications.

4222. By Gaurav Bhutani on 2013-12-06

PBM: Added a sum_aggregation kernel as implemented by Abbasi et. al. (2013). Also added the population balance statistic called MeanDia10 (d_10). DQMOM.F90 and fluidity_oprions.rnc modified for the two changes.

4221. By Gaurav Bhutani on 2013-11-25

Slight change to Schiller-Naumann (1933) drag implementation - multiplied by the volume fraction of continuous phase in drag force term as done by Abbasi et. al. (2013).

4220. By Gaurav Bhutani on 2013-11-20

Added modified Schiller & Naumann (1933) interphase drag coefficient to Multiphase.F90 as implemented in ANSYS CFX. Changes made to diamond tree for the option for this new drag coefficient. Slight change in Multiphase.F90 to change defualt to default in a switch-case statement.

4219. By Gaurav Bhutani on 2013-11-15

Added saling factor to scale Sauter Mean Diameter in case of non-dimensional implementation of PBE.

4218. By Gaurav Bhutani on 2013-11-06

Added a test case for homogeneous population balance with breakage and aggregation. The analytical solution of McCoy and Madras (2003) was used to compare against the numerical results.

4217. By Gaurav Bhutani on 2013-11-05

Included a breakage kernel as used by McCoy and Madras (2003) for the daughter distribution function. Also corrected the aggregation birth term contribution into source to handle negative abscissas and their cube roots.

4216. By Gaurav Bhutani on 2013-10-31

Corrected the PD algorithm implementation in DQMOM.F90

4215. By Gaurav Bhutani on 2013-10-28

Changes to Fluids.F90 to call dqmom_calculate_moments() and dqmom_calculate_statistics inside the nonlinear loop. This updates the SauterMeanDia after every nonlinear iteration to be used in the drag force term in the momentum equation in next nonlinear iteration.

4214. By Gaurav Bhutani on 2013-10-28

Included the scalar field -- SauterMeanDia -- in the list of statistics. Modified the comparison conditions for other statistics to exact string comparison. Modified the rnc file to include the scalar field called SauterMeanDia in the list of statistics. SauterMeanDia has been added to the statistics so that it can be calculated after every non-linear iteration in Fluids.F90 and the new sauter mean diameter be used in the momentum calculations.

=== modified file 'Makefile.in'
--- Makefile.in	2013-11-13 00:28:27 +0000
+++ Makefile.in	2013-12-12 22:35:40 +0000
@@ -310,6 +310,8 @@
 	@cd horizontal_adaptivity;   	$(MAKE)
 	@echo "    MAKE preprocessor"
 	@cd preprocessor; $(MAKE)
+	@echo "    MAKE population_balance"
+	@cd population_balance; $(MAKE)
 	@echo "    MAKE error_measures"
 	@cd error_measures;   	$(MAKE)
 	@echo "    MAKE assemble"
@@ -340,6 +342,7 @@
 forward_interfaces/*.o \
 horizontal_adaptivity/*.o \
 preprocessor/*.o \
+population_balance/*.o \
 error_measures/*.o \
 assemble/*.o \
 diagnostics/*.o \
@@ -387,6 +390,8 @@
 	@cd diagnostics;   $(MAKE) clean
 	@echo "    CLEAN preprocessor"
 	@cd preprocessor; $(MAKE) clean
+	@echo "    CLEAN population_balance"
+	@cd population_balance; $(MAKE) clean
 	@echo "    CLEAN error_measures/tests"
 	@cd error_measures/tests;   $(MAKE) clean
 	@echo "    CLEAN error_measures"
@@ -588,6 +593,8 @@
 	@echo "    Generating preprocessor dependencies"
 	@cd preprocessor; ../bin/create_makefile --exclude \
 	"register_diagnostics.F90 check_options.F90" $(TESTOPTS)
+	@echo "    Generating population_balance dependencies"
+	@cd population_balance; ../bin/create_makefile $(TESTOPTS)
 	@echo "    Generating error_measures dependencies"
 	@cd error_measures;   	../bin/create_makefile $(TESTOPTS)
 	@echo "    Generating assemble dependencies"

=== modified file 'assemble/Adapt_State.F90'
--- assemble/Adapt_State.F90	2013-11-13 01:51:23 +0000
+++ assemble/Adapt_State.F90	2013-12-12 22:35:40 +0000
@@ -74,6 +74,7 @@
   use diagnostic_fields_wrapper_new, only : calculate_diagnostic_variables_new => calculate_diagnostic_variables
   use pickers
   use write_gmsh
+  use dqmom
 #ifdef HAVE_ZOLTAN
   use zoltan_integration
   use mpi_interfaces
@@ -878,7 +879,7 @@
     type(state_type), dimension(:), intent(inout) :: states
 
     character(len = *), parameter :: base_path = "/mesh_adaptivity/hr_adaptivity/adapt_at_first_timestep"
-    integer :: adapt_iterations, i
+    integer :: adapt_iterations, i, i_state
     type(mesh_type), pointer :: old_mesh
     type(tensor_field) :: metric
     type(vector_field), pointer :: output_positions
@@ -896,7 +897,6 @@
       call copy_to_stored_values(states,"Old")
       call copy_to_stored_values(states,"Iterated")
       call relax_to_nonlinear(states)
-
       call calculate_diagnostic_variables(states)
       call calculate_diagnostic_variables_new(states)
 
@@ -916,6 +916,12 @@
       ! Adapt state, initialising fields from the options tree rather than
       ! interpolating them
       call adapt_state(states, metric, initialise_fields = .true.)
+      
+      ! dqmom_init helps to recalculate the abscissas and weights based on moment initial conditions (if provided)
+      call dqmom_init(states)
+!      do i_state = 1, option_count("/material_phase")
+!          call dqmom_calculate_moments(states(i_state))
+!      end do
     end do
 
     if(have_option(trim(base_path) // "/output_adapted_mesh")) then

=== modified file 'assemble/Diagnostic_fields_wrapper.F90'
--- assemble/Diagnostic_fields_wrapper.F90	2013-04-16 11:38:41 +0000
+++ assemble/Diagnostic_fields_wrapper.F90	2013-12-12 22:35:40 +0000
@@ -54,6 +54,7 @@
   use momentum_diagnostic_fields
   use spontaneous_potentials, only: calculate_formation_conductivity
   use sediment_diagnostics
+  use dqmom
   use geostrophic_pressure
   use multiphase_module
   
@@ -80,7 +81,7 @@
     type(state_type), dimension(:), pointer :: submaterials
     
     ewrite(1, *) "In calculate_diagnostic_variables"
- 
+
     do i = 1, size(state)
 
        ! start of fields that can be called through the generic calculate_diagnostic_variable
@@ -562,6 +563,11 @@
           call calculate_sediment_active_layer_volume_fractions(state(i))
        end if
        ! End of sediment diagnostics.
+       
+       ! Start of population balance diagnostics.
+       call dqmom_calculate_moments(state(i))
+       call dqmom_calculate_statistics(state(i))
+       ! End of population balance diagnostics.
 
        ! Multiphase-related diagnostic fields
        s_field => extract_scalar_field(state(i), "PhaseVolumeFraction", stat)

=== modified file 'assemble/Multiphase.F90'
--- assemble/Multiphase.F90	2013-03-03 01:36:16 +0000
+++ assemble/Multiphase.F90	2013-12-12 22:35:40 +0000
@@ -284,9 +284,10 @@
    
          ! Types of drag correlation
          integer, parameter :: DRAG_CORRELATION_TYPE_STOKES = 1, DRAG_CORRELATION_TYPE_WEN_YU = 2, DRAG_CORRELATION_TYPE_ERGUN = 3
+         integer, parameter :: DRAG_CORRELATION_TYPE_SCHILLER_NAUMANN = 4
          
          ewrite(1, *) "Entering add_fluid_particle_drag"
-         
+ 
          ! Let's check whether we actually have at least one particle phase.
          if(option_count("/material_phase/multiphase_properties/particle_diameter") == 0) then
             FLExit("Fluid-particle drag enabled but no particle_diameter has been specified for the particle phase(s).")
@@ -357,6 +358,8 @@
 
                type(scalar_field), pointer :: vfrac_fluid, vfrac_particle
                type(scalar_field), pointer :: density_fluid, density_particle
+               type(scalar_field), pointer :: pd_scalar_field ! scalar field defining particle diameter 
+
                type(vector_field), pointer :: velocity_fluid, velocity_particle
                type(vector_field), pointer :: oldu_fluid, oldu_particle
                type(vector_field), pointer :: nu_fluid, nu_particle ! Non-linear approximation to the Velocities
@@ -364,7 +367,9 @@
                type(scalar_field) :: nvfrac_fluid, nvfrac_particle
                real :: d ! Particle diameter
                character(len=OPTION_PATH_LEN) :: drag_correlation_name
+               character(len=OPTION_PATH_LEN) :: pd_field_name ! name of scalar field that defines particle diameter (can be the sauter mean dia)
                integer :: drag_correlation
+               logical :: have_constant_pd ! checks if the particle diameter is a constant or not 
 
                ! Get the necessary fields to calculate the drag force
                velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
@@ -376,8 +381,15 @@
                   density_fluid => extract_scalar_field(state(istate_fluid), "Density")
                   density_particle => extract_scalar_field(state(istate_particle), "Density")
                   viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
-         
-                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
+
+                  if(have_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/constant")) then
+                     call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/constant", d)
+                     have_constant_pd = .true.
+                  else if(have_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/use_scalar_field")) then
+                     call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/use_scalar_field", pd_field_name)
+                     pd_scalar_field => extract_scalar_field(state(istate_particle), pd_field_name)
+                     have_constant_pd = .false.
+                  end if
 
                   ! Calculate the non-linear approximation to the PhaseVolumeFractions
                   call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction")
@@ -401,7 +413,9 @@
                         drag_correlation = DRAG_CORRELATION_TYPE_WEN_YU
                      case("ergun")
                         drag_correlation = DRAG_CORRELATION_TYPE_ERGUN
-                     case("default")
+                     case("schiller_naumann")
+                        drag_correlation = DRAG_CORRELATION_TYPE_SCHILLER_NAUMANN
+                     case default
                         FLAbort("Unknown correlation for fluid-particle drag")
                   end select
                   
@@ -420,7 +434,9 @@
                                                             density_fluid, density_particle, &
                                                             nu_fluid, nu_particle, &
                                                             oldu_fluid, oldu_particle, &
-                                                            viscosity_fluid, d, drag_correlation)
+                                                            viscosity_fluid, &
+                                                            have_constant_pd, d, pd_scalar_field, &
+                                                            drag_correlation)
                      end if
 
                   end do element_loop
@@ -438,7 +454,9 @@
                                                       density_fluid, density_particle, &
                                                       nu_fluid, nu_particle, &
                                                       oldu_fluid, oldu_particle, &
-                                                      viscosity_fluid, d, drag_correlation)
+                                                      viscosity_fluid, &
+                                                      have_constant_pd, dia, pd_scalar_field, &
+                                                      drag_correlation)
                                                          
                integer, intent(in) :: ele
                type(element_type), intent(in) :: test_function
@@ -449,12 +467,14 @@
                     
                type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle
                type(scalar_field), intent(in) :: density_fluid, density_particle
+               type(scalar_field), intent(in) :: pd_scalar_field ! Scalar field representing particle diameter
                type(vector_field), intent(in) :: nu_fluid, nu_particle
                type(vector_field), intent(in) :: oldu_fluid, oldu_particle
                type(tensor_field), intent(in) :: viscosity_fluid    
-               real, intent(in) :: d ! Particle diameter 
+               real, intent(in) :: dia ! Constant particle diameter 
                integer, intent(in) :: drag_correlation
-               
+               logical, intent(in) :: have_constant_pd ! is true if particle diameter is a constant. is false if it is a scalar field (e.g. sauter mean dia)
+
                ! Local variables
                real, dimension(ele_ngi(u,ele)) :: vfrac_fluid_gi, vfrac_particle_gi
                real, dimension(ele_ngi(u,ele)) :: density_fluid_gi, density_particle_gi
@@ -473,6 +493,7 @@
                real, dimension(ele_ngi(u,ele)) :: particle_re_gi ! Particle Reynolds number
                real, dimension(ele_ngi(u,ele)) :: drag_coefficient_gi
                real, dimension(ele_ngi(u,ele)) :: magnitude_gi ! |v_f - v_p|
+               real, dimension(ele_ngi(u,ele)) :: d_gi ! particle diameter at the Gauss points
                
                real, dimension(ele_ngi(u,ele)) :: K
                real, dimension(ele_ngi(u,ele)) :: drag_force_big_m
@@ -497,9 +518,16 @@
                   magnitude_gi(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
                end do
 
+               ! Compute the particle diameter at the Gauss points
+               if(have_constant_pd) then
+                  d_gi = dia
+               else 
+                  d_gi = ele_val_at_quad(pd_scalar_field, ele)
+               end if
+
                ! Compute the particle Reynolds number
                ! (Assumes isotropic viscosity for now)
-               particle_re_gi = (vfrac_fluid_gi*density_fluid_gi*magnitude_gi*d) / viscosity_fluid_gi(1,1,:)
+               particle_re_gi = (vfrac_fluid_gi*density_fluid_gi*magnitude_gi*d_gi) / viscosity_fluid_gi(1,1,:)
            
                ! Compute the drag coefficient
                select case(drag_correlation)
@@ -525,6 +553,18 @@
                      
                   case(DRAG_CORRELATION_TYPE_ERGUN)
                      ! No drag coefficient is needed here.                  
+
+                  case(DRAG_CORRELATION_TYPE_SCHILLER_NAUMANN)
+                     ! Schiller & Naumann (1933) drag correlation
+                     ! Since the particle Reynolds number definition currently contains vfrac_fluid in numerator,
+                     ! we need to take that out by dividing as done below. Must be changed when Reynolds number defn is changed - gb812
+                     do gi = 1, ele_ngi(u,ele)
+                        if((particle_re_gi(gi)/vfrac_fluid_gi(gi)) < 1000) then
+                           drag_coefficient_gi(gi) = (24.0/(particle_re_gi(gi)/vfrac_fluid_gi(gi)))*(1.0+0.15*(particle_re_gi(gi)/vfrac_fluid_gi(gi))**0.687)
+                        else
+                           drag_coefficient_gi(gi) = 0.44
+                        end if
+                     end do
                end select
                       
                ! Don't let the drag_coefficient_gi be NaN
@@ -536,12 +576,15 @@
            
                select case(drag_correlation)
                   case(DRAG_CORRELATION_TYPE_STOKES)
-                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d)
+                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d_gi)
                   case(DRAG_CORRELATION_TYPE_WEN_YU)
                      ! Wen & Yu (1966) drag term
-                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d*vfrac_fluid_gi**2.7)
+                     K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d_gi*vfrac_fluid_gi**2.7)
                   case(DRAG_CORRELATION_TYPE_ERGUN)
-                     K = 150.0*((vfrac_particle_gi**2)*viscosity_fluid_gi(1,1,:))/(vfrac_fluid_gi*(d**2)) + 1.75*(vfrac_particle_gi*density_fluid_gi*magnitude_gi/d)
+                     K = 150.0*((vfrac_particle_gi**2)*viscosity_fluid_gi(1,1,:))/(vfrac_fluid_gi*(d_gi**2)) + &
+                                1.75*(vfrac_particle_gi*density_fluid_gi*magnitude_gi/d_gi)
+                  case(DRAG_CORRELATION_TYPE_SCHILLER_NAUMANN)
+                     K = vfrac_fluid_gi*vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(density_fluid_gi*magnitude_gi)/(d_gi)
                end select               
                
                if(is_particle_phase) then
@@ -680,6 +723,7 @@
 
                type(scalar_field), pointer :: vfrac_fluid, vfrac_particle
                type(scalar_field), pointer :: density_fluid, density_particle
+               type(scalar_field), pointer :: pd_scalar_field ! scalar field defining particle diameter
                type(vector_field), pointer :: velocity_fluid, velocity_particle
                type(scalar_field), pointer :: internal_energy_fluid, internal_energy_particle
                type(scalar_field), pointer :: old_internal_energy_fluid, old_internal_energy_particle
@@ -691,6 +735,8 @@
                real :: C_fluid, C_particle ! Specific heat of the fluid and particle phases at constant volume
                real :: gamma ! Ratio of specific heats for compressible phase
                integer :: kstat, cstat_fluid, cstat_particle, gstat
+               character(len=OPTION_PATH_LEN) :: pd_field_name ! name of scalar field that defines particle diameter (can be the sauter mean dia)
+               logical :: have_constant_pd ! checks if the particle diameter is a constant or not 
 
                ! Get the necessary fields to calculate the heat transfer term
                velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity")
@@ -703,7 +749,14 @@
                   density_particle => extract_scalar_field(state(istate_particle), "Density")
                   viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity")
          
-                  call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d)
+                  if(have_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/constant")) then
+                     call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/constant", d)
+                     have_constant_pd = .true.
+                  else if(have_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/use_scalar_field")) then
+                     call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter/use_scalar_field", pd_field_name)
+                     pd_scalar_field => extract_scalar_field(state(istate_particle), pd_field_name)
+                     have_constant_pd = .false.
+                  end if
                            
                   call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/effective_conductivity", k, kstat)
                            
@@ -759,7 +812,9 @@
                                                       internal_energy_particle, &
                                                       old_internal_energy_fluid, &
                                                       old_internal_energy_particle, &
-                                                      viscosity_fluid, d, k, C_fluid, &
+                                                      viscosity_fluid, &
+                                                      have_constant_pd, d, pd_scalar_field, &
+                                                      k, C_fluid, &
                                                       C_particle, gamma)
                      end if
 
@@ -781,7 +836,9 @@
                                                 internal_energy_particle, &
                                                 old_internal_energy_fluid, &
                                                 old_internal_energy_particle, &
-                                                viscosity_fluid, d, k, C_fluid, &
+                                                viscosity_fluid, &
+                                                have_constant_pd, dia, pd_scalar_field, &
+                                                k, C_fluid, &
                                                 C_particle, gamma)
                                                          
                integer, intent(in) :: ele
@@ -794,11 +851,14 @@
                     
                type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle
                type(scalar_field), intent(in) :: density_fluid, density_particle
+               type(scalar_field), intent(in) :: pd_scalar_field ! Scalar field representing particle diameter
                type(vector_field), intent(in) :: nu_fluid, nu_particle
                type(scalar_field), intent(in) :: internal_energy_fluid, internal_energy_particle
                type(scalar_field), intent(in) :: old_internal_energy_fluid, old_internal_energy_particle
                type(tensor_field), intent(in) :: viscosity_fluid    
-               real, intent(in) :: d, k, C_fluid, C_particle, gamma
+               real, intent(in) :: dia ! Constant particle diameter 
+               real, intent(in) :: k, C_fluid, C_particle, gamma
+               logical, intent(in) :: have_constant_pd ! is true if particle diameter is a constant. is false if it is a scalar field (e.g. sauter mean dia)
                
                ! Local variables
                real, dimension(ele_ngi(x,ele)) :: internal_energy_fluid_gi, internal_energy_particle_gi
@@ -814,6 +874,7 @@
                real, dimension(ele_loc(internal_energy,ele)) :: rhs_addto
                real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: matrix_addto
                
+               real, dimension(ele_ngi(x,ele)) :: d_gi ! particle diameter at quadrature points
                real, dimension(ele_ngi(x,ele)) :: particle_Re ! Particle Reynolds number
                real, dimension(ele_ngi(x,ele)) :: Pr ! Prandtl number
                real, dimension(ele_ngi(x,ele)) :: particle_Nu ! Particle Nusselt number
@@ -844,9 +905,16 @@
                   velocity_magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi))
                end do
 
+               ! Compute the particle diameter at quadrature points
+               if(have_constant_pd) then
+                  d_gi = dia
+               else 
+                  d_gi = ele_val_at_quad(pd_scalar_field, ele)
+               end if
+
                ! Compute the particle Reynolds number
                ! (Assumes isotropic viscosity for now)
-               particle_Re = (density_fluid_gi*velocity_magnitude*d) / viscosity_fluid_gi(1,1,:)
+               particle_Re = (density_fluid_gi*velocity_magnitude*d_gi) / viscosity_fluid_gi(1,1,:)
            
                ! Compute the Prandtl number
                ! (Assumes isotropic viscosity for now)
@@ -856,7 +924,7 @@
                particle_Nu = (7.0 - 10.0*vfrac_fluid_gi + 5.0*vfrac_fluid_gi**2)*(1.0 + 0.7*(particle_Re**0.2)*(Pr**(1.0/3.0))) + &
                               (1.33 - 2.4*vfrac_fluid_gi + 1.2*vfrac_fluid_gi**2)*(particle_Re**0.7)*(Pr**(1.0/3.0))
 
-               Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d**2)
+               Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d_gi**2)
                
                ! Note that the transfer term is defined in terms of temperatures (T_fluid and T_particle)
                ! Let's convert the temperatures to internal energy (E) using E_i = C_i*T_i,

=== modified file 'configure'
--- configure	2013-10-17 17:13:25 +0000
+++ configure	2013-12-12 22:35:40 +0000
@@ -14918,7 +14918,7 @@
 
 #*******************
 
-ac_config_files="$ac_config_files Makefile debug/Makefile bathymetry/Makefile ocean_forcing/Makefile ocean_forcing/tests/Makefile sediments/Makefile hyperlight/Makefile femtools/Makefile femtools/tests/Makefile forward_interfaces/Makefile horizontal_adaptivity/Makefile horizontal_adaptivity/tests/Makefile preprocessor/Makefile error_measures/Makefile error_measures/tests/Makefile parameterisation/Makefile parameterisation/tests/Makefile fldecomp/Makefile assemble/Makefile assemble/tests/Makefile diagnostics/Makefile main/Makefile tools/Makefile tools/version-info python/setup.py adjoint/Makefile adjoint/tests/Makefile climatology/Makefile libmba2d/Makefile libmba3d/Makefile libjudy/Makefile libjudy/src/Makefile libjudy/src/JudyCommon/Makefile libjudy/src/Judy1/Makefile libjudy/src/JudyL/Makefile libjudy/src/JudySL/Makefile libjudy/src/JudyHS/Makefile libwm/Makefile libvtkfortran/Makefile schemas/flml reduced_modelling/Makefile"
+ac_config_files="$ac_config_files Makefile debug/Makefile bathymetry/Makefile ocean_forcing/Makefile ocean_forcing/tests/Makefile sediments/Makefile population_balance/Makefile hyperlight/Makefile femtools/Makefile femtools/tests/Makefile forward_interfaces/Makefile horizontal_adaptivity/Makefile horizontal_adaptivity/tests/Makefile preprocessor/Makefile error_measures/Makefile error_measures/tests/Makefile parameterisation/Makefile parameterisation/tests/Makefile fldecomp/Makefile assemble/Makefile assemble/tests/Makefile diagnostics/Makefile main/Makefile tools/Makefile tools/version-info python/setup.py adjoint/Makefile adjoint/tests/Makefile climatology/Makefile libmba2d/Makefile libmba3d/Makefile libjudy/Makefile libjudy/src/Makefile libjudy/src/JudyCommon/Makefile libjudy/src/Judy1/Makefile libjudy/src/JudyL/Makefile libjudy/src/JudySL/Makefile libjudy/src/JudyHS/Makefile libwm/Makefile libvtkfortran/Makefile schemas/flml reduced_modelling/Makefile"
 
 cat >confcache <<\_ACEOF
 # This file is a shell script that caches the results of configure
@@ -15629,6 +15629,7 @@
     "ocean_forcing/Makefile") CONFIG_FILES="$CONFIG_FILES ocean_forcing/Makefile" ;;
     "ocean_forcing/tests/Makefile") CONFIG_FILES="$CONFIG_FILES ocean_forcing/tests/Makefile" ;;
     "sediments/Makefile") CONFIG_FILES="$CONFIG_FILES sediments/Makefile" ;;
+    "population_balance/Makefile") CONFIG_FILES="$CONFIG_FILES population_balance/Makefile" ;;
     "hyperlight/Makefile") CONFIG_FILES="$CONFIG_FILES hyperlight/Makefile" ;;
     "femtools/Makefile") CONFIG_FILES="$CONFIG_FILES femtools/Makefile" ;;
     "femtools/tests/Makefile") CONFIG_FILES="$CONFIG_FILES femtools/tests/Makefile" ;;

=== modified file 'configure.in'
--- configure.in	2013-09-23 20:07:28 +0000
+++ configure.in	2013-12-12 22:35:40 +0000
@@ -1503,6 +1503,7 @@
           bathymetry/Makefile
           ocean_forcing/Makefile ocean_forcing/tests/Makefile
           sediments/Makefile
+          population_balance/Makefile
           hyperlight/Makefile
           femtools/Makefile femtools/tests/Makefile
           forward_interfaces/Makefile

=== modified file 'diagnostics/Multiphase_Diagnostics.F90'
--- diagnostics/Multiphase_Diagnostics.F90	2011-10-24 02:20:27 +0000
+++ diagnostics/Multiphase_Diagnostics.F90	2013-12-12 22:35:40 +0000
@@ -56,7 +56,10 @@
          type(scalar_field), intent(inout) :: s_field ! Particle Reynolds number field
 
          !! Sub-options of the diagnostic field
-         real :: d ! Particle diameter
+         real :: dia ! Particle diameter
+         character(len = OPTION_PATH_LEN) :: pd_field_name !name of scalar field that defines particle diameter (can be the sauter mean dia)
+         type(scalar_field), pointer :: pd_scalar_field ! scalar field referring to particle diameter 
+         logical :: have_constant_pd ! checks if the particle diameter is a constant or not
          character(len = OPTION_PATH_LEN) :: continuous_phase_name
 
          !! Other local variables
@@ -78,6 +81,7 @@
          real, dimension(mesh_dim(s_field), ele_ngi(s_field, 1)) :: u_continuous_gi, u_particle_gi
          real, dimension(mesh_dim(s_field), mesh_dim(s_field), ele_ngi(s_field, 1)) :: viscosity_gi
          real, dimension(ele_ngi(s_field, 1)) :: density_gi, vfrac_gi
+         real, dimension(ele_ngi(s_field, 1)) :: d_gi ! particle diameter at the Gauss points
 
          real, dimension(:), allocatable :: magnitude ! |v_f - v_p|
 
@@ -93,7 +97,14 @@
          ewrite(1,*) 'Entering calculate_particle_reynolds_number'
 
          ! Get sub-options from under the diagnostic field in the options tree
-         call get_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/particle_diameter', d)
+         if(have_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/particle_diameter/constant')) then
+            call get_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/particle_diameter/constant', dia)
+            have_constant_pd = .true.
+         else if(have_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/particle_diameter/use_scalar_field')) then
+            call get_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/particle_diameter/use_scalar_field', pd_field_name)
+            pd_scalar_field => extract_scalar_field(states(state_index), pd_field_name)
+            have_constant_pd = .false.
+         end if
          call get_option(trim(s_field%option_path)//'/diagnostic/algorithm::particle_reynolds_number/continuous_phase_name', continuous_phase_name)
 
          ! Find the index of the continuous phase in states
@@ -130,13 +141,20 @@
             density_gi = ele_val_at_quad(density, ele)
             vfrac_gi = ele_val_at_quad(vfrac, ele)
 
+            ! Compute the particle diameter at the Gauss points
+            if(have_constant_pd) then
+               d_gi = dia
+            else
+               d_gi = ele_val_at_quad(pd_scalar_field, ele)
+            end if
+
             do gi = 1, ele_ngi(u_continuous, ele)
                magnitude(gi) = norm2(u_continuous_gi(:,gi) - u_particle_gi(:,gi))
             end do
 
             ! Compute the particle Reynolds number
             ! (Assumes isotropic viscosity for now)
-            particle_re_gi = (vfrac_gi*density_gi*magnitude*d) / viscosity_gi(1,1,:)
+            particle_re_gi = (vfrac_gi*density_gi*magnitude*d_gi) / viscosity_gi(1,1,:)
 
             ! Invert the mass matrix to get the particle_re value at each node
             particle_re_mat = matmul(inverse(shape_shape(particle_re_shape, particle_re_shape, detwei)), &

=== modified file 'examples/tephra_settling/tephra_settling.flml'
--- examples/tephra_settling/tephra_settling.flml	2013-05-30 11:47:42 +0000
+++ examples/tephra_settling/tephra_settling.flml	2013-12-12 22:35:40 +0000
@@ -584,7 +584,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">48e-6</real_value>
+        <constant>
+          <real_value rank="0">48e-6</real_value>
+        </constant>
       </particle_diameter>
     </multiphase_properties>
   </material_phase>

=== modified file 'main/Fluids.F90'
--- main/Fluids.F90	2013-10-07 16:23:35 +0000
+++ main/Fluids.F90	2013-12-12 22:35:40 +0000
@@ -103,6 +103,7 @@
   use detector_parallel, only: sync_detector_coordinates, deallocate_detector_list_array
   use momentum_diagnostic_fields, only: calculate_densities
   use sediment_diagnostics, only: calculate_sediment_flux
+  use DQMOM
 
   implicit none
 
@@ -248,13 +249,16 @@
          & default=1)
     call get_option("/timestepping/nonlinear_iterations/tolerance", &
          & nonlinear_iteration_tolerance, default=0.0)
-    
+
     if(have_option("/mesh_adaptivity/hr_adaptivity/adapt_at_first_timestep")) then
 
        if(have_option("/timestepping/nonlinear_iterations/nonlinear_iterations_at_adapt")) then
          call get_option('/timestepping/nonlinear_iterations/nonlinear_iterations_at_adapt',nonlinear_iterations_adapt)
          nonlinear_iterations = nonlinear_iterations_adapt
        end if
+      
+       ! set population balance initial conditions - for first adaptivity
+       call dqmom_init(state)
 
        call adapt_state_first_timestep(state)
 
@@ -299,6 +303,9 @@
        call populate_sub_state(state,sub_state)
     end if
 
+    ! set population balance initial conditions
+    call dqmom_init(state)
+
     ! Calculate the number of scalar fields to solve for and their correct
     ! solve order taking into account dependencies.
     call get_ntsol(ntsol)
@@ -571,6 +578,9 @@
 
           call compute_goals(state)
 
+          ! Calculate source terms for populations ------
+          call dqmom_calculate_source_terms(state, ITS)
+
           !------------------------------------------------
           ! Addition for calculating drag force ------ jem 05-06-2008
           if (have_option("/imported_solids/calculate_drag_on_surface")) then
@@ -750,6 +760,7 @@
             end if
           end do 
 
+
           ! This is where the non-legacy momentum stuff happens
           ! a loop over state (hence over phases) is incorporated into this subroutine call
           ! hence this lives outside the phase_loop
@@ -762,6 +773,14 @@
              call solve_momentum(state,at_first_timestep=((timestep==1).and.(its==1)),timestep=timestep, POD_state=POD_state)
           end if
 
+          ! calculate abscissa for populations -----
+          ! this must be done at the end of each non-linear iteration
+          call dqmom_calculate_abscissa(state)
+          do i = 1, size(state)
+             call dqmom_calculate_moments(state(i))
+             call dqmom_calculate_statistics(state(i))
+          end do
+
           if(nonlinear_iterations > 1) then
              ! Check for convergence between non linear iteration loops
              call test_and_write_convergence(state, current_time + dt, dt, its, change)

=== added directory 'population_balance'
=== added file 'population_balance/DQMOM.F90'
--- population_balance/DQMOM.F90	1970-01-01 00:00:00 +0000
+++ population_balance/DQMOM.F90	2013-12-12 22:35:40 +0000
@@ -0,0 +1,999 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!    
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!    
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied arranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+
+#include "fdebug.h"
+
+module dqmom
+
+  use vector_tools
+  use spud
+  use fields
+  use state_module
+  use state_fields_module
+  use initialise_fields_module
+  use global_parameters, only: OPTION_PATH_LEN, FIELD_NAME_LEN
+!  use sparsity_patterns_meshes
+  use sparse_tools
+  use solvers
+
+  implicit none
+
+  public dqmom_init, dqmom_calculate_source_terms, dqmom_calculate_abscissa,&
+       & dqmom_check_options, dqmom_calculate_moments, dqmom_calculate_statistics
+
+  private
+
+  !! TODO:
+  !! 1. Make field_finding algorithm work for multiple phases and make sure it works for
+  !!    several pop_balances
+  !! 3. Make check_options run
+  !! 4. Check all prognostic fields are identical (possible bar initial conditions)
+  !! 5. Check prognostic Source terms are set to diagnostic, Internal
+  !! 6. Minimum weight
+  !! 7. Check diagnostic fields are set to Internal
+  !! 8. Entire set up of fields is a horrible hack - could do with sorting out - suggest 
+  !!    new blueprint for how this could work - will only work for one pop_balance per phase. IMPORTANT
+  !! 9. Make changes to function dqmom_calculate_abscissa so that it works for a state insteaad of all states. 
+  !!    Make necessary changes to the place this is called in the code. 
+  !! 10.In Fluids.F90 use conditions on presence of population balance tree and call the functions inside the if condition.
+  !!    Loop on states inside Fluids.F90 (or any other file where DQMOM functions are called) instead of doing it inside DQMOM functions.  
+
+  !! Changes made by gb812
+  !! 1. in 'Construct C matrix (rhs pt. 2)' -> took grad_D out or squared term
+  !! 2. in 'Construct C matrix (rhs pt. 2)' -> while summing grad_D in the next line, changed from sum(grad_D) to sum(grad_D,1) 
+contains
+
+  subroutine get_pop_field(state, i_pop, i_field, type, item, stat, iterated)
+
+    ! Returns the required population balance field
+
+    type(state_type), intent(in) :: state
+    integer, intent(in) :: i_pop, i_field
+    character(len=FIELD_NAME_LEN) :: type ! moments, weights, ascissa, or weighted_abscissa
+    type(scalar_field), pointer, intent(out) :: item
+    integer, intent(out), optional :: stat
+    logical, intent(in), optional :: iterated
+    
+    character(len=FIELD_NAME_LEN) :: name
+    character(len=OPTION_PATH_LEN) :: option_path 
+
+    call get_pop_option_path(state, i_pop, option_path)
+    call get_option(trim(option_path)//'/'//trim(type)//'/scalar_field['//int2str(i_field -&
+         1)//']/name', name) 
+    if (present(iterated)) then
+       if (iterated) then
+          name = 'Iterated'//trim(name)
+       end if
+    end if
+    item => extract_scalar_field(state, trim(name), stat)
+
+  end subroutine get_pop_field
+
+  subroutine get_pop_option_path(state, i_pop, option_path)
+
+    type(state_type), intent(in) :: state
+    integer, intent(in) :: i_pop
+    character(len=OPTION_PATH_LEN), intent(out) :: option_path 
+
+    option_path = trim(state%option_path)//'/population_balance'
+    if (option_count(option_path) > 1) then
+       option_path = trim(state%option_path)//'/population_balance['&
+            &//int2str(i_pop)//']'
+    end if
+
+  end subroutine get_pop_option_path
+
+  subroutine dqmom_init(states)
+
+    type(state_type), intent(in), dimension(:) :: states
+
+    integer :: i_state, i_pop
+    character(len=OPTION_PATH_LEN) :: option_path 
+    
+    do i_state = 1, option_count("/material_phase")
+       do i_pop = 1, option_count(trim(states(i_state)%option_path)//&
+            '/population_balance')
+          call get_pop_option_path(states(i_state), i_pop, option_path)
+          if (have_option(trim(option_path)// &
+               '/calculate_initial_conditions_from_moments')) then
+             call dqmom_PD_algorithm(states(i_state), i_pop)
+          end if
+       end do
+    end do
+
+    call dqmom_calculate_abscissa(states)
+
+  end subroutine dqmom_init
+
+  subroutine dqmom_PD_algorithm(state, i_pop)
+
+    ! Algorithm for calculating abscissa and weights from the moment of a distribution
+    ! See Gordon 1968 and Mcgraw 1997
+   
+    type(state_type), intent(in) :: state
+    integer, intent(in) :: i_pop
+
+    type(scalar_field_pointer), dimension(:), allocatable :: moments
+    type(scalar_field_pointer), dimension(:), allocatable :: weighted_abscissa
+    type(scalar_field_pointer), dimension(:), allocatable :: weights
+    type(vector_field), pointer:: position
+
+    ! J - Jacobi matrix
+    ! P and alpha are required working arrays
+    ! e_values - eigenvalues
+    ! e_vectors - corresponding eigenvectors
+    real, dimension(:,:), allocatable :: P, Jac, e_vectors
+    real, dimension(:), allocatable :: alpha, e_values
+    character(len=OPTION_PATH_LEN) :: option_path
+    character(len=FIELD_NAME_LEN) :: type
+    integer :: i_field, n_abscissa, i_node, i, j, stat
+
+    call get_pop_option_path(state, i_pop, option_path)
+    n_abscissa = option_count(trim(option_path)//'/abscissa/scalar_field')
+
+    allocate(moments(2*n_abscissa))
+    allocate(weighted_abscissa(n_abscissa))
+    allocate(weights(n_abscissa))
+    allocate(P(2*n_abscissa + 1,2*n_abscissa + 1))
+    allocate(Jac(n_abscissa,n_abscissa))
+    allocate(alpha(2*n_abscissa))
+    allocate(e_vectors(n_abscissa,n_abscissa))
+    allocate(e_values(n_abscissa))
+
+    position => extract_vector_field(state, "Coordinate")
+    do i_field = 1, 2*n_abscissa
+       ! collect moment fields and initialise from initial conditions
+       ! will not be done automatically as this is a diagnostic field
+       type = 'moments'
+       call get_pop_field(state, i_pop, i_field, type, moments(i_field)%ptr)
+       call zero(moments(i_field)%ptr)
+       call initialise_field_over_regions(moments(i_field)%ptr, &
+            trim(moments(i_field)%ptr%option_path)//'/diagnostic/initial_condition', position)
+    end do
+    
+    do i_field = 1, n_abscissa
+       ! collect weighted_abscissa and weight fields and zero
+       type = 'weights'
+       call get_pop_field(state, i_pop, i_field, type, weights(i_field)%ptr)
+       type = 'weighted_abscissa'
+       call get_pop_field(state, i_pop, i_field, type, &
+            weighted_abscissa(i_field)%ptr)
+       call zero(weights(i_field)%ptr)
+       call zero(weighted_abscissa(i_field)%ptr)
+    end do
+
+    do i_node = 1, node_count(moments(1)%ptr)
+       ! zero working arrays
+       P = 0.0
+       Jac = 0.0
+       alpha = 0.0
+       e_vectors = 0.0
+       e_values = 0.0
+       
+       ! Construct P matrix
+       P(1,1) = 1.0
+       ! set the zero'th moment to 1.0 for the purposes of finding the abscissa
+       ! weights will be multiplied by the zero'th moment later to obtain their correct
+       ! values
+       P(1,2) = 1.0       
+       do i = 2, 2*n_abscissa
+          P(i,2) = (-1)**(i-1)*(node_val(moments(i)%ptr, i_node)/node_val(moments(1)%ptr, i_node))
+       end do
+       do j = 3, 2*n_abscissa + 1
+          do i = 1, (2*n_abscissa + 1) - (j - 1)
+             P(i,j) = P(1,j-1)*P(i+1,j-2) - P(1,j-2)*P(i+1,j-1)
+          end do
+       end do
+       ! Construct alpha array
+       alpha(1) = 0.0
+       do i = 2, 2*n_abscissa
+          alpha(i) = P(1,i+1)/(P(1,i)*P(1,i-1))
+       end do
+       ! Construct jacobi matrix
+       do i = 1, n_abscissa
+          Jac(i,i) = alpha(2*i) + alpha(2*i - 1)
+       end do
+       do i = 1, n_abscissa - 1
+          Jac(i, i+1) = ((alpha(2*i + 1)*alpha(2*i))**2.0)**0.25
+          Jac(i+1, i) = ((alpha(2*i + 1)*alpha(2*i))**2.0)**0.25
+       end do
+       
+       ! Calculate eigenvalues and eigenvectors
+       call eigendecomposition_symmetric(Jac, e_vectors, e_values, stat)
+       if (stat /= 0) then
+          FLExit('Cannot compute abscissa and weights using PD algorithm')
+       end if
+
+       do i_field = 1, n_abscissa
+          ! set prognostic field values
+          call set(weights(i_field)%ptr, i_node, &
+               node_val(moments(1)%ptr, i_node) * e_vectors(1,i_field)**2)
+          call set(weighted_abscissa(i_field)%ptr, i_node, &
+               node_val(weights(i_field)%ptr, i_node) * e_values(i_field))
+       end do
+       
+    end do
+
+    deallocate(moments, weighted_abscissa, weights, P, Jac, alpha, e_vectors, e_values)
+    
+  end subroutine dqmom_PD_algorithm
+
+  subroutine dqmom_calculate_abscissa(states)
+
+    type(state_type), dimension(:), intent(in) :: states
+
+    type(scalar_field), pointer :: abscissa, weight, weighted_abscissa
+    type(scalar_field) :: inv_weight
+    integer :: i_state, i_pop, i_abscissa
+    character(len=OPTION_PATH_LEN) :: option_path 
+    character(len=FIELD_NAME_LEN) :: type
+    
+    ewrite(1, *) "In dqmom_calculate_abscissa" 
+    do i_state = 1, option_count("/material_phase")
+       do i_pop = 1, option_count(trim(states(i_state)%option_path)//'/population_balance')
+          call get_pop_option_path(states(i_state), i_pop, option_path)
+          do i_abscissa = 1, option_count(trim(option_path)//'/abscissa/scalar_field')
+             ! get required fields
+             type = 'abscissa'
+             call get_pop_field(states(i_state), i_pop, i_abscissa, type, abscissa)
+             type = 'weights'
+             call get_pop_field(states(i_state), i_pop, i_abscissa, type, weight)
+             type = 'weighted_abscissa'
+             call get_pop_field(states(i_state), i_pop, i_abscissa, type, weighted_abscissa)
+
+             ! calculate the inverse of the weights
+             call allocate(inv_weight, weight%mesh, 'InverseWeight')
+             call set(inv_weight, weight)
+             where (inv_weight%val/=0.0)
+                inv_weight%val=1./inv_weight%val
+             end where
+
+             ! calculate abscissa
+             call set(abscissa, weighted_abscissa)
+             call scale(abscissa, inv_weight)
+
+             call deallocate(inv_weight)
+          end do
+       end do
+    end do
+    ewrite(1, *) "Exiting dqmom_calculate_abscissa"
+  end subroutine dqmom_calculate_abscissa
+
+  subroutine dqmom_calculate_source_terms(states, it)
+
+    type(state_type), dimension(:), intent(inout) :: states
+    integer, intent(in) :: it    
+
+    integer :: i_state, i_pop
+
+    do i_state = 1, option_count("/material_phase")
+       do i_pop = 1, option_count(trim(states(i_state)%option_path)//'/population_balance')
+          call dqmom_calculate_source_term_pop(states(i_state), it, i_pop)
+       end do
+    end do
+
+  end subroutine dqmom_calculate_source_terms
+
+  subroutine dqmom_calculate_source_term_pop(state, it, i_pop)
+
+    type(state_type), intent(inout) :: state
+    integer, intent(in) :: it
+
+    type(scalar_field_pointer), dimension(:), allocatable :: abscissa,&
+         weight, it_abscissa, it_weight, s_weighted_abscissa, s_weight
+    type(scalar_field), pointer :: lumped_mass
+    type(csr_matrix), pointer :: mass_matrix
+    type(scalar_field), dimension(:), allocatable :: r_abscissa, r_weight
+    type(tensor_field), pointer :: D
+    type(vector_field), pointer :: X
+    real :: theta, cond, growth_r, internal_dispersion_coeff, aggregation_freq_const, breakage_freq_const, breakage_freq_degree, perturb_val
+    integer :: i_pop, N, i, j, stat
+    character(len=OPTION_PATH_LEN) :: option_path 
+    character(len=FIELD_NAME_LEN) :: type, field_name, growth_type, aggregation_freq_type, breakage_freq_type, breakage_dist_type, singular_option
+    logical :: have_D = .false.
+    logical :: have_growth = .FALSE.
+    logical :: have_internal_dispersion = .FALSE.
+    logical :: have_aggregation = .FALSE.
+    logical :: have_breakage = .FALSE.
+
+    call get_pop_option_path(state, i_pop, option_path)
+    N = option_count(trim(option_path)//'/abscissa/scalar_field')
+    allocate(abscissa(N), weight(N), it_abscissa(N), it_weight(N), &
+         r_abscissa(N), r_weight(N), s_weighted_abscissa(N), s_weight(N))
+
+    do i = 1, N
+       ! collect abscissa and weight fields
+       type = 'weights'
+       call get_pop_field(state, i_pop, i, type, weight(i)%ptr)
+       call get_pop_field(state, i_pop, i, type, it_weight(i)%ptr, iterated=.true.)
+       type = 'abscissa'
+       call get_pop_field(state, i_pop, i, type, abscissa(i)%ptr)
+       call get_pop_field(state, i_pop, i, type, it_abscissa(i)%ptr, iterated=.true.)
+
+       ! get source fields (note this is the weighted abscissa source not the abscissa source)
+       s_weight(i)%ptr => extract_scalar_field(state, trim(weight(i)%ptr%name)//'Source')
+       call get_option(trim(option_path)//'/weighted_abscissa/scalar_field['// &
+            int2str(i - 1)//']/name', field_name) 
+       s_weighted_abscissa(i)%ptr => extract_scalar_field(state, trim(field_name)//'Source')
+       call zero(s_weight(i)%ptr)
+       call zero(s_weighted_abscissa(i)%ptr)
+
+       ! relax non-linear values
+       ! all temporal relaxations must be the same
+       call get_option(trim(weight(i)%ptr%option_path)// &
+            '/prognostic/temporal_discretisation/theta', theta)
+       ! do not recalculate source terms if theta = 0.0 after first non-linear iteration
+       if ((theta == 0.0) .and. (it /= 1)) then 
+          return
+       end if
+       call allocate(r_abscissa(i), abscissa(i)%ptr%mesh, &
+            "RelaxedAbscissa"//int2str(i))
+       call allocate(r_weight(i), weight(i)%ptr%mesh, "RelaxedWeight"//int2str(i))
+       call set(r_abscissa(i), it_abscissa(i)%ptr)
+       call scale(r_abscissa(i), theta)
+       call addto(r_abscissa(i), abscissa(i)%ptr, 1.0-theta)
+       call set(r_weight(i), it_weight(i)%ptr)
+       call scale(r_weight(i), theta)
+       call addto(r_weight(i), weight(i)%ptr, 1.0-theta)       
+    end do    
+
+    ! get diffusion
+    D => extract_tensor_field(state, trim(weight(1)%ptr%name)//'Diffusivity',stat)
+    if (stat == 0) then
+       have_D = .true.
+    end if
+
+    ! check for growth term
+    if (have_option(trim(option_path)//'/population_balance_source_terms/growth')) then
+       have_growth = .TRUE.
+       if (have_option(trim(option_path)//'/population_balance_source_terms/growth/power_law_growth')) then
+          growth_type = 'power_law_growth';
+          call get_option(trim(option_path)//'/population_balance_source_terms/growth/power_law_growth', growth_r)
+       end if
+    else
+       have_growth = .FALSE.
+    end if
+
+    ! check for internal dispersion term
+    if (have_option(trim(option_path)//'/population_balance_source_terms/internal_dispersion')) then
+       have_internal_dispersion = .TRUE.
+       call get_option(trim(option_path)//'/population_balance_source_terms/internal_dispersion', internal_dispersion_coeff)
+    else
+       have_internal_dispersion = .FALSE.
+    end if
+
+    ! check for aggregation term
+    if (have_option(trim(option_path)//'/population_balance_source_terms/aggregation')) then
+       have_aggregation = .TRUE.
+       if (have_option(trim(option_path)//'/population_balance_source_terms/aggregation/aggregation_frequency/constant_aggregation')) then
+          aggregation_freq_type = 'constant_aggregation'
+          call get_option(trim(option_path)//'/population_balance_source_terms/aggregation/aggregation_frequency/constant_aggregation', aggregation_freq_const)
+       else if (have_option(trim(option_path)//'/population_balance_source_terms/aggregation/aggregation_frequency/hydrodynamic_aggregation')) then
+          aggregation_freq_type = 'hydrodynamic_aggregation'
+       else if (have_option(trim(option_path)//'/population_balance_source_terms/aggregation/aggregation_frequency/sum_aggregation')) then
+          aggregation_freq_type = 'sum_aggregation'
+          call get_option(trim(option_path)//'/population_balance_source_terms/aggregation/aggregation_frequency/sum_aggregation', aggregation_freq_const)
+       end if
+    else
+       have_aggregation = .FALSE.
+    end if
+
+    ! check for breakage term
+    if (have_option(trim(option_path)//'/population_balance_source_terms/breakage')) then
+       have_breakage = .TRUE.
+       if (have_option(trim(option_path)//'/population_balance_source_terms/breakage/breakage_frequency/constant_breakage')) then
+          breakage_freq_type = 'constant_breakage'
+          call get_option(trim(option_path)//'/population_balance_source_terms/breakage/breakage_frequency/constant_breakage', breakage_freq_const)
+       else if (have_option(trim(option_path)//'/population_balance_source_terms/breakage/breakage_frequency/power_law_breakage')) then
+          breakage_freq_type = 'power_law_breakage'
+          call get_option(trim(option_path)//'/population_balance_source_terms/breakage/breakage_frequency/power_law_breakage/coefficient', breakage_freq_const)
+          call get_option(trim(option_path)//'/population_balance_source_terms/breakage/breakage_frequency/power_law_breakage/degree', breakage_freq_degree)
+       end if
+
+       if (have_option(trim(option_path)//'/population_balance_source_terms/breakage/distribution_function/symmetric_fragmentation')) then
+          breakage_dist_type = 'symmetric_fragmentation'
+       else if (have_option(trim(option_path)//'/population_balance_source_terms/breakage/distribution_function/mcCoy_madras_2003')) then
+          breakage_dist_type = 'mcCoy_madras_2003'
+       end if
+   
+    else
+       have_breakage = .FALSE.
+    end if
+
+    ! get ill-conditioned matrix settings
+    call get_option(trim(option_path)//'/ill_conditioned_matrices/required_condition_number', cond)
+    if (have_option(trim(option_path)//'/ill_conditioned_matrices/set_source_to_zero')) then
+       singular_option = 'set_source_to_zero';
+    else if (have_option(trim(option_path)//'/ill_conditioned_matrices/perturbate')) then 
+       singular_option = 'perturbate';
+       call get_option(trim(option_path)//'/ill_conditioned_matrices/perturbate/perturbation', perturb_val)
+!    else if (have_option(trim(option_path)//'/ill_conditioned_matrices/average_surrounding_nodes')) then        
+!       singular_option = 'average_surrounding_nodes';
+    end if
+
+    X => extract_vector_field(state, 'Coordinate')
+    ! assembly loop
+    do i = 1, ele_count(r_abscissa(1))
+       call dqmom_calculate_source_term_ele(r_abscissa, r_weight, s_weighted_abscissa, s_weight, &
+                &D, have_D, have_growth, growth_type, growth_r, have_internal_dispersion, internal_dispersion_coeff, &
+                &have_aggregation, aggregation_freq_type, aggregation_freq_const, &
+                &have_breakage, breakage_freq_type, breakage_freq_const, breakage_freq_degree, breakage_dist_type, X, singular_option, perturb_val, cond, i)       
+    end do
+
+    ! for non-DG we apply inverse mass globally
+    if(continuity(r_abscissa(1))>=0) then
+       if(have_option(trim(option_path)//'/adv_diff_source_term_interpolation/use_full_mass_matrix')) then
+          mass_matrix => get_mass_matrix(state, r_abscissa(1)%mesh)
+          ! r_abscissa(1) is used as a dummy scalar field here since it is not needed anymore. 
+          do j = 1, N
+             call zero(r_abscissa(1))
+             call petsc_solve(r_abscissa(1), mass_matrix, s_weight(j)%ptr, trim(option_path)//'/adv_diff_source_term_interpolation/use_full_mass_matrix')
+             call set(s_weight(j)%ptr, r_abscissa(1))
+             call zero(r_abscissa(1))
+             call petsc_solve(r_abscissa(1), mass_matrix, s_weighted_abscissa(j)%ptr, trim(option_path)//'/adv_diff_source_term_interpolation/use_full_mass_matrix')
+             call set(s_weighted_abscissa(j)%ptr, r_abscissa(1))
+          end do
+       else if(have_option(trim(option_path)//'/adv_diff_source_term_interpolation/use_mass_lumping')) then
+          lumped_mass => get_lumped_mass(state, r_abscissa(1)%mesh)
+          do j = 1, N
+             do i = 1, node_count(r_abscissa(1))
+                call set(s_weighted_abscissa(j)%ptr, i, node_val(s_weighted_abscissa(j)%ptr,i)&
+                     &/node_val(lumped_mass,i))
+                call set(s_weight(j)%ptr, i, node_val(s_weight(j)%ptr,i)&
+                     &/node_val(lumped_mass,i))
+             end do
+          end do
+       else 
+          FLAbort("Check the .flml file. You must specify an option under 'population_balance/adv_diff_source_term_interpolation'")
+       end if
+    end if
+
+    do i = 1, N
+       call deallocate(r_abscissa(i))
+       call deallocate(r_weight(i))
+    end do
+    deallocate(abscissa, weight, it_abscissa, it_weight, &
+         r_abscissa, r_weight)
+
+  end subroutine dqmom_calculate_source_term_pop
+
+  subroutine dqmom_calculate_source_term_ele(abscissa, weight, s_weighted_abscissa, s_weight, &
+                 &D, have_D, have_growth, growth_type, growth_r, have_internal_dispersion, internal_dispersion_coeff, &
+                 &have_aggregation, aggregation_freq_type, aggregation_freq_const, &
+                 &have_breakage, breakage_freq_type, breakage_freq_const, breakage_freq_degree, breakage_dist_type, X, singular_option, perturb_val, cond, ele)
+
+    type(scalar_field), dimension(:), intent(in) :: abscissa, weight
+    type(scalar_field_pointer), dimension(:), intent(inout) :: s_weighted_abscissa, s_weight
+    type(tensor_field), pointer, intent(in) :: D
+    type(vector_field), pointer, intent(in) :: X
+    integer, intent(in) :: ele
+    real, intent(in) :: cond, growth_r, internal_dispersion_coeff, aggregation_freq_const, breakage_freq_const, breakage_freq_degree, perturb_val
+    logical, intent(in) :: have_D, have_growth, have_internal_dispersion, have_aggregation, have_breakage
+    character(len=FIELD_NAME_LEN), intent(in) :: growth_type, aggregation_freq_type, breakage_freq_type, breakage_dist_type, singular_option
+
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)) :: abscissa_val_at_quad
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)*2, size(abscissa)*2) :: A
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)*2, size(abscissa)) :: A_3
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)) :: C
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)*2) :: S_rhs  ! source term (includes growth, breakage and coalescence term): gb 15-11-2012
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)*2, size(abscissa)) :: moment_daughter_dist_func
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)) :: break_freq
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa), size(abscissa)) :: aggregation_freq   ! at present it is not dependent on space coordinate, but can be dependent and will have to be a scalar fields
+    real, dimension(size(abscissa)*2, 1) :: b
+    real, dimension(ele_ngi(abscissa(1), ele), size(abscissa)) :: abscissa_S_at_quad
+    real, dimension(ele_ngi(abscissa(1), ele), size(weight)) :: weight_S_at_quad
+    real, dimension(ele_loc(abscissa(1), ele)) :: abscissa_S_at_nodes
+    real, dimension(ele_loc(abscissa(1), ele)) :: weight_S_at_nodes
+    real, dimension(ele_loc(abscissa(1), ele), ele_loc(abscissa(1), ele)) :: invmass
+    real, dimension(ele_ngi(abscissa(1), ele)) :: detwei
+    real, dimension(X%dim, ele_ngi(abscissa(1), ele)) :: grad_D
+    real, dimension(X%dim, X%dim, ele_ngi(abscissa(1), ele)) :: D_at_quad
+    type(element_type), pointer :: shape
+    integer, dimension(:), pointer :: nodes
+    real, dimension(ele_loc(abscissa(1), ele), ele_ngi(abscissa(1), ele), X%dim) :: dshape
+    real, dimension(size(abscissa)*2, size(abscissa)*2) :: svd_tmp1, svd_tmp2
+    real, dimension(size(abscissa)*2) :: SV
+    integer :: stat, N, i, j, k
+    real :: curr_time
+    real, dimension(size(abscissa)) :: perturb
+    integer :: p_num
+    real, dimension(size(abscissa)) :: absc_init
+
+    N = size(abscissa)
+    
+    nodes => ele_nodes(abscissa(1), ele)
+    shape => ele_shape(abscissa(1), ele)
+
+    call transform_to_physical(X, ele, shape, dshape=dshape, detwei=detwei)
+
+    ! construct A matrices (lhs knowns)
+    do i = 1, N
+       abscissa_val_at_quad(:,i) = ele_val_at_quad(abscissa(i), ele)       
+!!       print*, "abscissa = ", abscissa_val_at_quad(1,i)
+    end do
+    A = A_matrix(abscissa_val_at_quad)
+
+    ! construct A_3 matrix (rhs pt.1)
+    do i = 1, 2*N
+       do j = 1, N
+          A_3(:,i,j) = (i-1)*(i-2)*ele_val_at_quad(abscissa(j), ele)**(i-3)
+       end do
+    end do
+
+    ! construct C matrix (rhs pt.2)
+    if (have_D) then
+       do i = 1, N
+          D_at_quad = ele_val_at_quad(D, ele)
+          do j = 1, X%dim
+             grad_D(j,:) = D_at_quad(j,j,:)
+          end do
+          grad_D = ((ele_grad_at_quad(abscissa(i), ele, dshape))**2)*grad_D
+          C(:,i) = ele_val_at_quad(weight(i), ele)*sum(grad_D,1)
+       end do
+    else
+       C = 0.0
+    end if
+
+
+    !! gb nov-2012---
+    ! initialize dqmom source term to zero 
+    S_rhs = 0.0
+
+    ! construct S vector (rhs pt.3) for GROWTH term
+    if (have_growth) then
+       if (growth_type=='power_law_growth') then
+          do i = 1, 2*N
+             do j = 1, N
+                S_rhs(:,i) = S_rhs(:,i) + (i-1)*ele_val_at_quad(weight(j), ele)*(ele_val_at_quad(abscissa(j), ele)**(i-2))*(ele_val_at_quad(abscissa(j), ele)**growth_r)
+             end do 
+          end do 
+       end if
+    end if
+    
+    if (have_internal_dispersion) then
+       do i = 1, 2*N
+          do j = 1, N
+             ! S_rhs(:,i) = S_rhs(:,i) + (i-1)*(i-2)*ele_val_at_quad(weight(j), ele)*(ele_val_at_quad(abscissa(j), ele)**(i-3))*Diffusion_internal(j)
+             ! abscissa_val_at_quad() takes care of the perturbed abscissas
+             S_rhs(:,i) = S_rhs(:,i) + (i-1)*(i-2)*ele_val_at_quad(weight(j), ele)*(abscissa_val_at_quad(:,j)**(i-3))*internal_dispersion_coeff
+          end do
+       end do
+    end if
+
+    !!! construct S vector for BREAKAGE
+!    do i = 1, N
+!       break_freq(:,i) = 0.02*abscissa_val_at_quad(:,i)**3
+!    end do
+    
+    if (have_breakage) then
+       if (breakage_freq_type=='constant_breakage') then
+          break_freq = breakage_freq_const
+       else if (breakage_freq_type=='power_law_breakage') then
+          do i = 1, N
+             break_freq(:,i) = breakage_freq_const*abscissa_val_at_quad(:,i)**breakage_freq_degree
+          end do
+       end if
+
+       if (breakage_dist_type=='symmetric_fragmentation') then
+          do i = 1, 2*N
+             do j = 1, N
+                moment_daughter_dist_func(:,i,j) = (2.0**(((3-(i-1))/3.0)))*(abscissa_val_at_quad(:,j)**(i-1))   
+!                moment_daughter_dist_func(:,i,j) = (2.0**(((3-(i-1))/3.0)))*(ele_val_at_quad(abscissa(j), ele)**(i-1))
+             end do
+          end do
+       else if (breakage_dist_type=='mcCoy_madras_2003') then
+          do i = 1, 2*N
+             do j = 1, N
+                moment_daughter_dist_func(:,i,j) = (6.0/((i-1)+3.0))*(abscissa_val_at_quad(:,j)**(i-1))
+!                moment_daughter_dist_func(:,i,j) = (6.0/((i-1)+3.0))*(ele_val_at_quad(abscissa(j), ele)**(i-1))
+             end do
+          end do
+       end if
+
+       do i = 1, 2*N 
+          do j = 1, N
+             ! birth term due to breakage
+             S_rhs(:,i) = S_rhs(:,i) + break_freq(:,j)*ele_val_at_quad(weight(j), ele)*moment_daughter_dist_func(:,i,j)   ! daughter distribution function already includes the factor for number of particles formed after breakage
+             ! death term due to breakage
+             S_rhs(:,i) = S_rhs(:,i) - break_freq(:,j)*ele_val_at_quad(weight(j), ele)*(abscissa_val_at_quad(:,j)**(i-1))
+!             S_rhs(:,i) = S_rhs(:,i) - break_freq(:,j)*ele_val_at_quad(weight(j), ele)*(ele_val_at_quad(abscissa(j), ele)**(i-1))
+          end do
+       end do
+    endif
+
+    
+    !!! construct S vector for AGGREGATION
+    if (have_aggregation) then
+       if (aggregation_freq_type=='constant_aggregation') then
+          aggregation_freq = aggregation_freq_const
+       else if (aggregation_freq_type=='hydrodynamic_aggregation') then
+          do i = 1, N
+             do j = 1, N
+                aggregation_freq(:,i,j) = abscissa_val_at_quad(:,i)**3 + abscissa_val_at_quad(:,j)**3
+             end do
+          end do
+       else if (aggregation_freq_type=='sum_aggregation') then
+          do i = 1, N
+             do j = 1, N
+                aggregation_freq(:,i,j) = (abscissa_val_at_quad(:,i) + abscissa_val_at_quad(:,j))*aggregation_freq_const
+             end do
+          end do
+       end if
+
+       do i = 1, 2*N
+          do j = 1, N
+             do k = 1, N
+                ! birth term due to aggregation
+                S_rhs(:,i) = S_rhs(:,i) + 0.5 * aggregation_freq(:,j,k) * ele_val_at_quad(weight(j), ele) * ele_val_at_quad(weight(k), ele) * &
+                                           ((abs(abscissa_val_at_quad(:,j)**3 + abscissa_val_at_quad(:,k)**3)**(1.0/3.0)) * &
+                                           sign(1.0,(abscissa_val_at_quad(:,j)**3 + abscissa_val_at_quad(:,k)**3)))**(i-1)
+!                S_rhs(:,i) = S_rhs(:,i) + 0.5 * aggregation_freq(:,j,k) * ele_val_at_quad(weight(j), ele) * ele_val_at_quad(weight(k), ele) * &
+!                                           &(ele_val_at_quad(abscissa(j), ele)**3 + ele_val_at_quad(abscissa(k), ele)**3)**((i-1)/3.0)
+                ! death term due to aggregation
+                S_rhs(:,i) = S_rhs(:,i) - aggregation_freq(:,j,k) * ele_val_at_quad(weight(j), ele) * ele_val_at_quad(weight(k), ele) * abscissa_val_at_quad(:,j)**(i-1)
+!                S_rhs(:,i) = S_rhs(:,i) - aggregation_freq(:,j,k) * ele_val_at_quad(weight(j), ele)*ele_val_at_quad(weight(k), ele)*ele_val_at_quad(abscissa(j), ele)**(i-1)
+             end do
+          end do
+       end do
+    endif
+
+    !! ----------------
+
+!    perturb(1) = 0.00
+!    perturb(2) = 0.01
+!    perturb(3) = -0.01
+    ! check for ill-conditioned matrices
+    if (singular_option=='perturbate') then
+       do i = 1, ele_ngi(abscissa(1), ele)
+          call svd(A(i,:,:), svd_tmp1, SV, svd_tmp2)
+          p_num=1
+          do while (SV(size(SV))/SV(1) < cond)
+          !if (SV(size(SV))/SV(1) < cond) then
+             print*, "SINGULAR MATRIX"
+             print*, "condition number=", (SV(size(SV))/SV(1))
+             do j = 1, N
+                !! perturbation coefficient should be taken from diamond options with a check for default value
+                call random_number(perturb(1))
+                perturb(1)=-1.0+2.0*perturb(1)
+                print*, "perturb(1)=", N, perturb(1)
+                abscissa_val_at_quad(i,j) = abscissa_val_at_quad(i,j)*(1.0+perturb(1)*perturb_val)
+                print*, "abscissa", j, " = ", abscissa_val_at_quad(i,j)
+             end do
+             A = A_matrix(abscissa_val_at_quad)
+             call svd(A(i,:,:), svd_tmp1, SV, svd_tmp2)
+             print*, "perturbed condition number=", (SV(size(SV))/SV(1))
+             p_num=p_num+1
+             print*, "p_num=", p_num
+!             if (SV(size(SV))/SV(1) < cond) then
+!                print*, "perturbing does not help"
+!             end if
+!          end if
+          end do
+       end do
+    else if (singular_option=='set_source_to_zero') then
+       do i = 1, ele_ngi(abscissa(1), ele)
+          call svd(A(i,:,:), svd_tmp1, SV, svd_tmp2)
+          if (SV(size(SV))/SV(1) < cond) then
+             ewrite(2,*) 'ill-conditioned matrix found', SV(size(SV))/SV(1)
+             A(i,:,:) = 0.0
+             A_3(i,:,:) = 0.0
+             C(i,:) = 0.0
+             do j = 1, 2*N
+                A(i,j,j) = 1.0
+             end do
+          end if
+       end do
+    end if
+!!       p_num = 1
+
+
+!!       do j = 1, N
+!!          absc_init(j) = abscissa_val_at_quad(i,j)
+!!       end do
+
+!!       do while (SV(size(SV))/SV(1) < cond)
+!          call get_option("/timestepping/current_time", curr_time)   ! get the current simulation time
+          !! a better check for the first timestep instead of small value of the current time
+ !         if (curr_time < 1e-7) then   ! this block perturbates the abscissas in A matrix if current time is 0
+!!             print*, "perturbing %g \n", p_num
+!!             do j = 1, N
+                !! perturbation coefficient should be taken from diamond options with a check for default value
+!!                call random_number(perturb(1))
+!!                abscissa_val_at_quad(i,j) = absc_init(j) + perturb(1)
+!                abscissa_val_at_quad(:,j) = abscissa_val_at_quad(:,j) + (2-j)*(0.1) 
+!!                print*, "abscissa", j, " = ", abscissa_val_at_quad(i,j)
+!!             end do
+!!             A = A_matrix(abscissa_val_at_quad)
+!!             call svd(A(i,:,:), svd_tmp1, SV, svd_tmp2)
+!!             print*, "condition number=", (SV(size(SV))/SV(1))
+!!             p_num = p_num+1
+!             print*, "perturbed condition number=", (SV(size(SV))/SV(1))
+!             if (SV(size(SV))/SV(1) < cond) then
+!                print*, "perturbing does not help"
+!             end if
+
+!          else          
+!             ewrite(2,*) 'ill-conditioned matrix found'
+!             A(i,:,:) = 0.0
+!             A_3(i,:,:) = 0.0
+!             C(i,:) = 0.0
+!             do j = 1, 2*N
+!                A(i,j,j) = 1.0
+!             end do
+!          end if
+!       end if
+!!       end do
+!    end do
+           
+
+
+
+    ! solve linear system to find source values
+    do i = 1, ele_ngi(abscissa(1), ele)
+       b(:,1) = matmul(A_3(i,:,:), C(i,:)) + S_rhs(i,:)   !! gb 15-11-2012! added S_rhs term
+       call dqmom_solve(A(i,:,:), b, stat)
+       weight_S_at_quad(i,:) = b(:N,1)
+       abscissa_S_at_quad(i,:) = b(N+1:,1)
+    end do    
+
+    ! In the DG case we apply the inverse mass locally.
+    invmass = inverse(shape_shape(shape, shape, detwei))
+
+    ! integrate and add to source fields
+    do i = 1, N
+       weight_S_at_nodes = shape_rhs(shape, detwei* weight_S_at_quad(:,i))
+       abscissa_S_at_nodes = shape_rhs(shape, detwei* abscissa_S_at_quad(:,i))
+       if(continuity(abscissa(1))<0) then
+          weight_S_at_nodes = matmul(weight_S_at_nodes, invmass)
+          abscissa_S_at_nodes = matmul(abscissa_S_at_nodes, invmass)
+       end if
+       call addto(s_weight(i)%ptr, nodes, weight_S_at_nodes)
+       call addto(s_weighted_abscissa(i)%ptr, nodes, abscissa_S_at_nodes)
+    end do
+
+  end subroutine dqmom_calculate_source_term_ele
+
+  function A_matrix(abscissa)
+
+    real, dimension(:,:), intent(in) :: abscissa
+
+    real, dimension(size(abscissa,1), size(abscissa,2)*2, size(abscissa,2)*2) :: A_matrix
+    integer :: i, j, N
+
+    N = size(abscissa,2)        
+    do i = 1, 2*N
+       do j = 1, N
+          A_matrix(:,i,j) = (2-i)*abscissa(:,j)**(i-1)
+          A_matrix(:,i,j+N) = (i-1)*abscissa(:,j)**(i-2)
+       end do
+    end do
+
+  end function A_matrix
+
+  subroutine dqmom_solve(A, b, stat)
+    !!< Solve Ax=b for right hand sides b putting the result in b.
+    !!<
+    !!< This is simply a wrapper for lapack.
+    real, dimension(:,:), intent(in) :: A
+    real, dimension(:,:), intent(inout) :: b
+    integer, optional, intent(out) :: stat
+
+    real, dimension(size(A,1), size(A,2)) :: Atmp
+    integer, dimension(size(A,1)) :: ipiv
+    integer :: info
+
+    interface 
+#ifdef DOUBLEP
+       SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         INTEGER :: INFO, LDA, LDB, N, NRHS
+         INTEGER :: IPIV( * )
+         REAL ::  A( LDA, * ), B( LDB, * )
+       END SUBROUTINE DGESV
+#else
+       SUBROUTINE SGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         INTEGER :: INFO, LDA, LDB, N, NRHS
+         INTEGER :: IPIV( * )
+         REAL ::  A( LDA, * ), B( LDB, * )
+       END SUBROUTINE SGESV
+#endif
+    end interface
+
+    if (present(stat)) stat = 0
+
+    Atmp=A
+
+#ifdef DOUBLEP
+    call dgesv(&
+#else
+    call sgesv(&
+#endif
+    size(A,1), size(b,2), Atmp, size(A,1), ipiv, b, size(b,1), info) 
+    
+    if (present(stat)) then
+       stat = info
+    end if
+
+  end subroutine dqmom_solve
+
+  subroutine dqmom_calculate_moments(state)
+
+    type(state_type), intent(in) :: state
+
+    type(scalar_field), pointer :: abscissa, weight, moment
+    type(scalar_field) :: work
+    integer :: i_pop, N, i_moment, i_abscissa, i
+    character(len=OPTION_PATH_LEN) :: option_path 
+    character(len=FIELD_NAME_LEN) :: type
+    
+    ewrite(1, *) "In dqmom_calculate_moments"
+    do i_pop = 1, option_count(trim(state%option_path)//'/population_balance')
+       call get_pop_option_path(state, i_pop, option_path)
+       N = option_count(trim(option_path)//'/abscissa/scalar_field')
+       do i_moment = 1, 2*N
+          type = 'moments'
+          call get_pop_field(state, i_pop, i_moment, type, moment)
+          call zero(moment)
+          call allocate(work, moment%mesh, 'work')
+          do i_abscissa = 1, N
+             ! get required fields
+             type = 'abscissa'
+             call get_pop_field(state, i_pop, i_abscissa, type, abscissa)
+             type = 'weights'
+             call get_pop_field(state, i_pop, i_abscissa, type, weight)
+
+             ! calculate moment -> m(i) = sum_j(w(j)*q(j)**i)
+             call zero(work)
+             call set(work, 1.0)
+             do i = 1, i_moment - 1
+                call scale(work, abscissa)
+             end do
+             call scale(work, weight)
+             call addto(moment, work)
+          end do
+          call deallocate(work)
+       end do
+    end do
+    ewrite(1, *) "Exiting dqmom_calculate_moments"
+  end subroutine dqmom_calculate_moments
+
+  subroutine dqmom_calculate_statistics(state)
+
+    type(state_type), intent(in) :: state
+
+    type(scalar_field_pointer), dimension(:), allocatable :: moments
+    type(scalar_field), pointer :: stats
+    integer :: i_pop, N, i, j, i_stat
+    real :: mean, std, scaling_factor_Dia
+    character(len=OPTION_PATH_LEN) :: option_path 
+    character(len=FIELD_NAME_LEN) :: type
+    
+    ewrite(1, *) "In dqmom_calculate_statistics"
+    do i_pop = 1, option_count(trim(state%option_path)//'/population_balance')
+       call get_pop_option_path(state, i_pop, option_path)
+
+       if (option_count(trim(option_path)//'/statistics/scalar_field') > 0) then
+          
+          N = option_count(trim(option_path)//'/moments/scalar_field')
+          allocate(moments(N))
+          do i = 1, N
+             ! get required fields
+             type = 'moments'
+             call get_pop_field(state, i_pop, i, type, moments(i)%ptr)
+          end do
+
+          type = 'statistics'
+          do i_stat = 1, option_count(trim(option_path)//'/statistics/scalar_field')
+             call get_pop_field(state, i_pop, i_stat, type, stats)
+             if (trim(stats%name) == "Mean") then
+                call zero(stats)
+                do j = 1, node_count(stats)
+                   call set(stats, j, node_val(moments(2)%ptr,j)/node_val(moments(1)%ptr,j))
+                end do
+             end if
+             if (trim(stats%name) == "StandardDeviation") then
+                call zero(stats)
+                do j = 1, node_count(stats)
+                   mean = node_val(moments(2)%ptr,j)/node_val(moments(1)%ptr,j)
+                   call set(stats, j, (node_val(moments(3)%ptr,j)/node_val(moments(1)%ptr,j) -&
+                        mean**2)**0.5)
+                end do
+             end if
+             if (trim(stats%name) == "Skew") then
+                call zero(stats)
+                do j = 1, node_count(stats)
+                   mean = node_val(moments(2)%ptr,j)/node_val(moments(1)%ptr,j)
+                   std = (node_val(moments(3)%ptr,j)/node_val(moments(1)%ptr,j) - mean**2)**0.5
+                   call set(stats, j, &
+                        (node_val(moments(4)%ptr,j)/node_val(moments(1)%ptr,j) - &
+                        3*mean*std**2.0 - mean**3.0) / std**3.0)
+                end do
+             end if
+             if (trim(stats%name) == "SauterMeanDia") then
+                call zero(stats)
+                if (have_option(trim(option_path)//'/scaling_factor_Dia')) then
+                   call get_option(trim(option_path)//'/scaling_factor_Dia', scaling_factor_Dia)
+                else
+                   scaling_factor_Dia=1.0 
+                end if
+                do j = 1, node_count(stats)
+                   call set(stats, j, scaling_factor_Dia*(node_val(moments(4)%ptr,j)/node_val(moments(3)%ptr,j)) )
+                end do
+             end if
+             if (trim(stats%name) == "MeanDia10") then
+                call zero(stats)
+                if (have_option(trim(option_path)//'/scaling_factor_Dia')) then
+                   call get_option(trim(option_path)//'/scaling_factor_Dia', scaling_factor_Dia)
+                else
+                   scaling_factor_Dia=1.0
+                end if
+                do j = 1, node_count(stats)
+                   call set(stats, j, scaling_factor_Dia*(node_val(moments(2)%ptr,j)/node_val(moments(1)%ptr,j)) )
+                end do
+             end if
+
+          end do
+
+          deallocate(moments)
+
+       end if
+    end do
+    ewrite(1, *) "Exiting dqmom_calculate_statistics"
+
+
+  end subroutine dqmom_calculate_statistics
+
+  subroutine dqmom_check_options(state)
+
+    type(state_type), intent(in) :: state
+    
+    integer :: i_pop
+    integer :: n_abscissa, n_weights, n_weighted_abscissa, n_moments
+    character(len=OPTION_PATH_LEN) :: option_path 
+
+    write(*,*) 'in dqmom_check_options'
+
+    do i_pop = 1, option_count(trim(state%option_path)//'/population_balance')
+       option_path = trim(state%option_path)//'/population_balance['//int2str(i_pop)//']'
+
+       ! Check there are the same number of abscissa, weight and weighted abscissa fields
+       n_abscissa = option_count(trim(option_path)//'/abscissa/scalar_field')
+       n_weights = option_count(trim(option_path)//'/weights/scalar_field')
+       n_weighted_abscissa = option_count(trim(option_path)// &
+            '/weighted_abscissa/scalar_field')
+       assert(n_weights == n_abscissa) 
+       assert(n_weights == n_weighted_abscissa)
+
+       ! Check there are sufficient abscissa to calculate requested moments
+       n_moments = option_count(trim(option_path)//'/moments/scalar_field')
+       assert((n_moments / 2) == n_abscissa)       
+
+       ! Need to check that all fields are on the same mesh
+    end do
+
+  end subroutine dqmom_check_options
+
+end module dqmom

=== added file 'population_balance/Makefile.dependencies'
--- population_balance/Makefile.dependencies	1970-01-01 00:00:00 +0000
+++ population_balance/Makefile.dependencies	2013-12-12 22:35:40 +0000
@@ -0,0 +1,10 @@
+# Dependencies generated by create_makefile.py. DO NOT EDIT
+../include/dqmom.mod: DQMOM.o
+	@true
+
+DQMOM.o ../include/dqmom.mod: DQMOM.F90 ../include/fdebug.h \
+   ../include/fields.mod ../include/global_parameters.mod \
+   ../include/initialise_fields_module.mod ../include/spud.mod \
+   ../include/state_fields_module.mod ../include/state_module.mod \
+   ../include/vector_tools.mod
+

=== added file 'population_balance/Makefile.in'
--- population_balance/Makefile.in	1970-01-01 00:00:00 +0000
+++ population_balance/Makefile.in	2013-12-12 22:35:40 +0000
@@ -0,0 +1,66 @@
+#    Copyright (C) 2006 Imperial College London and others.
+#
+#    Please see the AUTHORS file in the main source directory for a full list
+#    of copyright holders.
+#
+#    Prof. C Pain
+#    Applied Modelling and Computation Group
+#    Department of Earth Science and Engineering
+#    Imperial College London
+#
+#    [e-mail]
+#
+#    This library is free software; you can redistribute it and/or
+#    modify it under the terms of the GNU Lesser General Public
+#    License as published by the Free Software Foundation,
+#    version 2.1 of the License.
+#
+#    This library is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with this library; if not, write to the Free Software
+#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+#    USA
+SHELL = @SHELL@
+
+PACKAGE_NAME = @PACKAGE_NAME@
+FLUIDITY = $(PACKAGE_NAME)
+
+FC      = @FC@
+FCFLAGS = -I../include @MOD_FLAG@../include @FCFLAGS@ -I/usr/local/include -I./
+
+CC      = @CC@
+CFLAGS  = @CFLAGS@ -I../include -I/usr/local/include
+
+CXX     = @CXX@
+CXXFLAGS= @CXXFLAGS@ -I../include -I/usr/local/include
+
+MAKE    = @MAKE@
+AR  = @AR@
+ARFLAGS = @ARFLAGS@
+
+LIB = ../lib/lib$(FLUIDITY).a
+
+OBJS = DQMOM.o \
+
+.SUFFIXES: .f90 .F90 .c .o .a
+
+.f90.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) -c $<
+.F90.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) $(GENFLAGS) -c $< 
+
+$(LIB): $(OBJS)
+	@mkdir -p ../lib
+
+default: $(LIB)
+
+clean:
+	rm -f *.o *.d *.mod
+
+include Makefile.dependencies

=== modified file 'preprocessor/Field_Priority_Lists.F90'
--- preprocessor/Field_Priority_Lists.F90	2013-02-01 17:50:12 +0000
+++ preprocessor/Field_Priority_Lists.F90	2013-12-12 22:35:40 +0000
@@ -133,6 +133,37 @@
              end do
           end if
 
+          ! added as hack for now - but this whole set up of fields could be way better!
+          ! prognostic pop balance fields - very limited applicability
+          if (have_option('/material_phase['//int2str(p)//']/population_balance/')) then
+             do f = 0, option_count('/material_phase['//int2str(p)//&
+                  ']/population_balance/weights/scalar_field') - 1
+                nsol=nsol+1
+                temp_field_optionpath_list(nsol)='/material_phase['//int2str(p)// &
+                     ']/population_balance/weights/scalar_field['//int2str(f)//']'
+                call get_option('/material_phase['//int2str(p)//&
+                     ']/population_balance/weights/scalar_field['//int2str(f)//&
+                     ']/name',temp_field_name_list(nsol))
+                call get_option('/material_phase['//int2str(p)//&
+                     ']/population_balance/weights/scalar_field['//int2str(f)//&
+                     ']/prognostic/priority', priority(nsol), default=0)
+                temp_field_state_list(nsol) = p+1
+             end do
+             do f = 0, option_count('/material_phase['//int2str(p)//&
+                  ']/population_balance/weighted_abscissa/scalar_field') - 1
+                nsol=nsol+1
+                temp_field_optionpath_list(nsol)='/material_phase['//int2str(p)// &
+                     ']/population_balance/weighted_abscissa/scalar_field['//int2str(f)//']'
+                call get_option('/material_phase['//int2str(p)//&
+                     ']/population_balance/weighted_abscissa/scalar_field['//int2str(f)//&
+                     ']/name',temp_field_name_list(nsol))
+                call get_option('/material_phase['//int2str(p)//&
+                     ']/population_balance/weighted_abscissa/scalar_field['//int2str(f)//&
+                     ']/prognostic/priority', priority(nsol), default=0)
+                temp_field_state_list(nsol) = p+1
+             end do
+          end if
+
           ! prognostic Mellor Yamada fields:
           if (have_option('/material_phase[' &
                //int2str(p)//']/subgridscale_parameterisations/Mellor_Yamada/scalar_field::KineticEnergy/prognostic')) then
@@ -315,6 +346,17 @@
              ntsol = ntsol + 1
           end if
        end do
+    
+       ! added as hack for now - but this whole set up of fields could be way better!
+       ! prognostic pop balance fields - very limited applicability
+       if (have_option('/material_phase['//int2str(p)//']/population_balance/')) then
+          ntsol = ntsol + &
+               option_count('/material_phase['//int2str(p)//&
+               ']/population_balance/weights/scalar_field') + &
+               option_count('/material_phase['//int2str(p)//&
+               ']/population_balance/weighted_abscissa/scalar_field')
+       end if
+
        ! prognostic scalar fields for Mellor Yamada:
        if (have_option('/material_phase[' &
             //int2str(p)//']/subgridscale_parameterisations/Mellor_Yamada/scalar_field::KineticEnergy/prognostic')) then

=== modified file 'preprocessor/Populate_State.F90'
--- preprocessor/Populate_State.F90	2013-10-28 06:46:48 +0000
+++ preprocessor/Populate_State.F90	2013-12-12 22:35:40 +0000
@@ -97,7 +97,7 @@
        
   !! A list of relative paths under /material_phase[i]
   !! that are searched for additional fields to be added.
-  character(len=OPTION_PATH_LEN), dimension(13) :: additional_fields_relative=&
+  character(len=OPTION_PATH_LEN), dimension(18) :: additional_fields_relative=&
        (/ &
        "/subgridscale_parameterisations/Mellor_Yamada                                                       ", &
        "/subgridscale_parameterisations/prescribed_diffusivity                                              ", &
@@ -111,7 +111,12 @@
        "/vector_field::Velocity/prognostic/spatial_discretisation/continuous_galerkin/les_model/dynamic_les ", &
        "/vector_field::Velocity/prognostic/equation::ShallowWater                                           ", &
        "/vector_field::Velocity/prognostic/equation::ShallowWater/bottom_drag                               ", &
-       "/vector_field::BedShearStress/diagnostic/calculation_method/velocity_gradient                       " &
+       "/vector_field::BedShearStress/diagnostic/calculation_method/velocity_gradient                       ", &
+       "/population_balance[#]/abscissa/                                                                    ", &
+       "/population_balance[#]/weights/                                                                     ", &
+       "/population_balance[#]/weighted_abscissa/                                                           ", &
+       "/population_balance[#]/moments/                                                                     ", &
+       "/population_balance[#]/statistics/                                                                  "  &
        /)
 
   !! Relative paths under a field that are searched for grandchildren
@@ -121,6 +126,22 @@
          &    "/spatial_discretisation/inner_element" &
          /)
 
+
+  !! Dynamic paths that are searched for fields
+  !! This allows for searching for field within paths that may branch several times
+  !! The index of any particular path should be replaced with #
+  character(len=OPTION_PATH_LEN), dimension(6):: &
+         dynamic_paths = (/&
+         &    "/material_phase[#]/equation_of_state/fluids/linear/        ", &
+         &    "/material_phase[#]/population_balance[#]/abscissa/         ", &
+         &    "/material_phase[#]/population_balance[#]/weights/          ", &
+         &    "/material_phase[#]/population_balance[#]/weighted_abscissa/", &
+         &    "/material_phase[#]/population_balance[#]/moments/          ", &
+         &    "/material_phase[#]/population_balance[#]/statistics/       " &
+         /)
+
+  character(len=OPTION_PATH_LEN), dimension(:), allocatable :: field_locations
+
 contains
 
 

=== modified file 'schemas/diagnostic_algorithms.rnc'
--- schemas/diagnostic_algorithms.rnc	2013-02-21 15:25:41 +0000
+++ schemas/diagnostic_algorithms.rnc	2013-12-12 22:35:40 +0000
@@ -1289,7 +1289,8 @@
          },
          ## The particle diameter
          element particle_diameter {
-            real
+            element constant { real }|
+            element use_scalar_field { anystring }
          },
          ## The name of the state representing the continuous phase
          element continuous_phase_name {

=== modified file 'schemas/diagnostic_algorithms.rng'
--- schemas/diagnostic_algorithms.rng	2013-02-21 15:25:41 +0000
+++ schemas/diagnostic_algorithms.rng	2013-12-12 22:35:40 +0000
@@ -1754,7 +1754,14 @@
       </element>
       <element name="particle_diameter">
         <a:documentation>The particle diameter</a:documentation>
-        <ref name="real"/>
+        <choice>
+          <element name="constant">
+            <ref name="real"/>
+          </element>
+          <element name="use_scalar_field">
+            <ref name="anystring"/>
+          </element>
+        </choice>
       </element>
       <element name="continuous_phase_name">
         <a:documentation>The name of the state representing the continuous phase</a:documentation>

=== modified file 'schemas/fluidity_options.rnc'
--- schemas/fluidity_options.rnc	2013-10-11 10:48:07 +0000
+++ schemas/fluidity_options.rnc	2013-12-12 22:35:40 +0000
@@ -1472,6 +1472,7 @@
             scalar_field_choice*,
             vector_field_choice*,
             tensor_field_choice*,
+            population_balance*,
             ## Parameters required to model spontaneous electrical potentials in porous media.
             element electrical_properties {
                (
@@ -1613,8 +1614,9 @@
             ## Phase interaction options for Fluidity's multiphase flow model
             element multiphase_properties {
                comment,
-               element particle_diameter { 
-                     real 
+               element particle_diameter {
+	          element constant { real }|
+		  element use_scalar_field { anystring }
                }?,
                ## If this is the fluid phase, 
                ## the effective conductivity is required 
@@ -1814,6 +1816,450 @@
       }
    )      
 
+population_balance = 
+   (
+      ##
+      element population_balance {
+         (
+            ## Determine initial weighted_abscissa and weights using the PD algorithm and the 
+            ## initial conditions set in the diagnostic moment fields.
+            ## 
+            ## In this case it does not matter what is set in the weighted_abscissa and weights
+            ## fields initial conditions option
+            element calculate_initial_conditions_from_moments {
+               empty
+            }|
+            ## Use initial conditions set in the weighted_abscissa and weights prognostic
+            ## fields
+            ## 
+            ## In this case it does not matter what is set in the moment field initial conditions 
+            ## options
+            element use_prognostic_field_initial_conditions {
+               empty
+            }
+         ),
+         attribute name { xsd:string },
+         element abscissa {
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Abscissa_0" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Abscissa_1" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { xsd:string },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }*
+         },
+         element weights {
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Weight_0" },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Weight_1" },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { xsd:string },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            }*
+         },
+         element weighted_abscissa {
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "WeightedAbscissa_0" },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "WeightedAbscissa_1" },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { xsd:string },
+               ## Field type
+               element prognostic {
+                  velocity_mesh_choice,
+                  prognostic_scalar_field
+               }
+            }*
+         },
+         ## This option only works if you are using a continuous mesh for the population balance scalar fields.
+         ## A DG mesh will use inverse mass calculated at element level.
+         ## In the population balance DQMOM algorithm all sources for the 
+         ## advection-diffusion equations are calculated at the quadrature points first 
+         ## and then interpolated back to the nodal points. 
+         ## For a continuous mesh, this interpolation is made by inverting the mass matrix. 
+         ## Mass lumping is fairly quick and can be used with 1st order continuous meshes. 
+         ## But the error is high for higher order continuous meshes and full mass matrix should be used.
+         element adv_diff_source_term_interpolation {
+            ## The mass matrix is lumped to obtain source fields at nodal points.
+            element use_mass_lumping { empty } |
+            ## Full mass marix is used and the linear system is solved iteratively to obtain source fields at nodal points.
+            element use_full_mass_matrix {
+               ## This block specifies the solver options for solving the linear system to interpolate advection-diffusion source terms at nodal points.
+               element solver {
+                  linear_solver_options_asym_scalar
+               }
+            }
+         },
+         element moments {
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Moment_0" },
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  (
+                     ## Initial condition for WholeMesh
+                     ##
+                     ## Only specify one condition if not using mesh regions.
+                     ## Otherwise select other initial_condition option, specify region_ids
+                     ## and distinct names.  Then add extra intial conditions for other regions.
+                     element initial_condition {
+                        attribute name { "WholeMesh" },
+                        input_choice_initial_condition_real
+                     }|
+                     ## Multiple initial_conditions are allowed if specifying
+                     ## different values in different
+                     ## regions of the mesh (defined by region_ids).  In this case
+                     ## each initial_condition
+                     ## requires a distinct name for the options dictionary.
+                     element initial_condition {
+                        attribute name { string },
+                        region_ids?,
+                        input_choice_initial_condition_real
+                     }
+                  )+,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Moment_1" },
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  (
+                     ## Initial condition for WholeMesh
+                     ##
+                     ## Only specify one condition if not using mesh regions.
+                     ## Otherwise select other initial_condition option, specify region_ids
+                     ## and distinct names.  Then add extra intial conditions for other regions.
+                     element initial_condition {
+                        attribute name { "WholeMesh" },
+                        input_choice_initial_condition_real
+                     }|
+                     ## Multiple initial_conditions are allowed if specifying
+                     ## different values in different
+                     ## regions of the mesh (defined by region_ids).  In this case
+                     ## each initial_condition
+                     ## requires a distinct name for the options dictionary.
+                     element initial_condition {
+                        attribute name { string },
+                        region_ids?,
+                        input_choice_initial_condition_real
+                     }
+                  )+,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Moment_2" },
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  (
+                     ## Initial condition for WholeMesh
+                     ##
+                     ## Only specify one condition if not using mesh regions.
+                     ## Otherwise select other initial_condition option, specify region_ids
+                     ## and distinct names.  Then add extra intial conditions for other regions.
+                     element initial_condition {
+                        attribute name { "WholeMesh" },
+                        input_choice_initial_condition_real
+                     }|
+                     ## Multiple initial_conditions are allowed if specifying
+                     ## different values in different
+                     ## regions of the mesh (defined by region_ids).  In this case
+                     ## each initial_condition
+                     ## requires a distinct name for the options dictionary.
+                     element initial_condition {
+                        attribute name { string },
+                        region_ids?,
+                        input_choice_initial_condition_real
+                     }
+                  )+,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Moment_3" },
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  (
+                     ## Initial condition for WholeMesh
+                     ##
+                     ## Only specify one condition if not using mesh regions.
+                     ## Otherwise select other initial_condition option, specify region_ids
+                     ## and distinct names.  Then add extra intial conditions for other regions.
+                     element initial_condition {
+                        attribute name { "WholeMesh" },
+                        input_choice_initial_condition_real
+                     }|
+                     ## Multiple initial_conditions are allowed if specifying
+                     ## different values in different
+                     ## regions of the mesh (defined by region_ids).  In this case
+                     ## each initial_condition
+                     ## requires a distinct name for the options dictionary.
+                     element initial_condition {
+                        attribute name { string },
+                        region_ids?,
+                        input_choice_initial_condition_real
+                     }
+                  )+,
+                  diagnostic_scalar_field
+               }
+            },
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { xsd:string },
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  (
+                     ## Initial condition for WholeMesh
+                     ##
+                     ## Only specify one condition if not using mesh regions.
+                     ## Otherwise select other initial_condition option, specify region_ids
+                     ## and distinct names.  Then add extra intial conditions for other regions.
+                     element initial_condition {
+                        attribute name { "WholeMesh" },
+                        input_choice_initial_condition_real
+                     }|
+                     ## Multiple initial_conditions are allowed if specifying
+                     ## different values in different
+                     ## regions of the mesh (defined by region_ids).  In this case
+                     ## each initial_condition
+                     ## requires a distinct name for the options dictionary.
+                     element initial_condition {
+                        attribute name { string },
+                        region_ids?,
+                        input_choice_initial_condition_real
+                     }
+                  )+,
+                  diagnostic_scalar_field
+               }
+            }*
+         },
+         element statistics {
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Mean" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }?,
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "StandardDeviation" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }?,
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "Skew" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }?,
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "SauterMeanDia" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }?,
+            element scalar_field {
+               attribute rank { "0" },
+               attribute name { "MeanDia10" },
+               ## Field type
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }
+            }?
+         },   
+
+         ## Source terms for the population balance equation like growth, 
+         ## breakage, aggegation and internal dispersion
+         element population_balance_source_terms {
+            ## Growth term in the internal coordinate space
+            element growth{ 
+               ## Enter the growth degree (r) for the power-law growth G(\xi)=\xi^{r}
+               element power_law_growth{ real }
+            }?,
+            ## Dispersion coefficient D. Currently it only works for constant dispersion and can't handle dispersion varying with internal coordinate
+            element internal_dispersion{ real }?,
+            ##
+            element aggregation{
+               ##
+               element aggregation_frequency{
+                  ##
+                  element constant_aggregation{ real }|
+                  ## aggregation frequency is given by \beta_{alpha gamma} = \xi_{alpha}^3 + \xi_{gamma}^3
+                  element hydrodynamic_aggregation{ empty }|
+                  ## aggregation frequency is given by \beta_{alpha gamma} = \beta_o (\xi_{alpha} + \xi_{gamma}). Enter the value of constant \beta_o below
+                  element sum_aggregation{ real }
+               }
+            }?,
+            ##
+            element breakage{
+               ##
+               element breakage_frequency{
+                  ##
+                  element constant_breakage{ real }|
+                  ## breakage frequency is given by a_{alpha} = p * \xi_{alpha}^r
+                  element power_law_breakage{
+                     ## coefficient p in breakage frequency a_{alpha} = p * \xi_{alpha}^r
+                     element coefficient{ real },
+                     ## degree r in breakage frequency a_{alpha} = p * \xi_{alpha}^r
+                     element degree{ real }
+                  }
+               },
+               element distribution_function{
+                  ##
+                  element symmetric_fragmentation{ empty }|
+                  ## McCoy and Madras (2003) used P(v|v') = 2/v' for v<v' and 0 otherwise, where v stands for volume. 
+                  ## When written in terms of particle diameter, this converts to P(x|x')=(6x^2)/x'^3 for x<x' and 0 otherwise.
+                  element mcCoy_madras_2003{ empty }
+               }
+            }?
+         },
+
+         ## At each vertex a linear system is solved to determine source terms
+         ## for the abscissa and weight transport equations. This requires inversion
+         ## of a matrix that has terms that are a function of the weighted abscissa.
+         ## If any two abscissa happen to lie very close to one another this can cause
+         ## this matrix to be ill-conditioned which can cause large errors or solver
+         ## failures when attempting to invert the matrix. There several options for 
+         ## evading this problem:
+         ##
+         ## (a) Set the abscissa and weight source terms to zero at the offending vertex.
+         ## (b) Perturbate the weighted absciss used to construct the matrix until the 
+         ## matrix is no longer ill conditioned and can be inverted.
+         ## (c) Use the source terms of neighbouring cells (not implemented)
+         element ill_conditioned_matrices {
+            ## The required condition number of the matrix. Recommended value of 1e-12
+            element required_condition_number{ real },
+            (
+               ## This is a better option when much of the domain has a minimum 0'th moment
+               element set_source_to_zero { empty }|
+               element perturbate{
+                  ## The amount to perturbate the weighted abscissa. The values are 
+                  ## iteratively perturbated by rand(-1,1)*perturbation until the matrix
+                  ## exceeds the required condition number. Too low a value can cause
+                  ## very long run times but a higher value will lead to larger errors.
+                  ## After every ten perturbations at any one grid location a message is 
+                  ## printed to the log file noting the difficulties in solving the linear 
+                  ## system.
+                  ## Recommended value of 1e-7
+                  element perturbation{ real }
+               }|
+               ## Do not check for ill-conditioning of matrix A
+               element do_nothing{ empty }
+            )
+         },
+         ## weights of zero can cause instabilities. It is recommended that a minimum
+         ## weight of approx 1e-10 is used.
+         element minimum_weight{ real },
+         ## Scaling factor to be used for the Sauter Mean Diameter / Mean Diameter calculated in the 
+         ## statistics. It is useful in a non-dimensional implementation of the PBE.
+         element scaling_factor_Dia{ real }?,
+         ## This is the submerged specific gravity, R, of this population balance
+         ## as a function of internal coordinate, Z. 
+         ## It will be used in the equation of state.
+         ## 
+         ## R(Z) = (rho_c(Z) - rho_a)/(rho_a)
+         ## 
+         ## Where: rho_c(Z) is the density of the dispersed phase and rho_a is 
+         ## the ambient fluid density
+         ##
+         ## Must be given as a python function in the form:
+         ## 
+         ##  def R(Z): 
+         ##     # Function code
+         ##     return # Return value of R
+         element submerged_specific_gravity {
+            python_code
+         }?,
+         element output_abscissa_and_weights{ empty }?,
+         element output_sources{ empty }?
+      }
+   )
+
 sediment = 
    (
       ## A sediment model. See the manual for details on how to use this model.

=== modified file 'schemas/fluidity_options.rng'
--- schemas/fluidity_options.rng	2013-11-06 21:30:29 +0000
+++ schemas/fluidity_options.rng	2013-12-12 22:35:40 +0000
@@ -1957,6 +1957,9 @@
           <zeroOrMore>
             <ref name="tensor_field_choice"/>
           </zeroOrMore>
+          <zeroOrMore>
+            <ref name="population_balance"/>
+          </zeroOrMore>
           <optional>
             <element name="electrical_properties">
               <a:documentation>Parameters required to model spontaneous electrical potentials in porous media.</a:documentation>
@@ -2191,7 +2194,14 @@
               <ref name="comment"/>
               <optional>
                 <element name="particle_diameter">
-                  <ref name="real"/>
+                  <choice>
+                    <element name="constant">
+                      <ref name="real"/>
+                    </element>
+                    <element name="use_scalar_field">
+                      <ref name="anystring"/>
+                    </element>
+                  </choice>
                 </element>
               </optional>
               <optional>
@@ -2547,6 +2557,644 @@
       </optional>
     </element>
   </start>
+  <define name="population_balance">
+    <element name="population_balance">
+      <a:documentation/>
+      <choice>
+        <element name="calculate_initial_conditions_from_moments">
+          <a:documentation>Determine initial weighted_abscissa and weights using the PD algorithm and the 
+initial conditions set in the diagnostic moment fields.
+
+In this case it does not matter what is set in the weighted_abscissa and weights
+fields initial conditions option</a:documentation>
+          <empty/>
+        </element>
+        <element name="use_prognostic_field_initial_conditions">
+          <a:documentation>Use initial conditions set in the weighted_abscissa and weights prognostic
+fields
+
+In this case it does not matter what is set in the moment field initial conditions 
+options</a:documentation>
+          <empty/>
+        </element>
+      </choice>
+      <attribute name="name">
+        <data type="string"/>
+      </attribute>
+      <element name="abscissa">
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Abscissa_0</value>
+          </attribute>
+          <element name="diagnostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Abscissa_1</value>
+          </attribute>
+          <element name="diagnostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <zeroOrMore>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <data type="string"/>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </zeroOrMore>
+      </element>
+      <element name="weights">
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Weight_0</value>
+          </attribute>
+          <element name="prognostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="prognostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Weight_1</value>
+          </attribute>
+          <element name="prognostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="prognostic_scalar_field"/>
+          </element>
+        </element>
+        <zeroOrMore>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <data type="string"/>
+            </attribute>
+            <element name="prognostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="prognostic_scalar_field"/>
+            </element>
+          </element>
+        </zeroOrMore>
+      </element>
+      <element name="weighted_abscissa">
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>WeightedAbscissa_0</value>
+          </attribute>
+          <element name="prognostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="prognostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>WeightedAbscissa_1</value>
+          </attribute>
+          <element name="prognostic">
+            <a:documentation>Field type</a:documentation>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="prognostic_scalar_field"/>
+          </element>
+        </element>
+        <zeroOrMore>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <data type="string"/>
+            </attribute>
+            <element name="prognostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="prognostic_scalar_field"/>
+            </element>
+          </element>
+        </zeroOrMore>
+      </element>
+      <element name="adv_diff_source_term_interpolation">
+        <a:documentation>This option only works if you are using a continuous mesh for the population balance scalar fields.
+A DG mesh will use inverse mass calculated at element level.
+In the population balance DQMOM algorithm all sources for the 
+advection-diffusion equations are calculated at the quadrature points first 
+and then interpolated back to the nodal points. 
+For a continuous mesh, this interpolation is made by inverting the mass matrix. 
+Mass lumping is fairly quick and can be used with 1st order continuous meshes. 
+But the error is high for higher order continuous meshes and full mass matrix should be used.</a:documentation>
+        <choice>
+          <element name="use_mass_lumping">
+            <a:documentation>The mass matrix is lumped to obtain source fields at nodal points.</a:documentation>
+            <empty/>
+          </element>
+          <element name="use_full_mass_matrix">
+            <a:documentation>Full mass marix is used and the linear system is solved iteratively to obtain source fields at nodal points.</a:documentation>
+            <element name="solver">
+              <a:documentation>This block specifies the solver options for solving the linear system to interpolate advection-diffusion source terms at nodal points.</a:documentation>
+              <ref name="linear_solver_options_asym_scalar"/>
+            </element>
+          </element>
+        </choice>
+      </element>
+      <element name="moments">
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Moment_0</value>
+          </attribute>
+          <element name="diagnostic">
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <oneOrMore>
+              <choice>
+                <element name="initial_condition">
+                  <a:documentation>Initial condition for WholeMesh
+
+Only specify one condition if not using mesh regions.
+Otherwise select other initial_condition option, specify region_ids
+and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
+                  <attribute name="name">
+                    <value>WholeMesh</value>
+                  </attribute>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+                <element name="initial_condition">
+                  <a:documentation>Multiple initial_conditions are allowed if specifying
+different values in different
+regions of the mesh (defined by region_ids).  In this case
+each initial_condition
+requires a distinct name for the options dictionary.</a:documentation>
+                  <attribute name="name">
+                    <data type="string" datatypeLibrary=""/>
+                  </attribute>
+                  <optional>
+                    <ref name="region_ids"/>
+                  </optional>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+              </choice>
+            </oneOrMore>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Moment_1</value>
+          </attribute>
+          <element name="diagnostic">
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <oneOrMore>
+              <choice>
+                <element name="initial_condition">
+                  <a:documentation>Initial condition for WholeMesh
+
+Only specify one condition if not using mesh regions.
+Otherwise select other initial_condition option, specify region_ids
+and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
+                  <attribute name="name">
+                    <value>WholeMesh</value>
+                  </attribute>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+                <element name="initial_condition">
+                  <a:documentation>Multiple initial_conditions are allowed if specifying
+different values in different
+regions of the mesh (defined by region_ids).  In this case
+each initial_condition
+requires a distinct name for the options dictionary.</a:documentation>
+                  <attribute name="name">
+                    <data type="string" datatypeLibrary=""/>
+                  </attribute>
+                  <optional>
+                    <ref name="region_ids"/>
+                  </optional>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+              </choice>
+            </oneOrMore>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Moment_2</value>
+          </attribute>
+          <element name="diagnostic">
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <oneOrMore>
+              <choice>
+                <element name="initial_condition">
+                  <a:documentation>Initial condition for WholeMesh
+
+Only specify one condition if not using mesh regions.
+Otherwise select other initial_condition option, specify region_ids
+and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
+                  <attribute name="name">
+                    <value>WholeMesh</value>
+                  </attribute>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+                <element name="initial_condition">
+                  <a:documentation>Multiple initial_conditions are allowed if specifying
+different values in different
+regions of the mesh (defined by region_ids).  In this case
+each initial_condition
+requires a distinct name for the options dictionary.</a:documentation>
+                  <attribute name="name">
+                    <data type="string" datatypeLibrary=""/>
+                  </attribute>
+                  <optional>
+                    <ref name="region_ids"/>
+                  </optional>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+              </choice>
+            </oneOrMore>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <element name="scalar_field">
+          <attribute name="rank">
+            <value>0</value>
+          </attribute>
+          <attribute name="name">
+            <value>Moment_3</value>
+          </attribute>
+          <element name="diagnostic">
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <oneOrMore>
+              <choice>
+                <element name="initial_condition">
+                  <a:documentation>Initial condition for WholeMesh
+
+Only specify one condition if not using mesh regions.
+Otherwise select other initial_condition option, specify region_ids
+and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
+                  <attribute name="name">
+                    <value>WholeMesh</value>
+                  </attribute>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+                <element name="initial_condition">
+                  <a:documentation>Multiple initial_conditions are allowed if specifying
+different values in different
+regions of the mesh (defined by region_ids).  In this case
+each initial_condition
+requires a distinct name for the options dictionary.</a:documentation>
+                  <attribute name="name">
+                    <data type="string" datatypeLibrary=""/>
+                  </attribute>
+                  <optional>
+                    <ref name="region_ids"/>
+                  </optional>
+                  <ref name="input_choice_initial_condition_real"/>
+                </element>
+              </choice>
+            </oneOrMore>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+        </element>
+        <zeroOrMore>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <data type="string"/>
+            </attribute>
+            <element name="diagnostic">
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <oneOrMore>
+                <choice>
+                  <element name="initial_condition">
+                    <a:documentation>Initial condition for WholeMesh
+
+Only specify one condition if not using mesh regions.
+Otherwise select other initial_condition option, specify region_ids
+and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
+                    <attribute name="name">
+                      <value>WholeMesh</value>
+                    </attribute>
+                    <ref name="input_choice_initial_condition_real"/>
+                  </element>
+                  <element name="initial_condition">
+                    <a:documentation>Multiple initial_conditions are allowed if specifying
+different values in different
+regions of the mesh (defined by region_ids).  In this case
+each initial_condition
+requires a distinct name for the options dictionary.</a:documentation>
+                    <attribute name="name">
+                      <data type="string" datatypeLibrary=""/>
+                    </attribute>
+                    <optional>
+                      <ref name="region_ids"/>
+                    </optional>
+                    <ref name="input_choice_initial_condition_real"/>
+                  </element>
+                </choice>
+              </oneOrMore>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </zeroOrMore>
+      </element>
+      <element name="statistics">
+        <optional>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>Mean</value>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>StandardDeviation</value>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>Skew</value>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>SauterMeanDia</value>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="scalar_field">
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>MeanDia10</value>
+            </attribute>
+            <element name="diagnostic">
+              <a:documentation>Field type</a:documentation>
+              <ref name="internal_algorithm"/>
+              <ref name="velocity_mesh_choice"/>
+              <ref name="diagnostic_scalar_field"/>
+            </element>
+          </element>
+        </optional>
+      </element>
+      <element name="population_balance_source_terms">
+        <a:documentation>Source terms for the population balance equation like growth, 
+breakage, aggegation and internal dispersion</a:documentation>
+        <optional>
+          <element name="growth">
+            <a:documentation>Growth term in the internal coordinate space</a:documentation>
+            <element name="power_law_growth">
+              <a:documentation>Enter the growth degree (r) for the power-law growth G(\xi)=\xi^{r}</a:documentation>
+              <ref name="real"/>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="internal_dispersion">
+            <a:documentation>Dispersion coefficient D. Currently it only works for constant dispersion and can't handle dispersion varying with internal coordinate</a:documentation>
+            <ref name="real"/>
+          </element>
+        </optional>
+        <optional>
+          <element name="aggregation">
+            <a:documentation/>
+            <element name="aggregation_frequency">
+              <a:documentation/>
+              <choice>
+                <element name="constant_aggregation">
+                  <a:documentation/>
+                  <ref name="real"/>
+                </element>
+                <element name="hydrodynamic_aggregation">
+                  <a:documentation>aggregation frequency is given by \beta_{alpha gamma} = \xi_{alpha}^3 + \xi_{gamma}^3</a:documentation>
+                  <empty/>
+                </element>
+                <element name="sum_aggregation">
+                  <a:documentation>aggregation frequency is given by \beta_{alpha gamma} = \beta_o (\xi_{alpha} + \xi_{gamma})</a:documentation>
+                  <ref name="real"/>
+                </element>
+              </choice>
+            </element>
+          </element>
+        </optional>
+        <optional>
+          <element name="breakage">
+            <a:documentation/>
+            <element name="breakage_frequency">
+              <a:documentation/>
+              <choice>
+                <element name="constant_breakage">
+                  <a:documentation/>
+                  <ref name="real"/>
+                </element>
+                <element name="power_law_breakage">
+                  <a:documentation>breakage frequency is given by a_{alpha} = p * \xi_{alpha}^r</a:documentation>
+                  <element name="coefficient">
+                    <a:documentation>coefficient p in breakage frequency a_{alpha} = p * \xi_{alpha}^r</a:documentation>
+                    <ref name="real"/>
+                  </element>
+                  <element name="degree">
+                    <a:documentation>degree r in breakage frequency a_{alpha} = p * \xi_{alpha}^r</a:documentation>
+                    <ref name="real"/>
+                  </element>
+                </element>
+              </choice>
+            </element>
+            <element name="distribution_function">
+              <choice>
+                <element name="symmetric_fragmentation">
+                  <a:documentation/>
+                  <empty/>
+                </element>
+                <element name="mcCoy_madras_2003">
+                  <a:documentation>McCoy and Madras (2003) used P(v|v') = 2/v' for v&lt;v' and 0 otherwise, where v stands for volume. 
+When written in terms of particle diameter, this converts to P(x|x')=(6x^2)/x'^3 for x&lt;x' and 0 otherwise.</a:documentation>
+                  <empty/>
+                </element>
+              </choice>
+            </element>
+          </element>
+        </optional>
+      </element>
+      <element name="ill_conditioned_matrices">
+        <a:documentation>At each vertex a linear system is solved to determine source terms
+for the abscissa and weight transport equations. This requires inversion
+of a matrix that has terms that are a function of the weighted abscissa.
+If any two abscissa happen to lie very close to one another this can cause
+this matrix to be ill-conditioned which can cause large errors or solver
+failures when attempting to invert the matrix. There several options for 
+evading this problem:
+
+(a) Set the abscissa and weight source terms to zero at the offending vertex.
+(b) Perturbate the weighted absciss used to construct the matrix until the 
+matrix is no longer ill conditioned and can be inverted.
+(c) Use the source terms of neighbouring cells (not implemented)</a:documentation>
+        <element name="required_condition_number">
+          <a:documentation>The required condition number of the matrix. Recommended value of 1e-12</a:documentation>
+          <ref name="real"/>
+        </element>
+        <choice>
+          <element name="set_source_to_zero">
+            <a:documentation>This is a better option when much of the domain has a minimum 0'th moment</a:documentation>
+            <empty/>
+          </element>
+          <element name="perturbate">
+            <element name="perturbation">
+              <a:documentation>The amount to perturbate the weighted abscissa. The values are 
+iteratively perturbated by rand(-1,1)*perturbation until the matrix
+exceeds the required condition number. Too low a value can cause
+very long run times but a higher value will lead to larger errors.
+After every ten perturbations at any one grid location a message is 
+printed to the log file noting the difficulties in solving the linear 
+system.
+Recommended value of 1e-7</a:documentation>
+              <ref name="real"/>
+            </element>
+          </element>
+          <element name="do_nothing">
+            <a:documentation>Do not check for ill-conditioning of matrix A</a:documentation>
+            <empty/>
+          </element>
+        </choice>
+      </element>
+      <element name="minimum_weight">
+        <a:documentation>weights of zero can cause instabilities. It is recommended that a minimum
+weight of approx 1e-10 is used.</a:documentation>
+        <ref name="real"/>
+      </element>
+      <optional>
+        <element name="scaling_factor_Dia">
+          <a:documentation>Scaling factor to be used for the Sauter Mean Diameter / Mean Diameter calculated in the 
+statistics. It is useful in a non-dimensional implementation of the PBE.</a:documentation>
+          <ref name="real"/>
+        </element>
+      </optional>
+      <optional>
+        <element name="submerged_specific_gravity">
+          <a:documentation>This is the submerged specific gravity, R, of this population balance
+as a function of internal coordinate, Z. 
+It will be used in the equation of state.
+
+R(Z) = (rho_c(Z) - rho_a)/(rho_a)
+
+Where: rho_c(Z) is the density of the dispersed phase and rho_a is 
+the ambient fluid density
+
+Must be given as a python function in the form:
+
+ def R(Z): 
+    # Function code
+    return # Return value of R</a:documentation>
+          <ref name="python_code"/>
+        </element>
+      </optional>
+      <optional>
+        <element name="output_abscissa_and_weights">
+          <empty/>
+        </element>
+      </optional>
+      <optional>
+        <element name="output_sources">
+          <empty/>
+        </element>
+      </optional>
+    </element>
+  </define>
   <define name="sediment">
     <element name="sediment">
       <a:documentation>A sediment model. See the manual for details on how to use this model.</a:documentation>

=== modified file 'schemas/multiphase_interaction.rnc'
--- schemas/multiphase_interaction.rnc	2012-06-15 19:32:03 +0000
+++ schemas/multiphase_interaction.rnc	2013-12-12 22:35:40 +0000
@@ -35,6 +35,11 @@
       element drag_correlation {
          attribute name { "ergun" },
          comment
+      }|
+      ## Fluid-particle drag term by Schiller-Naumann (1933).
+      element drag_correlation {
+         attribute name { "schiller_naumann" },
+         comment
       }
    )
    
@@ -55,4 +60,4 @@
          attribute name { "gunn" },
          comment
       }
-   )
\ No newline at end of file
+   )

=== modified file 'schemas/multiphase_interaction.rng'
--- schemas/multiphase_interaction.rng	2012-06-15 19:32:03 +0000
+++ schemas/multiphase_interaction.rng	2013-12-12 22:35:40 +0000
@@ -43,6 +43,13 @@
         </attribute>
         <ref name="comment"/>
       </element>
+      <element name="drag_correlation">
+        <a:documentation>Fluid-particle drag term by Schiller-Naumann (1933).</a:documentation>
+        <attribute name="name">
+          <value>schiller_naumann</value>
+        </attribute>
+        <ref name="comment"/>
+      </element>
     </choice>
   </define>
   <define name="heat_transfer">

=== modified file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml'
--- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml	2013-12-12 22:35:40 +0000
@@ -624,7 +624,9 @@
             <string_value lines="1">Dust::Velocity,Air::Velocity,Air::PhaseVolumeFraction,Air::Density</string_value>
           </depends>
           <particle_diameter>
-            <real_value rank="0">10e-6</real_value>
+            <constant>
+              <real_value rank="0">10e-6</real_value>
+            </constant>
           </particle_diameter>
           <continuous_phase_name>Air</continuous_phase_name>
         </algorithm>
@@ -726,7 +728,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">10e-6</real_value>
+        <constant>
+          <real_value rank="0">10e-6</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">758</real_value>

=== modified file 'tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml'
--- tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_inlet_velocity_bc_compressible/mphase_inlet_velocity_bc_compressible.flml	2013-12-12 22:35:40 +0000
@@ -845,7 +845,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">2e-4</real_value>
+        <constant>
+          <real_value rank="0">2e-4</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">954</real_value>

=== modified file 'tests/mphase_rogue_shock_tube_dense_bed_glass/mphase_rogue_shock_tube_dense_bed_glass.flml'
--- tests/mphase_rogue_shock_tube_dense_bed_glass/mphase_rogue_shock_tube_dense_bed_glass.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_rogue_shock_tube_dense_bed_glass/mphase_rogue_shock_tube_dense_bed_glass.flml	2013-12-12 22:35:40 +0000
@@ -842,7 +842,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">1.5e-3</real_value>
+        <constant>
+          <real_value rank="0">1.5e-3</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">758</real_value>

=== modified file 'tests/mphase_rogue_shock_tube_dense_bed_nylon/mphase_rogue_shock_tube_dense_bed_nylon.flml'
--- tests/mphase_rogue_shock_tube_dense_bed_nylon/mphase_rogue_shock_tube_dense_bed_nylon.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_rogue_shock_tube_dense_bed_nylon/mphase_rogue_shock_tube_dense_bed_nylon.flml	2013-12-12 22:35:40 +0000
@@ -842,7 +842,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">2.0e-3</real_value>
+        <constant>
+          <real_value rank="0">2.0e-3</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">758</real_value>

=== modified file 'tests/mphase_stokes_law/mphase_stokes_law.flml'
--- tests/mphase_stokes_law/mphase_stokes_law.flml	2012-05-19 22:54:15 +0000
+++ tests/mphase_stokes_law/mphase_stokes_law.flml	2013-12-12 22:35:40 +0000
@@ -502,7 +502,9 @@
             <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::Density</string_value>
           </depends>
           <particle_diameter>
-            <real_value rank="0">48e-6</real_value>
+            <constant>
+              <real_value rank="0">48e-6</real_value>
+            </constant>
           </particle_diameter>
           <continuous_phase_name>Water</continuous_phase_name>
         </algorithm>
@@ -676,7 +678,9 @@
     </vector_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">48e-6</real_value>
+        <constant>
+          <real_value rank="0">48e-6</real_value>
+        </constant>
       </particle_diameter>
     </multiphase_properties>
   </material_phase>

=== modified file 'tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml'
--- tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p0p1.flml	2013-12-12 22:35:40 +0000
@@ -725,7 +725,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">2e-4</real_value>
+        <constant>
+          <real_value rank="0">2.0e-4</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">954</real_value>

=== modified file 'tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml'
--- tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_strong_pressure_bc_compressible/mphase_strong_pressure_bc_compressible_p2p1.flml	2013-12-12 22:35:40 +0000
@@ -732,7 +732,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">2e-4</real_value>
+        <constant>
+          <real_value rank="0">2.0e-4</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">954</real_value>

=== modified file 'tests/mphase_subtract_out_reference_profile/mphase_subtract_out_reference_profile_1d.flml'
--- tests/mphase_subtract_out_reference_profile/mphase_subtract_out_reference_profile_1d.flml	2013-10-28 06:46:48 +0000
+++ tests/mphase_subtract_out_reference_profile/mphase_subtract_out_reference_profile_1d.flml	2013-12-12 22:35:40 +0000
@@ -811,7 +811,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">2e-4</real_value>
+        <constant>
+          <real_value rank="0">2e-4</real_value>
+        </constant>
       </particle_diameter>
       <specific_heat>
         <real_value rank="0">954</real_value>

=== modified file 'tests/mphase_tephra_settling/mphase_tephra_settling.flml'
--- tests/mphase_tephra_settling/mphase_tephra_settling.flml	2013-05-30 11:47:42 +0000
+++ tests/mphase_tephra_settling/mphase_tephra_settling.flml	2013-12-12 22:35:40 +0000
@@ -636,7 +636,9 @@
             <string_value lines="1">Water::Density,   Water::Velocity,  Water::PhaseVolumeFraction,  Tephra::Velocity</string_value>
           </depends>
           <particle_diameter>
-            <real_value rank="0">48e-6</real_value>
+            <constant>
+              <real_value rank="0">48e-6</real_value>
+            </constant>
           </particle_diameter>
           <continuous_phase_name>Water</continuous_phase_name>
         </algorithm>
@@ -656,7 +658,9 @@
     </scalar_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">48e-6</real_value>
+        <constant>
+          <real_value rank="0">48e-6</real_value>
+        </constant>
       </particle_diameter>
     </multiphase_properties>
   </material_phase>

=== modified file 'tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml'
--- tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml	2012-10-17 05:47:11 +0000
+++ tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml	2013-12-12 22:35:40 +0000
@@ -499,7 +499,9 @@
             <string_value lines="1">Water::Velocity,Tephra::Velocity,Water::Density</string_value>
           </depends>
           <particle_diameter>
-            <real_value rank="0">1e-3</real_value>
+            <constant>
+              <real_value rank="0">1e-3</real_value>
+            </constant>
           </particle_diameter>
           <continuous_phase_name>Water</continuous_phase_name>
         </algorithm>
@@ -781,7 +783,9 @@
     </vector_field>
     <multiphase_properties>
       <particle_diameter>
-        <real_value rank="0">1e-3</real_value>
+        <constant>
+          <real_value rank="0">1e-3</real_value>
+        </constant>
       </particle_diameter>
     </multiphase_properties>
   </material_phase>

=== added directory 'tests/popbal_homog_brk_aggr'
=== added file 'tests/popbal_homog_brk_aggr/Makefile'
--- tests/popbal_homog_brk_aggr/Makefile	1970-01-01 00:00:00 +0000
+++ tests/popbal_homog_brk_aggr/Makefile	2013-12-12 22:35:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+	gmsh -2 src/minimal.geo -optimize -bin
+
+clean:
+	rm -f  *.vtu *.pvtu fluidity.* *.s *.d.1 *.stat src/minimal.msh

=== added file 'tests/popbal_homog_brk_aggr/popbal_homog_N2.flml'
--- tests/popbal_homog_brk_aggr/popbal_homog_N2.flml	1970-01-01 00:00:00 +0000
+++ tests/popbal_homog_brk_aggr/popbal_homog_N2.flml	2013-12-12 22:35:40 +0000
@@ -0,0 +1,723 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">popbal_homog_N2_brk_aggr</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="src/minimal">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">5</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.001</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">5</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <vector_field name="Velocity" rank="1">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0 0</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <scalar_field name="m0_anal" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import tanh
+   phi_inf=10.0
+   phi_t = phi_inf*( (1.0+phi_inf*tanh(phi_inf*(t/2.0)))/(phi_inf+tanh(phi_inf*(t/2.0))) )
+   return phi_t</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="m1_anal" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import tanh, gamma, pi
+   phi_inf=10.0
+   phi_t = phi_inf*( (1.0+phi_inf*tanh(phi_inf*(t/2.0)))/(phi_inf+tanh(phi_inf*(t/2.0))) )
+   k1 = (pi/2.)*phi_t**2
+   k2 = (pi/6.)*phi_t
+   mom1 = (k1/(3.0*k2**(4./3.)))*gamma(4./3.)
+   return mom1</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="m2_anal" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import tanh, gamma, pi
+   phi_inf=10.0
+   phi_t = phi_inf*( (1.0+phi_inf*tanh(phi_inf*(t/2.0)))/(phi_inf+tanh(phi_inf*(t/2.0))) )
+   k1 = (pi/2.)*phi_t**2
+   k2 = (pi/6.)*phi_t
+   mom2 = (k1/(3.0*k2**(5./3.)))*gamma(5./3.)
+   return mom2</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="m3_anal" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import tanh, gamma, pi
+   phi_inf=10.0
+   phi_t = phi_inf*( (1.0+phi_inf*tanh(phi_inf*(t/2.0)))/(phi_inf+tanh(phi_inf*(t/2.0))) )
+   k1 = (pi/2.)*phi_t**2
+   k2 = (pi/6.)*phi_t
+   mom3 = k1/(3.*k2**2)
+   return mom3</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="diff_m0" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value lines="20" type="code" language="python">mom0_anal=state.scalar_fields["m0_anal"]
+mom0_dqmom=state.scalar_fields["Moment_0"]
+for i in range(field.node_count):
+  field.set(i, abs(mom0_dqmom.node_val(i)-mom0_anal.node_val(i))/abs(mom0_anal.node_val(i)))</string_value>
+          <depends>
+            <string_value lines="1">fluid::m0_anal</string_value>
+          </depends>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="diff_m1" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value lines="20" type="code" language="python">mom1_anal=state.scalar_fields["m1_anal"]
+mom1_dqmom=state.scalar_fields["Moment_1"]
+for i in range(field.node_count):
+  field.set(i, abs(mom1_dqmom.node_val(i)-mom1_anal.node_val(i))/abs(mom1_anal.node_val(i)))</string_value>
+          <depends>
+            <string_value lines="1">fluid::m1_anal</string_value>
+          </depends>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="diff_m2" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value lines="20" type="code" language="python">mom2_anal=state.scalar_fields["m2_anal"]
+mom2_dqmom=state.scalar_fields["Moment_2"]
+for i in range(field.node_count):
+  field.set(i, abs(mom2_dqmom.node_val(i)-mom2_anal.node_val(i))/abs(mom2_anal.node_val(i)))</string_value>
+          <depends>
+            <string_value lines="1">fluid::m2_anal</string_value>
+          </depends>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="diff_m3" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value lines="20" type="code" language="python">mom3_anal=state.scalar_fields["m3_anal"]
+mom3_dqmom=state.scalar_fields["Moment_3"]
+for i in range(field.node_count):
+  field.set(i, abs(mom3_dqmom.node_val(i)-mom3_anal.node_val(i))/abs(mom3_anal.node_val(i)))</string_value>
+          <depends>
+            <string_value lines="1">fluid::m3_anal</string_value>
+          </depends>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <population_balance name="DQMOM">
+      <calculate_initial_conditions_from_moments/>
+      <abscissa>
+        <scalar_field name="Abscissa_0" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Abscissa_1" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+      </abscissa>
+      <weights>
+        <scalar_field name="Weight_0" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms>
+                  <exclude_advection_terms/>
+                </advection_terms>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">0.99480625</real_value>
+              </constant>
+            </initial_condition>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+        <scalar_field name="Weight_1" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms>
+                  <exclude_advection_terms/>
+                </advection_terms>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">0.00519377</real_value>
+              </constant>
+            </initial_condition>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+      </weights>
+      <weighted_abscissa>
+        <scalar_field name="WeightedAbscissa_0" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms>
+                  <exclude_advection_terms/>
+                </advection_terms>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   return (0.99480625*1.13701522)</string_value>
+              </python>
+            </initial_condition>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+        <scalar_field name="WeightedAbscissa_1" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms>
+                  <exclude_advection_terms/>
+                </advection_terms>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   return (0.00519377*-4.46492922)</string_value>
+              </python>
+            </initial_condition>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+      </weighted_abscissa>
+      <adv_diff_source_term_interpolation>
+        <use_mass_lumping/>
+      </adv_diff_source_term_interpolation>
+      <moments>
+        <scalar_field name="Moment_0" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">1.0</real_value>
+              </constant>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_1" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import gamma, pi
+   mom1 = (6./pi)**(1./3.) * gamma(4./3.)
+   return mom1</string_value>
+              </python>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_2" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import gamma, pi
+   mom2 = (6./pi)**(2./3.) * gamma(5./3.)
+   return mom2</string_value>
+              </python>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_3" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X,t):
+   from math import pi
+   mom3 = (6./pi)
+   return mom3</string_value>
+              </python>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+          </diagnostic>
+        </scalar_field>
+      </moments>
+      <statistics/>
+      <population_balance_source_terms>
+        <aggregation>
+          <aggregation_frequency>
+            <constant_aggregation>
+              <real_value rank="0">1.0</real_value>
+            </constant_aggregation>
+          </aggregation_frequency>
+        </aggregation>
+        <breakage>
+          <breakage_frequency>
+            <power_law_breakage>
+              <coefficient>
+                <real_value rank="0">26.17993878</real_value>
+              </coefficient>
+              <degree>
+                <real_value rank="0">3</real_value>
+              </degree>
+            </power_law_breakage>
+          </breakage_frequency>
+          <distribution_function>
+            <mcCoy_madras_2003/>
+          </distribution_function>
+        </breakage>
+      </population_balance_source_terms>
+      <ill_conditioned_matrices>
+        <required_condition_number>
+          <real_value rank="0">1.0e-12</real_value>
+        </required_condition_number>
+        <set_source_to_zero/>
+      </ill_conditioned_matrices>
+      <minimum_weight>
+        <real_value rank="0">1e-10</real_value>
+      </minimum_weight>
+    </population_balance>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/popbal_homog_brk_aggr/popbal_homog_N2_brk_aggr.xml'
--- tests/popbal_homog_brk_aggr/popbal_homog_N2_brk_aggr.xml	1970-01-01 00:00:00 +0000
+++ tests/popbal_homog_brk_aggr/popbal_homog_N2_brk_aggr.xml	2013-12-12 22:35:40 +0000
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>Transient population balance using direct quadrature method of moments (N=2) for a homogeneous case with breakage and aggregation (McCoy and Madras (2003))</name>
+  <owner userid="gb812"/>
+  <tags>flml</tags>
+  <problem_definition length="medium" nprocs="1">
+    <command_line>fluidity popbal_homog_N2.flml</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="solvers_converged" language="python">import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files</variable>
+    <variable name="moment0_maxerror_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_homog_N2_brk_aggr.stat")
+moment0_maxerror_t_5s = s["fluid"]["diff_m0"]["max"][-1]</variable>
+    <variable name="moment1_maxerror_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_homog_N2_brk_aggr.stat")
+moment1_maxerror_t_5s = s["fluid"]["diff_m1"]["max"][-1]</variable>
+    <variable name="moment2_maxerror_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_homog_N2_brk_aggr.stat")
+moment2_maxerror_t_5s = s["fluid"]["diff_m2"]["max"][-1]</variable>
+    <variable name="moment3_maxerror_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_homog_N2_brk_aggr.stat")
+moment3_maxerror_t_5s = s["fluid"]["diff_m3"]["max"][-1]</variable>
+  </variables>
+  <pass_tests>
+    <test name="Solvers converged" language="python">assert(solvers_converged)</test>
+    <test name="moment0 max percent error at 5s less than 0.006" language="python">assert moment0_maxerror_t_5s &lt; 0.006</test>
+    <test name="moment1 max percent error at 5s less than 0.003" language="python">assert moment1_maxerror_t_5s &lt; 0.003</test>
+    <test name="moment2 max percent error at 5s less than 0.003" language="python">assert moment2_maxerror_t_5s &lt; 0.003</test>
+    <test name="moment3 max percent error at 5s less than 0.01" language="python">assert moment3_maxerror_t_5s &lt; 0.01</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added directory 'tests/popbal_homog_brk_aggr/src'
=== added file 'tests/popbal_homog_brk_aggr/src/minimal.geo'
--- tests/popbal_homog_brk_aggr/src/minimal.geo	1970-01-01 00:00:00 +0000
+++ tests/popbal_homog_brk_aggr/src/minimal.geo	2013-12-12 22:35:40 +0000
@@ -0,0 +1,28 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+7
+1 0 0 0
+2 1 0 0
+3 1 1 0
+4 0.4999999999990216 0 0
+5 1 0.4999999999990216 0
+6 0.5000000000013878 0.5000000000013878 0
+7 0.7499999999999051 0.2499999999999052 0
+$EndNodes
+$Elements
+12
+1 1 2 6 1 1 4
+2 1 2 6 1 4 2
+3 1 2 6 2 2 5
+4 1 2 6 2 5 3
+5 1 2 6 3 3 6
+6 1 2 6 3 6 1
[phone]
[phone]
[phone]
[phone]
[phone]
[phone]
+$EndElements

=== added directory 'tests/popbal_nonhomog2_2d_adapt'
=== added file 'tests/popbal_nonhomog2_2d_adapt/Makefile'
--- tests/popbal_nonhomog2_2d_adapt/Makefile	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_adapt/Makefile	2013-12-12 22:35:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+	gmsh -2 src/channel.geo -optimize
+
+clean:
+	rm -f  *.vtu *.pvtu fluidity.* *.s *.d.1 *.stat src/channel.msh

=== added file 'tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.flml'
--- tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.flml	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.flml	2013-12-12 22:35:40 +0000
@@ -0,0 +1,670 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">popbal_nonhomog2_2d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="./src/channel">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>continuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.005</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">20</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <vector_field name="Velocity" rank="1">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0 0</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <scalar_field name="m3_analytical" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp, pi, cos
+  if (abs(t-round(t))&lt;=0.01):
+     nu=0.001
+     L=1.0
+     n=100
+     phi_n=0.0
+     i=1
+     while i&lt;=n:
+         j=2*i
+         val = (1.0/j**2)*exp(-(nu*j**2*pi**2*t)/(L**2))*cos((j*pi*X[0])/L)
+         phi_n = phi_n + val
+         i=i+1
+     phi_n = phi_n*(-4.0/pi**2) + 1.0/6.0
+     return phi_n
+  else:
+     return 0.0</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ScalarAbsoluteDifference" rank="0">
+      <diagnostic field_name_b="m3_analytical" field_name_a="Moment_3">
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <population_balance name="DQMOM">
+      <use_prognostic_field_initial_conditions/>
+      <abscissa>
+        <scalar_field name="Abscissa_0" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Abscissa_1" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+      </abscissa>
+      <weights>
+        <scalar_field name="Weight_0" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms/>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">0.5</real_value>
+              </constant>
+            </initial_condition>
+            <tensor_field name="Diffusivity" rank="2">
+              <prescribed>
+                <value name="WholeMesh">
+                  <isotropic>
+                    <constant>
+                      <real_value rank="0">1e-3</real_value>
+                    </constant>
+                  </isotropic>
+                </value>
+                <output/>
+              </prescribed>
+            </tensor_field>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+        <scalar_field name="Weight_1" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms/>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">0.5</real_value>
+              </constant>
+            </initial_condition>
+            <tensor_field name="Diffusivity" rank="2">
+              <prescribed>
+                <value name="WholeMesh">
+                  <isotropic>
+                    <constant>
+                      <real_value rank="0">1e-3</real_value>
+                    </constant>
+                  </isotropic>
+                </value>
+                <output/>
+              </prescribed>
+            </tensor_field>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+      </weights>
+      <weighted_abscissa>
+        <scalar_field name="WeightedAbscissa_0" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms/>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X, t):
+  zeta1=0.5*  (1.0 + 2.0*X[0]*(1-X[0]))**(1.0/3.0)
+  return zeta1</string_value>
+              </python>
+            </initial_condition>
+            <tensor_field name="Diffusivity" rank="2">
+              <prescribed>
+                <value name="WholeMesh">
+                  <isotropic>
+                    <constant>
+                      <real_value rank="0">1e-3</real_value>
+                    </constant>
+                  </isotropic>
+                </value>
+                <output/>
+              </prescribed>
+            </tensor_field>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+        <scalar_field name="WeightedAbscissa_1" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms/>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>
+              <never_ignore_solver_failures/>
+              <diagnostics>
+                <monitors/>
+              </diagnostics>
+            </solver>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X, t):
+  zeta2=0.5*(-1.0)
+  return zeta2</string_value>
+              </python>
+            </initial_condition>
+            <tensor_field name="Diffusivity" rank="2">
+              <prescribed>
+                <value name="WholeMesh">
+                  <isotropic>
+                    <constant>
+                      <real_value rank="0">1e-3</real_value>
+                    </constant>
+                  </isotropic>
+                </value>
+                <output/>
+              </prescribed>
+            </tensor_field>
+            <scalar_field name="Source" rank="0">
+              <diagnostic>
+                <algorithm name="Internal" material_phase_support="multiple"/>
+                <output/>
+                <stat/>
+                <detectors>
+                  <include_in_detectors/>
+                </detectors>
+              </diagnostic>
+            </scalar_field>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </prognostic>
+        </scalar_field>
+      </weighted_abscissa>
+      <adv_diff_source_term_interpolation>
+        <use_mass_lumping/>
+      </adv_diff_source_term_interpolation>
+      <moments>
+        <scalar_field name="Moment_0" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X, t):
+  if X[0] &lt; 0.5:
+    return 0.6
+  else:
+    return 0.4</string_value>
+              </python>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_1" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <python>
+                <string_value lines="20" type="code" language="python">def val(X, t):
+  if X[0] &lt; 0.5:
+    return 1.0
+  else:
+    return 1.5</string_value>
+              </python>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_2" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">4.0</real_value>
+              </constant>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Moment_3" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <initial_condition name="WholeMesh">
+              <constant>
+                <real_value rank="0">12.0</real_value>
+              </constant>
+            </initial_condition>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <adaptivity_options>
+              <absolute_measure>
+                <scalar_field name="InterpolationErrorBound" rank="0">
+                  <prescribed>
+                    <value name="WholeMesh">
+                      <constant>
+                        <real_value rank="0">0.001</real_value>
+                      </constant>
+                    </value>
+                    <output/>
+                    <stat/>
+                    <detectors>
+                      <exclude_from_detectors/>
+                    </detectors>
+                  </prescribed>
+                </scalar_field>
+                <p_norm>
+                  <integer_value rank="0">2</integer_value>
+                </p_norm>
+              </absolute_measure>
+            </adaptivity_options>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+      </moments>
+      <statistics/>
+      <population_balance_source_terms/>
+      <ill_conditioned_matrices>
+        <required_condition_number>
+          <real_value rank="0">1e-12</real_value>
+        </required_condition_number>
+        <set_source_to_zero/>
+      </ill_conditioned_matrices>
+      <minimum_weight>
+        <real_value rank="0">1e-10</real_value>
+      </minimum_weight>
+    </population_balance>
+  </material_phase>
+  <mesh_adaptivity>
+    <hr_adaptivity>
+      <period_in_timesteps>
+        <integer_value rank="0">20</integer_value>
+      </period_in_timesteps>
+      <maximum_number_of_nodes>
+        <integer_value rank="0">20000</integer_value>
+      </maximum_number_of_nodes>
+      <enable_gradation/>
+      <tensor_field name="MinimumEdgeLengths">
+        <anisotropic_symmetric>
+          <constant>
+            <real_value symmetric="true" dim2="dim" shape="2 2" dim1="dim" rank="2">0.001 0.00 0.00 0.001</real_value>
+          </constant>
+        </anisotropic_symmetric>
+      </tensor_field>
+      <tensor_field name="MaximumEdgeLengths">
+        <anisotropic_symmetric>
+          <constant>
+            <real_value symmetric="true" dim2="dim" shape="2 2" dim1="dim" rank="2">0.1 0.0 0.0 0.1</real_value>
+          </constant>
+        </anisotropic_symmetric>
+      </tensor_field>
+      <adapt_at_first_timestep>
+        <number_of_adapts>
+          <integer_value rank="0">4</integer_value>
+        </number_of_adapts>
+      </adapt_at_first_timestep>
+      <debug/>
+    </hr_adaptivity>
+  </mesh_adaptivity>
+</fluidity_options>

=== added file 'tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.xml'
--- tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.xml	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_adapt/popbal_nonhomog2_2d_adapt.xml	2013-12-12 22:35:40 +0000
@@ -0,0 +1,46 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>Transient population balance using direct quadrature method of moments (N=2) for 2d mesh using P1 CG using ADAPTIVITY</name>
+  <owner userid="gb812"/>
+  <tags>flml</tags>
+  <problem_definition length="medium" nprocs="1">
+    <command_line>fluidity popbal_nonhomog2_2d_adapt.flml</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="solvers_converged" language="python">import os
+files = os.listdir("./")
+solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files</variable>
+    <variable name="moment3_L2error_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_nonhomog2_2d.stat")
+timestep=0.005
+t=5.0
+index_time_for_norm=((1.0/timestep)*t)-1
+moment3_L2error_t_5s = s["fluid"]["ScalarAbsoluteDifference"]["l2norm"][index_time_for_norm]</variable>
+    <variable name="moment3_L2error_t_20s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_nonhomog2_2d.stat")
+timestep=0.005
+t=20.0
+index_time_for_norm=((1.0/timestep)*t)-1
+moment3_L2error_t_20s = s["fluid"]["ScalarAbsoluteDifference"]["l2norm"][index_time_for_norm]</variable>
+    <variable name="moment3_maxerror_t_5s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_nonhomog2_2d.stat")
+timestep=0.005
+t=5.0
+index_time_for_norm=((1.0/timestep)*t)-1
+moment3_maxerror_t_5s = s["fluid"]["ScalarAbsoluteDifference"]["max"][index_time_for_norm]</variable>
+    <variable name="moment3_maxerror_t_20s" language="python">from fluidity_tools import stat_parser
+s = stat_parser("popbal_nonhomog2_2d.stat")
+timestep=0.005
+t=20.0
+index_time_for_norm=((1.0/timestep)*t)-1
+moment3_maxerror_t_20s = s["fluid"]["ScalarAbsoluteDifference"]["max"][index_time_for_norm]</variable>
+  </variables>
+  <pass_tests>
+    <test name="Solvers converged" language="python">assert(solvers_converged)</test>
+    <test name="moment3 L2error at 5s less than 0.0001" language="python">assert moment3_L2error_t_5s &lt; 0.0001</test>
+    <test name="moment3 L2error at 20s less than 0.0001" language="python">assert moment3_L2error_t_20s &lt; 0.0001</test>
+    <test name="moment3 maxerror at 5s less than 0.001" language="python">assert moment3_maxerror_t_5s &lt; 0.001</test>
+    <test name="moment3 maxerror at 20s less than 0.001" language="python">assert moment3_maxerror_t_20s &lt; 0.001</test>
+  </pass_tests>
+  <warn_tests/>
+</testproblem>

=== added file 'tests/popbal_nonhomog2_2d_adapt/readme'
--- tests/popbal_nonhomog2_2d_adapt/readme	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_adapt/readme	2013-12-12 22:35:40 +0000
@@ -0,0 +1,13 @@
+quad=3
+2d
+CG
+p1
+unstructured grid cl=1/25
+lumped mass
+time=0:0.005:20
+vel=0, transient external diffusion only. Diffusion 2 with a parabolic initial condition
+only moment 3 compared to the analytical case
+analytical case derived using separation of variables
+N=2
+tests for L2error and max error at 5s and 20s
+adaptivity

=== added directory 'tests/popbal_nonhomog2_2d_adapt/src'
=== added file 'tests/popbal_nonhomog2_2d_adapt/src/channel.geo'
--- tests/popbal_nonhomog2_2d_adapt/src/channel.geo	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_adapt/src/channel.geo	2013-12-12 22:35:40 +0000
@@ -0,0 +1,16 @@
+cl1 = 1.0/25.0;
+Point(1) = {0, 0, 0, cl1};
+Point(2) = {1, 0, 0, cl1};
+Point(3) = {0, 1, 0, cl1};
+Point(4) = {1, 1, 0, cl1};
+Line(1) = {1, 2};
+Line(2) = {3, 4};
+Line(3) = {1, 3};
+Line(4) = {2, 4};
+Line Loop(5) = {1, 4, -2, -3};
+Ruled Surface(5) = {5};
+Physical Line(1) = {1};
+Physical Line(2) = {2};
+Physical Line(3) = {3};
+Physical Line(4) = {4};
+Physical Surface(6) = {5};

=== added directory 'tests/popbal_nonhomog2_2d_cg_p1'
=== added file 'tests/popbal_nonhomog2_2d_cg_p1/Makefile'
--- tests/popbal_nonhomog2_2d_cg_p1/Makefile	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_cg_p1/Makefile	2013-12-12 22:35:40 +0000
@@ -0,0 +1,5 @@
+input: clean
+	gmsh -2 src/channel.geo -optimize -bin
+
+clean:
+	rm -f  *.vtu *.pvtu fluidity.* *.s *.d.1 *.stat src/channel.msh

=== added file 'tests/popbal_nonhomog2_2d_cg_p1/popbal_nonhomog2_2d.flml'
--- tests/popbal_nonhomog2_2d_cg_p1/popbal_nonhomog2_2d.flml	1970-01-01 00:00:00 +0000
+++ tests/popbal_nonhomog2_2d_cg_p1/popbal_nonhomog2_2d.flml	2013-12-12 22:35:40 +0000
@@ -0,0 +1,614 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">popbal_nonhomog2_2d</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="./src/channel">
+        <format name="gmsh"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">1</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <mesh_continuity>
+          <string_value>continuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">3</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">0.005</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">20</real_value>
+    </finish_time>
+  </timestepping>
+  <material_phase name="fluid">
+    <vector_field name="Velocity" rank="1">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0 0</real_value>
+          </constant>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <scalar_field name="m3_analytical" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+  from math import exp, pi, cos
+  if (abs(t-round(t))&lt;=0.01):
+     nu=0.001
+     L=1.0
+     n=100
+     phi_n=0.0
+     i=1
+     while i&lt;=n:
+         j=2*i
+         val = (1.0/j**2)*exp(-(nu*j**2*pi**2*t)/(L**2))*cos((j*pi*X[0])/L)
+         phi_n = phi_n + val
+         i=i+1
+     phi_n = phi_n*(-4.0/pi**2) + 1.0/6.0
+     return phi_n
+  else:
+     return 0.0</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="ScalarAbsoluteDifference" rank="0">
+      <diagnostic field_name_b="m3_analytical" field_name_a="Moment_3">
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <exclude_from_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <population_balance name="DQMOM">
+      <use_prognostic_field_initial_conditions/>
+      <abscissa>
+        <scalar_field name="Abscissa_0" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+        <scalar_field name="Abscissa_1" rank="0">
+          <diagnostic>
+            <algorithm name="Internal" material_phase_support="multiple"/>
+            <mesh name="VelocityMesh"/>
+            <output/>
+            <stat/>
+            <convergence>
+              <include_in_convergence/>
+            </convergence>
+            <detectors>
+              <include_in_detectors/>
+            </detectors>
+            <steady_state>
+              <include_in_steady_state/>
+            </steady_state>
+            <consistent_interpolation/>
+          </diagnostic>
+        </scalar_field>
+      </abscissa>
+      <weights>
+        <scalar_field name="Weight_0" rank="0">
+          <prognostic>
+            <mesh name="VelocityMesh"/>
+            <equation name="AdvectionDiffusion"/>
+            <spatial_discretisation>
+              <continuous_galerkin>
+                <stabilisation>
+                  <no_stabilisation/>
+                </stabilisation>
+                <advection_terms/>
+                <mass_terms/>
+              </continuous_galerkin>
+              <conservative_advection>
+                <real_value rank="0">0.0</real_value>
+              </conservative_advection>
+            </spatial_discretisation>
+            <temporal_discretisation>
+              <theta>
+                <real_value rank="0">0.5</real_value>
+              </theta>
+            </temporal_discretisation>
+            <solver>
+              <iterative_method name="gmres">
+                <restart>
+                  <integer_value rank="0">30</integer_value>
+                </restart>
+              </iterative_method>
+              <preconditioner name="sor"/>
+              <relative_error>
+                <real_value rank="0">1e-7</real_value>
+              </relative_error>
+              <max_iterations>
+                <integer_value rank="0">1000</integer_value>
+              </max_iterations>