Merge lp:~michael-lange/fluidity/fix_detector_init into lp:fluidity

[Michael Lange]

Fixing the initialisation of detectors in parallel (lp:967093) and adding an according test case. Many thanks to M. Filipiak for providing the patch.

Also fixing detector initialisation from file in parallel, which would previously segfault, due to the legacy restriction to only one proc.

[lp:~michael-lange/fluidity/fix_detector_init]

4030. By Michael Lange on 2012-10-25

Merge from trunk

4031. By Michael Lange on 2012-10-26

The new test case should use flredecomp.

[Mark Filipiak (mjf-q)]

Not sure if I'm allowed to review the changes since I'm one of the authors but I'd like to get this fix into the trunk so I can then add the fldecomp-to-flredecomp change.

The change to get the detectors read in from file in parallel looks ok (and matches what is done for a checkpointed detector file in parallel). The new test fails on the old trunk and passes with this branch.

I looked again at the code I added to work out the owner of a detector. I still think it's OK. It's has a slight inefficiency as James Maddison noticed but I don't know if changing the picker construction to consider only owned elements will break other things (e.g. interpolate_boundary_values(Geostrophic_Pressure.F90)). Changing find_serial() to have the same ownership test as find_parallel() might *seem* cleaner but will certainly break interpolate_boundary_values(Geostrophic_Pressure.F90).

I think the names of the routine find_serial() and find_parallel() are perhaps not descriptive, the first really means 'find any element, that this process has in its mesh, that contains the point, even if this process doesn't own that element': I don't think it is worth changing the names.

It can be merged into the trunk (Michael, I have a version of fix_detector_init that I merged trunk r4125 into and tested up to medium tests locally, not in buildbot).

review: Approve

[lp:~michael-lange/fluidity/fix_detector_init]

4032. By Michael Lange on 2012-11-08

Merge from trunk

=== modified file 'femtools/Diagnostic_variables.F90'
--- femtools/Diagnostic_variables.F90	2012-01-20 15:54:44 +0000
+++ femtools/Diagnostic_variables.F90	2012-11-08 12:08:22 +0000
@@ -1369,7 +1369,10 @@
     assert(xfield%dim+1==local_coord_count(shape))
 
     ! Determine element and local_coords from position
-    call picker_inquire(xfield,position,element,local_coord=lcoords,global=.false.)
+    ! In parallel, global=.false. can often work because there will be
+    ! a halo of non-owned elements in your process and so you can work out
+    ! ownership without communication.  But in general it won't work.
+    call picker_inquire(xfield,position,element,local_coord=lcoords,global=.true.)
 
     ! If we're in parallel and don't own the element, skip this detector
     if (isparallel()) then
@@ -1568,26 +1571,24 @@
           else
 
              ! Reading from a binary file where the user has placed the detector positions
-             if (getprocno() == 1) then
-                 default_stat%detector_file_unit=free_unit()
-                 call get_option("/io/detectors/detector_array/from_file/file_name",detector_file_filename)
+             default_stat%detector_file_unit=free_unit()
+             call get_option("/io/detectors/detector_array/from_file/file_name",detector_file_filename)
 
 #ifdef STREAM_IO
-                 open(unit = default_stat%detector_file_unit, file = trim(detector_file_filename), &
-                      & action = "read", access = "stream", form = "unformatted")
+             open(unit = default_stat%detector_file_unit, file = trim(detector_file_filename), &
+                  & action = "read", access = "stream", form = "unformatted")
 #else
-                 FLAbort("No stream I/O support")
+             FLAbort("No stream I/O support")
 #endif
-       
-                 do j=1,ndete
-                    write(detector_name, fmt) trim(funcnam)//"_", j
-                    default_stat%detector_list%detector_names(k)=trim(detector_name)
-                    read(default_stat%detector_file_unit) detector_location
-                    call create_single_detector(default_stat%detector_list, xfield, &
-                          detector_location, k, type_det, trim(detector_name))
-                    k=k+1          
-                 end do
-              end if                
+   
+             do j=1,ndete
+                write(detector_name, fmt) trim(funcnam)//"_", j
+                default_stat%detector_list%detector_names(k)=trim(detector_name)
+                read(default_stat%detector_file_unit) detector_location
+                call create_single_detector(default_stat%detector_list, xfield, &
+                      detector_location, k, type_det, trim(detector_name))
+                k=k+1          
+             end do
           end if
        end do      
     else 
 
=== modified file 'femtools/Node_Owner_Finder_Fortran.F90'
--- femtools/Node_Owner_Finder_Fortran.F90	2010-11-09 14:08:19 +0000
+++ femtools/Node_Owner_Finder_Fortran.F90	2012-11-08 12:08:22 +0000
@@ -268,8 +268,16 @@
         call cnode_owner_finder_query_output(id, nele_ids)
 
         closest_ele_id = -1
-        do j = 1, nele_ids
+        possible_elements_loop: do j = 1, nele_ids
           call cnode_owner_finder_get_output(id, possible_ele_id, j)
+          ! If this process does not own this possible_ele_id element then
+          ! don't consider it.  This filter is needed to make this subroutine work in
+          ! parallel without having to use universal numbers, which aren't defined
+          ! for all the meshes that use this subroutine.
+          if(.not.element_owned(positions_a,possible_ele_id)) then
+             assert(isparallel())
+             cycle possible_elements_loop
+          end if
           ! Zero tolerance - we're not using an "epsilon-ball" approach here
           if(ownership_predicate(positions_a, possible_ele_id, positions(:, i), 0.0, miss = miss)) then
             ele_ids(i) = possible_ele_id
@@ -281,7 +289,7 @@
             closest_ele_id = possible_ele_id
             closest_misses(i) = miss
           end if
-        end do
+        end do possible_elements_loop
 
         ! We didn't find an owner, so choose the element with the closest miss
         ele_ids(i) = closest_ele_id
@@ -310,6 +318,9 @@
           ! Another processes has a smaller miss for this coordinate
           ele_ids(i) = -1
         end if
+        ! if no process has closest_misses(i) < out_of_bounds_tolerance
+        ! then ele_ids(i) is already set to -1 in positions_loop above
+        ! on all processes.  This matches the find_serial(0) behaviour.
       end do
 #endif
             
 
=== added directory 'tests/detectors_parallel_init'
=== added file 'tests/detectors_parallel_init/Makefile'
--- tests/detectors_parallel_init/Makefile	1970-01-01 00:00:00 +0000
+++ tests/detectors_parallel_init/Makefile	2012-11-08 12:08:22 +0000
@@ -0,0 +1,15 @@
+default: input
+
+input: clean
+	gmsh -2 -o src/square.msh src/square.geo
+	../../bin/gmsh2triangle --2d src/square.msh
+	mpiexec -n 2 ../../bin/flredecomp -i 1 -o 2 detectors detectors_flredecomp
+
+clean: clean-mesh clean-run
+clean-mesh:
+	rm -f *.ele *.edge *.face *.node *.halo *.pvtu *.vtu *.detectors *.detectors.dat *.stat
+	rm -f gmsh.log gmsh_err.log
+clean-run:
+	rm -f matrixdump matrixdump.info
+	rm -f fluidity.err-* fluidity.log-*
+	rm -rf detectors_[0-4]* *flredecomp*
 
=== added file 'tests/detectors_parallel_init/detectors.flml'
--- tests/detectors_parallel_init/detectors.flml	1970-01-01 00:00:00 +0000
+++ tests/detectors_parallel_init/detectors.flml	2012-11-08 12:08:22 +0000
@@ -0,0 +1,121 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">detectors</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="square">
+        <format name="triangle"/>
+        <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">6</integer_value>
+      </degree>
+    </quadrature>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">0.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+    <detectors>
+      <static_detector name="MidPoint">
+        <location>
+          <real_value shape="2" dim1="dim" rank="1">0.5 0.5</real_value>
+        </location>
+      </static_detector>
+      <detector_array name="DArray">
+        <number_of_detectors>
+          <integer_value rank="0">3</integer_value>
+        </number_of_detectors>
+        <static/>
+        <from_file file_name="positions.dat">
+          <format>
+            <string_value>binary</string_value>
+          </format>
+        </from_file>
+      </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">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <physical_parameters/>
+  <material_phase name="Water">
+    <vector_field name="Velocity" rank="1">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+	return(X[1],0.0)</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </prescribed>
+    </vector_field>
+    <scalar_field name="Tracer" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="code" language="python">def val(X,t):
+	return X[0]</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+  </material_phase>
+  <flredecomp>
+    <final_partitioner>
+      <zoltan>
+        <method>graph</method>
+      </zoltan>
+    </final_partitioner>
+  </flredecomp>
+</fluidity_options>
 
=== added file 'tests/detectors_parallel_init/detectors_parallel_init.xml'
--- tests/detectors_parallel_init/detectors_parallel_init.xml	1970-01-01 00:00:00 +0000
+++ tests/detectors_parallel_init/detectors_parallel_init.xml	2012-11-08 12:08:22 +0000
@@ -0,0 +1,61 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+
+<testproblem>
+  <name>detectors_parallel_init</name>
+  <owner userid="mlange"/>
+  <tags>flml</tags>
+  <problem_definition length="short" nprocs="1">
+    <command_line>mpiexec -n 2 ../../bin/fluidity detectors_flredecomp.flml </command_line>
+    <!-- Test that detectors are initialised correctly in parallel. -->
+  </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 = "time" language = "python">
+from fluidity_tools import stat_parser
+time = stat_parser("detectors.stat")["ElapsedTime"]["value"]
+    </variable>
+    <variable name = "det_midpoint" language = "python">
+from fluidity_tools import stat_parser
+det_midpoint = stat_parser("detectors.detectors")["MidPoint"]['position']
+    </variable>
+    <variable name = "det_arr1" language = "python">
+from fluidity_tools import stat_parser
+det_arr1 = stat_parser("detectors.detectors")["DArray_1"]['position']
+    </variable>
+    <variable name = "det_arr2" language = "python">
+from fluidity_tools import stat_parser
+det_arr2 = stat_parser("detectors.detectors")["DArray_2"]['position']
+    </variable>
+    <variable name = "det_arr3" language = "python">
+from fluidity_tools import stat_parser
+det_arr3 = stat_parser("detectors.detectors")["DArray_3"]['position']
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Solvers converged" language="python">
+      assert(solvers_converged)
+    </test>
+    <test name="Test finished" language="python">
+assert abs(time - 1.0) &lt; 1.0e-6
+    </test>
+    <test name="Detector output" language="python">
+assert abs(len(time) - len(det_midpoint[-1])) &lt; 1.0e-6
+    </test>
+    <test name="Detectors initialised correctly" language="python">
+assert abs(det_midpoint[0] - 0.5).all() &lt; 1.0e-6
+assert abs(det_midpoint[1] - 0.5).all() &lt; 1.0e-6
+assert abs(det_arr1[0]).all() &lt; 1.0e-6
+assert abs(det_arr1[1] - 1.0).all() &lt; 1.0e-6
+assert abs(det_arr2[0] - 0.5).all() &lt; 1.0e-6
+assert abs(det_arr2[1] - 0.5).all() &lt; 1.0e-6
+assert abs(det_arr3[0] - 1.0).all() &lt; 1.0e-6
+assert abs(det_arr3[1]).all() &lt; 1.0e-6
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
 
=== added file 'tests/detectors_parallel_init/positions.dat'
 Binary files tests/detectors_parallel_init/positions.dat	1970-01-01 00:00:00 +0000 and tests/detectors_parallel_init/positions.dat	2012-11-08 12:08:22 +0000 differ
=== added directory 'tests/detectors_parallel_init/src'
=== added file 'tests/detectors_parallel_init/src/square.geo'
--- tests/detectors_parallel_init/src/square.geo	1970-01-01 00:00:00 +0000
+++ tests/detectors_parallel_init/src/square.geo	2012-11-08 12:08:22 +0000
@@ -0,0 +1,7 @@
+Point (1) = {0, 0, 0};
+Extrude {1,0,0} {
+  Point{1}; Layers{2};
+}
+Extrude {0,1,0} {
+  Line{1}; Layers{2};
+}