Merge lp:~fluidity-core/fluidity/fieldweightedpartition into lp:fluidity

[Rhodri Davies]

This branch includes changes that allow one to weight mesh partitions based upon a user prescribed scalar field. Based upon this information, flredecomp will try to load balance in such a way that the sum of user prescribed weights in each partition is equal.

I have added a test for this new functionality. Results look good.

I have also updated the manual to describe this change.

The buildbot branch is here:
http://buildbot-ocean.ese.ic.ac.uk:8080/builders/cap_fields

Will update when green (hopefully).

As usual, a big fat thank you and man hug to my office hommie skramer!!!

[lp:~fluidity-core/fluidity/fieldweightedpartition]

4056. By Rhodri Davies on 2012-09-11

Correct flredecomp path

[Rhodri Davies (rhodri-davies)]

Update: Buildbot is as green as Wales' beautiful valleys!

[lp:~fluidity-core/fluidity/fieldweightedpartition]

4057. By Rhodri Davies on 2012-09-12

Merge trunk change into branch

[Jon Hill (jon-hill)]

Looks good. Go for it.

review: Approve

[lp:~fluidity-core/fluidity/fieldweightedpartition]

4058. By Rhodri Davies on 2012-09-18

Merge trunk changes into branch

=== modified file 'assemble/Zoltan_callbacks.F90'
--- assemble/Zoltan_callbacks.F90	2011-10-18 09:36:51 +0000
+++ assemble/Zoltan_callbacks.F90	2012-09-18 19:42:23 +0000
@@ -64,7 +64,6 @@
     ierr = ZOLTAN_OK
   end function zoltan_cb_owned_node_count
 
-
   subroutine zoltan_cb_get_owned_nodes(data, num_gid_entries, num_lid_entries, global_ids, local_ids, wgt_dim, obj_wgts, ierr)
     integer(zoltan_int), dimension(*), intent(in) :: data ! not used
     integer(zoltan_int), intent(in) :: num_gid_entries, num_lid_entries 
@@ -75,7 +74,7 @@
     integer(zoltan_int), intent(out) :: ierr
     
     integer :: count, i
-    real(zoltan_float) :: max_obj_wgt
+    real(zoltan_float) :: max_obj_wgt, min_obj_wgt
     
     ewrite(1,*) "In zoltan_cb_get_owned_nodes"
     
@@ -96,24 +95,38 @@
     if(zoltan_global_migrate_extruded_mesh) then
        ! weight the nodes according to the number of nodes in the column beneath it
        max_obj_wgt = 1.0
-       do i=1,count
+       do i = 1, count
           obj_wgts(i) = float(row_length(zoltan_global_columns_sparsity, i))
           max_obj_wgt = max(max_obj_wgt, obj_wgts(i))
        end do
        ! normalise according to the most nodes in a column
-       do i=1,count
+       do i = 1, count
           obj_wgts(i) = obj_wgts(i)/max_obj_wgt
        end do
     else
-       do i=1,count
+       do i = 1, count
           obj_wgts(i) = 1.0
        end do
     end if
+
+    if(zoltan_global_field_weighted_partitions) then
+       max_obj_wgt = 1.0
+       min_obj_wgt = 0.0
+       do i = 1, count
+          obj_wgts(i) = node_val(zoltan_global_field_weighted_partition_values,i)
+          max_obj_wgt = max(max_obj_wgt, obj_wgts(i))
+          min_obj_wgt = min(min_obj_wgt, obj_wgts(i))
+       end do
+
+       if((max_obj_wgt > 1.0) .OR. (min_obj_wgt < 0.0)) then
+        FLExit("0.0 <= FieldWeightedPartitionValues <= 1.0: condition not satisfied")
+       end if
+
+    end if
     
     ierr = ZOLTAN_OK
   end subroutine zoltan_cb_get_owned_nodes
 
-
   subroutine zoltan_cb_get_num_edges(data, num_gid_entries, num_lid_entries, num_obj, global_ids, local_ids, num_edges, ierr)  
     integer(zoltan_int), dimension(*), intent(in) :: data 
     integer(zoltan_int), intent(in) :: num_gid_entries, num_lid_entries, num_obj
@@ -150,7 +163,6 @@
     ierr = ZOLTAN_OK
   end subroutine zoltan_cb_get_num_edges
 
-
   subroutine zoltan_cb_get_edge_list(data, num_gid_entries, num_lid_entries, num_obj, global_ids, local_ids, &
        &  num_edges, nbor_global_id, nbor_procs, wgt_dim, ewgts, ierr)
     integer(zoltan_int), intent(in) :: data 
 
=== modified file 'assemble/Zoltan_global_variables.F90'
--- assemble/Zoltan_global_variables.F90	2011-09-26 17:04:31 +0000
+++ assemble/Zoltan_global_variables.F90	2012-09-18 19:42:23 +0000
@@ -34,11 +34,12 @@
   ! Needed for zoltan_cb_get_owned_nodes
   type(csr_sparsity), save :: zoltan_global_columns_sparsity
   logical, save :: zoltan_global_migrate_extruded_mesh
+  logical, save :: zoltan_global_field_weighted_partitions
+  type(scalar_field), save :: zoltan_global_field_weighted_partition_values
 
   ! Needed for zoltan_cb_get_num_edges
   type(csr_sparsity), save, pointer :: zoltan_global_zz_sparsity_one
 
-
   ! Needed for zoltan_cb_get_edge_list
   logical, save :: zoltan_global_calculate_edge_weights
   ! elements with quality greater than this value are ok
 
=== modified file 'assemble/Zoltan_integration.F90'
--- assemble/Zoltan_integration.F90	2011-12-07 11:33:45 +0000
+++ assemble/Zoltan_integration.F90	2012-09-18 19:42:23 +0000
@@ -136,12 +136,20 @@
     zoltan_global_migrate_extruded_mesh = option_count('/geometry/mesh/from_mesh/extrude') > 0 &
       .and. .not. present_and_true(ignore_extrusion)
 
+    zoltan_global_field_weighted_partitions = &
+     have_option(trim(zoltan_global_base_option_path) // "/field_weighted_partitions")
+
+    if(zoltan_global_migrate_extruded_mesh .AND. zoltan_global_field_weighted_partitions) then
+        ewrite(-1,*) "Cannot weight mesh partitions based upon extruded columns"// &
+                     "and a prescribed field. Select one option only or fix the code."
+        FLExit("Use Weighted mesh partitions for EITHER extruded meshes or prescribed fields")
+    end if
+    
     call setup_module_variables(states, final_adapt_iteration, zz)
-    
+
     call setup_quality_module_variables(states, metric) ! this needs to be called after setup_module_variables
                                         ! (but only on the 2d mesh with 2+1d adaptivity)
 
-
     load_imbalance_tolerance = get_load_imbalance_tolerance(final_adapt_iteration)
     call set_zoltan_parameters(final_adapt_iteration, flredecomp, flredecomp_target_procs, load_imbalance_tolerance, zz)
 
@@ -424,7 +432,22 @@
           call insert(zoltan_global_universal_element_number_to_region_id, universal_element_number, zoltan_global_zz_positions%mesh%region_ids(i))
        end do
     end if
-    
+
+    if(zoltan_global_field_weighted_partitions) then
+       zoltan_global_field_weighted_partition_values = extract_scalar_field(states, "FieldWeightedPartitionValues") 
+       assert(zoltan_global_field_weighted_partition_values%mesh == zoltan_global_zz_mesh)
+
+       if(zoltan_global_field_weighted_partition_values%mesh%name /= zoltan_global_zz_mesh%name) then
+          ewrite(-1,*) "FieldWeightedPartitionValues and Zoltan Global ZZ Mesh must be on the " // &
+                       "same mesh. 99.9% of the time, this means that FieldWeightedPartitionValues " // & 
+                       "must be on the external mesh."
+          FLExit("FieldWeightedPartitionValues must be on the external mesh")
+       end if
+
+       call incref(zoltan_global_field_weighted_partition_values)
+
+    end if
+
   end subroutine setup_module_variables
 
   subroutine setup_quality_module_variables(states, metric)
@@ -724,6 +747,10 @@
     if(zoltan_global_migrate_extruded_mesh) then
        call deallocate(zoltan_global_columns_sparsity)
     end if
+    if(zoltan_global_field_weighted_partitions) then
+       call deallocate(zoltan_global_field_weighted_partition_values)
+    end if
+
   end subroutine cleanup_quality_module_variables
 
   subroutine cleanup_other_module_variables
@@ -1787,7 +1814,6 @@
     
   end subroutine initialise_transfer
 
-
   subroutine update_detector_list_element(detector_list_array)
     ! Update the detector%element field for every detector left in our list
     ! and check that we did not miss any in the first send
 
=== modified file 'main/Fluids.F90'
--- main/Fluids.F90	2012-09-05 08:44:55 +0000
+++ main/Fluids.F90	2012-09-18 19:42:23 +0000
@@ -893,6 +893,7 @@
        end if
 
     end do timestep_loop
+
     ! ****************************
     ! *** END OF TIMESTEP LOOP ***
     ! ****************************
 
=== modified file 'manual/meshes.tex'
--- manual/meshes.tex	2011-11-07 12:52:36 +0000
+++ manual/meshes.tex	2012-09-18 19:42:23 +0000
@@ -238,7 +238,13 @@
 The output of running flredecomp is a series of mesh and vtu files as well
 as the new flml; in this case \lstinline+foo_flredecomp.flml+. 
 Note that \lstinline[language=bash]+flredecomp+ must be run on a number of processors equal to the larger number of processors between input and output.
-More information can be found in section~\ref{sec:flredecomp}.
+
+When using flredecomp, it is possible to partition the mesh based upon a user defined weighting. 
+This is achieved by prescribing a scalar field, bounded between values of 0 and 1, under 
+\option{/flredecomp/field\_weighted\_partitions}. Flredecomp will then try to ensure that the sum 
+of weights on each partition is approximately equal.
+
+Further information on flredecomp can be found in section~\ref{sec:flredecomp}.
 
 \subsection{Decomposing a periodic mesh}
 \index{mesh!meshing tools!periodise}
 
=== modified file 'preprocessor/Populate_State.F90'
--- preprocessor/Populate_State.F90	2012-09-03 11:38:44 +0000
+++ preprocessor/Populate_State.F90	2012-09-18 19:42:23 +0000
@@ -1211,6 +1211,15 @@
             states(1), dont_allocate_prognostic_value_spaces=dont_allocate_prognostic_value_spaces)
     end if
 
+    ! Field that controls the weighting of partitions:
+    if (have_option('/flredecomp/field_weighted_partitions')) then
+       call allocate_and_insert_scalar_field('/flredecomp/field_weighted_partitions/scalar_field::FieldWeightedPartitionValues', states(1))
+    end if
+    
+    if (have_option('/mesh_adaptivity/hr_adaptivity/zoltan_options/field_weighted_partitions')) then
+       call allocate_and_insert_scalar_field('/mesh_adaptivity/hr_adaptivity/zoltan_options/field_weighted_partitions/scalar_field::FieldWeightedPartitionValues', states(1))
+    end if
+    
     ! grid velocity
     if (have_option('/mesh_adaptivity/mesh_movement/vector_field::GridVelocity')) then
        call allocate_and_insert_vector_field('/mesh_adaptivity/mesh_movement/vector_field::GridVelocity', &
 
=== modified file 'schemas/adaptivity_options.rnc'
--- schemas/adaptivity_options.rnc	2012-03-21 10:13:54 +0000
+++ schemas/adaptivity_options.rnc	2012-09-18 19:42:23 +0000
@@ -974,6 +974,20 @@
             element additional_adapt_iterations {
                 integer
             }?,
+            element field_weighted_partitions {
+            ## Field weighted partitions: turning on this option allows one to 
+            ## weight mesh partitions, based upon a prescribed scalar field. 
+            ## Note that the field should have a minimum value of 0 and a maximum
+            ## value of 1 (no normalisation is done within the code).
+	        element scalar_field {
+                  attribute rank { "0" },
+                  attribute name { "FieldWeightedPartitionValues" },
+                  element prescribed {
+                     coordinate_mesh_choice,
+                     prescribed_scalar_field
+                  }
+               }
+            }?,
             ## Zoltan Debugging
             ## 
             ## Turn on more verbose output for use when debugging Zoltan.
 
=== modified file 'schemas/adaptivity_options.rng'
--- schemas/adaptivity_options.rng	2012-03-21 10:13:54 +0000
+++ schemas/adaptivity_options.rng	2012-09-18 19:42:23 +0000
@@ -1143,6 +1143,26 @@
             </element>
           </optional>
           <optional>
+            <element name="field_weighted_partitions">
+              <element name="scalar_field">
+                <a:documentation>Field weighted partitions: turning on this option allows one to 
+weight mesh partitions, based upon a prescribed scalar field. 
+Note that the field should have a minimum value of 0 and a maximum
+value of 1 (no normalisation is done within the code).</a:documentation>
+                <attribute name="rank">
+                  <value>0</value>
+                </attribute>
+                <attribute name="name">
+                  <value>FieldWeightedPartitionValues</value>
+                </attribute>
+                <element name="prescribed">
+                  <ref name="coordinate_mesh_choice"/>
+                  <ref name="prescribed_scalar_field"/>
+                </element>
+              </element>
+            </element>
+          </optional>
+          <optional>
             <element name="zoltan_debug">
               <a:documentation>Zoltan Debugging
 
 
=== modified file 'schemas/flredecomp.rnc'
--- schemas/flredecomp.rnc	2011-09-12 11:17:09 +0000
+++ schemas/flredecomp.rnc	2012-09-18 19:42:23 +0000
@@ -30,6 +30,20 @@
                 }
             }
         }?,
+        element field_weighted_partitions {
+        ## Field weighted partitions: turning on this option allows one to 
+        ## weight mesh partitions, based upon a prescribed scalar field.
+        ## Note that the field should have a minimum value of 0 and a maximum
+        ## value of 1 (no normalisation is done within the code).
+	    element scalar_field {
+               attribute rank { "0" },
+               attribute name { "FieldWeightedPartitionValues" },
+               element prescribed {
+                  coordinate_mesh_choice,
+                  prescribed_scalar_field
+               }
+            }
+        }?,
         ## Zoltan Debugging
         ## 
         ## Turn on more verbose output for use when debugging Zoltan.
 
=== modified file 'schemas/flredecomp.rng'
--- schemas/flredecomp.rng	2011-09-12 11:17:09 +0000
+++ schemas/flredecomp.rng	2012-09-18 19:42:23 +0000
@@ -38,6 +38,26 @@
         </element>
       </optional>
       <optional>
+        <element name="field_weighted_partitions">
+          <element name="scalar_field">
+            <a:documentation>Field weighted partitions: turning on this option allows one to 
+weight mesh partitions, based upon a prescribed scalar field.
+Note that the field should have a minimum value of 0 and a maximum
+value of 1 (no normalisation is done within the code).</a:documentation>
+            <attribute name="rank">
+              <value>0</value>
+            </attribute>
+            <attribute name="name">
+              <value>FieldWeightedPartitionValues</value>
+            </attribute>
+            <element name="prescribed">
+              <ref name="coordinate_mesh_choice"/>
+              <ref name="prescribed_scalar_field"/>
+            </element>
+          </element>
+        </element>
+      </optional>
+      <optional>
         <element name="zoltan_debug">
           <a:documentation>Zoltan Debugging
 
 
=== added directory 'tests/flredecomp_2d_fieldweighted'
=== added file 'tests/flredecomp_2d_fieldweighted/Makefile'
--- tests/flredecomp_2d_fieldweighted/Makefile	1970-01-01 00:00:00 +0000
+++ tests/flredecomp_2d_fieldweighted/Makefile	2012-09-18 19:42:23 +0000
@@ -0,0 +1,13 @@
+MODEL=flredecomp-2d-fieldweighted
+
+default: input
+
+input:  clean
+	gmsh -2 -optimize src/Subduction_Mesh.geo
+	cp src/Subduction_Mesh.msh .
+	../../bin/gmsh2triangle --2d -i Subduction_Mesh.msh
+
+clean:
+	rm -f  *.ele *.edge *.node *.vtu *.stat *.msh *.detectors fluidity.* Parallel*
+	rm -R -f $(MODEL)_*
+
 
=== added file 'tests/flredecomp_2d_fieldweighted/check_partitions.py'
--- tests/flredecomp_2d_fieldweighted/check_partitions.py	1970-01-01 00:00:00 +0000
+++ tests/flredecomp_2d_fieldweighted/check_partitions.py	2012-09-18 19:42:23 +0000
@@ -0,0 +1,28 @@
+#! /usr/bin/env python
+
+import sys
+import numpy
+
+num_partitions = 16
+
+def calculate_partition_weights():
+    partition_weights = numpy.zeros(num_partitions,dtype=float)
+    
+    for i in range(num_partitions):
+        filename = 'Parallel-NP16-flredecomp-2d-fieldweighted_CoordinateMesh_%d.ele'%i
+        elements = numpy.loadtxt(filename, skiprows=1)
+        # Sum weights across each partition.
+        # Note prognostic region ID == 27: weighting (in .flml) == 1.0
+        # Prescribed regions ID == 24, 30: weighting (in .flml) == 0.01
+        for j in range(elements.shape[0]):
+            if(elements[j,4] == 27):
+                partition_weights[i] = partition_weights[i] + 1.0
+            else:
+                partition_weights[i] = partition_weights[i] + 0.005
+
+    min_weight = min(partition_weights)
+    max_weight = max(partition_weights)
+    max_load_imbalance = max_weight / min_weight
+
+    return max_load_imbalance
+
 
=== added file 'tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml'
--- tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml	1970-01-01 00:00:00 +0000
+++ tests/flredecomp_2d_fieldweighted/flredecomp-2d-fieldweighted.flml	2012-09-18 19:42:23 +0000
@@ -0,0 +1,410 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">flredecomp_field_weighted_partition_test</string_value>
+    <comment>Checks for field_weighted_partitions in flredecomp and gmsh2triangle -i. </comment>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">stokes</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">2</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="Subduction_Mesh">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">5</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">2e12</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="CoordinateMesh"/>
+    <stat>
+      <output_at_start/>
+    </stat>
+    <detectors>
+      <static_detector name="Temperature6060">
+        <location>
+          <real_value shape="2" dim1="dim" rank="1">60000. -60000.</real_value>
+        </location>
+      </static_detector>
+      <detector_array name="TSlab">
+        <number_of_detectors>
+          <integer_value rank="0">36</integer_value>
+        </number_of_detectors>
+        <static/>
+        <python>
+          <string_value lines="20" type="code" language="python">def val(t):
+  import numpy
+  start = numpy.array([0, 0])
+  end = numpy.array([210000, -210000])
+  output = []
+  for i in range(36):
+    output.append(start + (float(i)/35) * end)
+  return output</string_value>
+        </python>
+      </detector_array>
+      <detector_array name="TWedge">
+        <number_of_detectors>
+          <integer_value rank="0">78</integer_value>
+        </number_of_detectors>
+        <static/>
+        <python>
+          <string_value lines="20" type="code" language="python">def val(t):
+  output1 = []
+  for i in range(12):
+    for j in range(12):
+        if i &gt;= j:
+           output1.append([54000 + (float(i)/20) * 120000, -54000 + (float(j)/20) * (-120000)])
+  return output1</string_value>
+        </python>
+      </detector_array>
+      <fail_outside_domain/>
+    </detectors>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1e11</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0e15</real_value>
+    </finish_time>
+    <adaptive_timestep>
+      <requested_cfl>
+        <real_value rank="0">3.0</real_value>
+      </requested_cfl>
+      <courant_number name="CFLNumber">
+        <mesh name="CoordinateMesh"/>
+      </courant_number>
+      <increase_tolerance>
+        <real_value rank="0">5</real_value>
+      </increase_tolerance>
+    </adaptive_timestep>
+    <steady_state>
+      <tolerance>
+        <real_value rank="0">1e-3</real_value>
+        <infinity_norm/>
+      </tolerance>
+    </steady_state>
+  </timestepping>
+  <material_phase name="Fluid">
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="CoordinateMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <remove_stabilisation_term/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method>
+            <full_schur_complement>
+              <inner_matrix name="FullMomentumMatrix">
+                <solver>
+                  <iterative_method name="preonly"/>
+                  <preconditioner name="lu"/>
+                  <relative_error>
+                    <real_value rank="0">1.0e-6</real_value>
+                  </relative_error>
+                  <max_iterations>
+                    <integer_value rank="0">10000</integer_value>
+                  </max_iterations>
+                  <start_from_zero/>
+                  <never_ignore_solver_failures/>
+                  <diagnostics>
+                    <monitors/>
+                  </diagnostics>
+                </solver>
+              </inner_matrix>
+              <preconditioner_matrix name="DiagonalSchurComplement"/>
+            </full_schur_complement>
+          </use_projection_method>
+        </scheme>
+        <solver>
+          <iterative_method name="fgmres"/>
+          <preconditioner name="jacobi"/>
+          <relative_error>
+            <real_value rank="0">1.0e-4</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <exclude_from_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <exclude_mass_terms/>
+            </mass_terms>
+            <advection_terms>
+              <exclude_advection_terms/>
+            </advection_terms>
+            <stress_terms>
+              <tensor_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">1.0</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">1.0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="preonly"/>
+          <preconditioner name="lu"/>
+          <relative_error>
+            <real_value rank="0">1.0e-6</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="Batchelor">
+          <region_ids>
+            <integer_value shape="1" rank="1">27</integer_value>
+          </region_ids>
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X, t):
+  
+   	from math import sqrt, atan, cos, sin, pi
+
+   	Xcoord = X[0]
+        Ycoord = X[1]
+  
+   	Depth = (X[1] * -1.) -50000.
+   	Xdist  = X[0] -50000.
+
+   	if(Xdist == 0.) :
+	   	Xdist = 0.000000000000001
+
+   	if Depth == 0 :
+ 
+      		Vx = 0.
+      		Vy = 0.
+
+   	else :
+ 
+      		ALFA = atan(1)
+      		THETA = atan(Depth/Xdist)
+
+		VTHETA = (((ALFA - THETA) * sin(THETA) * sin(ALFA)) - (ALFA * THETA * sin(ALFA-THETA))) / (ALFA**2 - (sin(ALFA)**2))
+		VTHETA = VTHETA * (0.05/(3600.*24.*365.)) *-1.
+
+		VR = (((ALFA - THETA) * cos(THETA) * sin(ALFA)) - (sin(ALFA) * sin(THETA)) - (ALFA * sin(ALFA-THETA)) + (ALFA*THETA*cos(ALFA-THETA))) / ((ALFA**2) - (sin(ALFA))**2)
+      		VR = VR * (0.05/(3600.*24.*365.))
+
+      		Vx = -( VTHETA*sin(THETA) - VR*cos(THETA) )
+      		Vy = - (VTHETA*cos(THETA) + VR*sin(THETA) )
+
+   	return [Vx, Vy]</string_value>
+          </python>
+        </initial_condition>
+        <prescribed_region name="Crust">
+          <region_ids>
+            <integer_value shape="1" rank="1">24</integer_value>
+          </region_ids>
+          <constant>
+            <real_value shape="2" dim1="dim" rank="1">0. 0.</real_value>
+          </constant>
+        </prescribed_region>
+        <prescribed_region name="Slab">
+          <region_ids>
+            <integer_value shape="1" rank="1">30</integer_value>
+          </region_ids>
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X, t):
+        from math import sqrt
+        vx = sqrt(((0.05/(3600.*24.*365.))**2.)/2.)
+        vy = -sqrt(((0.05/(3600.*24.*365.))**2.)/2.)
+        return [vx,vy]</string_value>
+          </python>
+        </prescribed_region>
+        <boundary_conditions name="Crust">
+          <surface_ids>
+            <integer_value shape="1" rank="1">32</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </x_component>
+              <y_component>
+                <constant>
+                  <real_value rank="0">0.0</real_value>
+                </constant>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="Slab">
+          <surface_ids>
+            <integer_value shape="1" rank="1">31</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <align_bc_with_cartesian>
+              <x_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X, t):
+        from math import sqrt
+        vx = sqrt(((0.05/(3600.*24.*365.))**2.)/2.)
+        return vx</string_value>
+                </python>
+              </x_component>
+              <y_component>
+                <python>
+                  <string_value lines="20" type="code" language="python">def val(X, t):
+        from math import sqrt
+        vy = -sqrt(((0.05/(3600.*24.*365.))**2.)/2.)
+        return vy</string_value>
+                </python>
+              </y_component>
+            </align_bc_with_cartesian>
+          </type>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <isotropic>
+                <constant>
+                  <real_value rank="0">1.0e21</real_value>
+                </constant>
+              </isotropic>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <surface_integral type="normal" name="Outflow">
+            <surface_ids>
+              <integer_value shape="2" rank="1">13 16</integer_value>
+            </surface_ids>
+          </surface_integral>
+          <surface_integral type="normal" name="All">
+            <surface_ids>
+              <integer_value shape="3" rank="1">12 13 16</integer_value>
+            </surface_ids>
+          </surface_integral>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+  </material_phase>
+  <flredecomp>
+    <field_weighted_partitions>
+      <scalar_field name="FieldWeightedPartitionValues" rank="0">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="Prognostic_Regions">
+            <region_ids>
+              <integer_value shape="1" rank="1">27</integer_value>
+            </region_ids>
+            <constant>
+              <real_value rank="0">1.</real_value>
+            </constant>
+          </value>
+          <value name="Prescribed_Regions">
+            <region_ids>
+              <integer_value shape="2" rank="1">24 30</integer_value>
+            </region_ids>
+            <constant>
+              <real_value rank="0">0.005</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat/>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </scalar_field>
+    </field_weighted_partitions>
+  </flredecomp>
+</fluidity_options>
 
=== added file 'tests/flredecomp_2d_fieldweighted/flredecomp_2d_fieldweighted.xml'
--- tests/flredecomp_2d_fieldweighted/flredecomp_2d_fieldweighted.xml	1970-01-01 00:00:00 +0000
+++ tests/flredecomp_2d_fieldweighted/flredecomp_2d_fieldweighted.xml	2012-09-18 19:42:23 +0000
@@ -0,0 +1,25 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>flredecomp-2d-fieldweighted</name>
+  <owner userid="rhodrid"/>
+  <tags>flml</tags>
+
+  <problem_definition length="medium" nprocs="1">
+    <command_line>mpiexec -n 16 flredecomp -i 1 -o 16 flredecomp-2d-fieldweighted Parallel-NP16-flredecomp-2d-fieldweighted</command_line>
+  </problem_definition>
+
+  <variables>
+    <variable name="max_load_imbalance" language="python">
+from check_partitions import *
+max_load_imbalance = calculate_partition_weights()
+      <comment>Calculates the weight of each partition (i.e. sums weights of elements in each partition). The weights are prescribed in the flml to be 1.0 for region 27, and 0.005 elsewhere.</comment>
+    </variable>
+  </variables>
+
+  <pass_tests>
+    <test name="maximum load imbalance .lt. 1.5 " language="python">
+      assert ( max_load_imbalance &lt; 1.5 )
+    </test>
+  </pass_tests>
+
+</testproblem>
 
=== added directory 'tests/flredecomp_2d_fieldweighted/src'
=== added file 'tests/flredecomp_2d_fieldweighted/src/Subduction_Mesh.geo'
--- tests/flredecomp_2d_fieldweighted/src/Subduction_Mesh.geo	1970-01-01 00:00:00 +0000
+++ tests/flredecomp_2d_fieldweighted/src/Subduction_Mesh.geo	2012-09-18 19:42:23 +0000
@@ -0,0 +1,52 @@
+// Gmsh project created on Fri Dec  4 14:59:47 2009
+// We reproduce a subduction zone on a (660x600) km domain. 
+// The kinematically moving slab (see image, on the left) subducts under a 50km-thick plate (up, right) and hence dynamically inducing a flow in the upper-mantle wedge (down, on the right). 
+// We therefore are trying to match the benchmark results established by vanKeken2008 for this problem, computing the thermal field throughout the domain by use of finite element methods. 
+// The structured, non-adaptive mesh shown in the images contains for instance the sought information at the mesh nodes, spaced out 2.5 km from each other near the critical point of contact between the three// different regions of our domain, and 7.5km everywhere else. In this particular case, the test uses a prescribed wedge flow (Batchelor1967) and is mainly intended to check for the validity of initial condi// tions at the various boundaries. 
+
+minsize=4000;
+maxsize=4000;
+
+Point(1) = {0,0,0,minsize};
+Point(2) = {660000,0,0,maxsize};
+Point(3) = {660000,-50000,0,maxsize};
+Point(4) = {660000,-600000,0,maxsize};
+Point(5) = {600000,-600000,0,maxsize};
+Point(6) = {0,-600000,0,maxsize};
+Point(7) = {50000,-50000,0,minsize};
+Point(8) = {660000,-420000,0,maxsize};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,8};
+Line(4) = {4,5};
+Line(5) = {5,6};
+Line(7) = {1,7};
+Line(9) = {7,3};
+Line(10) = {8,4};
+Line(11) = {7,5};
+Line(14) = {6,1};
+
+Physical Line(10) = {1};
+Physical Line(11) = {2};
+Physical Line(12) = {3};
+Physical Line(13) = {4};
+Physical Line(19) = {5};
+Physical Line(15) = {14};
+Physical Line(16) = {10};
+// Physical Line(17) = {9};
+// Physical Line(18) = {11};
+// Physical Line(20) = {7};
+
+Line Loop(22) = {1,2,-9,-7};
+Plane Surface(23) = {22};
+Physical Surface(24) = {23};
+Line Loop(25) = {9,3,10,4,-11};
+Plane Surface(26) = {25};
+Physical Surface(27) = {26};
+Line Loop(28) = {5,14,7,11};
+Plane Surface(29) = {28};
+Physical Surface(30) = {29};
+
+Physical Line(31) = {11};
+Physical Line(32) = {9};