Improve allocator to support arbitary domain sizes.

src/allocator.f90

@@ -31,10 +31,14 @@ module m_allocator
      !! [[m_allocator(module):release_block(subroutine)]].  The
      !! released block is then pushed in front of the block list.
 
-      integer :: dims(3)
+      integer :: nx_padded, ny_padded, nz_padded, sz
       !> The id for the next allocated block.  This counter is
       !> incremented each time a new block is allocated.
       integer :: next_id = 0
+      !> Padded dimensions for x, y, and z oriented fields
+      integer :: xdims(3), ydims(3), zdims(3)
+      !> Padded dimensions for natural Cartesian ordering
+      integer :: cdims(3)
       !> The pointer to the first block on the list.  Non associated if
       !> the list is empty
       ! TODO: Rename first to head
@@ -47,20 +51,25 @@ module m_allocator
       procedure :: destroy
    end type allocator_t
 
+   interface allocator_t
+      module procedure allocator_init
+   end interface allocator_t
+
    type :: field_t
      !! Memory block type holding both a data field and a pointer
      !! to the next block.  The `field_t` type also holds a integer
      !! `refcount` that counts the number of references to this
      !! field.  User code is currently responsible for incrementing
      !! the reference count.
       class(field_t), pointer :: next
-      real(dp), allocatable :: data(:)
+      real(dp), pointer, private :: p_data(:)
+      real(dp), pointer, contiguous :: data(:, :, :)
       integer :: refcount = 0
       integer :: id !! An integer identifying the memory block.
    end type field_t
 
    interface field_t
-      module procedure field_constructor
+      module procedure field_init
    end interface field_t
 
    type :: flist_t
@@ -69,16 +78,40 @@ module m_allocator
 
 contains
 
-   function field_constructor(dims, next, id) result(m)
-      integer, intent(in) :: dims(3), id
+   function field_init(nx, ny, nz, sz, next, id) result(f)
+      integer, intent(in) :: nx, ny, nz, sz, id
       type(field_t), pointer, intent(in) :: next
-      type(field_t) :: m
-
-      allocate (m%data(dims(1)*dims(2)*dims(3)))
-      m%refcount = 0
-      m%next => next
-      m%id = id
-   end function field_constructor
+      type(field_t) :: f
+
+      allocate (f%p_data(nx*ny*nz))
+      ! will be removed, bounds remapping will be carried out by get_block.
+      f%data(1:sz, 1:nx, 1:ny*nz/sz) => f%p_data
+      f%refcount = 0
+      f%next => next
+      f%id = id
+   end function field_init
+
+   function allocator_init(nx, ny, nz, sz) result(allocator)
+      integer, intent(in) :: nx, ny, nz, sz
+      type(allocator_t) :: allocator
+
+      integer :: nx_padded, ny_padded, nz_padded
+
+      ! Apply padding based on sz
+      nx_padded = nx
+      ny_padded = ny
+      nz_padded = nz
+
+      allocator%nx_padded = nx_padded
+      allocator%ny_padded = ny_padded
+      allocator%nz_padded = nz_padded
+      allocator%sz = sz
+
+      allocator%xdims = [sz, nx_padded, ny_padded*nz_padded/sz]
+      allocator%ydims = [sz, ny_padded, nx_padded*nz_padded/sz]
+      allocator%zdims = [sz, nz_padded, nx_padded*ny_padded/sz]
+      allocator%cdims = [nx_padded, ny_padded, nz_padded]
+   end function allocator_init
 
    function create_block(self, next) result(ptr)
     !! Allocate memory for a new block and return a pointer to a new
@@ -89,7 +122,8 @@ function create_block(self, next) result(ptr)
       class(field_t), pointer :: ptr
       self%next_id = self%next_id + 1
       allocate (newblock)
-      newblock = field_t(self%dims, next, id=self%next_id)
+      newblock = field_t(self%nx_padded, self%ny_padded, self%nz_padded, &
+                         self%sz, next, id=self%next_id)
       ptr => newblock
    end function create_block
 
@@ -115,6 +149,12 @@ function get_block(self) result(handle)
       handle => self%first
       self%first => self%first%next ! 2nd block becomes head block
       handle%next => null() ! Detach ex-head block from the block list
+
+      ! Bounds remapping will be carried out by a dedicated function
+      ! here based on the optional ordering get_block is passed
+      ! something like:
+      !handle%data(1:self%xdims(1), 1:self%xdims(2), 1:self%xdims(3)) &
+      !   => handle%p_data
    end function get_block
 
    subroutine release_block(self, handle)

src/cuda/allocator.f90

@@ -9,26 +9,40 @@ module m_cuda_allocator
       procedure :: create_block => create_cuda_block
    end type cuda_allocator_t
 
+   interface cuda_allocator_t
+      module procedure cuda_allocator_init
+   end interface cuda_allocator_t
+
    type, extends(field_t) :: cuda_field_t
-      real(dp), allocatable, device :: data_d(:)
+      real(dp), device, pointer, private :: p_data_d(:)
+      real(dp), device, pointer, contiguous :: data_d(:, :, :)
    end type cuda_field_t
 
    interface cuda_field_t
-      module procedure cuda_field_constructor
+      module procedure cuda_field_init
    end interface cuda_field_t
 
 contains
 
-   function cuda_field_constructor(dims, next, id) result(m)
-      integer, intent(in) :: dims(3), id
+   function cuda_field_init(nx, ny, nz, sz, next, id) result(f)
+      integer, intent(in) :: nx, ny, nz, sz, id
       type(cuda_field_t), pointer, intent(in) :: next
-      type(cuda_field_t) :: m
+      type(cuda_field_t) :: f
+
+      allocate (f%p_data_d(nx*ny*nz))
+      ! will be removed, bounds remapping will be carried out by get_block.
+      f%data_d(1:sz, 1:nx, 1:ny*nz/sz) => f%p_data_d
+      f%refcount = 0
+      f%next => next
+      f%id = id
+   end function cuda_field_init
+
+   function cuda_allocator_init(nx, ny, nz, sz) result(allocator)
+      integer, intent(in) :: nx, ny, nz, sz
+      type(cuda_allocator_t) :: allocator
 
-      allocate (m%data_d(dims(1)*dims(2)*dims(3)))
-      m%refcount = 0
-      m%next => next
-      m%id = id
-   end function cuda_field_constructor
+      allocator%allocator_t = allocator_t(nx, ny, nz, sz)
+   end function cuda_allocator_init
 
    function create_cuda_block(self, next) result(ptr)
       class(cuda_allocator_t), intent(inout) :: self
@@ -37,7 +51,8 @@ function create_cuda_block(self, next) result(ptr)
       class(field_t), pointer :: ptr
       allocate (newblock)
       self%next_id = self%next_id + 1
-      newblock = cuda_field_t(self%dims, next, id=self%next_id)
+      newblock = cuda_field_t(self%nx_padded, self%ny_padded, self%nz_padded, &
+                              self%sz, next, id=self%next_id)
       ptr => newblock
    end function create_cuda_block
 

src/solver.f90

@@ -101,7 +101,7 @@ function init(backend, time_integrator, xdirps, ydirps, zdirps, globs) &
       solver%w => solver%backend%allocator%get_block()
 
       ! Set initial conditions
-      dims(:) = solver%backend%allocator%dims(:)
+      dims(:) = solver%backend%allocator%xdims(:)
       allocate(u_init(dims(1), dims(2), dims(3)))
       allocate(v_init(dims(1), dims(2), dims(3)))
       allocate(w_init(dims(1), dims(2), dims(3)))

src/xcompact.f90

@@ -111,15 +111,16 @@ program xcompact
    zdirps%n_blocks = globs%n_groups_z
 
 #ifdef CUDA
-   cuda_allocator = cuda_allocator_t([SZ, globs%nx_loc, globs%n_groups_x])
+   cuda_allocator = cuda_allocator_t(globs%nx_loc, globs%ny_loc, &
+                                     globs%nz_loc, SZ)
    allocator => cuda_allocator
    print*, 'CUDA allocator instantiated'
 
    cuda_backend = cuda_backend_t(globs, allocator)
    backend => cuda_backend
    print*, 'CUDA backend instantiated'
 #else
-   omp_allocator = allocator_t([SZ, globs%nx_loc, globs%n_groups_x])
+   omp_allocator = allocator_t(globs%nx_loc, globs%ny_loc, globs%nz_loc, SZ)
    allocator => omp_allocator
    print*, 'OpenMP allocator instantiated'
 

tests/cuda/test_cuda_allocator.f90

@@ -12,7 +12,7 @@ program test_allocator_cuda
   class(field_t), pointer :: ptr1, ptr2, ptr3
   integer, allocatable :: l(:)
 
-  allocator = cuda_allocator_t(dims)
+  allocator = cuda_allocator_t(dims(1), dims(2), dims(3), 8)
 
   allpass = .true.
 

tests/omp/test_omp_transeq.f90

@@ -68,7 +68,7 @@ program test_omp_transeq
    ydirps%n_blocks = globs%n_groups_y
    zdirps%n_blocks = globs%n_groups_z
 
-   omp_allocator = allocator_t([SZ, globs%nx_loc, globs%n_groups_x])
+   omp_allocator = allocator_t(xdirps%n, ydirps%n, zdirps%n, SZ)
    allocator => omp_allocator
    print*, 'OpenMP allocator instantiated'
 

tests/test_allocator.f90

@@ -11,7 +11,7 @@ program test_allocator
   class(field_t), pointer :: ptr1, ptr2, ptr3
   integer, allocatable :: l(:)
 
-  allocator = allocator_t(dims)
+  allocator = allocator_t(dims(1), dims(2), dims(3), 8)
 
   allpass = .true.