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

[Adam Candy]

Introduces a diagnostic for the buoyancy adjustment by vertical mixing scheme.
Logging for the scheme has been leashed and generally improved.

Introduces a diagnostic for the buoyancy adjustment by vertical mixing scheme.
Logging for the scheme has been leashed and generally improved.

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

4054. By Adam Candy on 2012-09-18

Tidied up commented-out code.

[Alexandros Avdis (alexandros-avdis)]

Looks good, it helps with some of my work too.

review: Approve

=== modified file 'assemble/Advection_Diffusion_DG.F90'
--- assemble/Advection_Diffusion_DG.F90	2012-09-17 19:22:44 +0000
+++ assemble/Advection_Diffusion_DG.F90	2012-09-18 16:34:19 +0000
@@ -143,6 +143,7 @@
   logical :: on_sphere
   ! Vertical diffusion by mixing option
   logical :: have_buoyancy_adjustment_by_vertical_diffusion
+  logical :: have_buoyancy_adjustment_diffusivity
 
 contains
 
@@ -751,6 +752,9 @@
     logical :: cache_valid
     integer :: num_threads
 
+    !! Diffusivity to add due to the buoyancy adjustment by vertical mixing scheme
+    type(scalar_field) :: buoyancy_adjustment_diffusivity
+
     ewrite(1,*) "Writing advection-diffusion equation for "&
          &//trim(field_name)
 
@@ -945,6 +949,7 @@
     if (present(diffusion_m)) call zero(diffusion_m)
     if (present(diffusion_RHS)) call zero(diffusion_RHS)
     if (have_buoyancy_adjustment_by_vertical_diffusion) then
+      ewrite(3,*) "Buoyancy adjustment by vertical mixing: enabled"
       if (have_option(trim(T%option_path)//"/prognostic/buoyancy_adjustment"//&
           &"/by_vertical_diffusion/project_buoyancy_to_continuous_space")) then    
         buoyancy_from_state = extract_scalar_field(state, "VelocityBuoyancyDensity", stat)
@@ -956,7 +961,9 @@
         ! Grab an extra reference to cause the deallocate below to be safe.
         ! Check this is OK
         call lumped_mass_galerkin_projection_scalar(state, buoyancy, buoyancy_from_state)
+        ewrite(3,*) "Buoyancy adjustment by vertical mixing: projecting to continuous space"
       else
+        ewrite(3,*) "Buoyancy adjustment by vertical mixing: no projection"
         buoyancy = extract_scalar_field(state, "VelocityBuoyancyDensity", stat)
         if (stat/=0) FLAbort('Error extracting buoyancy field.')
         call incref(buoyancy)
@@ -976,6 +983,16 @@
       end if
       ! Set direction of mixing diffusion, default is in the y- and z-direction for 2- and 3-d spaces respectively
       ! TODO: Align this direction with gravity local to an element
+      
+      ! Check if the diagnostic associated with the buoyancy adjustment by vertical mixing scheme is required
+      buoyancy_adjustment_diffusivity = extract_scalar_field(state, "BuoyancyAdjustmentDiffusivity", stat)
+      if (stat==0) then
+        have_buoyancy_adjustment_diffusivity = .true.
+        ewrite(3,*) "Buoynacy adjustment by vertical mixing: Updating BuoyancyAdjustmentDiffusivity field."
+      else
+        have_buoyancy_adjustment_diffusivity = .false.
+      end if
+
     end if
 
     if (include_diffusion) then
@@ -1011,6 +1028,7 @@
             & rhs_diff, X, X_old, X_new, T, U_nl, U_mesh, Source, &
             & Absorption, Diffusivity, bc_value, bc_type, q_mesh, mass, &
             & buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude, &
+            & buoyancy_adjustment_diffusivity, &
             & add_src_directly_to_rhs, porosity_theta) 
        
       end do element_loop
@@ -1127,6 +1145,7 @@
        & X, X_old, X_new, T, U_nl, U_mesh, Source, Absorption, Diffusivity,&
        & bc_value, bc_type, &
        & q_mesh, mass, buoyancy, gravity, gravity_magnitude, mixing_diffusion_amplitude, &
+       & buoyancy_adjustment_diffusivity, &
        & add_src_directly_to_rhs, porosity_theta)
     !!< Construct the advection_diffusion equation for discontinuous elements in
     !!< acceleration form.
@@ -1154,6 +1173,9 @@
     type(scalar_field), intent(in) :: T, Source, Absorption
     !! Diffusivity
     type(tensor_field), intent(in) :: Diffusivity
+
+    !! Diffusivity to add due to the buoyancy adjustment by vertical mixing scheme
+    type(scalar_field), intent(inout) :: buoyancy_adjustment_diffusivity
     
     !! If adding Source directly to rhs then
     !! do nothing with it here
@@ -1423,12 +1445,21 @@
                  &* gravity_magnitude * dr**2 * gravity_at_node(i) * drho_dz(:)
          end do
        end if
-       
-       !ewrite(3,*) "mixing_density_values", buoyancysample
-       ewrite(3,*) "mixing_grad_rho", minval(grad_rho(:,:)), maxval(grad_rho(:,:))
-       ewrite(3,*) "mixing_drho_dz", minval(drho_dz(:)), maxval(drho_dz(:))
-       ewrite(3,*) "mixing_coeffs amp dt g dr", mixing_diffusion_amplitude, dt, gravity_magnitude, dr**2
-       ewrite(3,*) "mixing_diffusion", minval(mixing_diffusion(2,2,:)), maxval(mixing_diffusion(2,2,:))
+        
+       if(have_buoyancy_adjustment_diffusivity) then
+         ewrite (-1,*) 'asc nodes', T_ele
+         do i = 1,size(T_ele)
+           call set(buoyancy_adjustment_diffusivity, T_ele(i), mixing_diffusion_amplitude * dt&
+                 &* gravity_magnitude * dr**2 * maxval(drho_dz(:)))
+         end do
+         ewrite(4,*) "Buoynacy adjustment diffusivity, ele:", ele, "diffusivity:", mixing_diffusion_amplitude * dt * gravity_magnitude * dr**2 * maxval(drho_dz(:))
+       end if 
+        
+       !! Buoyancy adjustment by vertical mixing scheme debugging statements
+       ewrite(4,*) "mixing_grad_rho", minval(grad_rho(:,:)), maxval(grad_rho(:,:))
+       ewrite(4,*) "mixing_drho_dz", minval(drho_dz(:)), maxval(drho_dz(:))
+       ewrite(4,*) "mixing_coeffs amp dt g dr", mixing_diffusion_amplitude, dt, gravity_magnitude, dr**2
+       ewrite(4,*) "mixing_diffusion", minval(mixing_diffusion(2,2,:)), maxval(mixing_diffusion(2,2,:))
 
        mixing_diffusion_rhs=shape_tensor_rhs(T%mesh%shape, mixing_diffusion, detwei_rho)
        t_mass=shape_shape(T%mesh%shape, T%mesh%shape, detwei_rho)

=== modified file 'schemas/fluidity_options.rnc'
--- schemas/fluidity_options.rnc	2012-09-03 09:42:37 +0000
+++ schemas/fluidity_options.rnc	2012-09-18 16:34:19 +0000
@@ -4185,6 +4185,20 @@
                }
             )
          }|
+         element scalar_field {
+            attribute rank { "0" },
+            attribute name { "BuoyancyAdjustmentDiffusivity" },
+            (
+               element diagnostic {
+                  internal_algorithm,
+                  velocity_mesh_choice,
+                  diagnostic_scalar_field
+               }|
+               element aliased {
+                  generic_aliased_field
+               }
+            )
+         }|
          ## The equilibrium pressure perturbation 
          ## 
          ## Further details to be added

=== modified file 'schemas/fluidity_options.rng'
--- schemas/fluidity_options.rng	2012-09-03 09:42:37 +0000
+++ schemas/fluidity_options.rng	2012-09-18 16:34:19 +0000
@@ -5368,6 +5368,24 @@
         </choice>
       </element>
       <element name="scalar_field">
+        <attribute name="rank">
+          <value>0</value>
+        </attribute>
+        <attribute name="name">
+          <value>BuoyancyAdjustmentDiffusivity</value>
+        </attribute>
+        <choice>
+          <element name="diagnostic">
+            <ref name="internal_algorithm"/>
+            <ref name="velocity_mesh_choice"/>
+            <ref name="diagnostic_scalar_field"/>
+          </element>
+          <element name="aliased">
+            <ref name="generic_aliased_field"/>
+          </element>
+        </choice>
+      </element>
+      <element name="scalar_field">
         <a:documentation>The equilibrium pressure perturbation 
 
 Further details to be added