Merge lp:~s-parkinson11/fluidity/sinking_velocity_schema_fix into lp:fluidity

[Samuel Parkinson]

Edited the schema to allow specification of the mesh used for sinking velocity fields. This is needed when the velocity field is on a different mesh to the scalar field with the sinking velocity or when using P1DGP2 with a DG field and projecting velocity to continuous.

[Samuel Parkinson (s-parkinson11)]

I realised I hadn't actually pushed the changes to this branch before requesting a review. I have now :)

[lp:~s-parkinson11/fluidity/sinking_velocity_schema_fix]

4068. By Samuel Parkinson on 2012-09-29

changes to schema so that sinking velocity can have any mesh

[Jon Hill (jon-hill)]

Add some documentation to the schema to say that this should be the velocity mesh?

Also, an options check somewhere wouldn't go amiss to make sure this is the case.

[lp:~s-parkinson11/fluidity/sinking_velocity_schema_fix]

4069. By Samuel Parkinson on 2012-10-04

Added SinkingVelocity mesh checking to check_options routines of:
Advection_Diffusion_CG.F90
Advection_Diffusion_FV.F90
Field_Equations_CV.F90

Added check_options routine in Advection_Diffusion_DG.F90 with check for SinkingVelocity mesh.

not tested yet.

4070. By Samuel Parkinson on 2012-10-04

merge with trunk

4071. By Samuel Parkinson on 2012-10-04

fixed bugs in options checking routines - these now work as planned

4072. By Samuel Parkinson on 2012-10-04

updated test flml files

[Samuel Parkinson (s-parkinson11)]

Thanks Jon,

I have made some changes following your review. Do you think these meet the requirements?

I did notice that all prognostic field advection-diffusion discretisation options checking only happens on scalar fields that are in the root node of a material phase in the schema file. This could do with fixing at some point.

Cheers,

Sam

[Jon Hill (jon-hill)]

Perfect.

review: Approve

[lp:~s-parkinson11/fluidity/sinking_velocity_schema_fix]

4073. By Samuel Parkinson on 2012-10-04

fixed bug that made all runs with scalar fields, but no sinking velocity, fail

4074. By Samuel Parkinson on 2012-10-25

merge with broken branch

4075. By Samuel Parkinson on 2012-11-27

merge with trunk

=== modified file 'assemble/Advection_Diffusion_CG.F90'
--- assemble/Advection_Diffusion_CG.F90	2012-08-12 22:16:26 +0000
+++ assemble/Advection_Diffusion_CG.F90	2012-11-27 18:25:25 +0000
@@ -1431,7 +1431,7 @@
   subroutine advection_diffusion_cg_check_options
     !!< Check CG advection-diffusion specific options
     
-    character(len = FIELD_NAME_LEN) :: field_name, state_name
+    character(len = FIELD_NAME_LEN) :: field_name, state_name, mesh_0, mesh_1
     character(len = OPTION_PATH_LEN) :: path
     integer :: i, j, stat
     real :: beta, l_theta
@@ -1504,7 +1504,20 @@
               call field_error(state_name, field_name, &
                 & "Weak Dirichlet boundary conditions with diffusivity are not supported by CG advection-diffusion")
             end if
-            
+                 
+            if (have_option(trim(path) // "/scalar_field::SinkingVelocity")) then
+               call get_option(trim(complete_field_path(trim(path) // &
+                    "/scalar_field::SinkingVelocity"))//"/mesh[0]/name", &
+                    mesh_0, stat)
+               if(stat == SPUD_NO_ERROR) then
+                  call get_option(trim(complete_field_path("/material_phase[" // int2str(i) // &
+                       "]/vector_field::Velocity")) // "/mesh[0]/name", mesh_1)
+                  if(trim(mesh_0) /= trim(mesh_1)) then
+                     call field_warning(state_name, field_name, &
+                          & "SinkingVelocity is on a different mesh to the Velocity field. This could cause problems")
+                  end if
+               end if
+            end if
             if(have_option(trim(path) // "/spatial_discretisation/continuous_galerkin/advection_terms/exclude_advection_terms")) then
               if(have_option(trim(path) // "/scalar_field::SinkingVelocity")) then
                 call field_warning(state_name, field_name, &
 
=== modified file 'assemble/Advection_Diffusion_DG.F90'
--- assemble/Advection_Diffusion_DG.F90	2012-10-11 17:52:01 +0000
+++ assemble/Advection_Diffusion_DG.F90	2012-11-27 18:25:25 +0000
@@ -45,6 +45,7 @@
   use petsc_solve_state_module
   use boundary_conditions
   use boundary_conditions_from_options
+  use field_options
   use spud
   use upwind_stabilisation
   use slope_limiters_dg
@@ -67,7 +68,8 @@
   character(len=255), private :: message
 
   private
-  public solve_advection_diffusion_dg, construct_advection_diffusion_dg
+  public solve_advection_diffusion_dg, construct_advection_diffusion_dg, &
+       advection_diffusion_dg_check_options
 
   ! Local private control parameters. These are module-global parameters
   ! because it would be expensive and/or inconvenient to re-evaluate them
@@ -3120,5 +3122,84 @@
     end subroutine initialise_aijxy
      
   end function masslumped_rt0_aij
+    
+  subroutine advection_diffusion_dg_check_options
+    !!< Check DG advection-diffusion specific options
+    
+    character(len = FIELD_NAME_LEN) :: field_name, mesh_0, mesh_1, state_name
+    character(len = OPTION_PATH_LEN) :: path
+    integer :: i, j, stat
+    real :: beta, l_theta
+
+    if(option_count("/material_phase/scalar_field/prognostic/spatial_discretisation/discontinuous_galerkin") == 0) then
+       ! Nothing to check
+       return
+    end if
+
+    ewrite(2, *) "Checking DG advection-diffusion options"
+
+    do i = 0, option_count("/material_phase") - 1
+       path = "/material_phase[" // int2str(i) // "]"
+       call get_option(trim(path) // "/name", state_name)
+
+       do j = 0, option_count(trim(path) // "/scalar_field") - 1
+          path = "/material_phase[" // int2str(i) // "]/scalar_field[" // int2str(j) // "]"
+          call get_option(trim(path) // "/name", field_name)
+
+          if(field_name /= "Pressure") then
+
+             path = trim(path) // "/prognostic"
+
+             if(have_option(trim(path) // "/spatial_discretisation/discontinuous_galerkin").and.&
+                  have_option(trim(path) // "/equation[0]")) then   
+                 
+                if (have_option(trim(path) // "/scalar_field::SinkingVelocity")) then
+                   call get_option(trim(complete_field_path(trim(path) // &
+                        "/scalar_field::SinkingVelocity"))//"/mesh[0]/name", &
+                        mesh_0, stat)
+                   if(stat == SPUD_NO_ERROR) then
+                      call get_option(trim(complete_field_path("/material_phase[" // int2str(i) // &
+                           "]/vector_field::Velocity")) // "/mesh[0]/name", mesh_1)
+                      if(trim(mesh_0) /= trim(mesh_1)) then
+                         call field_warning(state_name, field_name, &
+                              "SinkingVelocity is on a different mesh to the Velocity field this could "//&
+                              "cause problems. If using advection_scheme/project_velocity_to_continuous "//&
+                              "and a discontinuous Velocity field, then SinkingVelocity must be on a "//&
+                              "continuous mesh and hence should not be on the same mesh as the Velocity")
+                      end if
+                   end if
+                end if
+
+             end if
+
+          end if
+       end do
+    end do
+
+    ewrite(2, *) "Finished checking CG advection-diffusion options"
+
+  contains
+
+    subroutine field_warning(state_name, field_name, msg)
+      character(len = *), intent(in) :: state_name
+      character(len = *), intent(in) :: field_name
+      character(len = *), intent(in) :: msg
+
+      ewrite(0, *) "Warning: For field " // trim(field_name) // " in state " // trim(state_name)
+      ewrite(0, *) trim(msg)
+
+    end subroutine field_warning
+
+    subroutine field_error(state_name, field_name, msg)
+      character(len = *), intent(in) :: state_name
+      character(len = *), intent(in) :: field_name
+      character(len = *), intent(in) :: msg
+
+      ewrite(-1, *) "For field " // trim(field_name) // " in state " // trim(state_name)
+      FLExit(trim(msg))
+
+    end subroutine field_error
+
+  end subroutine advection_diffusion_dg_check_options
 
 end module advection_diffusion_DG
 
=== modified file 'assemble/Advection_Diffusion_FV.F90'
--- assemble/Advection_Diffusion_FV.F90	2012-08-09 12:50:11 +0000
+++ assemble/Advection_Diffusion_FV.F90	2012-11-27 18:25:25 +0000
@@ -598,7 +598,7 @@
   
   subroutine advection_diffusion_fv_check_options
   
-    character(len = FIELD_NAME_LEN) :: field_name, state_name
+    character(len = FIELD_NAME_LEN) :: field_name, mesh_0, mesh_1, state_name
     character(len = OPTION_PATH_LEN) :: path
     integer :: i, j, stat
     real :: beta, l_theta
@@ -652,7 +652,20 @@
               call field_warning(state_name, field_name, &
                 & "Implicitness factor (theta) should = 1.0 when excluding mass")
             end if
-  
+                 
+            if (have_option(trim(path) // "/scalar_field::SinkingVelocity")) then
+               call get_option(trim(complete_field_path(trim(path) // &
+                    "/scalar_field::SinkingVelocity"))//"/mesh[0]/name", &
+                    mesh_0, stat)
+               if(stat == SPUD_NO_ERROR) then
+                  call get_option(trim(complete_field_path("/material_phase[" // int2str(i) // &
+                       "]/vector_field::Velocity")) // "/mesh[0]/name", mesh_1)
+                  if(trim(mesh_0) /= trim(mesh_1)) then
+                     call field_warning(state_name, field_name, &
+                          & "SinkingVelocity is on a different mesh to the Velocity field this could cause problems")
+                  end if
+               end if
+            end if
             if(have_option(trim(path) // "/spatial_discretisation/finite_volume/advection_terms/exclude_advection_terms")) then
               if(have_option(trim(path) // "/scalar_field::SinkingVelocity")) then
                 call field_warning(state_name, field_name, &
 
=== modified file 'assemble/Field_Equations_CV.F90'
--- assemble/Field_Equations_CV.F90	2012-08-18 21:34:37 +0000
+++ assemble/Field_Equations_CV.F90	2012-11-27 18:25:25 +0000
@@ -3544,8 +3544,8 @@
     !************************************************************************
     ! control volume options checking
     subroutine field_equations_cv_check_options
-      integer :: nmat, nfield, m, f
-      character(len=OPTION_PATH_LEN) :: mat_name, field_name, diff_scheme
+      integer :: nmat, nfield, m, f, stat
+      character(len=OPTION_PATH_LEN) :: mat_name, field_name, mesh_0, mesh_1, diff_scheme
       integer :: weakdirichlet_count
       logical :: cv_disc, mmat_cv_disc, diff, conv_file, subcycle, cv_temp_disc, tolerance, explicit
       real :: theta, p_theta
@@ -3675,6 +3675,23 @@
                 FLExit("Only pure control volume or coupled_cv discretisations can solve explicitly")
               end if
             end if
+                 
+            if(mmat_cv_disc .or. cv_disc) then
+               if (have_option("/material_phase["//int2str(m)//"]/scalar_field["//int2str(f)//&
+                    "]/prognostic/scalar_field::SinkingVelocity")) then
+                  call get_option(trim(complete_field_path("/material_phase["//&
+                       int2str(m)//"]/scalar_field["//int2str(f)//&
+                       "]/prognostic/scalar_field::SinkingVelocity"))//"/mesh[0]/name", mesh_0, stat)
+                  if(stat == 0) then
+                     call get_option(trim(complete_field_path("/material_phase[" // int2str(m) // &
+                          "]/vector_field::Velocity")) // "/mesh[0]/name", mesh_1)
+                     if(trim(mesh_0) /= trim(mesh_1)) then
+                        ewrite(0, *) "SinkingVelocity for "//trim(field_name)//&
+                             " is on a different mesh to the Velocity field this could cause problems"
+                     end if
+                  end if
+               end if
+            end if
 
          end do
        end do
 
=== modified file 'schemas/fluidity_options.rnc'
--- schemas/fluidity_options.rnc	2012-10-04 16:29:40 +0000
+++ schemas/fluidity_options.rnc	2012-11-27 18:25:25 +0000
@@ -1904,10 +1904,13 @@
                ##
                ## As with the sinking velocity, this velocity is in the direction of gravity so if the substance
                ## floats upwards, this field should be negative.
+               ##
+               ## This should be on the same mesh as the associated SinkingVelocity field.
                element scalar_field {
                   attribute name { "UnhinderedSinkingVelocity" },
                   attribute rank { "0" },
                   element prescribed {
+                     velocity_mesh_choice,
                      prescribed_scalar_field
                   }
                }?,
 
=== modified file 'schemas/fluidity_options.rng'
--- schemas/fluidity_options.rng	2012-11-05 14:25:55 +0000
+++ schemas/fluidity_options.rng	2012-11-27 18:25:25 +0000
@@ -2662,7 +2662,9 @@
 If the 'SinkingVelocity' is not set to diagnostic, this field will not be used.
 
 As with the sinking velocity, this velocity is in the direction of gravity so if the substance
-floats upwards, this field should be negative.</a:documentation>
+floats upwards, this field should be negative.
+
+This should be on the same mesh as the associated SinkingVelocity field.</a:documentation>
                 <attribute name="name">
                   <value>UnhinderedSinkingVelocity</value>
                 </attribute>
@@ -2670,6 +2672,7 @@
                   <value>0</value>
                 </attribute>
                 <element name="prescribed">
+                  <ref name="velocity_mesh_choice"/>
                   <ref name="prescribed_scalar_field"/>
                 </element>
               </element>
 
=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2012-09-29 11:32:16 +0000
+++ schemas/prognostic_field_options.rnc	2012-11-27 18:25:25 +0000
@@ -295,14 +295,20 @@
       ## 
       ## This velocity is in the direction of gravity so if the substance
       ## floats or swims upwards, this field should be negative.
+      ##
+      ## This should be on the same mesh as the velocity field unless you are 
+      ## using a DG discretisation and projecting the velocity to contiuous, in 
+      ## which case it should be on a continuous mesh.
       element scalar_field {
          attribute name { "SinkingVelocity" },
          attribute rank { "0" },
          (
             element prescribed {
+               velocity_mesh_choice,
                prescribed_scalar_field_no_adapt
             }|
             element diagnostic {
+               velocity_mesh_choice,
                (
                   scalar_python_diagnostic_algorithm |
                   internal_algorithm
 
=== modified file 'schemas/prognostic_field_options.rng'
--- schemas/prognostic_field_options.rng	2012-09-29 11:32:16 +0000
+++ schemas/prognostic_field_options.rng	2012-11-27 18:25:25 +0000
@@ -381,7 +381,11 @@
         <a:documentation>Velocity at which this substance sinks through the water column.
 
 This velocity is in the direction of gravity so if the substance
-floats or swims upwards, this field should be negative.</a:documentation>
+floats or swims upwards, this field should be negative.
+
+This should be on the same mesh as the velocity field unless you are 
+using a DG discretisation and projecting the velocity to contiuous, in 
+which case it should be on a continuous mesh.</a:documentation>
         <attribute name="name">
           <value>SinkingVelocity</value>
         </attribute>
@@ -390,9 +394,11 @@
         </attribute>
         <choice>
           <element name="prescribed">
+            <ref name="velocity_mesh_choice"/>
             <ref name="prescribed_scalar_field_no_adapt"/>
           </element>
           <element name="diagnostic">
+            <ref name="velocity_mesh_choice"/>
             <choice>
               <ref name="scalar_python_diagnostic_algorithm"/>
               <ref name="internal_algorithm"/>
 
=== modified file 'tests/biology_6_component_parallel/biology_conservation.flml'
--- tests/biology_6_component_parallel/biology_conservation.flml	2012-03-11 02:25:56 +0000
+++ tests/biology_6_component_parallel/biology_conservation.flml	2012-11-27 18:25:25 +0000
@@ -517,6 +517,7 @@
         </scalar_field>
         <scalar_field name="SinkingVelocity" rank="0">
           <prescribed>
+            <mesh name="VelocityMesh"/>
             <value name="WholeMesh">
               <python>
                 <string_value lines="20" type="code" language="python">def val(X,t):
 
=== modified file 'tests/biology_6_component_parallel/biology_conservation_parallel.flml'
--- tests/biology_6_component_parallel/biology_conservation_parallel.flml	2012-03-11 02:25:56 +0000
+++ tests/biology_6_component_parallel/biology_conservation_parallel.flml	2012-11-27 18:25:25 +0000
@@ -517,6 +517,7 @@
         </scalar_field>
         <scalar_field name="SinkingVelocity" rank="0">
           <prescribed>
+            <mesh name="VelocityMesh"/>
             <value name="WholeMesh">
               <python>
                 <string_value lines="20" type="code" language="python">def val(X,t):
 
=== modified file 'tests/biology_conservation_adapt/biology_conservation.flml'
--- tests/biology_conservation_adapt/biology_conservation.flml	2012-03-11 02:25:56 +0000
+++ tests/biology_conservation_adapt/biology_conservation.flml	2012-11-27 18:25:25 +0000
@@ -558,6 +558,7 @@
         </scalar_field>
         <scalar_field name="SinkingVelocity" rank="0">
           <prescribed>
+            <mesh name="VelocityMesh"/>
             <value name="WholeMesh">
               <python>
                 <string_value lines="20" type="code" language="python">def val(X,t):
 
=== modified file 'tests/mms_sediment/MMS_X.flml'
--- tests/mms_sediment/MMS_X.flml	2012-10-05 13:00:44 +0000
+++ tests/mms_sediment/MMS_X.flml	2012-11-27 18:25:25 +0000
@@ -1194,6 +1194,7 @@
           </scalar_field>
           <scalar_field name="SinkingVelocity" rank="0">
             <diagnostic>
+              <mesh name="VelocityMesh"/>
               <algorithm name="Internal" material_phase_support="multiple"/>
               <output/>
               <stat/>
@@ -1327,6 +1328,7 @@
           </scalar_field>
           <scalar_field name="UnhinderedSinkingVelocity" rank="0">
             <prescribed>
+              <mesh name="VelocityMesh"/>
               <value name="WholeMesh">
                 <constant>
                   <real_value rank="0">0.33</real_value>
@@ -1450,6 +1452,7 @@
           </scalar_field>
           <scalar_field name="SinkingVelocity" rank="0">
             <diagnostic>
+              <mesh name="VelocityMesh"/>
               <algorithm name="Internal" material_phase_support="multiple"/>
               <output/>
               <stat/>
@@ -1583,6 +1586,7 @@
           </scalar_field>
           <scalar_field name="UnhinderedSinkingVelocity" rank="0">
             <prescribed>
+              <mesh name="VelocityMesh"/>
               <value name="WholeMesh">
                 <constant>
                   <real_value rank="0">0.66</real_value>
 
=== modified file 'tests/sinking_velocity/sinking_velocity.flml'
--- tests/sinking_velocity/sinking_velocity.flml	2012-03-11 02:25:56 +0000
+++ tests/sinking_velocity/sinking_velocity.flml	2012-11-27 18:25:25 +0000
@@ -176,6 +176,7 @@
         </boundary_conditions>
         <scalar_field name="SinkingVelocity" rank="0">
           <prescribed>
+            <mesh name="VelocityMesh"/>
             <value name="WholeMesh">
               <python>
                 <string_value lines="20" type="code" language="python">def val(X,t):
 
=== modified file 'tests/sinking_velocity_cv/sinking_velocity.flml'
--- tests/sinking_velocity_cv/sinking_velocity.flml	2012-03-11 02:25:56 +0000
+++ tests/sinking_velocity_cv/sinking_velocity.flml	2012-11-27 18:25:25 +0000
@@ -178,6 +178,7 @@
         </boundary_conditions>
         <scalar_field name="SinkingVelocity" rank="0">
           <prescribed>
+            <mesh name="VelocityMesh"/>
             <value name="WholeMesh">
               <python>
                 <string_value lines="20" type="code" language="python">def val(X,t):