feat(subdomain): make fully asynchronous

example/time-paradigm.f90

@@ -12,7 +12,7 @@ program time_paradigm_m
   character(len=:), allocatable :: steps_string, resolution_string
   type(command_line_t) command_line
   integer(int64) counter_start, counter_end, clock_rate
-  integer :: steps=200, resolution=64
+  integer :: steps=300, resolution=128
 
   associate(me => this_image())
     if (command_line%argument_present(["--help"])) then
@@ -56,16 +56,18 @@ function functional_programming_time() result(system_time)
 
     call T%define(side=1., boundary_val=T_boundary, internal_val=T_internal_initial, n=resolution)
 
-    call system_clock(t_start_functional)
-
     associate(dt => T%dx()*T%dy()/(4*alpha))
+      call system_clock(t_start_functional)
+
       functional_programming: &
       do step = 1, steps
+        sync all
         T =  T + dt * alpha * .laplacian. T
       end do functional_programming
+
+      call system_clock(t_end_functional, clock_rate)
     end associate
 
-    call system_clock(t_end_functional, clock_rate)
     system_time = real(t_end_functional - t_start_functional)/real(clock_rate)
 
     associate(L_infinity_norm => maxval(abs(T%values() - T_steady)))
@@ -80,16 +82,20 @@ function procedural_programming_time() result(system_time)
     real system_time
     type(subdomain_t) T
 
+    call T%define(side=1., boundary_val=0., internal_val=1., n=resolution)
+
     associate(dt => T%dx()*T%dy()/(4*alpha))
-      call T%define(side=1., boundary_val=0., internal_val=1., n=resolution)
       call system_clock(t_start_procedural)
+
       procedural_programming: &
       do step = 1, steps
+        sync all
         call T%step(alpha*dt)
       end do procedural_programming
+
+      call system_clock(t_end_procedural, clock_rate)
     end associate
 
-    call system_clock(t_end_procedural, clock_rate)
     system_time = real(t_end_procedural - t_start_procedural)/real(clock_rate)
 
     associate(L_infinity_norm => maxval(abs(T%values() - T_steady)))

src/matcha/subdomain_m.f90

@@ -14,14 +14,14 @@ module subdomain_m
     generic :: operator(.laplacian.) => laplacian
     generic :: operator(*) => multiply
     generic :: operator(+) => add
-    generic :: assignment(=) => assign_and_sync
+    generic :: assignment(=) => assign_
     procedure dx
     procedure dy
     procedure dz
     procedure values
     procedure, private :: laplacian
     procedure, private :: add
-    procedure, private :: assign_and_sync
+    procedure, private :: assign_
   end type
 
   interface
@@ -83,7 +83,7 @@ pure module function add(lhs, rhs) result(total)
       type(subdomain_t) total
     end function
 
-    module subroutine assign_and_sync(lhs, rhs)
+    module subroutine assign_(lhs, rhs)
       implicit none
       class(subdomain_t), intent(out) :: lhs
       type(subdomain_t), intent(in) :: rhs

src/matcha/subdomain_s.f90

@@ -55,7 +55,6 @@
     allocate(halo_x(west:east, ny, nz)[*])
     if (me>1) halo_x(east,:,:)[me-1] = self%s_(1,:,:)
     if (me<num_subdomains) halo_x(west,:,:)[me+1] = self%s_(my_nx,:,:)
-    sync all
   end procedure
 
   module procedure dx
@@ -123,20 +122,23 @@
     total%s_ =  lhs%s_ + rhs%s_
   end procedure
 
-  module procedure assign_and_sync
-    call assert(allocated(rhs%s_), "subdomain_t%assign_and_sync: allocated(rhs%s_)")
-    sync all
+  module procedure assign_
+    call assert(allocated(rhs%s_), "subdomain_t%assign_: allocated(rhs%s_)")
     lhs%s_ =  rhs%s_
-    if (me>1) halo_x(east,:,:)[me-1] = rhs%s_(1,:,:)
-    if (me<num_subdomains) halo_x(west,:,:)[me+1] = rhs%s_(my_nx,:,:)
-    sync all
+    call exchange_halo(rhs%s_)
   end procedure
 
   module procedure values
     call assert(allocated(self%s_), "subdomain_t%values: allocated(self%s_)")
     my_values =  self%s_
   end procedure
 
+  subroutine exchange_halo(s)
+    real, intent(in) :: s(:,:,:)
+    if (me>1) halo_x(east,:,:)[me-1] = s(1,:,:)
+    if (me<num_subdomains) halo_x(west,:,:)[me+1] = s(my_nx,:,:)
+  end subroutine
+
   module procedure step
 
     call assert(allocated(self%s_), "subdomain_t%laplacian: allocated(rhs%s_)")
@@ -149,10 +151,7 @@
     call internal_points(increment)
     call edge_points(increment)
     call apply_boundary_condition(increment)
-
-    sync all
     self%s_ = self%s_ + increment
-    sync all
     call exchange_halo(self%s_)
 
   contains
@@ -207,11 +206,6 @@ subroutine apply_boundary_condition(ds)
       if (me==num_subdomains) ds(my_nx,:,:) = 0.
     end subroutine
 
-    subroutine exchange_halo(s)
-      real, intent(in) :: s(:,:,:)
-      if (me>1) halo_x(east,:,:)[me-1] = s(1,:,:)
-      if (me<num_subdomains) halo_x(west,:,:)[me+1] = s(my_nx,:,:)
-    end subroutine
 
   end procedure
 

test/subdomain_test_m.f90

@@ -53,7 +53,8 @@ subroutine output(v)
     critical
       do j = 1, size(v,2)
         do k = 1, size(v,3)
-          print *,"image ",this_image(),": ",j,k,v(:,j,k)
+          !print *,"image ",this_image(),": ",j,k,v(:,j,k)
+          print *,j,k,v(:,j,k)
         end do
       end do
     end critical
@@ -65,7 +66,8 @@ function concave_laplacian() result(test_passes)
     type(subdomain_t) f, laplacian_f
     real, allocatable :: lap_f_vals(:,:,:)
 
-    call f%define(side=1., boundary_val=1., internal_val=2., n=21) ! internally constant subdomain with a step down at all surfaces
+    call f%define(side=1., boundary_val=1., internal_val=2., n=32) ! internally constant subdomain with a step down at all surfaces
+    sync all
     laplacian_f = .laplacian. f
     lap_f_vals = laplacian_f%values()
 
@@ -154,14 +156,15 @@ function concave_laplacian() result(test_passes)
   function correct_steady_state() result(test_passes)
     logical test_passes
     type(subdomain_t) T
-    real, parameter :: T_boundary = 1., T_initial = 2., tolerance = 0.01, T_steady = T_boundary, alpha = 1.
-    integer, parameter :: steps = 6000
+    real, parameter :: T_boundary = 1., T_initial = 2., tolerance = 5.E-03, T_steady = T_boundary, alpha = 1.
+    integer, parameter :: steps = 25000
     integer step
 
-    call T%define(side=1., boundary_val=T_boundary, internal_val=T_initial, n=21) ! const. internally with a step down at boundaries
+    call T%define(side=1., boundary_val=T_boundary, internal_val=T_initial, n=32) ! const. internally with a step down at boundaries
 
     associate(dt => T%dx()*T%dy()*T%dz()/(4*alpha))
       do step = 1, steps
+        sync all
         T =  T + dt * alpha * .laplacian. T
       end do
     end associate
@@ -174,7 +177,7 @@ function correct_steady_state() result(test_passes)
   function functional_matches_procedural() result(test_passes)
     logical test_passes
     real, parameter :: tolerance = 1.E-06
-    integer, parameter :: steps = 6000, n=21
+    integer, parameter :: steps = 6000, n=32
     real, parameter :: alpha = 1.
     real, parameter :: side=1., boundary_val=1., internal_val=2.
 
@@ -195,6 +198,7 @@ function T_functional()
 
       associate(dt => T%dx()*T%dy()/(4*alpha))
         do step = 1, steps
+          sync all
           T =  T + dt * alpha * .laplacian. T
         end do
       end associate
@@ -211,6 +215,7 @@ function T_procedural()
 
       associate(dt => T%dx()*T%dy()/(4*alpha))
         do step = 1, steps
+          sync all
           call T%step(alpha*dt)
         end do
       end associate