Merge lp:~fluidity-core/fluidity/moving-mesh-checkpoint into lp:fluidity

[Stephan Kramer]

Fix checkpointing for free surface moving mesh simulations by using OriginalCoordinate for the dumping the mesh and the field vtus. The mesh movement will be reapplied at the start of the restart run. Checkpointing with a moving mesh not controlled by free surface movement, now gives a warning (not supported). Switching on mesh movement in spherical patch test case to test this checkpointing.

[Stephan Kramer (s-kramer)]

Thanks for reviewing, Sam. I'll merge in the trunk (as this is a while ago) and see if it goes green again before merging.

[lp:~fluidity-core/fluidity/moving-mesh-checkpoint]

3992. By Stephan Kramer on 2013-06-10

Merging in trunk.

3993. By Stephan Kramer on 2013-06-12

Relax those tolerances a bit.

=== modified file 'femtools/Checkpoint.F90'
--- femtools/Checkpoint.F90	2012-07-27 11:01:30 +0000
+++ femtools/Checkpoint.F90	2013-06-12 14:29:25 +0000
@@ -486,22 +486,35 @@
           end if
         end if
 
+        ! Write out the mesh using a suitable coordinate field
+        if (mesh%name=="CoordinateMesh") then
+          if (have_option("/mesh_adaptivity/mesh_movement/free_surface")) then
+            ! we don't want/need to checkpoint the moved mesh, as the mesh movement will again be applied after the restart
+            ! based on the checkpointed FreeSurface field.
+            position => extract_vector_field(state(1), "OriginalCoordinate", stat=stat1)
+            if (stat1/=0) then
+              ! some cases (e.g. flredecomp) OriginalCoordinate doesn't exist
+              position => extract_vector_field(state(1), "Coordinate")
+            end if
+          else
+            position => extract_vector_field(state(1), "Coordinate")
+          end if
+        else
+          position => extract_vector_field(state(1), trim(mesh%name)//"Coordinate")
+        end if
+
         if(get_active_nparts(ele_count(mesh)) > 1) then
-          ! Write out the mesh
-          if (mesh%name=="CoordinateMesh") then
-            position => extract_vector_field(state(1), "Coordinate")
-          else
-            position => extract_vector_field(state(1), trim(mesh%name)//"Coordinate")
-          end if
           call write_mesh_files(parallel_filename(mesh_filename), mesh_format, position)
           ! Write out the halos
           ewrite(2, *) "Checkpointing halos"
           call write_halos(mesh_filename, mesh)
         else
           ! Write out the mesh
-          call write_mesh_files(mesh_filename, mesh_format, state(1), mesh)
+          call write_mesh_files(mesh_filename, mesh_format, position)
         end if
-     end if
+
+      end if
+
    end do
 
   end subroutine checkpoint_meshes
@@ -534,7 +547,24 @@
     assert(len_trim(prefix) > 0)
 
     do i = 1, size(state)
-      positions => extract_vector_field(state(i), "Coordinate")
+      if (have_option("/mesh_adaptivity/mesh_movement/free_surface")) then
+        ! we don't want/need to checkpoint the moved mesh, as the mesh movement will again be applied after the restart
+        ! based on the checkpointed FreeSurface field.
+        positions => extract_vector_field(state(i), "OriginalCoordinate", stat=stat)
+        if (stat/=0) then
+          ! some cases (e.g. flredecomp) OriginalCoordinate doesn't exist
+          positions => extract_vector_field(state(1), "Coordinate")
+        end if
+      else
+        if (have_option("/mesh_adaptivity/mesh_movement")) then
+          ! for other mesh movement schemes we should probably write out the "moved" mesh to retain that information
+          ! However if the checkpointed mesh is not the "CoordinateMesh", it will have its own MeshNameCoordinate field
+          ! that hasn't moved; leading to a failure to restart from the checkpoint. This is only one of the possible problems
+          ! generally this functionality is untested.
+          ewrite(0,*) "WARNING: using mesh_movement with checkpointing is untested and likely broken."
+        end if
+        positions => extract_vector_field(state(i), "Coordinate")
+      end if
       do j = 1, size(state(i)%meshes)
         mesh => state(i)%meshes(j)%ptr
         nparts = get_active_nparts(ele_count(mesh))

=== modified file 'tests/spherical_patch/spherical_patch.xml'
--- tests/spherical_patch/spherical_patch.xml	2012-07-24 16:34:10 +0000
+++ tests/spherical_patch/spherical_patch.xml	2013-06-12 14:29:25 +0000
@@ -11,16 +11,10 @@
   <variables>
     <variable name="volume_before" language="python">from fluidity_tools import stat_parser
 stat=stat_parser('spherical_segment.stat')
-# this only works because we don't move the free surface
-# we actually compute H \int_top \eta, i.e. the volume difference with
-# equilibrium state multiplied by depth H
-volume_before=stat['Fields']['FreeSurface']['integral']</variable>
+volume_before=stat['Fields']['Density']['integral']</variable>
     <variable name="volume_after" language="python">from fluidity_tools import stat_parser
 stat=stat_parser('spherical_segment_checkpoint.stat')
-# this only works because we don't move the free surface
-# we actually compute H \int_top \eta, i.e. the volume difference with
-# equilibrium state multiplied by depth H
-volume_after=stat['Fields']['FreeSurface']['integral']</variable>
+volume_after=stat['Fields']['Density']['integral']</variable>
     <variable name="solvers_converged" language="python">import os
 files = os.listdir("./")
 solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files</variable>
@@ -30,19 +24,33 @@
     <variable name="max_velocity_after" language="python">from fluidity_tools import stat_parser
 stat=stat_parser('spherical_segment_checkpoint.stat')
 max_velocity_after=stat['Fields']['Velocity%magnitude']['max']</variable>
+    <variable name="fs_integral_before" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment.stat')
+fs_integral_before=stat['Fields']['FreeSurface']['surface_integral%Bottom']</variable>
+    <variable name="fs_integral_after" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('spherical_segment_checkpoint.stat')
+fs_integral_after=stat['Fields']['FreeSurface']['surface_integral%Bottom']</variable>
   </variables>
   <pass_tests>
     <test name="volume_conservation_run1" language="python">from fluidity_tools import compare_variable
-compare_variable(volume_before[0], volume_before[-1], 1e-5)</test>
+compare_variable(volume_before[0], volume_before[-1], 1e-9)</test>
     <test name="volume_conservation_run12" language="python">from fluidity_tools import compare_variable
 # compares volume at the end of first run with that of end first timestep second run
-compare_variable(volume_before[-1], volume_after[0], 1e-5)</test>
+compare_variable(volume_before[-1], volume_after[0], 1e-9)</test>
     <test name="volume_conservation_run2" language="python">from fluidity_tools import compare_variable
-compare_variable(volume_after[0], volume_after[-1], 1e-5)</test>
+compare_variable(volume_after[0], volume_after[-1], 1e-9)</test>
     <test name="solvers_converged" language="python">assert solvers_converged</test>
     <test name="max_velocity_run1" language="python"># make sure we're actually doing something and not blowing up
 assert all(max_velocity_before&gt;0.01) and all(max_velocity_before&lt;0.02)</test>
     <test name="max_velocity_run2" language="python"># make sure we're actually doing something and not blowing up
 assert all(max_velocity_after&gt;0.01) and all(max_velocity_after&lt;0.02)</test>
+    <test name="fs_integral_run1" language="python">from fluidity_tools import compare_variable
+compare_variable(fs_integral_before[0], fs_integral_before[-1], 1e-5)</test>
+    <test name="fs_integral_run2" language="python">from fluidity_tools import compare_variable
+compare_variable(fs_integral_after[0], fs_integral_after[-1], 1e-5)</test>
+    <test name="fs_integral_run12" language="python">from fluidity_tools import compare_variable
+# compares freesurface integral at the end of first run with that of end first timestep second run
+compare_variable(fs_integral_before[-1], fs_integral_after[0], 1e-5)
+</test>
   </pass_tests>
 </testproblem>

=== modified file 'tests/spherical_patch/spherical_segment.flml'
--- tests/spherical_patch/spherical_segment.flml	2013-05-30 11:47:42 +0000
+++ tests/spherical_patch/spherical_segment.flml	2013-06-12 14:29:25 +0000
@@ -101,14 +101,14 @@
     </mesh>
     <quadrature>
       <degree>
-        <integer_value rank="0">4</integer_value>
+        <integer_value rank="0">6</integer_value>
       </degree>
       <surface_degree>
-        <integer_value rank="0">3</integer_value>
+        <integer_value rank="0">6</integer_value>
       </surface_degree>
     </quadrature>
     <spherical_earth>
-      <superparametric_mapping/>
+      <linear_mapping/>
     </spherical_earth>
     <ocean_boundaries>
       <top_surface_ids>
@@ -366,7 +366,7 @@
             <real_value rank="0">1.0e-8</real_value>
           </relative_error>
           <max_iterations>
-            <integer_value rank="0">10000</integer_value>
+            <integer_value rank="0">1000</integer_value>
           </max_iterations>
           <never_ignore_solver_failures/>
           <diagnostics>
@@ -453,7 +453,13 @@
         <algorithm name="Internal" material_phase_support="multiple"/>
         <mesh name="PressureMesh"/>
         <output/>
-        <stat/>
+        <stat>
+          <surface_integral type="value" name="Bottom">
+            <surface_ids>
+              <integer_value shape="1" rank="1">2</integer_value>
+            </surface_ids>
+          </surface_integral>
+        </stat>
         <convergence>
           <include_in_convergence/>
         </convergence>
@@ -500,4 +506,30 @@
       </diagnostic>
     </scalar_field>
   </material_phase>
+  <mesh_adaptivity>
+    <mesh_movement>
+      <free_surface>
+        <move_surface_nodes/>
+      </free_surface>
+      <vector_field name="GridVelocity" rank="1">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </vector_field>
+    </mesh_movement>
+  </mesh_adaptivity>
 </fluidity_options>