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

[Samuel Parkinson]

Changed bedload field to prognostic. I think it is a prognostic field really despite being trivial to solve and this way:
A) it can be checkpointed
B) we can start with an abritrary initial condition

I know you said you didn't like this solution but it seems the simplest option. What do you think?

[Jon Hill (jon-hill)]

Approve, but grudgingly. The checkpointing *should* work, but I agree that an abritrary initial condition is required here.

There might be a better way in the long term.

An options check for a CV field would be useful though. Can this be added?

[Samuel Parkinson (s-parkinson11)]

Thanks again Jon. Yes I agree that perhaps this isn't the nicest way to do this. I am not totally sure what you mean by an options check for a CV field, can you elaborate at all?

[Jon Hill (jon-hill)]

The schema says:
No solver settings are required as the sediment model only works
 for p1 CG or DG fields where the mass matrix can be easily inverted.

So if you have a CV field this won't work. You should check which discretisation has been set and error if CV.

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

3488. By Samuel Parkinson on 2012-10-04

fixed error in schema description

[Samuel Parkinson (s-parkinson11)]

This was actually an error in the description. It will work for CV just not for >P1 discretisations

[Jon Hill (jon-hill)]

OK, merge it then.

review: Approve

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

3489. By Samuel Parkinson on 2012-10-05

changed mms test schema so that bedload is not calculated

3490. By Samuel Parkinson on 2012-10-08

merge with trunk

3491. By Samuel Parkinson on 2012-10-11

merge with trunk

3492. By Samuel Parkinson on 2012-10-22

removed some multiphase compatibility due to bug in flredecomp

3493. By Samuel Parkinson on 2012-10-24

merge with trunk

=== modified file 'schemas/fluidity_options.rnc'
--- schemas/fluidity_options.rnc	2012-09-25 06:55:02 +0000
+++ schemas/fluidity_options.rnc	2012-10-24 15:56:21 +0000
@@ -1829,24 +1829,33 @@
                ## deposited through the prescribed boundary.
                ##
                ## The units are in unit distances as a depth of sediment.
+               ##
+               ## - The equation type must be set to SedimentBedload
+               ##
+               ## - Spatial discretisation should be as the parent sediment class
+               ##
+               ## - Initial conditions can be set as normal
+               ##
+               ## - Most other settings are irrelevant
+               ##     - No solver settings are required as the sediment model only works
+               ##       for p1 discretisations where the mass matrix can be easily inverted.
+               ##     - The bedload is always calculated completely explicitly
+               ##     - You cannot apply Source or Absorbtion terms to the bedload
+               ##     - There is no diffusion of the bedload
                element scalar_field {
                   attribute rank { "0" },
                   attribute name { "Bedload" },
                   ## Field type
                   (
-                     element diagnostic {
+                     element prognostic {
                         ## Surface ids over which to calculate the bedload:
                         element surface_ids {
                            integer_vector
                         },
-                        diagnostic_scalar_field
-                     }|
-                     element prescribed {
-                        ## Surface ids over which to calculate the bedload deposition:
-                        element surface_ids {
-                           integer_vector
-                        },
-                        prescribed_scalar_field
+                        prognostic_scalar_field,
+                        ## Disables calculation of change to bedload due to 
+                        ## erosion/deposition (generally only used for mms testing):
+                        element disable_calculation { empty }?
                      }
                   )
                },

=== modified file 'schemas/fluidity_options.rng'
--- schemas/fluidity_options.rng	2012-09-25 06:55:02 +0000
+++ schemas/fluidity_options.rng	2012-10-24 15:56:21 +0000
@@ -2562,30 +2562,43 @@
               <a:documentation>Sediment bedload diagnostic which records the sediment 
 deposited through the prescribed boundary.
 
-The units are in unit distances as a depth of sediment.</a:documentation>
+The units are in unit distances as a depth of sediment.
+
+- The equation type must be set to SedimentBedload
+
+- Spatial discretisation should be as the parent sediment class
+
+- Initial conditions can be set as normal
+
+- Most other settings are irrelevant
+    - No solver settings are required as the sediment model only works
+      for p1 CG or DG fields where the mass matrix can be easily inverted.
+    - The bedload is always calculated completely explicitly
+    - You cannot apply Source or Absorbtion terms to the bedload
+    - There is no diffusion of the bedload</a:documentation>
               <attribute name="rank">
                 <value>0</value>
               </attribute>
               <attribute name="name">
                 <value>Bedload</value>
               </attribute>
-              <choice>
+              <group>
                 <a:documentation>Field type</a:documentation>
-                <element name="diagnostic">
+                <element name="prognostic">
                   <element name="surface_ids">
                     <a:documentation>Surface ids over which to calculate the bedload:</a:documentation>
                     <ref name="integer_vector"/>
                   </element>
-                  <ref name="diagnostic_scalar_field"/>
-                </element>
-                <element name="prescribed">
-                  <element name="surface_ids">
-                    <a:documentation>Surface ids over which to calculate the bedload deposition:</a:documentation>
-                    <ref name="integer_vector"/>
-                  </element>
-                  <ref name="prescribed_scalar_field"/>
-                </element>
-              </choice>
+                  <ref name="prognostic_scalar_field"/>
+                  <optional>
+                    <element name="disable_calculation">
+                      <a:documentation>Disables calculation of change to bedload due to 
+erosion/deposition (generally only used for mms testing):</a:documentation>
+                      <empty/>
+                    </element>
+                  </optional>
+                </element>
+              </group>
             </element>
             <element name="scalar_field">
               <a:documentation>Sediment bedload diagnostic which records the sediment 

=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2012-09-18 08:03:27 +0000
+++ schemas/prognostic_field_options.rnc	2012-10-24 15:56:21 +0000
@@ -2068,6 +2068,11 @@
          ## Note: Only works for continuous galerkin or control volume discretisations.
          element equation {
             attribute name { "KEpsilon" }
+         }|
+         ## Equation type for the sediment bedload fields.
+         ## Note: Only works for P1 discretisations.
+         element equation {
+            attribute name { "SedimentBedload" }
          }
       )
    )

=== modified file 'schemas/prognostic_field_options.rng'
--- schemas/prognostic_field_options.rng	2012-09-18 08:03:27 +0000
+++ schemas/prognostic_field_options.rng	2012-10-24 15:56:21 +0000
@@ -2553,6 +2553,13 @@
           <value>KEpsilon</value>
         </attribute>
       </element>
+      <element name="equation">
+        <a:documentation>Equation type for the sediment bedload fields.
+Note: Only works for P1 discretisations.</a:documentation>
+        <attribute name="name">
+          <value>SedimentBedload</value>
+        </attribute>
+      </element>
     </choice>
   </define>
   <define name="equation_coefficients">

=== modified file 'sediments/Sediment.F90'
--- sediments/Sediment.F90	2012-09-17 13:35:07 +0000
+++ sediments/Sediment.F90	2012-10-24 15:56:21 +0000
@@ -81,7 +81,8 @@
     integer, intent(out), optional              :: stat
     character(len=FIELD_NAME_LEN)               :: name
 
-    call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
+	! had to remove trim(state%option_path)// as this didn't work with flredecomp
+    call get_option('/material_phase[0]/sediment/scalar_field['//int2str(i_field -&
          & 1)//']/name', name) 
     item => extract_scalar_field(state, trim(name), stat)
 
@@ -95,7 +96,8 @@
     character(len=FIELD_NAME_LEN), intent(out)  :: item
     integer, intent(out), optional              :: stat
 
-    call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
+	! had to remove trim(state%option_path)// as this didn't work with flredecomp
+    call get_option('/material_phase[0]/sediment/scalar_field['//int2str(i_field -&
          & 1)//']/name', item, stat=stat) 
 
   end subroutine get_sediment_field_name
@@ -109,7 +111,8 @@
     character(len=FIELD_NAME_LEN), intent(out)  :: item
     integer, intent(out), optional              :: stat
 
-    call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
+	! had to remove trim(state%option_path)// as this didn't work with flredecomp
+    call get_option('/material_phase[0]/sediment/scalar_field['//int2str(i_field -&
          & 1)//']/prognostic/'//option, item, stat=stat) 
 
   end subroutine get_sediment_option_string
@@ -124,7 +127,8 @@
     integer, intent(out), optional              :: stat
     real, intent(in), optional                  :: default
     
-    call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
+	! had to remove trim(state%option_path)// as this didn't work with flredecomp
+    call get_option('/material_phase[0]/sediment/scalar_field['//int2str(i_field -&
          & 1)//']/prognostic/'//option, item, stat = stat, default = default) 
 
   end subroutine get_sediment_option_real
@@ -140,7 +144,8 @@
     
     character(len=FIELD_NAME_LEN)               :: field_name
 
-    call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
+	! had to remove trim(state%option_path)// as this didn't work with flredecomp
+    call get_option('/material_phase[0]/sediment/scalar_field['//int2str(i_field -&
          & 1)//']/name', field_name) 
     item => extract_scalar_field(state, trim(field_name)//option, stat)
 
@@ -499,7 +504,6 @@
 
     character(len=FIELD_NAME_LEN)           :: field_mesh, sediment_mesh, bc_type
     character(len=OPTION_PATH_LEN)          :: field_option_path 
-    character(len=10)                       :: type
     integer                                 :: i_field, i_bc, i_bc_surf, i_bedload_surf,&
          & n_sediment_fields, nbcs
     integer, dimension(2)                   :: bc_surface_id_count, bedload_surface_id_count
@@ -536,10 +540,6 @@
           ! Loop over boundary conditions for field
           boundary_conditions: do i_bc=0, nbcs-1
 
-             ! get name and type of boundary condition
-             call get_option(trim(field_option_path)//'/boundary_conditions['//int2str(i_bc)//&
-                  &']/type[0]/name', bc_type)
-
              ! check whether this is a reentrainment boundary
              if (.not. (trim(bc_type) .eq. "sediment_reentrainment")) then
                 cycle boundary_conditions
@@ -549,22 +549,15 @@
              if (.not.(have_option('/material_phase[0]/vector_field::BedShearStress'))) then
                 FLExit("Reentrainment boundary condition requires a BedShearStress field")
              end if
-
-             ! warn if bedload field is prescribed
-             type = 'diagnostic'
-             if (have_option(trim(field_option_path)//'/scalar_field::Bedload/prescribed')) then
-                ewrite(0,*) 'WARNING: Bedload field is prescribed'
-                type = 'prescribed'
-             end if
              
              ! check boundary id's are the same for re-entrainment and bedload
              
              ! get bedload surface ids
              bedload_surface_id_count=option_shape(trim(field_option_path)// &
-                  '/scalar_field::Bedload/'//type//'/surface_ids')
+                  '/scalar_field::Bedload/prognostic/surface_ids')
              allocate(bedload_surface_ids(bedload_surface_id_count(1)))
              call get_option(trim(field_option_path)// &
-                  '/scalar_field::Bedload/'//type//'/surface_ids', bedload_surface_ids) 
+                  '/scalar_field::Bedload/prognostic/surface_ids', bedload_surface_ids) 
 
              ! get reentrainment surface ids
              bc_surface_id_count=option_shape(trim(field_option_path)//'/boundary_conditions['&

=== modified file 'sediments/Sediment_Diagnostics.F90'
--- sediments/Sediment_Diagnostics.F90	2012-10-22 19:47:09 +0000
+++ sediments/Sediment_Diagnostics.F90	2012-10-24 15:56:21 +0000
@@ -173,19 +173,11 @@
        end if
        
        ! obtain surface ids over which to record deposition
-       if (have_option(trim(bedload_field%option_path)//"/diagnostic")) then
-          surface_id_count=option_shape(trim(bedload_field%option_path)//&
-               &"/diagnostic/surface_ids")
-          allocate(surface_ids(surface_id_count(1)))
-          call get_option(trim(bedload_field%option_path)//"/diagnostic/surface_ids", &
-               & surface_ids)
-       else
-          surface_id_count=option_shape(trim(bedload_field%option_path)//&
-               &"/prescribed/surface_ids")
-          allocate(surface_ids(surface_id_count(1)))
-          call get_option(trim(bedload_field%option_path)//"/prescribed/surface_ids", &
-               & surface_ids)
-       end if
+       surface_id_count=option_shape(trim(bedload_field%option_path)//&
+            &"/prognostic/surface_ids")
+       allocate(surface_ids(surface_id_count(1)))
+       call get_option(trim(bedload_field%option_path)//"/prognostic/surface_ids", &
+            & surface_ids)
 
        ! loop through elements in surface field
        elements: do ele=1,element_count(deposited_sediment(i_field))
@@ -243,7 +235,7 @@
           call scale(diagnostic_field, 1./dt)
        end if
        
-       if (.not. have_option(trim(bedload_field%option_path)//'/prescribed')) then
+       if (.not. have_option(trim(bedload_field%option_path)//'/prognostic/disable_calculation')) then
           ! Add on sediment falling in and subtract sediment coming out
           do i_node = 1, node_count(surface_mesh(i_field))
              ! add deposited sediment
@@ -687,19 +679,11 @@
        end if
 
        ! obtain sediment bedload surface ids
-       if (have_option(trim(bedload%option_path)//"/diagnostic")) then
-          surface_id_count=option_shape(trim(bedload%option_path)//&
-               &"/diagnostic/surface_ids")
-          allocate(surface_ids(surface_id_count(1)))
-          call get_option(trim(bedload%option_path)//"/diagnostic/surface_ids", &
-               & surface_ids)
-       else
-          surface_id_count=option_shape(trim(bedload%option_path)//&
-               &"/prescribed/surface_ids")
-          allocate(surface_ids(surface_id_count(1)))
-          call get_option(trim(bedload%option_path)//"/prescribed/surface_ids", &
-               & surface_ids)
-       end if
+       surface_id_count=option_shape(trim(bedload%option_path)//&
+            &"/prognostic/surface_ids")
+       allocate(surface_ids(surface_id_count(1)))
+       call get_option(trim(bedload%option_path)//"/prognostic/surface_ids", &
+            & surface_ids)
        
        ! loop through elements in surface field
        elements: do i_ele=1, element_count(volume_fraction_surface)

=== modified file 'tests/mms_sediment/MMS_X.flml'
--- tests/mms_sediment/MMS_X.flml	2012-08-30 10:34:23 +0000
+++ tests/mms_sediment/MMS_X.flml	2012-10-24 15:56:21 +0000
@@ -1215,23 +1215,70 @@
           </steady_state>
           <consistent_interpolation/>
           <scalar_field name="Bedload" rank="0">
-            <prescribed>
+            <prognostic>
               <surface_ids>
                 <integer_value shape="1" rank="1">2</integer_value>
               </surface_ids>
-              <value name="WholeMesh">
+              <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-6</real_value>
+                </relative_error>
+                <absolute_error>
+                  <real_value rank="0">1e-15</real_value>
+                </absolute_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):
    import mms_sediment_tools as k
    return k.s1_d(X)</string_value>
                 </python>
-              </value>
+              </initial_condition>
               <output/>
               <stat/>
+              <convergence>
+                <include_in_convergence/>
+              </convergence>
               <detectors>
-                <exclude_from_detectors/>
+                <include_in_detectors/>
               </detectors>
-            </prescribed>
+              <steady_state>
+                <include_in_steady_state/>
+              </steady_state>
+              <consistent_interpolation/>
+              <disable_calculation/>
+            </prognostic>
           </scalar_field>
           <scalar_field name="BedloadVolumeFraction" rank="0">
             <diagnostic>
@@ -1424,23 +1471,70 @@
           </steady_state>
           <consistent_interpolation/>
           <scalar_field name="Bedload" rank="0">
-            <prescribed>
+            <prognostic>
               <surface_ids>
                 <integer_value shape="1" rank="1">2</integer_value>
               </surface_ids>
-              <value name="WholeMesh">
+              <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-6</real_value>
+                </relative_error>
+                <absolute_error>
+                  <real_value rank="0">1e-15</real_value>
+                </absolute_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):
    import mms_sediment_tools as k
    return k.s2_d(X)</string_value>
                 </python>
-              </value>
+              </initial_condition>
               <output/>
               <stat/>
+              <convergence>
+                <include_in_convergence/>
+              </convergence>
               <detectors>
-                <exclude_from_detectors/>
+                <include_in_detectors/>
               </detectors>
-            </prescribed>
+              <steady_state>
+                <include_in_steady_state/>
+              </steady_state>
+              <consistent_interpolation/>
+              <disable_calculation/>
+            </prognostic>
           </scalar_field>
           <scalar_field name="BedloadVolumeFraction" rank="0">
             <diagnostic>