Merge pull request #334 from FluidityProject/parallel-non-extruded-ocean-boundaries

Parallel non extruded ocean boundaries

Author: stephankramer

assemble/Adapt_State.F90

@@ -1235,19 +1235,17 @@ subroutine adapt_state_internal(states, metric, initialise_fields)
         call profiler_tic("find_particles_mesh_adapt")
         do j = 1, size(detector_list_array)
            call search_for_detectors(detector_list_array(j)%ptr, new_positions)
-        end do
-        call profiler_toc("find_particles_mesh_adapt")
-#ifdef DDEBUG
-        ! Sanity check that all local detectors are owned
-        call get_registered_detector_lists(detector_list_array)
-        do j = 1, size(detector_list_array)
+
+           ! Sanity check that all local detectors are owned
            detector=>detector_list_array(j)%ptr%first
            do k = 1, detector_list_array(j)%ptr%length
-              assert(detector%element>0)
+              if (detector%element<=0) then
+                FLAbort("Lost one of the detectors during an adapt")
+              end if
               detector=>detector%next
            end do
         end do
-#endif
+        call profiler_toc("find_particles_mesh_adapt")
       end if
 
       ! Then reallocate all fields

assemble/Discrete_Properties.F90

@@ -58,7 +58,6 @@ subroutine enforce_discrete_properties(states, only_prescribed, exclude_interpol
     type(state_type) :: alg_state
 
     integer :: state, state_cnt, mesh_i
-    integer :: mesh_cnt
 
     character(len=255), dimension(1), parameter :: algorithms = (/&
                        & "solenoidal" /)
@@ -71,7 +70,6 @@ subroutine enforce_discrete_properties(states, only_prescribed, exclude_interpol
 
     ewrite(1, *) "In enforce_discrete_properties"
 
-    mesh_cnt = option_count("/geometry/mesh")
     alg_cnt = size(algorithms)
     state_cnt = size(states)
 
@@ -85,7 +83,7 @@ subroutine enforce_discrete_properties(states, only_prescribed, exclude_interpol
         ! only actually care about one algorithm at the moment but for futureproofing we'll do this in a loop
         ewrite(2, *) "  Considering algorithm " // trim(algorithms(alg))
 
-        do mesh_i = 1, mesh_cnt
+        do mesh_i = 1, mesh_count(states(state))
           mesh => extract_mesh(states(state), mesh_i)
           call insert(alg_state, mesh, name=trim(mesh%name))
         end do

assemble/Free_Surface.F90

@@ -45,9 +45,9 @@ module free_surface_module
   use state_module
   use sparse_matrices_fields
   use boundary_conditions
-  use vertical_extrapolation_module
   use halos
   use field_options
+  use vertical_extrapolation_module
   use physics_from_options
   use coordinates
   use tidal_module, only: calculate_diagnostic_equilibrium_pressure
@@ -223,25 +223,28 @@ subroutine add_free_surface_to_cmc_projection(state, cmc, dt, &
   !!< that weakly enforces the kinematic boundary condition.
   !!<
   !!< The free surface is combined in with the pressure field such that
-  !!< *at* the free surface p=g rho0 \eta. This has the advantage of combining
+  !!< *at* the free surface p=g \Delta rho \eta. This has the advantage of combining
   !!< the pressure gradient and free surface gradient terms in the momentum
   !!< equation. It solves the continuity equation directly coupled with
   !!< the free surface.
-  !!< With this approach all pressures are considered at the integer time 
+  !!< With this approach all pressures are considered at the integer time
   !!< levels, i.e. we apply a theta weighting for the pressure gradient term
   !!< in the momentum equation:
-  !!<   M (u^n+1-u^n) - dt C p^{n+theta_pressure_gradient} + ... = 0
+  !!<   M (u^n+1-u^n) - dt C p^{n+theta_pg} + ... = 0
   !!< We're solving the continuity equation:
-  !!<   C^T u^{n+theta_divergence}+ M_fs p^{n+1}-p^n = 0
+  !!<   C^T u^{n+theta_div}+ alpha M_fs p^{n+1}-p^n = 0
   !!< which leads to a projection equation of:
-  !!<   ( C^T M^-1 C dt dp + coef M_fs ) phi = 
-  !!<     theta_divergence C^T u* + (1-theta_divergence) C^T u^n - 
+  !!<   ( C^T M^-1 C dt dp + coef M_fs ) phi =
+  !!<     theta_div C^T u* + (1-theta_div) C^T u^n -
   !!<     alpha M_fs (p*-p^n)
-  !!< where M_fs is the free surface integral of M_i M_j, 
-  !!< alpha=1/(g dt), coef=alpha/(theta_divergence theta_pressure_gradient dt)
-  !!< and phi=dp theta_divergence theta_pressure_gradient dt. Note however,
-  !!< that dp in the routine stands for phi, see correct_pressure in 
-  !!< Momentum_Equation.F90
+  !!< where M_fs is the free surface integral of M_i M_j,
+  !!< alpha=1/(g dt), coef=alpha/(theta_div theta_pg dt)
+  !!< and phi=dp theta_div theta_pg dt and dp=p^n+1-p^*
+  !!< Notes:
+  !!< - the above is slightly misleading because with a prognostic f.s. we only apply theta weighting
+  !!< to the free surface (not interior pressure), but for the purposes of this routine we can ignore that
+  !!< - the above assumes \Delta rho=1 but should be included in alpha M_fs (divided by). For 
+  !!< this will be included in the mass matrix, otherwise we apply it as an overall factor 1/\Delta rho
 
     type(state_type), intent(inout) :: state
     type(csr_matrix), intent(inout) :: cmc
@@ -413,10 +416,6 @@ subroutine add_free_surface_to_cmc_projection(state, cmc, dt, &
           else
             delta_rho = rho0
           end if
-          if (move_mesh .and. delta_rho/=1.0) then
-            ! Someone needs to go through the maths to see where we should divide by delta_rho
-            FLExit("Free surface with a density difference that is not 1.0 and mesh movement not supported")
-          end if
 
           alpha=1.0/g/dt                        ! delta_rho included in alpha and coeff within element loop
           coef = alpha/(theta_pressure_gradient*theta_divergence*dt)
@@ -570,9 +569,10 @@ subroutine add_free_surface_element(surface_mesh_ele, sele)
         end if
 
         if(present(rhs)) then
+           assert(.not. variable_density)
            call addto(rhs, nodes, &
                 -(matmul(mass_ele, top_pressures) &
-                -matmul(mass_ele_old, old_top_pressures))*alpha)
+                -matmul(mass_ele_old, old_top_pressures))*alpha/delta_rho)
         end if
         if (have_wd .and. present(rhs)) then
            call addto(rhs, nodes, &

assemble/Geostrophic_Pressure.F90

@@ -2428,7 +2428,7 @@ subroutine initialise_geostrophic_interpolation_states(old_states, new_states)
     type(vector_field), pointer :: new_velocity, old_velocity
 
     do i = 1, size(new_states)
-      do j = 1, vector_field_count(old_states(i))
+      do j = 1, vector_field_count(new_states(i))
         new_velocity => extract_vector_field(new_states(i), j)
         if(have_option(trim(complete_field_path(new_velocity%option_path, stat = stat)) // "/geostrophic_interpolation")) then
           old_velocity => extract_vector_field(old_states(i), new_velocity%name)

femtools/Boundary_Conditions.F90

@@ -83,9 +83,16 @@ module boundary_conditions
   end interface extract_scalar_surface_field
 
   interface has_surface_field
-     module procedure has_scalar_surface_field_specific, &
-       has_vector_surface_field_specific
+     module procedure has_scalar_surface_field_by_bc_number, &
+       has_vector_surface_field_by_bc_number, &
+       has_scalar_surface_field_by_bc_name, &
+       has_vector_surface_field_by_bc_name
   end interface
+
+  interface has_scalar_surface_field
+    module procedure vector_has_scalar_surface_field_by_bc_number, &
+      vector_has_scalar_surface_field_by_bc_name
+  end interface has_scalar_surface_field
 
   interface get_boundary_condition_count
      module procedure get_scalar_boundary_condition_count, &
@@ -792,9 +799,9 @@ function extract_vector_scalar_surface_field_by_name(field, bc_name, name, stat)
 
   end function extract_vector_scalar_surface_field_by_name
 
-  function has_scalar_surface_field_specific(field, n, name)
+  function has_scalar_surface_field_by_bc_number(field, n, name)
   !!< Tells whether a surface_field with the given is present
-  logical :: has_scalar_surface_field_specific
+  logical :: has_scalar_surface_field_by_bc_number
   type(scalar_field), intent(in):: field
   integer, intent(in):: n
   character(len=*), intent(in):: name
@@ -808,19 +815,19 @@ function has_scalar_surface_field_specific(field, n, name)
     if (associated(bc%surface_fields)) then
       do i=1, size(bc%surface_fields)
         if (bc%surface_fields(i)%name==name) then
-          has_scalar_surface_field_specific=.true.
+          has_scalar_surface_field_by_bc_number=.true.
           return
         end if
       end do
     end if
 
-    has_scalar_surface_field_specific=.false.
+    has_scalar_surface_field_by_bc_number=.false.
 
-  end function has_scalar_surface_field_specific
+  end function has_scalar_surface_field_by_bc_number
 
-  function has_vector_surface_field_specific(field, n, name)
+  function has_vector_surface_field_by_bc_number(field, n, name)
   !!< Tells whether a surface_field with the given is present
-  logical :: has_vector_surface_field_specific
+  logical :: has_vector_surface_field_by_bc_number
   type(vector_field), intent(in):: field
   integer, intent(in):: n
   character(len=*), intent(in):: name
@@ -834,19 +841,19 @@ function has_vector_surface_field_specific(field, n, name)
     if (associated(bc%surface_fields)) then
       do i=1, size(bc%surface_fields)
         if (bc%surface_fields(i)%name==name) then
-          has_vector_surface_field_specific=.true.
+          has_vector_surface_field_by_bc_number=.true.
           return
         end if
       end do
     end if
 
-    has_vector_surface_field_specific=.false.
+    has_vector_surface_field_by_bc_number=.false.
 
-  end function has_vector_surface_field_specific
+  end function has_vector_surface_field_by_bc_number
 
-  function has_scalar_surface_field(field, n, name)
+  function vector_has_scalar_surface_field_by_bc_number(field, n, name)
   !!< Tells whether a surface_field with the given is present
-  logical :: has_scalar_surface_field
+  logical :: vector_has_scalar_surface_field_by_bc_number
   type(vector_field), intent(in):: field
   integer, intent(in):: n
   character(len=*), intent(in):: name
@@ -860,16 +867,79 @@ function has_scalar_surface_field(field, n, name)
     if (associated(bc%scalar_surface_fields)) then
       do i=1, size(bc%scalar_surface_fields)
         if (bc%scalar_surface_fields(i)%name==name) then
-          has_scalar_surface_field=.true.
+          vector_has_scalar_surface_field_by_bc_number=.true.
           return
         end if
       end do
     end if
 
-    has_scalar_surface_field=.false.
+    vector_has_scalar_surface_field_by_bc_number=.false.
 
-  end function has_scalar_surface_field
+  end function vector_has_scalar_surface_field_by_bc_number
+  
+  function has_scalar_surface_field_by_bc_name(field, bc_name, name)
+  !!< Tells whether a surface_field with the given name is present
+  !!< If the bc_name does not exist an error is given
+  logical :: has_scalar_surface_field_by_bc_name
+  type(scalar_field), intent(in):: field
+  character(len=*), intent(in):: bc_name, name
 
+    integer i
+
+    do i=1, size(field%bc%boundary_condition)
+      if (field%bc%boundary_condition(i)%name==bc_name) then
+        has_scalar_surface_field_by_bc_name = has_scalar_surface_field_by_bc_number(field, i, name)
+        return
+      end if
+    end do
+
+    ewrite(-1,*) 'Unknown boundary condition: ', name
+    FLAbort("Sorry!")
+
+  end function has_scalar_surface_field_by_bc_name
+
+  function has_vector_surface_field_by_bc_name(field, bc_name, name)
+  !!< Tells whether a surface_field with the given name is present
+  !!< If the bc_name does not exist an error is given
+  logical :: has_vector_surface_field_by_bc_name
+  type(vector_field), intent(in):: field
+  character(len=*), intent(in):: bc_name, name
+
+    integer i
+
+    do i=1, size(field%bc%boundary_condition)
+      if (field%bc%boundary_condition(i)%name==bc_name) then
+        has_vector_surface_field_by_bc_name = has_vector_surface_field_by_bc_number(field, i, name)
+        return
+      end if
+    end do
+
+    ewrite(-1,*) 'Unknown boundary condition: ', name
+    FLAbort("Sorry!")
+
+  end function has_vector_surface_field_by_bc_name
+
+  function vector_has_scalar_surface_field_by_bc_name(field, bc_name, name)
+  !!< Tells whether a scalar surface_field with the given name is present under a vector field bc
+  !!< If the bc_name does not exist an error is given
+  logical :: vector_has_scalar_surface_field_by_bc_name
+  type(vector_field), intent(in):: field
+  character(len=*), intent(in):: bc_name, name
+
+    integer i
+
+    do i=1, size(field%bc%boundary_condition)
+      if (field%bc%boundary_condition(i)%name==bc_name) then
+        vector_has_scalar_surface_field_by_bc_name = vector_has_scalar_surface_field_by_bc_number(field, i, name)
+        return
+      end if
+    end do
+
+    ewrite(-1,*) 'Unknown boundary condition: ', name
+    FLAbort("Sorry!")
+
+  end function vector_has_scalar_surface_field_by_bc_name
+
   integer function get_scalar_boundary_condition_count(field)
   !!< Get number of boundary conditions of a scalar field
   type(scalar_field), intent(in):: field

femtools/Detector_Move_Lagrangian.F90

@@ -231,7 +231,8 @@ subroutine move_lagrangian_detectors(state, detector_list, dt)
 
              !This call serialises send_list_array, sends it, 
              !receives serialised receive_list_array, and unserialises that.
-             call exchange_detectors(state(1),detector_list, send_list_array)
+             call exchange_detectors(state(1),detector_list, send_list_array, &
+               include_update_vector=.true.)
 
           end do detector_timestepping_loop
        end do RKstages_loop