Add BC support in mesh class.

src/common.f90

@@ -19,6 +19,7 @@ module m_common
                         Y_EDGE = 102, & ! Data on edges along Y
                         Z_EDGE = 103, & ! Data on edges along Z
                         none = -1 ! The location of data isn't specified
+  integer, parameter :: BC_PERIODIC = 0, BC_NEUMANN = 1, BC_DIRICHLET = 2
   integer, protected :: &
     rdr_map(4, 4) = reshape([0, RDR_Y2X, RDR_Z2X, RDR_C2X, &
                              RDR_X2Y, 0, RDR_Z2Y, RDR_C2Y, &

src/mesh.f90

@@ -4,7 +4,8 @@ module m_mesh
   use mpi
   use m_common, only: dp, DIR_X, DIR_Y, DIR_Z, DIR_C, &
                       CELL, VERT, none, X_FACE, Y_FACE, Z_FACE, &
-                      X_EDGE, Y_EDGE, Z_EDGE
+                      X_EDGE, Y_EDGE, Z_EDGE, &
+                      BC_PERIODIC, BC_NEUMANN, BC_DIRICHLET
   use m_field, only: field_t
 
   implicit none
@@ -38,6 +39,7 @@ module m_mesh
     integer, dimension(3), private :: vert_dims ! local number of vertices in each direction without padding (cartesian structure)
     integer, dimension(3), private :: cell_dims ! local number of cells in each direction without padding (cartesian structure)
     logical, dimension(3), private :: periodic_BC ! Whether or not a direction has a periodic BC
+    integer, dimension(3, 2), private :: BCs_global
     integer, private :: sz
     type(geo_t), allocatable :: geo ! object containing geometry information
     type(parallel_t), allocatable :: par ! object containing parallel domain decomposition information
@@ -77,30 +79,46 @@ module m_mesh
 
 contains
 
-  function mesh_init(dims_global, nproc_dir, L_global, &
-                     periodic_BC) result(mesh)
+  function mesh_init(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z) &
+    result(mesh)
     !! Completely initialise the mesh object.
     !! Upon initialisation the mesh object can be read-only and shouldn't be edited
     !! Takes as argument global information about the mesh like its length, number of cells and decomposition in each direction
     integer, dimension(3), intent(in) :: dims_global
     integer, dimension(3), intent(in) :: nproc_dir ! Number of proc in each direction
     real(dp), dimension(3), intent(in) :: L_global
-    logical, dimension(3), optional, intent(in) :: periodic_BC
+    character(len=*), dimension(2), intent(in) :: BC_x, BC_y, BC_z
     type(mesh_t) :: mesh
 
-    integer :: dir
+    character(len=20), dimension(3, 2) :: BC_all
+    integer :: dir, j
     integer :: ierr
 
     allocate (mesh%geo)
     allocate (mesh%par)
     mesh%global_vert_dims(:) = dims_global
 
-    if (present(periodic_BC)) then
-      mesh%periodic_BC(:) = periodic_BC
-    else
-      ! Default to periodic BC
-      mesh%periodic_BC(:) = .true.
-    end if
+    BC_all(1, 1) = BC_x(1); BC_all(1, 2) = BC_x(2)
+    BC_all(2, 1) = BC_y(1); BC_all(2, 2) = BC_y(2)
+    BC_all(3, 1) = BC_z(1); BC_all(3, 2) = BC_z(2)
+    do dir = 1, 3
+      do j = 1, 2
+        select case (trim(BC_all(dir, j)))
+        case ('periodic')
+          mesh%BCs_global(dir, j) = BC_PERIODIC
+        case ('neumann')
+          mesh%BCs_global(dir, j) = BC_NEUMANN
+        case ('dirichlet')
+          mesh%BCs_global(dir, j) = BC_DIRICHLET
+        case default
+          error stop 'Unknown BC'
+        end select
+      end do
+    end do
+
+    do dir = 1, 3
+      mesh%periodic_BC(dir) = all(mesh%BCs_global(dir, :) == BC_PERIODIC)
+    end do
 
     do dir = 1, 3
       if (mesh%periodic_BC(dir)) then

src/xcompact.f90

@@ -42,12 +42,13 @@ program xcompact
 
   character(len=200) :: input_file
   character(len=3) :: method
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   integer, dimension(3) :: dims_global
   integer, dimension(3) :: nproc_dir = 0
   real(dp), dimension(3) :: L_global
   integer :: nrank, nproc, ierr
 
-  namelist /domain_params/ L_global, dims_global, nproc_dir
+  namelist /domain_params/ L_global, dims_global, nproc_dir, BC_x, BC_y, BC_z
   namelist /time_intg_params/ method
 
   call MPI_Init(ierr)
@@ -63,9 +64,12 @@ program xcompact
 #endif
 
   ! set defaults
-  L_global = [2*pi, 2*pi, 2*pi]
   dims_global = [256, 256, 256]
+  L_global = [2*pi, 2*pi, 2*pi]
   nproc_dir = [1, 1, nproc]
+  BC_x = [character(len=20) :: 'periodic', 'periodic']
+  BC_y = [character(len=20) :: 'periodic', 'periodic']
+  BC_z = [character(len=20) :: 'periodic', 'periodic']
   method = 'AB3'
 
   if (command_argument_count() >= 1) then
@@ -83,7 +87,7 @@ program xcompact
     nproc_dir = [1, 1, nproc]
   end if
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global)
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
 #ifdef CUDA
   cuda_allocator = cuda_allocator_t(mesh, SZ)

tests/cuda/test_cuda_allocator.f90

@@ -12,6 +12,7 @@ program test_allocator_cuda
   logical :: allpass
   integer, dimension(3) :: dims, nproc_dir
   real(dp) :: L_global(3)
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   class(allocator_t), allocatable :: allocator
   class(mesh_t), allocatable :: mesh
   class(field_t), pointer :: ptr1, ptr2, ptr3
@@ -24,7 +25,11 @@ program test_allocator_cuda
   nproc_dir = [1, 1, 1]
   L_global = [2*pi, 2*pi, 2*pi]
 
-  mesh = mesh_t(dims, nproc_dir, L_global)
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
+  mesh = mesh_t(dims, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
   allocator = cuda_allocator_t(mesh, 8)
 

tests/omp/test_omp_transeq.f90

@@ -20,6 +20,7 @@ program test_omp_transeq
   class(mesh_t), allocatable :: mesh
   integer, dimension(3) :: dims_global, nproc_dir
   real(dp), dimension(3) :: L_global
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
 
   integer :: n, n_groups, i, j, k
   integer :: nrank, nproc
@@ -48,7 +49,11 @@ program test_omp_transeq
   ! Domain decomposition in each direction
   nproc_dir = [nproc, 1, 1]
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global)
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
   xdirps%dir = DIR_X; ydirps%dir = DIR_Y; zdirps%dir = DIR_Z
 

tests/test_allocator.f90

@@ -12,6 +12,7 @@ program test_allocator
   logical :: allpass
   integer, dimension(3) :: nproc_dir, dims_global
   real(dp), dimension(3) :: L_global
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   class(allocator_t), allocatable :: allocator
   class(mesh_t), allocatable :: mesh
   class(field_t), pointer :: ptr1, ptr2, ptr3
@@ -28,7 +29,11 @@ program test_allocator
   ! Domain decomposition in each direction
   nproc_dir = [1, 1, 1]
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global)
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
   allocator = allocator_t(mesh, 8)
 

tests/test_mesh.f90

@@ -11,8 +11,8 @@ program test_allocator
 
   logical :: allpass
   integer, dimension(3) :: nproc_dir, dims_global
-  logical, dimension(3) :: periodic_BC
   real(dp), dimension(3) :: L_global
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   class(allocator_t), allocatable :: allocator
   class(mesh_t), allocatable :: mesh
   class(field_t), pointer :: ptr1, ptr2, ptr3
@@ -34,9 +34,11 @@ program test_allocator
   ! Domain decomposition in each direction
   nproc_dir = [4, 1, 1]
 
-  periodic_BC = [.false., .true., .false.]
+  BC_x = [character(len=20) :: 'dirichlet', 'neumann']
+  BC_y = [character(len=20) :: 'periodic', 'periodic']
+  BC_z = [character(len=20) :: 'neumann', 'dirichlet']
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global, periodic_BC)
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
   if (mesh%par%nproc /= 4) then
     allpass = .false.

tests/test_reordering.f90

@@ -45,6 +45,7 @@ program test_reorder
   type(dirps_t), target :: xdirps, ydirps, zdirps
   integer, dimension(3) :: dims_padded, dims_global, nproc_dir
   real(dp), dimension(3) :: L_global
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   logical :: pass_X, pass_Y, pass_Z
 
 #ifdef CUDA
@@ -76,7 +77,11 @@ program test_reorder
   ! Domain decomposition in each direction
   nproc_dir = [nproc, 1, 1]
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global)
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
 #ifdef CUDA
   cuda_allocator = cuda_allocator_t(mesh, SZ)

tests/test_scalar_product.f90

@@ -29,6 +29,8 @@ program test_scalar_product
 
   type(mesh_t) :: mesh
 
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
+
   character(len=5), dimension(4), parameter :: test = &
     ["DIR_X", "DIR_Y", "DIR_Z", "DIR_C"]
   integer, dimension(4), parameter :: dir = [DIR_X, DIR_Y, DIR_Z, DIR_C]
@@ -43,9 +45,14 @@ program test_scalar_product
   call MPI_Comm_rank(MPI_COMM_WORLD, nrank, ierr)
   call MPI_Comm_size(MPI_COMM_WORLD, nproc, ierr)
 
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
   mesh = mesh_t([nx, ny, nz], &
                 [1, 1, nproc], &
-                [lx, ly, lz])
+                [lx, ly, lz], &
+                BC_x, BC_y, BC_z)
 
 #ifdef CUDA
 #else

tests/test_sum_intox.f90

@@ -27,6 +27,7 @@ program test_sum_intox
 
   type(mesh_t) :: mesh
 
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
   integer :: nrank, nproc
   integer :: ierr
 
@@ -36,9 +37,14 @@ program test_sum_intox
   call MPI_Comm_rank(MPI_COMM_WORLD, nrank, ierr)
   call MPI_Comm_size(MPI_COMM_WORLD, nproc, ierr)
 
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
   mesh = mesh_t([nx, ny, nz], &
                 [1, 1, nproc], &
-                [lx, ly, lz])
+                [lx, ly, lz], &
+                BC_x, BC_y, BC_z)
 
 #ifdef CUDA
 #else

tests/test_time_integrator.f90

@@ -30,6 +30,7 @@ program test_omp_adamsbashforth
   class(field_t), pointer :: du, dv, dw
   integer, dimension(3) :: dims_global, nproc_dir
   real(dp), dimension(3) :: L_global
+  character(len=20) :: BC_x(2), BC_y(2), BC_z(2)
 
   class(base_backend_t), pointer :: backend
   class(allocator_t), pointer :: allocator
@@ -64,7 +65,11 @@ program test_omp_adamsbashforth
   ! Domain decomposition in each direction
   nproc_dir = [1, 1, 1]
 
-  mesh = mesh_t(dims_global, nproc_dir, L_global)
+  BC_x = ['periodic', 'periodic']
+  BC_y = ['periodic', 'periodic']
+  BC_z = ['periodic', 'periodic']
+
+  mesh = mesh_t(dims_global, nproc_dir, L_global, BC_x, BC_y, BC_z)
 
   ! allocate object
 #ifdef CUDA