Merge pull request #220 from edoddridge/halos

Halos - making better use of halos to improve parallel performance

aronnax.conf

@@ -103,7 +103,8 @@ rho0 = 1035.
 # nx is the number of grid points in the x direction
 # ny is the number of grid points in the y direction
 # layers is the number of active layers
-# OL is the width of the halo region around tiles (must be set to 1)
+# OL is the width of the halo region around tiles
+#   must be set to >= 3 if using multiple tiles
 # dx is the x grid spacing in metres
 # dy is the y grid spacing in metres
 # fUfile defines the Coriolis parameter on the u grid points in 1/s
@@ -116,7 +117,7 @@ rho0 = 1035.
 nx = 10
 ny = 10
 layers = 2
-OL = 1
+OL = 3
 dx = 2e4
 dy = 2e4
 fUfile = :beta_plane_f_u:1e-5,2e-11

src/advection_schemes.f90

@@ -25,8 +25,8 @@ subroutine h_advec_1_centered(dhdt_advec, h, u, v, dx, dy, &
     dhdt_advec = 0d0
 
     do k = 1, layers
-      do j = ylow, yhigh
-        do i = xlow, xhigh
+      do j = ylow-OL+1, yhigh+OL-1
+        do i = xlow-OL+1, xhigh+OL-1
           dhdt_advec(i,j,k) = & 
               - ((h(i,j,k)+h(i+1,j,k))*u(i+1,j,k) &
                - (h(i-1,j,k)+h(i,j,k))*u(i,j,k))/(dx*2d0) & ! d(hu)/dx
@@ -60,8 +60,8 @@ subroutine h_advec_1_upwind(dhdt_advec, h, u, v, dx, dy, &
     dhdt_advec = 0d0
 
     do k = 1, layers
-      do j = ylow, yhigh
-        do i = xlow, xhigh
+      do j = ylow-OL+1, yhigh+OL-1
+        do i = xlow-OL+1, xhigh+OL-1
           dhdt_advec(i,j,k) = &
               ! d(hu)/dx
               ! u(i+1,j,k) point

src/aronnax.f90

@@ -97,7 +97,7 @@ program aronnax
   kv = 0d0
 
   ! size of halo region
-  OL = 1
+  OL = 3
 
   open(unit=8, file="parameters.in", status='OLD', recl=80)
   read(unit=8, nml=NUMERICS)

src/barotropic_mode.f90

@@ -30,9 +30,9 @@ subroutine calc_baro_u(ub, u, h, eta, freesurfFac, xlow, xhigh, ylow, yhigh, lay
     ! add free surface elevation to the upper layer
     h_temp(:, :, 1) = h(:, :, 1) + eta*freesurfFac
 
-    do i = xlow, xhigh+1
-      do j = ylow, yhigh
-        do k = 1, layers
+    do k = 1, layers
+      do j = ylow-OL+1, yhigh+OL-1
+        do i = xlow-OL+1, xhigh+OL
           ub(i,j) = ub(i,j) + u(i,j,k)*(h_temp(i,j,k)+h_temp(i-1,j,k))/2d0
         end do
       end do
@@ -63,9 +63,9 @@ subroutine calc_baro_v(vb, v, h, eta, freesurfFac, xlow, xhigh, ylow, yhigh, lay
     ! add free surface elevation to the upper layer
     h_temp(:, :, 1) = h(:, :, 1) + eta*freesurfFac
 
-    do i = xlow, xhigh
-      do j = ylow, yhigh+1
-        do k = 1, layers
+    do k = 1, layers
+      do j = ylow-OL+1, yhigh+OL
+        do i = xlow-OL+1, xhigh+OL-1
           vb(i,j) = vb(i,j) + v(i,j,k)*(h_temp(i,j,k)+h_temp(i,j-1,k))/2d0
         end do
       end do
@@ -95,8 +95,8 @@ subroutine calc_eta_star(ub, vb, eta, etastar, &
 
     etastar = 0d0
 
-    do i = xlow, xhigh
-      do j = ylow, yhigh
+    do j = ylow-OL+1, yhigh+OL-1
+      do i = xlow-OL+1, xhigh+OL-1
         etastar(i,j) = freesurfFac*eta(i,j) - &
             dt*((ub(i+1,j) - ub(i,j))/dx + (vb(i,j+1) - vb(i,j))/dy)
       end do
@@ -118,7 +118,7 @@ subroutine SOR_solver(a, etanew, etastar, nx, ny, OL, dt, &
 
     double precision, intent(in)  :: a(5, nx, ny)
     double precision, intent(out) :: etanew(1-OL:nx+OL, 1-OL:ny+OL)
-    double precision, intent(in)  :: etastar(1-OL:nx+OL, 1-OL:ny+OL)
+    double precision, intent(inout)  :: etastar(1-OL:nx+OL, 1-OL:ny+OL)
     integer, intent(in) :: nx, ny, OL
     double precision, intent(in) :: dt
     double precision, intent(in) :: rjac, eps
@@ -130,6 +130,8 @@ subroutine SOR_solver(a, etanew, etastar, nx, ny, OL, dt, &
     double precision norm, norm0
     double precision relax_param
 
+    call wrap_fields_2D(etastar, nx, ny, OL)
+
     rhs = -etastar(1:nx,1:ny)/dt**2
     ! first guess for etanew
     etanew = etastar
@@ -492,9 +494,10 @@ subroutine update_velocities_for_barotropic_tendency(array, etanew, g_vec, &
 
     ! TODO Assert that xstep and ystep are either 1, 0 or 0, 1.
 
-    do i = xlow-OL+xstep, xhigh
-      do j = ylow-OL+ystep, yhigh
-        do k = 1, layers
+    do k = 1, layers
+      do j = ylow-OL+ystep, yhigh+OL-1
+        do i = xlow-OL+xstep, xhigh+OL-1
+
           baro_contrib = &
               -g_vec(1)*(etanew(i,j) - etanew(i-xstep,j-ystep))/(dspace)
           array(i,j,k) = array(i,j,k) + dt*baro_contrib
@@ -646,13 +649,14 @@ subroutine barotropic_correction(hnew, unew, vnew, eta, etanew, depth, a, &
       etanew, nx, ny, xlow, xhigh, ylow, yhigh, OL, dt, maxits, eps, ierr)
 #endif
 
+    etanew = etanew*wetmask
+
     if (debug_level .ge. 4) then
       call write_output_3d(etanew, nx, ny, 1, ilower, iupper, &
                           xlow, xhigh, ylow, yhigh, OL, 0, 0, &
                           n, 'output/snap.eta_new.', num_procs, myid)
     end if
 
-    etanew = etanew*wetmask
 
     call update_halos(etanew, nx, ny, 1, ilower, iupper, &
                           xlow, xhigh, ylow, yhigh, OL, num_procs, myid)

src/bernoulli.f90

@@ -62,17 +62,15 @@ subroutine evaluate_b_iso(b, h, u, v, g_vec, depth, &
 
     ! For the rest of the layers we get a baroclinic pressure contribution
     do k = 1, layers ! move through the different layers of the model
-      do j = ylow, yhigh ! move through longitude
-        do i = xlow, xhigh ! move through latitude
+      do j = ylow-OL+1, yhigh+OL-1 ! move through longitude
+        do i = xlow-OL+1, xhigh+OL-1 ! move through latitude
           b(i,j,k) = M(i,j,k) &
               + (u(i,j,k)**2+u(i+1,j,k)**2+v(i,j,k)**2+v(i,j+1,k)**2)/4.0d0
           ! Add the (u^2 + v^2)/2 term to the Montgomery Potential
         end do
       end do
     end do
 
-    call update_halos(b, nx, ny, layers, ilower, iupper, &
-                          xlow, xhigh, ylow, yhigh, OL, num_procs, myid)
 
     return
   end subroutine evaluate_b_iso
@@ -101,8 +99,8 @@ subroutine evaluate_b_RedGrav(b, h, u, v, gr, nx, ny, layers, xlow, xhigh, ylow,
     b = 0d0
 
     do k = 1, layers ! move through the different layers of the model
-      do j = ylow, yhigh ! move through longitude
-        do i = xlow, xhigh ! move through latitude
+      do j = ylow-OL+1, yhigh+OL-1 ! move through longitude
+        do i = xlow-OL+1, xhigh+OL-1 ! move through latitude
           ! The following loops are to get the pressure term in the
           ! Bernoulli Potential
           b_proto = 0d0
@@ -124,8 +122,6 @@ subroutine evaluate_b_RedGrav(b, h, u, v, gr, nx, ny, layers, xlow, xhigh, ylow,
       end do
     end do
 
-    call update_halos(b, nx, ny, layers, ilower, iupper, &
-                          xlow, xhigh, ylow, yhigh, OL, num_procs, myid)
     return
   end subroutine evaluate_b_RedGrav
 

src/boundaries.f90

@@ -176,11 +176,15 @@ subroutine wrap_fields_3D(array, nx, ny, layers, OL)
     double precision, intent(inout) :: array(1-OL:nx+OL, 1-OL:ny+OL, layers)
     integer, intent(in) :: nx, ny, layers, OL
 
-    ! wrap array around for periodicity
-    array(0, :, :) = array(nx, :, :)
-    array(nx+1, :, :) = array(1, :, :)
-    array(:, 0, :) = array(:, ny, :)
-    array(:, ny+1, :) = array(:, 1, :)
+    integer :: k
+
+    do k = 1,layers
+      ! wrap array around for periodicity
+      array(1-OL:0, :, k) = array(nx-OL+1:nx, :, k)
+      array(nx+1:nx+OL, :, k) = array(1:OL, :, k)
+      array(:, 1-OL:0, k) = array(:, ny-OL+1:ny, k)
+      array(:, ny+1:ny+OL, k) = array(:, 1:OL, k)
+    end do
 
     return
   end subroutine wrap_fields_3D
@@ -195,10 +199,10 @@ subroutine wrap_fields_2D(array, nx, ny, OL)
     integer, intent(in) :: nx, ny, OL
 
     ! wrap array around for periodicity
-    array(0, :) = array(nx, :)
-    array(nx+1, :) = array(1, :)
-    array(:, 0) = array(:, ny)
-    array(:, ny+1) = array(:, 1)
+    array(1-OL:0, :) = array(nx-OL+1:nx, :)
+    array(nx+1:nx+OL, :) = array(1:OL, :)
+    array(:, 1-OL:0) = array(:, ny-OL+1:ny)
+    array(:, ny+1:ny+OL) = array(:, 1:OL)
 
     return
   end subroutine wrap_fields_2D

src/io.f90

@@ -651,11 +651,11 @@ subroutine write_diag_output(array, nx, ny, layers, ilower, iupper, &
     if (myid .eq. 0) then
       ! prepare data for file
       do k = 1, layers
-        diag_out(1+(4*(k-1))) = sum(global_array(:,:,k))/dble(size(global_array(:,:,k))) ! mean
-        diag_out(2+(4*(k-1))) = maxval(global_array(:,:,k))
-        diag_out(3+(4*(k-1))) = minval(global_array(:,:,k))
-        diag_out(4+(4*(k-1))) = sqrt( sum( (global_array(:,:,k) - diag_out(1+(4*(k-1))))**2)/ &
-                              dble(size(global_array(:,:,k))))
+        diag_out(1+(4*(k-1))) = sum(global_array(1:nx,1:ny,k))/dble(size(global_array(1:nx,1:ny,k))) ! mean
+        diag_out(2+(4*(k-1))) = maxval(global_array(1:nx,1:ny,k))
+        diag_out(3+(4*(k-1))) = minval(global_array(1:nx,1:ny,k))
+        diag_out(4+(4*(k-1))) = sqrt( sum( (global_array(1:nx,1:ny,k) - diag_out(1+(4*(k-1))))**2)/ &
+                              dble(size(global_array(1:nx,1:ny,k))))
       end do
 
       ! Output the data to a file

src/model_main.f90

@@ -331,10 +331,7 @@ subroutine model_run(h, u, v, eta, depth, dx, dy, wetmask, &
                           xlow, xhigh, ylow, yhigh, OL, num_procs, myid)
       call update_halos(v_new, nx, ny, layers, ilower, iupper, &
                           xlow, xhigh, ylow, yhigh, OL, num_procs, myid)
-      if (.not. RedGrav) then
-        call update_halos(etanew, nx, ny, 1, ilower, iupper, &
-                          xlow, xhigh, ylow, yhigh, OL, num_procs, myid)
-      end if   
+      ! eta_new halo is updated in barotropic_correction
 
       ! Accumulate average fields
       if (avwrite .ne. 0) then