diff --git a/CMakeLists.txt b/CMakeLists.txt index 0a1510e..802bfd1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,7 @@ project( HOMEPAGE_URL "https://github.com/smillerc/cato" LANGUAGES Fortran C) -set(CMAKE_VERBOSE_MAKEFILE OFF) + set(CMAKE_VERBOSE_MAKEFILE OFF) set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) diff --git a/include/version.h.in b/include/version.h.in index ab956db..942d107 100755 --- a/include/version.h.in +++ b/include/version.h.in @@ -2,7 +2,7 @@ GIT_HASH = "@GIT_SHA1@" GIT_REF = "@GIT_REFSPEC@" GIT_LOCAL_CHANGES = "@GIT_LOCAL_CHANGES@" -CATO_VERSION = "@CATO_VERSION_MAJOR@.@CATO_VERSION_MINOR@.@CATO_PATCH_VERSION@" +CATO_VERSION = "@CATO_VERSION_MAJOR@.@CATO_VERSION_MINOR@.@CATO_VERSION_PATCH@" COMPILE_HOST = "@HOST_NAME@" ! COMPILE_OS = "@CMAKE_HOST_SYSTEM@" BUILD_TYPE = "@CMAKE_BUILD_TYPE@" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d7dae9b..0360d8e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -138,9 +138,9 @@ add_library( lib/io/mod_hdf5.f90 lib/io/mod_input.f90 lib/grid/grid_block.f90 - lib/grid/grid_block_1d.f90 + # lib/grid/grid_block_1d.f90 lib/grid/grid_block_2d.f90 - lib/grid/grid_block_3d.f90 + # lib/grid/grid_block_3d.f90 # lib/grid/grid_type.f90 # lib/grid/regular_2d_grid.f90 lib/grid/grid_factory.f90 diff --git a/src/lib/boundary_conditions/bc_factory.f90 b/src/lib/boundary_conditions/bc_factory.f90 index a47a800..28c47d3 100644 --- a/src/lib/boundary_conditions/bc_factory.f90 +++ b/src/lib/boundary_conditions/bc_factory.f90 @@ -21,7 +21,7 @@ module mod_bc_factory use iso_fortran_env, only: ik => int32 use mod_input, only: input_t - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_boundary_conditions, only: boundary_condition_t use mod_periodic_bc, only: periodic_bc_t, periodic_bc_constructor use mod_symmetry_bc, only: symmetry_bc_t, symmetry_bc_constructor @@ -40,7 +40,7 @@ function bc_factory(bc_type, location, input, grid) result(bc) character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y) class(boundary_condition_t), pointer :: bc class(input_t), intent(in) :: input - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid select case(trim(bc_type)) case('periodic') diff --git a/src/lib/boundary_conditions/boundary_conditions.f90 b/src/lib/boundary_conditions/boundary_conditions.f90 index bb1cd41..5e6b270 100644 --- a/src/lib/boundary_conditions/boundary_conditions.f90 +++ b/src/lib/boundary_conditions/boundary_conditions.f90 @@ -25,7 +25,7 @@ module mod_boundary_conditions use mod_functional, only: operator(.reverse.) use mod_field, only: field_2d_t use mod_input, only: input_t - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t implicit none @@ -85,7 +85,7 @@ subroutine set_indices(self, grid) !< by the other procedures that apply the boundary conditions class(boundary_condition_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik) :: i ! integer(ik), dimension(:, :), intent(in) :: ghost_layers !< (ilo_layers(n), ihi_layers(n), jlo_layers(n), jhi_layers(n)); indices to the ghost layers. @@ -106,31 +106,31 @@ subroutine set_indices(self, grid) allocate(self%jhi_ghost(self%n_ghost_layers)) if(self%n_ghost_layers == 1) then - self%ilo_ghost = grid%cell_lbounds_halo(1) - self%ihi_ghost = grid%cell_ubounds_halo(1) - self%jlo_ghost = grid%cell_lbounds_halo(2) - self%jhi_ghost = grid%cell_ubounds_halo(2) - self%ilo = grid%cell_lbounds(1) - self%ihi = grid%cell_ubounds(1) - self%jlo = grid%cell_lbounds(2) - self%jhi = grid%cell_ubounds(2) + self%ilo_ghost = grid%lbounds_halo(1) + self%ihi_ghost = grid%ubounds_halo(1) + self%jlo_ghost = grid%lbounds_halo(2) + self%jhi_ghost = grid%ubounds_halo(2) + self%ilo = grid%lbounds(1) + self%ihi = grid%ubounds(1) + self%jlo = grid%lbounds(2) + self%jhi = grid%ubounds(2) else if(self%n_ghost_layers > 1) then do i = 1, self%n_ghost_layers - self%ilo_ghost(i) = grid%cell_lbounds_halo(1) + i - 1 - self%ihi_ghost(i) = grid%cell_ubounds_halo(1) - i + 1 + self%ilo_ghost(i) = grid%lbounds_halo(1) + i - 1 + self%ihi_ghost(i) = grid%ubounds_halo(1) - i + 1 - self%jlo_ghost(i) = grid%cell_lbounds_halo(2) + i - 1 - self%jhi_ghost(i) = grid%cell_ubounds_halo(2) - i + 1 + self%jlo_ghost(i) = grid%lbounds_halo(2) + i - 1 + self%jhi_ghost(i) = grid%ubounds_halo(2) - i + 1 end do self%ilo_ghost = .reverse.self%ilo_ghost self%jlo_ghost = .reverse.self%jlo_ghost self%ihi_ghost = .reverse.self%ihi_ghost self%jhi_ghost = .reverse.self%jhi_ghost - self%ilo = grid%cell_lbounds(1) - self%ihi = grid%cell_ubounds(1) - self%jlo = grid%cell_lbounds(2) - self%jhi = grid%cell_ubounds(2) + self%ilo = grid%lbounds(1) + self%ihi = grid%ubounds(1) + self%jlo = grid%lbounds(2) + self%jhi = grid%ubounds(2) else error stop "Error in boundary_condition_t%set_indicies(), n_ghost_layers <= 0" end if diff --git a/src/lib/boundary_conditions/periodic_bc.fypp b/src/lib/boundary_conditions/periodic_bc.fypp index 63addca..3a25fc2 100644 --- a/src/lib/boundary_conditions/periodic_bc.fypp +++ b/src/lib/boundary_conditions/periodic_bc.fypp @@ -24,7 +24,7 @@ module mod_periodic_bc use mod_functional, only: reverse, operator(.reverse.) use mod_field, only: field_2d_t use mod_globals, only: enable_debug_print, debug_print - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_boundary_conditions, only: boundary_condition_t use mod_input, only: input_t @@ -46,7 +46,7 @@ contains type(periodic_bc_t), pointer :: bc character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y) class(input_t), intent(in) :: input - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid allocate(bc) bc%name = 'periodic' diff --git a/src/lib/boundary_conditions/pressure_input_bc.f90 b/src/lib/boundary_conditions/pressure_input_bc.f90 index 4bdaa88..4b1ac91 100644 --- a/src/lib/boundary_conditions/pressure_input_bc.f90 +++ b/src/lib/boundary_conditions/pressure_input_bc.f90 @@ -22,7 +22,7 @@ module mod_pressure_input_bc use, intrinsic :: iso_fortran_env, only: ik => int32, rk => real64, std_err => error_unit use mod_globals, only: enable_debug_print, debug_print use mod_field, only: field_2d_t - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_boundary_conditions, only: boundary_condition_t use mod_nondimensionalization, only: p_0, rho_0, t_0 use mod_input, only: input_t @@ -62,7 +62,7 @@ function pressure_input_bc_constructor(location, input, grid) result(bc) type(pressure_input_bc_t), pointer :: bc character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y) class(input_t), intent(in) :: input - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid allocate(bc) bc%name = 'pressure_input' diff --git a/src/lib/boundary_conditions/symmetry_bc.f90 b/src/lib/boundary_conditions/symmetry_bc.f90 index 9ba2acc..6ecf018 100644 --- a/src/lib/boundary_conditions/symmetry_bc.f90 +++ b/src/lib/boundary_conditions/symmetry_bc.f90 @@ -23,7 +23,7 @@ module mod_symmetry_bc use mod_globals, only: enable_debug_print, debug_print use mod_field, only: field_2d_t use mod_error, only: error_msg - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_boundary_conditions, only: boundary_condition_t use mod_input, only: input_t @@ -183,7 +183,7 @@ function symmetry_bc_constructor(location, input, grid) result(bc) type(symmetry_bc_t), pointer :: bc character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y) class(input_t), intent(in) :: input - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid allocate(bc) bc%name = 'symmetry' diff --git a/src/lib/boundary_conditions/zero_gradient_bc.f90 b/src/lib/boundary_conditions/zero_gradient_bc.f90 index a0b0546..d4172ae 100644 --- a/src/lib/boundary_conditions/zero_gradient_bc.f90 +++ b/src/lib/boundary_conditions/zero_gradient_bc.f90 @@ -22,7 +22,7 @@ module mod_zero_gradient_bc use, intrinsic :: iso_fortran_env, only: ik => int32, rk => real64 use mod_globals, only: enable_debug_print, debug_print use mod_field, only: field_2d_t - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_boundary_conditions, only: boundary_condition_t use mod_input, only: input_t @@ -42,7 +42,7 @@ function zero_gradient_bc_constructor(location, input, grid) result(bc) type(zero_gradient_bc_t), pointer :: bc character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y) class(input_t), intent(in) :: input - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid allocate(bc) bc%name = 'zero_gradient' diff --git a/src/lib/error.f90 b/src/lib/error.f90 index 05235f3..7ba12d1 100644 --- a/src/lib/error.f90 +++ b/src/lib/error.f90 @@ -20,20 +20,22 @@ subroutine error_msg(message, module_name, class_name, procedure_name, file_name logical, intent(in), optional :: error_stop if(present(class_name)) then - write(std_err, '(a, i0)') "Error: "//trim(message)//"; in "//module_name & + write(std_out, '(a, i0)') "Error: "//trim(message)//"; in "//module_name & //"::"//class_name//"%"//procedure_name// & "() in "//file_name//":", line_number else - write(std_err, '(a, i0)') "Error: "//trim(message)//"; in "//module_name & + write(std_out, '(a, i0)') "Error: "//trim(message)//"; in "//module_name & //"::"//procedure_name//"() in "//file_name//":", line_number end if + ! Catch-all error code if(present(error_stop)) then if(error_stop) error stop 1 else error stop 1 end if + flush(std_err) end subroutine error_msg end module mod_error diff --git a/src/lib/evolution_operator/local_evo_operator.f90 b/src/lib/evolution_operator/local_evo_operator.f90 index cb54de3..3358e37 100644 --- a/src/lib/evolution_operator/local_evo_operator.f90 +++ b/src/lib/evolution_operator/local_evo_operator.f90 @@ -33,7 +33,7 @@ module mod_local_evo_operator use mod_input, only: input_t use mod_mach_cone_collection, only: mach_cone_collection_t - use mod_grid_block, only: grid_block_t, C1, M1, C2, M2, C3, M3, C4, M4 + use mod_grid_block_2d, only: grid_block_2d_t, C1, M1, C2, M2, C3, M3, C4, M4 use mod_abstract_reconstruction, only: abstract_reconstruction_t implicit none @@ -47,7 +47,7 @@ module mod_local_evo_operator !< operator will inherit from this class. Most of the attributes are the same, just that the !< implementation varies between those who inherit this class - class(grid_block_t), pointer :: grid => null() + class(grid_block_2d_t), pointer :: grid => null() !< pointer to the grid object real(rk), dimension(:, :, :), pointer :: reconstructed_rho => null() @@ -79,7 +79,7 @@ module mod_local_evo_operator subroutine initialize(self, dt, grid_target, recon_operator_target) !< Constructor for the FVLEG operator class(local_evo_operator_t), intent(inout) :: self - class(grid_block_t), intent(in), target :: grid_target + class(grid_block_2d_t), intent(in), target :: grid_target real(rk), intent(in) :: dt !< timestep class(abstract_reconstruction_t), intent(in), target :: recon_operator_target diff --git a/src/lib/fluid/mod_fluid.f90 b/src/lib/fluid/mod_fluid.f90 index c6395c6..324c0be 100644 --- a/src/lib/fluid/mod_fluid.f90 +++ b/src/lib/fluid/mod_fluid.f90 @@ -135,7 +135,7 @@ function new_fluid(input, grid) result(fluid) type(fluid_t), pointer :: fluid select type(grid) - class is (grid_block_2d_t) + class is(grid_block_2d_t) allocate(fluid) call fluid%initialize(input, grid) end select @@ -162,39 +162,39 @@ subroutine initialize(self, input, grid) self%rho = field_2d(name='rho', long_name='Density', & descrip='Cell Density', units='g/cm^3', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%rho_u = field_2d(name='rhou', long_name='rhou', descrip='Cell Conserved quantity (Density * X-Velocity)', & units='g cm/cm^2 s', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%rho_v = field_2d(name='rhov', long_name='rhov', descrip='Cell Conserved quantity (Density * Y-Velocity)', & units='g cm/cm^2 s', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%rho_E = field_2d(name='rhoE', long_name='rhoE', descrip='Cell Conserved quantity (Density * Total Energy)', & units='g erg / cm^3', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%u = field_2d(name='u', long_name='X Velocity', descrip='Cell X-Velocity', units='cm/s', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%v = field_2d(name='v', long_name='Y Velocity', descrip='Cell Y-Velocity', units='cm/s', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%p = field_2d(name='p', long_name='Pressure', descrip='Cell Pressure', units='barye', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%cs = field_2d(name='cs', long_name='Sound Speed', descrip='Cell Sound Speed', units='cm/s', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%mach_u = field_2d(name='mach_u', long_name='Mach X', descrip='Cell Mach number in x-direction', & units='dimensionless', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%mach_v = field_2d(name='mach_v', long_name='Mach Y', descrip='Cell Mach number in y-direction', & units='dimensionless', & - global_dims=grid%cell_dim_global, n_halo_cells=input%n_ghost_layers) + global_dims=grid%global_dims, n_halo_cells=input%n_ghost_layers) self%smooth_residuals = input%smooth_residuals @@ -371,19 +371,21 @@ subroutine initialize_from_hdf5(self, input) call eos%primitive_to_conserved(rho=self%rho, u=self%u, v=self%v, p=self%p, & rho_u=self%rho_u, rho_v=self%rho_v, rho_E=self%rho_E) - write(*, '(a)') 'Initial fluid stats' - write(*, '(a)') '===================================================' - write(*, '(a, f0.4)') 'EOS Gamma: ', eos%get_gamma() - write(*, '(a, 2(es16.6, 1x))') 'Min/Max density [non-dim]: ', minval(density), maxval(density) - write(*, '(a, 2(es16.6, 1x))') 'Min/Max x-velocity [non-dim]: ', minval(x_velocity), maxval(x_velocity) - write(*, '(a, 2(es16.6, 1x))') 'Min/Max y-velocity [non-dim]: ', minval(x_velocity), maxval(y_velocity) - write(*, '(a, 2(es16.6, 1x))') 'Min/Max pressure [non-dim]: ', minval(pressure), maxval(pressure) - write(*, '(a, 2(es16.6, 1x))') 'Min/Max density [dim]: ', minval(density) * rho_0, maxval(density) * rho_0 - write(*, '(a, 2(es16.6, 1x))') 'Min/Max x-velocity [dim]: ', minval(x_velocity) * v_0, maxval(x_velocity) * v_0 - write(*, '(a, 2(es16.6, 1x))') 'Min/Max y-velocity [dim]: ', minval(y_velocity) * v_0, maxval(y_velocity) * v_0 - write(*, '(a, 2(es16.6, 1x))') 'Min/Max pressure [dim]: ', minval(pressure) * p_0, maxval(pressure) * p_0 - write(*, '(a)') '===================================================' - write(*, *) + if (this_image() == 1) then + write(*, '(a)') 'Initial fluid stats' + write(*, '(a)') '===================================================' + write(*, '(a, f0.4)') 'EOS Gamma: ', eos%get_gamma() + write(*, '(a, 2(es16.6, 1x))') 'Min/Max density [non-dim]: ', minval(density), maxval(density) + write(*, '(a, 2(es16.6, 1x))') 'Min/Max x-velocity [non-dim]: ', minval(x_velocity), maxval(x_velocity) + write(*, '(a, 2(es16.6, 1x))') 'Min/Max y-velocity [non-dim]: ', minval(x_velocity), maxval(y_velocity) + write(*, '(a, 2(es16.6, 1x))') 'Min/Max pressure [non-dim]: ', minval(pressure), maxval(pressure) + write(*, '(a, 2(es16.6, 1x))') 'Min/Max density [dim]: ', minval(density) * rho_0, maxval(density) * rho_0 + write(*, '(a, 2(es16.6, 1x))') 'Min/Max x-velocity [dim]: ', minval(x_velocity) * v_0, maxval(x_velocity) * v_0 + write(*, '(a, 2(es16.6, 1x))') 'Min/Max y-velocity [dim]: ', minval(y_velocity) * v_0, maxval(y_velocity) * v_0 + write(*, '(a, 2(es16.6, 1x))') 'Min/Max pressure [dim]: ', minval(pressure) * p_0, maxval(pressure) * p_0 + write(*, '(a)') '===================================================' + write(*, *) + endif end subroutine initialize_from_hdf5 @@ -475,22 +477,19 @@ subroutine integrate(self, dt, grid, error_code) real(rk), intent(in) :: dt !< time step class(grid_block_t), intent(in) :: grid !< grid class - the solver needs grid topology integer(ik), intent(out) :: error_code - + if(enable_debug_print) call debug_print('Running fluid_t%integerate()', __FILE__, __LINE__) self%time = self%time + dt self%dt = dt - select type(grid) - class is (grid_block_2d_t) - select case(trim(self%time_integration_scheme)) - case('ssp_rk2') - call self%ssp_rk2(grid, error_code) - case('ssp_rk3') - call self%ssp_rk3(grid, error_code) - case default - call error_msg(module_name='mod_fluid', class_name='fluid_t', procedure_name='assign_fluid', & - message="Unknown time integration scheme", file_name=__FILE__, line_number=__LINE__) - end select + select case(trim(self%time_integration_scheme)) + case('ssp_rk2') + call self%ssp_rk2(grid, error_code) + case('ssp_rk3') + call self%ssp_rk3(grid, error_code) + case default + call error_msg(module_name='mod_fluid', class_name='fluid_t', procedure_name='assign_fluid', & + message="Unknown time integration scheme", file_name=__FILE__, line_number=__LINE__) end select end subroutine integrate @@ -526,19 +525,21 @@ real(rk) function get_timestep(self, grid) result(delta_t) real(rk), allocatable, save, dimension(:, :) :: dx, dy - g_ilo = grid%cell_lbounds(1) - g_ihi = grid%cell_ubounds(1) + ! This seems silly to have to do, but the master class requires + ! that the grid is a grid_block_t type, even though this fluid class + ! will always use a grid_block_2d_t type. + select type(grid) + class is(grid_block_2d_t) + g_ilo = grid%lbounds(1) + g_ihi = grid%ubounds(1) + g_jlo = grid%lbounds(2) + g_jhi = grid%ubounds(2) - g_jlo = grid%cell_lbounds(2) - g_jhi = grid%cell_ubounds(2) + if(.not. allocated(dx)) allocate(dx(g_ilo:g_ihi, g_jlo:g_jhi)) + if(.not. allocated(dy)) allocate(dy(g_ilo:g_ihi, g_jlo:g_jhi)) - if (.not. allocated(dx)) allocate(dx(g_ilo:g_ihi, g_jlo:g_jhi)) - if (.not. allocated(dy)) allocate(dy(g_ilo:g_ihi, g_jlo:g_jhi)) - - select type(grid) - class is (grid_block_2d_t) + dx = grid%dx(g_ilo:g_ihi, g_jlo:g_jhi) dy = grid%dy(g_ilo:g_ihi, g_jlo:g_jhi) - dx = grid%dy(g_ilo:g_ihi, g_jlo:g_jhi) end select err_msg = '' @@ -548,21 +549,20 @@ real(rk) function get_timestep(self, grid) result(delta_t) jlo = self%u%lbounds(2) jhi = self%u%ubounds(2) - print*,'fluid: ', ilo, ihi, jlo, jhi - print*,'grid : ', g_ilo, g_ihi, g_jlo, g_jhi + ! I would have put this in a cleaner associate block, but GFortran+OpenCoarrays bugs out on this coarray_max_delta_t = minval(self%cfl / & - (((abs(self%u%data(ilo:ihi, jlo:jhi)) + & - self%cs%data(ilo:ihi, jlo:jhi)) / dx) + & - ((abs(self%v%data(ilo:ihi, jlo:jhi)) + & - self%cs%data(ilo:ihi, jlo:jhi)) / dy))) - + (((abs(self%u%data(ilo:ihi, jlo:jhi)) + & + self%cs%data(ilo:ihi, jlo:jhi)) / dx) + & + ((abs(self%v%data(ilo:ihi, jlo:jhi)) + & + self%cs%data(ilo:ihi, jlo:jhi)) / dy))) + sync all ! Get the minimum timestep across all the images and save it on each image call co_min(coarray_max_delta_t, stat=ierr) - if (ierr /= 0) then + if(ierr /= 0) then call error_msg(module_name='mod_fluid', class_name='fluid_t', procedure_name='get_timestep', & - message="Unable to run the co_min() to get min timestep across images: err_msg=" // trim(err_msg), & + message="Unable to run the co_min() to get min timestep across images: err_msg="//trim(err_msg), & file_name=__FILE__, line_number=__LINE__) end if delta_t = coarray_max_delta_t @@ -797,7 +797,7 @@ end subroutine sanity_check subroutine ssp_rk3(U, grid, error_code) !< Strong-stability preserving Runge-Kutta 3rd order class(fluid_t), intent(inout) :: U - class(grid_block_2d_t), intent(in) :: grid + class(grid_block_t), intent(in) :: grid type(fluid_t), allocatable :: U_1 !< first stage type(fluid_t), allocatable :: U_2 !< second stage @@ -810,25 +810,29 @@ subroutine ssp_rk3(U, grid, error_code) allocate(U_1, source=U) allocate(U_2, source=U) - ! 1st stage + select type(grid) + class is(grid_block_2d_t) + ! 1st stage if(enable_debug_print) call debug_print(new_line('a')//'Running fluid_t%ssp_rk3() 1st stage'//new_line('a'), __FILE__, __LINE__) - call U%apply_bc() - U_1 = U + U%t(grid, stage=1) * dt + call U%apply_bc() + U_1 = U + U%t(grid, stage=1) * dt - ! 2nd stage + ! 2nd stage if(enable_debug_print) call debug_print(new_line('a')//'Running fluid_t%ssp_rk3() 2nd stage'//new_line('a'), __FILE__, __LINE__) - call U_1%apply_bc() - U_2 = U * (3.0_rk / 4.0_rk) & - + U_1 * (1.0_rk / 4.0_rk) & - + U_1%t(grid, stage=2) * ((1.0_rk / 4.0_rk) * dt) + call U_1%apply_bc() + U_2 = U * (3.0_rk / 4.0_rk) & + + U_1 * (1.0_rk / 4.0_rk) & + + U_1%t(grid, stage=2) * ((1.0_rk / 4.0_rk) * dt) - ! Final stage + ! Final stage if(enable_debug_print) call debug_print(new_line('a')//'Running fluid_t%ssp_rk2() 3rd stage'//new_line('a'), __FILE__, __LINE__) - call U_2%apply_bc() - U = U * (1.0_rk / 3.0_rk) & - + U_2 * (2.0_rk / 3.0_rk) & - + U_2%t(grid, stage=3) * ((2.0_rk / 3.0_rk) * dt) + call U_2%apply_bc() + U = U * (1.0_rk / 3.0_rk) & + + U_2 * (2.0_rk / 3.0_rk) & + + U_2%t(grid, stage=3) * ((2.0_rk / 3.0_rk) * dt) + end select + ! Check for NaNs and improper negatives (e.g. in pressure and density) call U%sanity_check(error_code) ! Convergence history @@ -842,7 +846,7 @@ end subroutine ssp_rk3 subroutine ssp_rk2(U, grid, error_code) !< Strong-stability preserving Runge-Kutta 2nd order class(fluid_t), intent(inout) :: U - class(grid_block_2d_t), intent(in) :: grid + class(grid_block_t), intent(in) :: grid type(fluid_t), allocatable :: U_1 !< first stage integer(ik), intent(out) :: error_code @@ -855,20 +859,25 @@ subroutine ssp_rk2(U, grid, error_code) ! 1st stage if(enable_debug_print) then - call debug_print(new_line('a')//'Running fluid_t%ssp_rk2_t() 1st stage'//new_line('a'), & + call debug_print('Running fluid_t%ssp_rk2_t() 1st stage', & __FILE__, __LINE__) end if - call U%apply_bc() - U_1 = U + U%t(grid, stage=1) * dt - ! Final stage - if(enable_debug_print) then - call debug_print(new_line('a')//'Running fluid_t%ssp_rk2_t() 2nd stage'//new_line('a'), & - __FILE__, __LINE__) - end if - call U_1%apply_bc() - U = U * 0.5_rk + U_1 * 0.5_rk + U_1%t(grid, stage=2) * (0.5_rk * dt) + select type(grid) + class is(grid_block_2d_t) + call U%apply_bc() + U_1 = U + U%t(grid, stage=1) * dt + + ! Final stage + if(enable_debug_print) then + call debug_print('Running fluid_t%ssp_rk2_t() 2nd stage', & + __FILE__, __LINE__) + end if + call U_1%apply_bc() + U = U * 0.5_rk + U_1 * 0.5_rk + U_1%t(grid, stage=2) * (0.5_rk * dt) + end select + ! Check for NaNs and improper negatives (e.g. in pressure and density) call U%sanity_check(error_code) ! Convergence history diff --git a/src/lib/flux_solvers/ausm_plus_solver.f90 b/src/lib/flux_solvers/ausm_plus_solver.f90 index 1c714ae..630432b 100644 --- a/src/lib/flux_solvers/ausm_plus_solver.f90 +++ b/src/lib/flux_solvers/ausm_plus_solver.f90 @@ -36,7 +36,7 @@ module mod_ausm_plus_solver use mod_muscl_interpolation, only: muscl_interpolation_t use mod_flux_solver, only: edge_split_flux_solver_t use mod_eos, only: eos - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_input, only: input_t implicit none @@ -142,7 +142,7 @@ end subroutine copy_ausm_plus subroutine solve_ausm_plus(self, dt, grid, lbounds, rho, u, v, p, d_rho_dt, d_rho_u_dt, d_rho_v_dt, d_rho_E_dt) !< Solve and flux the edges class(ausm_plus_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(2), intent(in) :: lbounds real(rk), intent(in) :: dt !< timestep (not really used in this solver, but needed for others) real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(inout) :: rho !< density diff --git a/src/lib/flux_solvers/fvleg_solver.f90 b/src/lib/flux_solvers/fvleg_solver.f90 index 8b36623..b3118c7 100644 --- a/src/lib/flux_solvers/fvleg_solver.f90 +++ b/src/lib/flux_solvers/fvleg_solver.f90 @@ -35,7 +35,7 @@ module mod_fvleg_solver use mod_abstract_reconstruction, only: abstract_reconstruction_t use mod_reconstruction_factory, only: reconstruction_factory use mod_flux_array, only: get_fluxes, flux_array_t - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_input, only: input_t implicit none @@ -96,7 +96,7 @@ subroutine solve_fvleg(self, dt, grid, lbounds, rho, u, v, p, & d_rho_dt, d_rho_u_dt, d_rho_v_dt, d_rho_E_dt) !< Solve and flux the edges class(fvleg_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(2), intent(in) :: lbounds real(rk), intent(in) :: dt real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(inout) :: rho !< density @@ -165,8 +165,8 @@ subroutine solve_fvleg(self, dt, grid, lbounds, rho, u, v, p, & call debug_print('Running fvleg_t%solve_fvleg()', __FILE__, __LINE__) - associate(imin => grid%cell_lbounds_halo(1), imax => grid%cell_ubounds_halo(1), & - jmin => grid%cell_lbounds_halo(2), jmax => grid%cell_ubounds_halo(2)) + associate(imin => grid%lbounds_halo(1), imax => grid%ubounds_halo(1), & + jmin => grid%lbounds_halo(2), jmax => grid%ubounds_halo(2)) allocate(rho_recon_state(1:8, imin:imax, jmin:jmax)) allocate(u_recon_state(1:8, imin:imax, jmin:jmax)) @@ -325,7 +325,7 @@ subroutine flux_edges(grid, & evolved_du_mid_rho, evolved_du_mid_u, evolved_du_mid_v, evolved_du_mid_p, & d_rho_dt, d_rhou_dt, d_rhov_dt, d_rhoE_dt) !< Evaluate the fluxes along the edges. This is equation 13 in the paper - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid real(rk), dimension(grid%ilo_node:, grid%jlo_node:), contiguous, intent(in) :: evolved_corner_rho !< (i,j); Reconstructed rho at the corners real(rk), dimension(grid%ilo_node:, grid%jlo_node:), contiguous, intent(in) :: evolved_corner_u !< (i,j); Reconstructed u at the corners @@ -339,10 +339,10 @@ subroutine flux_edges(grid, & real(rk), dimension(grid%ilo_node:, grid%cell_lbounds(2):), contiguous, intent(in) :: evolved_du_mid_u !< (i,j); Reconstructed u at the down/up midpoints real(rk), dimension(grid%ilo_node:, grid%cell_lbounds(2):), contiguous, intent(in) :: evolved_du_mid_v !< (i,j); Reconstructed v at the down/up midpoints real(rk), dimension(grid%ilo_node:, grid%cell_lbounds(2):), contiguous, intent(in) :: evolved_du_mid_p !< (i,j); Reconstructed p at the down/up midpoints - real(rk), dimension(grid%cell_lbounds_halo(1):, grid%cell_lbounds_halo(2):), contiguous, intent(inout) :: d_rho_dt - real(rk), dimension(grid%cell_lbounds_halo(1):, grid%cell_lbounds_halo(2):), contiguous, intent(inout) :: d_rhou_dt - real(rk), dimension(grid%cell_lbounds_halo(1):, grid%cell_lbounds_halo(2):), contiguous, intent(inout) :: d_rhov_dt - real(rk), dimension(grid%cell_lbounds_halo(1):, grid%cell_lbounds_halo(2):), contiguous, intent(inout) :: d_rhoE_dt + real(rk), dimension(grid%lbounds_halo(1):, grid%lbounds_halo(2):), contiguous, intent(inout) :: d_rho_dt + real(rk), dimension(grid%lbounds_halo(1):, grid%lbounds_halo(2):), contiguous, intent(inout) :: d_rhou_dt + real(rk), dimension(grid%lbounds_halo(1):, grid%lbounds_halo(2):), contiguous, intent(inout) :: d_rhov_dt + real(rk), dimension(grid%lbounds_halo(1):, grid%lbounds_halo(2):), contiguous, intent(inout) :: d_rhoE_dt integer(ik) :: ilo, ihi, jlo, jhi integer(ik) :: i, j, k, edge, xy @@ -507,8 +507,8 @@ subroutine flux_edges(grid, & !$omp end do simd !$omp end parallel - ilo = grid%cell_lbounds_halo(1); ihi = grid%cell_ubounds_halo(1) - jlo = grid%cell_lbounds_halo(2); jhi = grid%cell_ubounds_halo(2) + ilo = grid%lbounds_halo(1); ihi = grid%ubounds_halo(1) + jlo = grid%lbounds_halo(2); jhi = grid%ubounds_halo(2) d_rho_dt(ilo, :) = 0.0_rk d_rho_dt(ihi, :) = 0.0_rk d_rho_dt(:, jlo) = 0.0_rk diff --git a/src/lib/flux_solvers/m_ausmpw_plus_solver.fypp b/src/lib/flux_solvers/m_ausmpw_plus_solver.fypp index 06152ff..65ca82b 100644 --- a/src/lib/flux_solvers/m_ausmpw_plus_solver.fypp +++ b/src/lib/flux_solvers/m_ausmpw_plus_solver.fypp @@ -54,7 +54,7 @@ module mod_m_ausmpw_plus_solver use mod_muscl_interpolation, only: muscl_interpolation_t use mod_flux_solver, only: edge_split_flux_solver_t use mod_eos, only: eos - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_input, only: input_t implicit none @@ -112,7 +112,7 @@ contains subroutine solve_m_ausmpw_plus(self, dt, grid, lbounds, rho, u, v, p, d_rho_dt, d_rho_u_dt, d_rho_v_dt, d_rho_E_dt) !< Solve and flux the edges class(m_ausmpw_plus_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(2), intent(in) :: lbounds real(rk), intent(in) :: dt !< timestep (not really used in this solver, but needed for others) real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(inout) :: rho !< density @@ -220,10 +220,10 @@ contains ! ! This is the numbering convention that this module uses - ilo = grid%ilo_cell - ihi = grid%ihi_cell - jlo = grid%jlo_cell - jhi = grid%jhi_cell + ilo = grid%lbounds(1) + ihi = grid%ubounds(1) + jlo = grid%lbounds(2) + jhi = grid%ubounds(2) if(allocated(self%iflux)) deallocate(self%iflux) if(allocated(self%jflux)) deallocate(self%jflux) @@ -231,8 +231,8 @@ contains allocate(self%iflux(4, ilo - 1:ihi, jlo:jhi)) allocate(self%jflux(4, ilo:ihi, jlo - 1:jhi)) - ilo_bc = grid%cell_lbounds_halo(1) - ihi_bc = grid%cell_ubounds_halo(1) + ilo_bc = grid%lbounds_halo(1) + ihi_bc = grid%ubounds_halo(1) jlo_bc = grid%jlo_bc_cell jhi_bc = grid%jhi_bc_cell @@ -339,7 +339,7 @@ contains rho_ave, u_ave, v_ave, p_ave, w2, eflux_lbounds, edge_flux) !< Construct the fluxes for each edge in the ${DIR}$ direction. This is templated via the Fypp pre-processor class(m_ausmpw_plus_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(3), intent(in) :: lbounds !< bounds of the primitive variable arrays real(rk), dimension(lbounds(2):, lbounds(3):), contiguous, intent(in) :: rho_ave !< (i,j); cell averaged value ofdensity; needed for critical Mach number calcs real(rk), dimension(lbounds(2):, lbounds(3):), contiguous, intent(in) :: u_ave !< (i,j); cell averaged value of x-velocity; needed for critical Mach number calcs diff --git a/src/lib/flux_solvers/slau_solver.fypp b/src/lib/flux_solvers/slau_solver.fypp index 72b83bb..1eeea08 100644 --- a/src/lib/flux_solvers/slau_solver.fypp +++ b/src/lib/flux_solvers/slau_solver.fypp @@ -51,7 +51,7 @@ module mod_slau_solver use mod_muscl_interpolation, only: muscl_interpolation_t use mod_flux_solver, only: edge_split_flux_solver_t use mod_eos, only: eos - use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t use mod_input, only: input_t implicit none @@ -120,7 +120,7 @@ contains subroutine solve_slau(self, dt, grid, lbounds, rho, u, v, p, d_rho_dt, d_rho_u_dt, d_rho_v_dt, d_rho_E_dt) !< Solve and flux the edges class(slau_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(2), intent(in) :: lbounds real(rk), intent(in) :: dt !< timestep (not really used in this solver, but needed for others) real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(inout) :: rho !< density @@ -196,10 +196,10 @@ contains ! ! This is the numbering convention that this module uses - ilo = grid%ilo_cell - ihi = grid%ihi_cell - jlo = grid%jlo_cell - jhi = grid%jhi_cell + ilo = grid%lbounds(1) + ihi = grid%ubounds(1) + jlo = grid%lbounds(2) + jhi = grid%ubounds(2) if(allocated(self%iflux)) deallocate(self%iflux) if(allocated(self%jflux)) deallocate(self%jflux) @@ -207,8 +207,8 @@ contains allocate(self%iflux(4, ilo - 1:ihi, jlo:jhi)) allocate(self%jflux(4, ilo:ihi, jlo - 1:jhi)) - ilo_bc = grid%cell_lbounds_halo(1) - ihi_bc = grid%cell_ubounds_halo(1) + ilo_bc = grid%lbounds_halo(1) + ihi_bc = grid%ubounds_halo(1) jlo_bc = grid%jlo_bc_cell jhi_bc = grid%jhi_bc_cell @@ -259,7 +259,7 @@ contains rho_ave, u_ave, v_ave, p_ave, eflux_lbounds, edge_flux) !< Construct the fluxes for each edge in the ${DIR}$ direction. This is templated via the Fypp pre-processor class(slau_solver_t), intent(inout) :: self - class(grid_block_t), intent(in) :: grid + class(grid_block_2d_t), intent(in) :: grid integer(ik), dimension(3), intent(in) :: lbounds !< bounds of the primitive variable arrays real(rk), dimension(lbounds(2):, lbounds(3):), contiguous, intent(in) :: rho_ave !< (i,j); cell averaged value ofdensity; needed for critical Mach number calcs real(rk), dimension(lbounds(2):, lbounds(3):), contiguous, intent(in) :: u_ave !< (i,j); cell averaged value of x-velocity; needed for critical Mach number calcs diff --git a/src/lib/globals.f90 b/src/lib/globals.f90 index f89750e..8514e06 100644 --- a/src/lib/globals.f90 +++ b/src/lib/globals.f90 @@ -24,7 +24,7 @@ module mod_globals implicit none - logical, parameter :: enable_debug_print = .true. + logical, parameter :: enable_debug_print = .false. logical, parameter :: enable_file_and_line_stats = .false. logical, parameter :: plot_gradients = .false. logical, parameter :: plot_limiters = .false. @@ -130,19 +130,15 @@ subroutine debug_print(str, file, line_number) integer, intent(in), optional :: line_number if(enable_debug_print) then - ! if(this_image() == 1) then - if(enable_file_and_line_stats) then if(present(file) .and. present(line_number)) then - write(std_out, '(a,a,i0,3(a))') file, ':', line_number, ' "', str, '"' + write(std_out, '(a,i0,4(a),i0,3(a))') "Image: ", this_image(), ", ", file, ':', line_number, ' "', str, '"' else - write(std_out, '(a)') str + write(std_out, '(a,i0,2(a))') "Image: ", this_image(), ", ", str end if else - write(std_out, '(a)') str + write(std_out, '(a,i0,2(a))') "Image: ", this_image(), ", ", str end if - - ! end if end if end subroutine diff --git a/src/lib/grid/grid_block.f90 b/src/lib/grid/grid_block.f90 index 2cc6373..13353aa 100644 --- a/src/lib/grid/grid_block.f90 +++ b/src/lib/grid/grid_block.f90 @@ -23,40 +23,28 @@ module mod_grid_block end type type, abstract :: grid_block_t - - integer(ik) :: host_image_id = 0 !< # of halo cells used in this block; typically 2 or 3 depending on spatial order - integer(ik) :: n_halo_cells = 0 !< # of halo cells used in this block; typically 2 or 3 depending on spatial order logical :: is_axisymmetric = .false. !< Is this an axisymmetric grid? (need for volume computation) logical :: is_uniform = .false. !< Are all the cells the same size/shape? - - ! Bounds information - integer(ik), dimension(:), allocatable :: global_node_dims !< (i, j, k); # of array indices in each dimension - integer(ik), dimension(:), allocatable :: global_cell_dims !< (i, j, k); # of array indices in each dimension - - integer(ik), dimension(:), allocatable :: domain_node_shape !< (i, j, k); # of array indices in each dimension - integer(ik), dimension(:), allocatable :: domain_cell_shape !< (i, j, k); # of array indices in each dimension - - integer(ik), dimension(:), allocatable :: global_node_dims_halo !< (i, j, k); # of array indices in each dimension - integer(ik), dimension(:), allocatable :: global_cell_dims_halo !< (i, j, k); # of array indices in each dimension - - integer(ik), dimension(:), allocatable :: domain_node_shape_halo !< (i, j, k); # of array indices in each dimension - integer(ik), dimension(:), allocatable :: domain_cell_shape_halo !< (i, j, k); # of array indices in each dimension - - integer(ik), dimension(:), allocatable :: node_lbounds !< (i, j, k); lower bounds (not including halo cells) - integer(ik), dimension(:), allocatable :: node_ubounds !< (i, j, k); upper bounds (not including halo cells) - - integer(ik), dimension(:), allocatable :: node_lbounds_halo !< (i, j, k); lower bounds (including halo cells) - integer(ik), dimension(:), allocatable :: node_ubounds_halo !< (i, j, k); upper bounds (including halo cells) - - integer(ik), dimension(:), allocatable :: cell_lbounds !< (i, j, k); lower bounds (not including halo cells) - integer(ik), dimension(:), allocatable :: cell_ubounds !< (i, j, k); upper bounds (not including halo cells) - - integer(ik), dimension(:), allocatable :: cell_lbounds_halo !< (i, j, k); lower bounds (including halo cells) - integer(ik), dimension(:), allocatable :: cell_ubounds_halo !< (i, j, k); upper bounds (including halo cells) - - ! Boundary condition checks - logical, dimension(:,:), allocatable :: on_bc !< ((lo,hi),(i,j,k)) does this grid live on the ihi boundary? - + integer(ik) :: host_image_id = 0 !< # of halo cells used in this block; typically 2 or 3 depending on spatial order + integer(ik) :: n_halo_cells = 0 !< # of halo cells used in this block; typically 2 or 3 depending on spatial order + integer(ik) :: total_cells = 0 !< + integer(ik), dimension(3) :: block_dims = 0 !< (ni,nj,nk); local block cell dimensions + integer(ik), dimension(3) :: global_dims = 0 !< (ni,nj,nk); global cell dimensions + integer(ik), dimension(3) :: lbounds = 0 !< (i,j,k); block lower cell bounds excluding halo cells + integer(ik), dimension(3) :: ubounds = 0 !< (i,j,k); block upper cell bounds excluding halo cells + integer(ik), dimension(3) :: lbounds_global = 0 !< (i,j,k); block lower cell bounds excluding halo cells + integer(ik), dimension(3) :: ubounds_global = 0 !< (i,j,k); block upper cell bounds excluding halo cells + integer(ik), dimension(3) :: lbounds_halo = 0 !< (i,j,k); block lower cell bounds including halo cells + integer(ik), dimension(3) :: ubounds_halo = 0 !< (i,j,k); block upper cell bounds including halo cells + + logical :: on_ilo_bc = .false. + logical :: on_ihi_bc = .false. + logical :: on_jlo_bc = .false. + logical :: on_jhi_bc = .false. + logical :: on_klo_bc = .false. + logical :: on_khi_bc = .false. + + integer(ik), dimension(:), allocatable :: neighbors !< (direction); neighbor image contains procedure(initialize), deferred :: initialize end type grid_block_t @@ -70,5 +58,4 @@ subroutine initialize(self, input) end interface contains - end module mod_grid_block diff --git a/src/lib/grid/grid_block_1d.f90 b/src/lib/grid/grid_block_1d.f90 index 4c362c7..4e1fa5b 100644 --- a/src/lib/grid/grid_block_1d.f90 +++ b/src/lib/grid/grid_block_1d.f90 @@ -35,32 +35,32 @@ subroutine init_1d_block(self, input) class(grid_block_1d_t), intent(inout) :: self class(input_t), intent(in) :: input - allocate(self%global_node_dims (1)) - self%global_node_dims = 0 - allocate(self%global_cell_dims (1)) - self%global_cell_dims = 0 + allocate(self%global_node_dims(1)) + self%global_node_dims = 0 + allocate(self%global_cell_dims(1)) + self%global_cell_dims = 0 allocate(self%domain_node_shape(1)) self%domain_node_shape = 0 allocate(self%domain_cell_shape(1)) self%domain_cell_shape = 0 - allocate(self%node_lbounds (1)) - self%node_lbounds = 0 - allocate(self%node_ubounds (1)) - self%node_ubounds = 0 + allocate(self%node_lbounds(1)) + self%node_lbounds = 0 + allocate(self%node_ubounds(1)) + self%node_ubounds = 0 allocate(self%node_lbounds_halo(1)) self%node_lbounds_halo = 0 allocate(self%node_ubounds_halo(1)) self%node_ubounds_halo = 0 - allocate(self%cell_lbounds (1)) - self%cell_lbounds = 0 - allocate(self%cell_ubounds (1)) - self%cell_ubounds = 0 + allocate(self%cell_lbounds(1)) + self%cell_lbounds = 0 + allocate(self%cell_ubounds(1)) + self%cell_ubounds = 0 allocate(self%cell_lbounds_halo(1)) self%cell_lbounds_halo = 0 allocate(self%cell_ubounds_halo(1)) self%cell_ubounds_halo = 0 - allocate(self%on_bc (2,1)) - self%on_bc = .false. + allocate(self%on_bc(2, 1)) + self%on_bc = .false. end subroutine subroutine finalize_1d_block(self) @@ -74,16 +74,16 @@ subroutine finalize_1d_block(self) if(allocated(self%centroid_y)) deallocate(self%centroid_y) if(allocated(self%edge_lengths)) deallocate(self%edge_lengths) if(allocated(self%edge_norm_vectors)) deallocate(self%edge_norm_vectors) - if(allocated(self%global_node_dims )) deallocate(self%global_node_dims ) - if(allocated(self%global_cell_dims )) deallocate(self%global_cell_dims ) + if(allocated(self%global_node_dims)) deallocate(self%global_node_dims) + if(allocated(self%global_cell_dims)) deallocate(self%global_cell_dims) if(allocated(self%domain_node_shape)) deallocate(self%domain_node_shape) if(allocated(self%domain_cell_shape)) deallocate(self%domain_cell_shape) - if(allocated(self%node_lbounds )) deallocate(self%node_lbounds ) - if(allocated(self%node_ubounds )) deallocate(self%node_ubounds ) + if(allocated(self%node_lbounds)) deallocate(self%node_lbounds) + if(allocated(self%node_ubounds)) deallocate(self%node_ubounds) if(allocated(self%node_lbounds_halo)) deallocate(self%node_lbounds_halo) if(allocated(self%node_ubounds_halo)) deallocate(self%node_ubounds_halo) - if(allocated(self%cell_lbounds )) deallocate(self%cell_lbounds ) - if(allocated(self%cell_ubounds )) deallocate(self%cell_ubounds ) + if(allocated(self%cell_lbounds)) deallocate(self%cell_lbounds) + if(allocated(self%cell_ubounds)) deallocate(self%cell_ubounds) if(allocated(self%cell_lbounds_halo)) deallocate(self%cell_lbounds_halo) if(allocated(self%cell_ubounds_halo)) deallocate(self%cell_ubounds_halo) if(allocated(self%on_bc)) deallocate(self%on_bc) diff --git a/src/lib/grid/grid_block_2d.f90 b/src/lib/grid/grid_block_2d.f90 index 02bb2b5..5fbfe81 100644 --- a/src/lib/grid/grid_block_2d.f90 +++ b/src/lib/grid/grid_block_2d.f90 @@ -15,8 +15,7 @@ module mod_grid_block_2d private public :: grid_block_2d_t, new_2d_grid_block - type :: grid_block_2d_t - logical :: is_uniform = .false. + type, extends(grid_block_t) :: grid_block_2d_t real(rk), dimension(:, :), allocatable :: volume !< (i, j); volume of each cell real(rk), dimension(:, :), allocatable :: dx !< (i, j); dx spacing of each cell real(rk), dimension(:, :), allocatable :: dy !< (i, j); dy spacing of each cell @@ -26,56 +25,12 @@ module mod_grid_block_2d real(rk), dimension(:, :), allocatable :: centroid_y !< (i, j); y location of the cell centroid real(rk), dimension(:, :, :), allocatable :: edge_lengths !< ((edge_1:edge_n), i, j); length of each edge real(rk), dimension(:, :, :, :), allocatable :: edge_norm_vectors !< ((x,y), edge, i, j); normal direction vector of each face - - integer(ik) :: host_image_id = 0 - integer(ik) :: n_halo_cells = 0 - integer(ik) :: total_cells = 0 - integer(ik) :: total_cells_halo = 0 - integer(ik) :: ni_cells = 0 - integer(ik) :: nj_cells = 0 - integer(ik) :: ni_cells_halo = 0 - integer(ik) :: nj_cells_halo = 0 - - ! All quantities are with respect to the current local block, unless 'global' is specified. - ! The 'halo' suffix means that halo/ghost/bc cells are included - integer(ik) :: ilo_cell = 0 !< cell low i index (excluding halo) - integer(ik) :: jlo_cell = 0 !< cell low j index (excluding halo) - integer(ik) :: ihi_cell = 0 !< cell high i index (excluding halo) - integer(ik) :: jhi_cell = 0 !< cell high j index (excluding halo) - integer(ik) :: ilo_cell_halo = 0 !< cell low i index (including halo) - integer(ik) :: jlo_cell_halo = 0 !< cell low j index (including halo) - integer(ik) :: ihi_cell_halo = 0 !< cell high i index (including halo) - integer(ik) :: jhi_cell_halo = 0 !< cell high j index (including halo) - integer(ik) :: ilo_cell_global = 0 !< cell low i index (excluding halo) in the global scope - integer(ik) :: jlo_cell_global = 0 !< cell low j index (excluding halo) in the global scope - integer(ik) :: ihi_cell_global = 0 !< cell high i index (excluding halo) in the global scope - integer(ik) :: jhi_cell_global = 0 !< cell high j index (excluding halo) in the global scope - - integer(ik), dimension(2) :: cell_dim_global = 0 - !< (i, j); number of cells in each direction in the global scope (excluding halo) - - integer(ik) :: ilo_node = 0 - integer(ik) :: jlo_node = 0 - integer(ik) :: ihi_node = 0 - integer(ik) :: jhi_node = 0 - integer(ik) :: ilo_node_halo = 0 - integer(ik) :: jlo_node_halo = 0 - integer(ik) :: ihi_node_halo = 0 - integer(ik) :: jhi_node_halo = 0 - - logical :: on_ilo_bc = .false. - logical :: on_ihi_bc = .false. - logical :: on_jlo_bc = .false. - logical :: on_jhi_bc = .false. - - ! Parallel neighbor information - integer(ik), dimension(8) :: neighbors = 0 !< (N, S, E, W, NE, SE, SW, NW); parallel neighbor image indices contains ! Private methods private procedure :: populate_element_specifications procedure :: scale_and_nondimensionalize - + ! Public methods procedure, public :: initialize => init_2d_block procedure, public :: gather @@ -109,72 +64,53 @@ subroutine init_2d_block(self, input) call h5%initialize(trim(input%initial_condition_file), status='old', action='r') call h5%shape('/density', cell_shape) - call h5%read('/n_ghost_layers',self%n_halo_cells) + call h5%read('/n_ghost_layers', self%n_halo_cells) call h5%finalize() ! The grid includes the ghost layer information. We want to tile based on ! the real domain b/c it makes the book-keeping a bit easier. ! There are two sets of ghost/halo/boundary cells for each direction - + ! Get the total # of cells in each direction global_cell_dims = int(cell_shape, ik) ! Remove the halo cell count for now, b/c this makes partitioning easier - global_cell_dims = global_cell_dims - (2*self%n_halo_cells) - self%ilo_cell_global = 1 - self%jlo_cell_global = 1 - self%ihi_cell_global = global_cell_dims(1) - self%jhi_cell_global = global_cell_dims(2) - self%cell_dim_global = global_cell_dims + global_cell_dims = global_cell_dims - (2 * self%n_halo_cells) + + self%global_dims(1:2) = global_cell_dims + self%lbounds_global = [1, 1, 0] + self%ubounds_global = [self%global_dims(1), self%global_dims(2), 0] ! Now partition via the tile_indices function (this splits based on current image number) indices = tile_indices(global_cell_dims) lbounds = indices([1, 3]) ubounds = indices([2, 4]) - self%ilo_cell = lbounds(1) - self%jlo_cell = lbounds(2) - self%ihi_cell = ubounds(1) - self%jhi_cell = ubounds(2) - - self%ilo_cell_halo = self%ilo_cell - self%n_halo_cells - self%jlo_cell_halo = self%jlo_cell - self%n_halo_cells - self%ihi_cell_halo = self%ihi_cell + self%n_halo_cells - self%jhi_cell_halo = self%jhi_cell + self%n_halo_cells - - self%ilo_node = self%ilo_cell - self%jlo_node = self%jlo_cell - self%ihi_node = self%ihi_cell + 1 - self%jhi_node = self%jhi_cell + 1 - - self%ilo_node_halo = self%ilo_node - self%n_halo_cells - self%jlo_node_halo = self%jlo_node - self%n_halo_cells - self%ihi_node_halo = self%ihi_node + self%n_halo_cells - self%jhi_node_halo = self%jhi_node + self%n_halo_cells + self%lbounds(1:2) = lbounds + self%ubounds(1:2) = ubounds + + self%lbounds_halo(1:2) = self%lbounds(1:2) - self%n_halo_cells + self%ubounds_halo(1:2) = self%ubounds(1:2) + self%n_halo_cells self%host_image_id = this_image() - self%ni_cells = self%ihi_cell - self%ilo_cell + 1 - self%nj_cells = self%jhi_cell - self%jlo_cell + 1 - self%ni_cells_halo = self%ihi_cell_halo - self%ilo_cell_halo + 1 - self%nj_cells_halo = self%jhi_cell_halo - self%jlo_cell_halo + 1 - self%total_cells = self%ni_cells * self%nj_cells - self%total_cells_halo = self%ni_cells_halo * self%nj_cells_halo + self%block_dims(1:2) = self%ubounds(1:2) - self%lbounds(1:2) + 1 + self%total_cells = size(self%block_dims) ! Determine if this grid block has an edge on one of the global boundaries - if(self%ilo_cell == 1) self%on_ilo_bc = .true. - if(self%jlo_cell == 1) self%on_jlo_bc = .true. - if(self%ihi_cell == self%ihi_cell_global) self%on_ihi_bc = .true. - if(self%jhi_cell == self%jhi_cell_global) self%on_jhi_bc = .true. - - ! Allocate the node-based arrays - associate(ilo => self%ilo_node_halo, ihi => self%ihi_node_halo, & - jlo => self%jlo_node_halo, jhi => self%jhi_node_halo) + if(self%lbounds(1) == 1) self%on_ilo_bc = .true. + if(self%lbounds(2) == 1) self%on_jlo_bc = .true. + if(self%ubounds(1) == self%ubounds_global(1)) self%on_ihi_bc = .true. + if(self%ubounds(2) == self%ubounds_global(2)) self%on_jhi_bc = .true. + + ! Allocate the node-based arrays (thus the + 1 in the ilo/ihi) + associate(ilo => self%lbounds_halo(1), ihi => self%ubounds_halo(1) + 1, & + jlo => self%lbounds_halo(2), jhi => self%ubounds_halo(2) + 1) allocate(self%node_x(ilo:ihi, jlo:jhi)) !< (i, j); x location of each node allocate(self%node_y(ilo:ihi, jlo:jhi)) !< (i, j); y location of each node end associate ! Allocate all of the cell-based arrays - associate(ilo => self%ilo_cell_halo, ihi => self%ihi_cell_halo, & - jlo => self%jlo_cell_halo, jhi => self%jhi_cell_halo) + associate(ilo => self%lbounds_halo(1), ihi => self%ubounds_halo(1), & + jlo => self%lbounds_halo(2), jhi => self%ubounds_halo(2)) allocate(self%volume(ilo:ihi, jlo:jhi)) !< (i, j); volume of each cell allocate(self%dx(ilo:ihi, jlo:jhi)) !< (i, j); dx spacing of each cell allocate(self%dy(ilo:ihi, jlo:jhi)) !< (i, j); dy spacing of each cell @@ -202,21 +138,19 @@ end subroutine init_2d_block subroutine print_grid_stats(self) class(grid_block_2d_t), intent(in) :: self integer :: ni, nj - if (this_image() == 1) then - ni = self%ihi_cell_global - self%ilo_cell_global + 1 - nj = self%jhi_cell_global - self%jlo_cell_global + 1 + if(this_image() == 1) then write(*, '(a)') "Grid stats:" write(*, '(a)') "=========================" write(*, '(a)') "Blocks" write(*, '(a)') "-------------" - write(*, '(a, i0)') "Number of blocks : ", num_images() - write(*, '(2(a,i0))') "Average block size (cells): ", self%ni_cells, " x " , self%nj_cells - write(*, '(a, i6)') "Number of halo cells : ", self%n_halo_cells + write(*, '(a, i0)') "Number of blocks : ", num_images() + write(*, '(2(a,i0))') "Average block size (cells): ", self%block_dims(1), " x ", self%block_dims(2) + write(*, '(a, i6)') "Number of halo cells : ", self%n_halo_cells write(*, '(a)') "Global" write(*, '(a)') "-------------" - write(*, '(a, i0)') "Number of i cells: ", ni - write(*, '(a, i0)') "Number of j cells: ", nj - write(*, '(a, i0)') "Total cells : ", ni*nj + write(*, '(a, i0)') "Number of i cells: ", self%global_dims(1) + write(*, '(a, i0)') "Number of j cells: ", self%global_dims(2) + write(*, '(a, i0)') "Total cells : ", size(self%global_dims) end if end subroutine print_grid_stats @@ -231,7 +165,7 @@ subroutine finalize_2d_block(self) if(allocated(self%centroid_x)) deallocate(self%centroid_x) if(allocated(self%centroid_y)) deallocate(self%centroid_y) if(allocated(self%edge_lengths)) deallocate(self%edge_lengths) - if(allocated(self%edge_norm_vectors)) deallocate(self%edge_norm_vectors) + if(allocated(self%edge_norm_vectors)) deallocate(self%edge_norm_vectors) end subroutine ! -------------------------------------------------------------------- @@ -251,7 +185,7 @@ subroutine read_from_h5(self, input) character(:), allocatable :: filename character(32) :: str_buff = '' character(300) :: msg = '' - real(rk), allocatable, dimension(:,:) :: x + real(rk), allocatable, dimension(:, :) :: x integer(hsize_t), allocatable :: dims(:) @@ -286,16 +220,16 @@ subroutine read_from_h5(self, input) ! I couldn't get the slice read to work, so we read the whole array ! in and extract the slice later - call h5%shape('/x',dims) + call h5%shape('/x', dims) allocate(x(dims(1), dims(2))) - - call h5%read(dname='/x', value=x) - self%node_x = x(self%ilo_cell_halo:self%ihi_cell_halo, & - self%jlo_cell_halo:self%jhi_cell_halo) - - call h5%read(dname='/y', value=x) - self%node_y = x(self%ilo_cell_halo:self%ihi_cell_halo, & - self%jlo_cell_halo:self%jhi_cell_halo) + + call h5%read(dname='/x', value=x) + self%node_x = x(self%lbounds_halo(1) + self%n_halo_cells:self%ubounds_halo(1) + self%n_halo_cells + 1, & + self%lbounds_halo(2) + self%n_halo_cells:self%ubounds_halo(2) + self%n_halo_cells + 1) + + call h5%read(dname='/y', value=x) + self%node_y = x(self%lbounds_halo(1) + self%n_halo_cells:self%ubounds_halo(1) + self%n_halo_cells + 1, & + self%lbounds_halo(2) + self%n_halo_cells:self%ubounds_halo(2) + self%n_halo_cells + 1) if(input%restart_from_file) then call h5%readattr('/x', 'units', str_buff) @@ -345,11 +279,12 @@ function gather(self, var, image) select case(var) case('x', 'y') - ni = self%ihi_cell_global - self%ilo_cell_global + 2 ! the +2 is b/c of nodes vs cells - nj = self%jhi_cell_global - self%jlo_cell_global + 2 ! the +2 is b/c of nodes vs cells + ni = self%global_dims(1) + 1 + (self%n_halo_cells * 2) ! the +1 is b/c of nodes vs cells + nj = self%global_dims(2) + 1 + (self%n_halo_cells * 2) ! the +1 is b/c of nodes vs cells allocate(gather_coarray(ni, nj)[*]) - associate(ilo => self%ilo_node_halo, ihi => self%ihi_node_halo, & - jlo => self%jlo_node_halo, jhi => self%jhi_node_halo) + ! Allocate the node-based arrays (thus the + 1 in the ilo/ihi) + associate(ilo => self%lbounds_halo(1), ihi => self%ubounds_halo(1) + 1, & + jlo => self%lbounds_halo(2), jhi => self%ubounds_halo(2) + 1) select case(var) case('x') @@ -364,11 +299,11 @@ function gather(self, var, image) case('volume') - ni = self%ihi_cell_global - self%ilo_cell_global + 1 - nj = self%jhi_cell_global - self%jlo_cell_global + 1 + ni = self%global_dims(1) + (self%n_halo_cells * 2) + nj = self%global_dims(2) + (self%n_halo_cells * 2) allocate(gather_coarray(ni, nj)[*]) - associate(ilo => self%ilo_cell_halo, ihi => self%ihi_cell_halo, & - jlo => self%jlo_cell_halo, jhi => self%jhi_cell_halo) + associate(ilo => self%lbounds_halo(1), ihi => self%ubounds_halo(1), & + jlo => self%lbounds_halo(2), jhi => self%ubounds_halo(2)) gather_coarray(ilo:ihi, jlo:jhi)[image] = self%volume(ilo:ihi, jlo:jhi) sync all if(this_image() == image) gather = gather_coarray @@ -396,8 +331,8 @@ subroutine populate_element_specifications(self) x_coords = 0.0_rk y_coords = 0.0_rk - do j = self%jlo_cell_halo, self%jhi_cell_halo - do i = self%ilo_cell_halo, self%ihi_cell_halo + do j = self%lbounds_halo(2), self%ubounds_halo(2) + do i = self%lbounds_halo(1), self%ubounds_halo(1) associate(x => self%node_x, y => self%node_y) x_coords = [x(i, j), x(i + 1, j), x(i + 1, j + 1), x(i, j + 1)] y_coords = [y(i, j), y(i + 1, j), y(i + 1, j + 1), y(i, j + 1)] diff --git a/src/lib/grid/grid_block_3d.f90 b/src/lib/grid/grid_block_3d.f90 index a89ff5d..7ce6e15 100644 --- a/src/lib/grid/grid_block_3d.f90 +++ b/src/lib/grid/grid_block_3d.f90 @@ -36,32 +36,32 @@ subroutine init_3d_block(self, input) class(grid_block_3d_t), intent(inout) :: self class(input_t), intent(in) :: input - allocate(self%global_node_dims (3)) - self%global_node_dims = 0 - allocate(self%global_cell_dims (3)) - self%global_cell_dims = 0 + allocate(self%global_node_dims(3)) + self%global_node_dims = 0 + allocate(self%global_cell_dims(3)) + self%global_cell_dims = 0 allocate(self%domain_node_shape(3)) self%domain_node_shape = 0 allocate(self%domain_cell_shape(3)) self%domain_cell_shape = 0 - allocate(self%node_lbounds (3)) - self%node_lbounds = 0 - allocate(self%node_ubounds (3)) - self%node_ubounds = 0 + allocate(self%node_lbounds(3)) + self%node_lbounds = 0 + allocate(self%node_ubounds(3)) + self%node_ubounds = 0 allocate(self%node_lbounds_halo(3)) self%node_lbounds_halo = 0 allocate(self%node_ubounds_halo(3)) self%node_ubounds_halo = 0 - allocate(self%cell_lbounds (3)) - self%cell_lbounds = 0 - allocate(self%cell_ubounds (3)) - self%cell_ubounds = 0 + allocate(self%cell_lbounds(3)) + self%cell_lbounds = 0 + allocate(self%cell_ubounds(3)) + self%cell_ubounds = 0 allocate(self%cell_lbounds_halo(3)) self%cell_lbounds_halo = 0 allocate(self%cell_ubounds_halo(3)) self%cell_ubounds_halo = 0 - allocate(self%on_bc (2,3)) - self%on_bc = .false. + allocate(self%on_bc(2, 3)) + self%on_bc = .false. end subroutine @@ -76,16 +76,16 @@ subroutine finalize_3d_block(self) if(allocated(self%centroid_y)) deallocate(self%centroid_y) if(allocated(self%edge_lengths)) deallocate(self%edge_lengths) if(allocated(self%edge_norm_vectors)) deallocate(self%edge_norm_vectors) - if(allocated(self%global_node_dims )) deallocate(self%global_node_dims ) - if(allocated(self%global_cell_dims )) deallocate(self%global_cell_dims ) + if(allocated(self%global_node_dims)) deallocate(self%global_node_dims) + if(allocated(self%global_cell_dims)) deallocate(self%global_cell_dims) if(allocated(self%domain_node_shape)) deallocate(self%domain_node_shape) if(allocated(self%domain_cell_shape)) deallocate(self%domain_cell_shape) - if(allocated(self%node_lbounds )) deallocate(self%node_lbounds ) - if(allocated(self%node_ubounds )) deallocate(self%node_ubounds ) + if(allocated(self%node_lbounds)) deallocate(self%node_lbounds) + if(allocated(self%node_ubounds)) deallocate(self%node_ubounds) if(allocated(self%node_lbounds_halo)) deallocate(self%node_lbounds_halo) if(allocated(self%node_ubounds_halo)) deallocate(self%node_ubounds_halo) - if(allocated(self%cell_lbounds )) deallocate(self%cell_lbounds ) - if(allocated(self%cell_ubounds )) deallocate(self%cell_ubounds ) + if(allocated(self%cell_lbounds)) deallocate(self%cell_lbounds) + if(allocated(self%cell_ubounds)) deallocate(self%cell_ubounds) if(allocated(self%cell_lbounds_halo)) deallocate(self%cell_lbounds_halo) if(allocated(self%cell_ubounds_halo)) deallocate(self%cell_ubounds_halo) if(allocated(self%on_bc)) deallocate(self%on_bc) diff --git a/src/lib/grid/grid_factory.f90 b/src/lib/grid/grid_factory.f90 index 66e4997..e7eacc0 100644 --- a/src/lib/grid/grid_factory.f90 +++ b/src/lib/grid/grid_factory.f90 @@ -24,9 +24,9 @@ module mod_grid_factory use mod_globals, only: debug_print, set_domain_dimensionality use mod_input, only: input_t use mod_grid_block, only: grid_block_t - use mod_grid_block_1d, only: new_1d_grid_block + ! use mod_grid_block_1d, only: new_1d_grid_block use mod_grid_block_2d, only: new_2d_grid_block - use mod_grid_block_3d, only: new_3d_grid_block + ! use mod_grid_block_3d, only: new_3d_grid_block implicit none private @@ -41,18 +41,18 @@ function grid_factory(input) result(grid) call debug_print('Creating a grid in grid_factory', __FILE__, __LINE__) - ! select case(trim(input%grid_type)) - ! case('X') - ! grid => new_1d_grid_block(input) - ! case('XY') - ! grid => new_2d_grid_block(input) - ! ! call set_domain_dimensionality(dimensionality='2D_XY', & - ! ! grid_orthogonality=.true., num_ghost_layers=grid%n_halo_cells) - ! ! case('RZ') - ! case('XYZ') - ! grid => new_3d_grid_block(input) - ! case default - ! error stop 'Unsupported grid type in grid_factory' - ! end select + select case(trim(input%grid_type)) + ! case('X') + ! grid => new_1d_grid_block(input) + case('XY') + grid => new_2d_grid_block(input) + ! call set_domain_dimensionality(dimensionality='2D_XY', & + ! grid_orthogonality=.true., num_ghost_layers=grid%n_halo_cells) + ! case('RZ') + ! case('XYZ') + ! grid => new_3d_grid_block(input) + case default + error stop 'Unsupported grid type in grid_factory' + end select end function grid_factory end module mod_grid_factory diff --git a/src/lib/io/contour_writer.f90 b/src/lib/io/contour_writer.f90 index 540b32a..6949afb 100644 --- a/src/lib/io/contour_writer.f90 +++ b/src/lib/io/contour_writer.f90 @@ -24,6 +24,9 @@ module mod_contour_writer use mod_master_puppeteer, only: master_puppeteer_t use mod_globals, only: n_ghost_layers use mod_fluid, only: fluid_t + use mod_grid_block, only: grid_block_t + use mod_grid_block_2d, only: grid_block_2d_t + use mod_error, only: error_msg use mod_units use mod_nondimensionalization, only: rho_0, v_0, p_0, t_0, l_0, e_0 use mod_eos, only: eos @@ -116,12 +119,14 @@ type(contour_writer_t) function constructor(input) result(writer) writer%plot_ghost_cells = input%plot_ghost_cells writer%plot_64bit = input%plot_64bit - write(*, '(a)') "Contour Writer Info" - write(*, '(a)') "===================" - write(*, '(a, a)') 'Output folder: ', writer%results_folder - write(*, '(a, l1)') 'plot_64bit: ', writer%plot_64bit - write(*, '(a, l1)') 'plot_ghost_cells: ', writer%plot_ghost_cells - write(*, *) + if(this_image() == 1) then + write(*, '(a)') "Contour Writer Info" + write(*, '(a)') "===================" + write(*, '(a, a)') 'Output folder: ', writer%results_folder + write(*, '(a, l1)') 'plot_64bit: ', writer%plot_64bit + write(*, '(a, l1)') 'plot_ghost_cells: ', writer%plot_ghost_cells + write(*, *) + endif end function @@ -144,28 +149,28 @@ subroutine write_contour(self, master, time, iteration) self%xdmf_filename = trim(char_buff)//'.xdmf' if(self%plot_ghost_cells) then - self%ilo_cell = master%grid%cell_lbounds_halo(1) - self%ihi_cell = master%grid%cell_ubounds_halo(1) - self%jlo_cell = master%grid%cell_lbounds_halo(2) - self%jhi_cell = master%grid%cell_ubounds_halo(2) - - self%ilo_node = master%grid%node_lbounds_halo(1) - self%ihi_node = master%grid%node_ubounds_halo(1) - self%jlo_node = master%grid%node_lbounds_halo(2) - self%jhi_node = master%grid%node_ubounds_halo(2) + self%ilo_cell = master%grid%lbounds_halo(1) + self%ihi_cell = master%grid%ubounds_halo(1) + self%jlo_cell = master%grid%lbounds_halo(2) + self%jhi_cell = master%grid%ubounds_halo(2) + + self%ilo_node = master%grid%lbounds_halo(1) + self%ihi_node = master%grid%ubounds_halo(1) + 1 + self%jlo_node = master%grid%lbounds_halo(2) + self%jhi_node = master%grid%ubounds_halo(2) + 1 else - self%ilo_cell = master%grid%cell_lbounds(1) - self%ihi_cell = master%grid%cell_ubounds(1) - self%jlo_cell = master%grid%cell_lbounds(2) - self%jhi_cell = master%grid%cell_ubounds(2) - - self%ilo_node = master%grid%node_lbounds(1) - self%ihi_node = master%grid%node_ubounds(1) - self%jlo_node = master%grid%node_lbounds(2) - self%jhi_node = master%grid%node_ubounds(2) + self%ilo_cell = master%grid%lbounds(1) + self%ihi_cell = master%grid%ubounds(1) + self%jlo_cell = master%grid%lbounds(2) + self%jhi_cell = master%grid%ubounds(2) + + self%ilo_node = master%grid%lbounds(1) + self%ihi_node = master%grid%ubounds(1) + 1 + self%jlo_node = master%grid%lbounds(2) + self%jhi_node = master%grid%ubounds(2) + 1 end if - write(*, '(a,a)') "Saving contour file: "//self%hdf5_filename + if (this_image() == 1) write(*, '(a,a)') "Saving contour file: "//self%hdf5_filename select case(self%format) case('xdmf') call self%write_hdf5(master, time, iteration) @@ -173,8 +178,9 @@ subroutine write_contour(self, master, time, iteration) case('hdf5', 'h5') call self%write_hdf5(master, time, iteration) case default - print *, 'Contour format:', self%format - error stop "Unsupported I/O contour format" + call error_msg(module_name='mod_contour_writer', class_name='contour_writer_t', procedure_name='write_contour', & + message="Unsupported I/O contour format'"// self%format // "' (must be .xdmf .h5 or .hdf5)", & + file_name=__FILE__, line_number=__LINE__) end select end subroutine write_contour @@ -194,9 +200,9 @@ subroutine write_hdf5(self, master, time, iteration) time_w_dims = time * io_time_units * t_0 delta_t_w_dims = master%dt * t_0 - if (this_image() == 1) then + if(this_image() == 1) then call self%hdf5_file%initialize(filename=self%results_folder//'/'//self%hdf5_filename, & - status='new', action='w', comp_lvl=self%compression_level) + status='new', action='w', comp_lvl=self%compression_level) ! Header info call self%hdf5_file%add('/title', master%title) @@ -328,7 +334,7 @@ subroutine write_hdf5(self, master, time, iteration) ! if(allocated(int_data_buffer)) deallocate(int_data_buffer) ! if(allocated(io_data_buffer)) deallocate(io_data_buffer) - if (this_image() == 1) call self%hdf5_file%finalize() + if(this_image() == 1) call self%hdf5_file%finalize() end subroutine write_hdf5 subroutine write_xdmf(self, master, time, iteration) diff --git a/src/lib/io/mod_input.f90 b/src/lib/io/mod_input.f90 index 7435556..c83f626 100644 --- a/src/lib/io/mod_input.f90 +++ b/src/lib/io/mod_input.f90 @@ -396,7 +396,7 @@ end subroutine read_from_ini subroutine display_config(self) class(input_t), intent(in) :: self - if (this_image() == 1) then + if(this_image() == 1) then write(*, '(a)') "Input settings:" write(*, '(a)') "===============" write(*, *) diff --git a/src/lib/master_puppeteer/master_puppeteer.f90 b/src/lib/master_puppeteer/master_puppeteer.f90 index 7d9ea10..7b5eaa5 100644 --- a/src/lib/master_puppeteer/master_puppeteer.f90 +++ b/src/lib/master_puppeteer/master_puppeteer.f90 @@ -31,9 +31,9 @@ module mod_master_puppeteer use mod_input, only: input_t use mod_eos, only: eos use mod_grid_block, only: grid_block_t - use mod_grid_block_1d, only: grid_block_1d_t + ! use mod_grid_block_1d, only: grid_block_1d_t use mod_grid_block_2d, only: grid_block_2d_t - use mod_grid_block_3d, only: grid_block_3d_t + ! use mod_grid_block_3d, only: grid_block_3d_t use mod_grid_factory, only: grid_factory use hdf5_interface, only: hdf5_file use mod_nondimensionalization, only: set_scale_factors @@ -146,29 +146,33 @@ subroutine integrate(self, dt, error_code) real(rk) :: max_dt - if (present(dt)) then + if(present(dt)) then self%dt = dt else self%dt = 0.0_rk end if self%iteration = self%iteration + 1 - + ! Get the maximum allowable timestep from each physics package - max_dt = self%fluid%get_timestep(self%grid) - if (self%dt > max_dt) then - write(*, '(a,i0,a,2(es16.6,a))') "Warning on image ",this_image(), & - ", the input dt (",self%dt,& - ") is larger than the max allowable dt (", max_dt, & - ") based on the CFL condition" + select type(grid => self%grid) + class is(grid_block_2d_t) + max_dt = self%fluid%get_timestep(self%grid) + end select + + if(self%dt > max_dt) then + write(*, '(a,i0,a,2(es16.6,a))') "Warning on image ", this_image(), & + ", the input dt (", self%dt, & + ") is larger than the max allowable dt (", max_dt, & + ") based on the CFL condition" else self%dt = max_dt - end if + end if self%time = self%time + self%dt select type(grid => self%grid) - class is (grid_block_2d_t) + class is(grid_block_2d_t) call self%fluid%integrate(dt=self%dt, grid=self%grid, error_code=error_code) end select diff --git a/src/lib/non_dimensional.f90 b/src/lib/non_dimensional.f90 index 1d89443..cf7c263 100644 --- a/src/lib/non_dimensional.f90 +++ b/src/lib/non_dimensional.f90 @@ -81,17 +81,19 @@ subroutine set_scale_factors(density_scale, pressure_scale) e_0 = p_0 t_0 = l_0 / v_0 - print * - write(*, '(a)') "Scale factors (CATO is non-dimensional)" - write(*, '(a)') "========================================" - write(*, '(a, es10.3)') "Time scale factor (t_0): ", t_0 - write(*, '(a, es10.3)') "Density scale factor (rho_0): ", rho_0 - write(*, '(a, es10.3)') "Length scale factor (l_0): ", l_0 - write(*, '(a, es10.3)') "Pressure scale factor (p_0): ", p_0 - write(*, '(a, es10.3)') "Velocity scale factor (v_0): ", v_0 - write(*, '(a, es10.3)') "Energy scale factor (e_0): ", e_0 - write(*, '(a)') "========================================" - print * + if (this_image() == 1) then + print * + write(*, '(a)') "Scale factors (CATO is non-dimensional)" + write(*, '(a)') "========================================" + write(*, '(a, es10.3)') "Time scale factor (t_0): ", t_0 + write(*, '(a, es10.3)') "Density scale factor (rho_0): ", rho_0 + write(*, '(a, es10.3)') "Length scale factor (l_0): ", l_0 + write(*, '(a, es10.3)') "Pressure scale factor (p_0): ", p_0 + write(*, '(a, es10.3)') "Velocity scale factor (v_0): ", v_0 + write(*, '(a, es10.3)') "Energy scale factor (e_0): ", e_0 + write(*, '(a)') "========================================" + print * + end if scale_factors_set = .true. diff --git a/src/main.f90 b/src/main.f90 index 9908e5f..2258eb5 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -55,9 +55,8 @@ program cato open(std_error, file='std.err') - ! ascii art for the heck of it :) - if (this_image() == 1) then + if(this_image() == 1) then write(std_out, '(a)') write(std_out, '(a)') " ______ ___ .___________. ______ " write(std_out, '(a)') " / | / \ | | / __ \ " @@ -73,7 +72,6 @@ program cato write(std_out, '(a, i0, a)') "Running with ", omp_get_max_threads(), " OpenMP threads" #endif /* USE_OPENMP */ - call open_debug_files() call print_version_stats() call get_command_argument(1, command_line_arg) @@ -108,7 +106,7 @@ program cato call contour_writer%write_contour(master, time, iteration) end if - if (this_image() == 1) then + if(this_image() == 1) then print * write(std_out, '(a)') '--------------------------------------------' write(std_out, '(a)') ' Starting time loop:' @@ -130,19 +128,19 @@ program cato call master%set_time(time, iteration) - if (this_image() == 1) then + if(this_image() == 1) then write(std_out, '(a, es10.3)') "Starting time:", time write(std_out, '(a, es10.3)') "Starting timestep:", delta_t endif do while(time < max_time .and. iteration < input%max_iterations) - if (this_image() == 1) then + if(this_image() == 1) then write(std_out, '(2(a, es10.3), a, i0)') 'Time =', time * io_time_units * t_0, & ' '//trim(io_time_label)//', Delta t =', delta_t * t_0, ' s, Iteration: ', iteration endif - + ! Integrate in time - if (input%use_constant_delta_t) then + if(input%use_constant_delta_t) then call master%integrate(error_code=error_code, dt=input%initial_delta_t / t_0) else call master%integrate(error_code=error_code) @@ -160,7 +158,7 @@ program cato ! I/O if(abs(time - next_output_time) < epsilon(1.0_rk)) then next_output_time = next_output_time + contour_interval_dt - if (this_image() == 1) then + if(this_image() == 1) then write(std_out, '(a, es10.3, a)') 'Saving Contour, Next Output Time: ', & next_output_time * t_0 * io_time_units, ' '//trim(io_time_label) endif diff --git a/src/third_party/h5fortran/src/reader.f90 b/src/third_party/h5fortran/src/reader.f90 index dbda93e..18563e6 100644 --- a/src/third_party/h5fortran/src/reader.f90 +++ b/src/third_party/h5fortran/src/reader.f90 @@ -134,50 +134,50 @@ end procedure hdf_read_1d module procedure hdf_read_2d - integer(HSIZE_T) :: dims(rank(value)) - integer(HID_T) :: did, sid, mem_sid - integer :: ier +integer(HSIZE_T) :: dims(rank(value)) +integer(HID_T) :: did, sid, mem_sid +integer :: ier - did = 0 !< sentinel - sid = H5S_ALL_F - mem_sid = H5S_ALL_F - dims = shape(value) +did = 0 !< sentinel +sid = H5S_ALL_F +mem_sid = H5S_ALL_F +dims = shape(value) - if(present(istart) .and. present(iend)) then - if(present(stride)) then +if(present(istart) .and. present(iend)) then + if(present(stride)) then !! necessary to use this present check for Intel and GCC - call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride) - else - call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend) - endif + call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend, stride) else - call hdf_shape_check(self, dname, dims) - call h5dopen_f(self%lid, dname, did, ier) - if(ier /= 0) then - write(stderr, *) 'h5fortran:ERROR:reader: could not setup read ', dname, ' from ', self%filename - error stop - endif + call hdf_get_slice(self, dname, did, sid, mem_sid, istart, iend) endif - - select type(value) - type is(real(real64)) - call h5dread_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid) - type is(real(real32)) - call h5dread_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid) - type is(integer(int32)) - call h5dread_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid) - class default - error stop 'h5fortran:reader: incorrect type' - end select +else + call hdf_shape_check(self, dname, dims) + call h5dopen_f(self%lid, dname, did, ier) if(ier /= 0) then - write(stderr, *) 'h5fortran:ERROR:reader: could not read ', dname, ' from ', self%filename + write(stderr, *) 'h5fortran:ERROR:reader: could not setup read ', dname, ' from ', self%filename error stop endif +endif - call hdf_wrapup(did, sid, ier) +select type(value) +type is(real(real64)) + call h5dread_f(did, H5T_NATIVE_DOUBLE, value, dims, ier, mem_sid, sid) +type is(real(real32)) + call h5dread_f(did, H5T_NATIVE_REAL, value, dims, ier, mem_sid, sid) +type is(integer(int32)) + call h5dread_f(did, H5T_NATIVE_INTEGER, value, dims, ier, mem_sid, sid) +class default + error stop 'h5fortran:reader: incorrect type' +end select +if(ier /= 0) then + write(stderr, *) 'h5fortran:ERROR:reader: could not read ', dname, ' from ', self%filename + error stop +endif - if(present(ierr)) ierr = ier - if(check(ier, self%filename, dname) .and. .not. present(ierr)) error stop +call hdf_wrapup(did, sid, ier) + +if(present(ierr)) ierr = ier +if(check(ier, self%filename, dname) .and. .not. present(ierr)) error stop end procedure hdf_read_2d diff --git a/tests/unit/test_bc/test_bc.pf b/tests/unit/test_bc/test_bc.pf index 0e320b2..2621737 100644 --- a/tests/unit/test_bc/test_bc.pf +++ b/tests/unit/test_bc/test_bc.pf @@ -61,8 +61,8 @@ contains print *, "Setting up the periodic grid" ! These are normally handled by the master puppeteer, but for now we make them ourselves - associate(left => grid%cell_lbounds_halo(1), right => grid%cell_ubounds_halo(1), & - bottom => grid%cell_lbounds_halo(2), top => grid%cell_ubounds_halo(2)) + associate(left => grid%lbounds_halo(1), right => grid%ubounds_halo(1), & + bottom => grid%lbounds_halo(2), top => grid%ubounds_halo(2)) if(allocated(rho)) deallocate(rho) allocate(rho(left:right, bottom:top), stat=alloc_status) @@ -85,7 +85,7 @@ contains ! Domain cartoon layout ! |---|---||---|---|---|---||---|---| - ! 6 | 0 | 6 || 8 | 0 | 0 | 6 || 8 | 0 | <- j=6; aka top_ghost (grid%cell_ubounds_halo(2)) + ! 6 | 0 | 6 || 8 | 0 | 0 | 6 || 8 | 0 | <- j=6; aka top_ghost (grid%ubounds_halo(2)) ! |---|---||---|---|---|---||---|---| ! 5 | 5 | 2 || 1 | 5 | 5 | 2 || 1 | 5 | ! |===|===||===|===|===|===||===|===| @@ -99,14 +99,14 @@ contains ! |===|===||===|===|===|===||===|===| ! 0 | 7 | 3 || 4 | 7 | 7 | 3 || 4 | 7 | ! |---|---||---|---|---|---||---|---| - ! -1 | 0 | 6 || 8 | 0 | 0 | 6 || 8 | 0 | <- j=-1; aka bottom_ghost (grid%cell_lbounds_halo(2)) + ! -1 | 0 | 6 || 8 | 0 | 0 | 6 || 8 | 0 | <- j=-1; aka bottom_ghost (grid%lbounds_halo(2)) ! |---|---||---|---|---|---||---|---| ! -1 0 1 2 3 4 5 6 - ! i0; aka left_ghost, grid%cell_lbounds_halo(1) + ! i0; aka left_ghost, grid%lbounds_halo(1) ! i1; aka left, grid%cell_lbounds(1) ! i4; aka right, grid%cell_ubounds(1) - ! i5; aka right_ghost grid%cell_ubounds_halo(1) + ! i5; aka right_ghost grid%ubounds_halo(1) ! Test the conserved var state rho = 0.0_rk @@ -127,8 +127,8 @@ contains associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) print *, "Initial conditions:" do j = top_ghost, bottom_ghost, -1 write(*, '(i3, a, 10(f6.0))') j, " : ", rho(:, j) @@ -162,8 +162,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" @@ -211,8 +211,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" @@ -260,8 +260,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" @@ -309,8 +309,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" @@ -369,8 +369,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" @@ -436,15 +436,15 @@ contains real(rk), dimension(4) :: right_sym_col = 0.0_rk real(rk), dimension(4) :: left_sym_col = 0.0_rk - bottom_ghost = grid%cell_lbounds_halo(2) - top_ghost = grid%cell_ubounds_halo(2) + bottom_ghost = grid%lbounds_halo(2) + top_ghost = grid%ubounds_halo(2) - left_ghost = grid%cell_lbounds_halo(1) - right_ghost = grid%cell_ubounds_halo(1) + left_ghost = grid%lbounds_halo(1) + right_ghost = grid%ubounds_halo(1) ! These are normally handled by the master puppeteer, but for now we make them ourselves - associate(left => grid%cell_lbounds_halo(1), right => grid%cell_ubounds_halo(1), & - bottom => grid%cell_lbounds_halo(2), top => grid%cell_ubounds_halo(2)) + associate(left => grid%lbounds_halo(1), right => grid%ubounds_halo(1), & + bottom => grid%lbounds_halo(2), top => grid%ubounds_halo(2)) if(allocated(rho)) deallocate(rho) allocate(rho(left:right, bottom:top), stat=alloc_status) @@ -483,7 +483,7 @@ contains ! Domain cartoon layout ! |---|---||---|---|---|---||---|---| - ! | - | - || 8 | 0 | 0 | 6 || - | - | <- j=6; aka top_ghost (grid%cell_ubounds_halo(2)) + ! | - | - || 8 | 0 | 0 | 6 || - | - | <- j=6; aka top_ghost (grid%ubounds_halo(2)) ! |---|---||---|---|---|---||---|---| ! | - | - || 4 | 7 | 7 | 3 || - | - | ! |===|===||===|===|===|===||===|===| @@ -497,14 +497,14 @@ contains ! |===|===||===|===|===|===||===|===| ! | - | - || 1 | 5 | 5 | 2 || - | - | ! |---|---||---|---|---|---||---|---| - ! | - | - || 8 | 0 | 0 | 6 || - | - | <- j=-1; aka bottom_ghost (grid%cell_lbounds_halo(2)) + ! | - | - || 8 | 0 | 0 | 6 || - | - | <- j=-1; aka bottom_ghost (grid%lbounds_halo(2)) ! |---|---||---|---|---|---||---|---| ! i0 i1 i2 i3 i4 i5 - ! i0; aka left_ghost, grid%cell_lbounds_halo(1) + ! i0; aka left_ghost, grid%lbounds_halo(1) ! i1; aka left, grid%cell_lbounds(1) ! i4; aka right, grid%cell_ubounds(1) - ! i5; aka right_ghost grid%cell_ubounds_halo(1) + ! i5; aka right_ghost grid%ubounds_halo(1) ! Set unique values along the (real) edges so I can test the bc routines associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & @@ -530,8 +530,8 @@ contains ! Now test the results associate(left => grid%cell_lbounds(1), right => grid%cell_ubounds(1), & bottom => grid%cell_lbounds(2), top => grid%cell_ubounds(2), & - left_ghost => grid%cell_lbounds_halo(1), right_ghost => grid%cell_ubounds_halo(1), & - bottom_ghost => grid%cell_lbounds_halo(2), top_ghost => grid%cell_ubounds_halo(2)) + left_ghost => grid%lbounds_halo(1), right_ghost => grid%ubounds_halo(1), & + bottom_ghost => grid%lbounds_halo(2), top_ghost => grid%ubounds_halo(2)) ! Print out in all it's glory print *, "After BC's applied" diff --git a/tests/unit/test_geometry/test_grid.pf b/tests/unit/test_geometry/test_grid.pf index 70dc6490..08e4fe8 100644 --- a/tests/unit/test_geometry/test_grid.pf +++ b/tests/unit/test_geometry/test_grid.pf @@ -210,13 +210,13 @@ contains @assertEqual(-1, grid%ilo_bc_node) @assertEqual(-1, grid%jlo_bc_node) - @assertEqual(-1, grid%cell_lbounds_halo(1)) - @assertEqual(-1, grid%cell_lbounds_halo(2)) + @assertEqual(-1, grid%lbounds_halo(1)) + @assertEqual(-1, grid%lbounds_halo(2)) @assertEqual(7, grid%ihi_bc_node) @assertEqual(5, grid%jhi_bc_node) - @assertEqual(6, grid%cell_ubounds_halo(1)) - @assertEqual(4, grid%cell_ubounds_halo(2)) + @assertEqual(6, grid%ubounds_halo(1)) + @assertEqual(4, grid%ubounds_halo(2)) @assertEqual(1.0_rk, grid%min_dx) @assertEqual(2.0_rk, grid%min_dy) @@ -232,8 +232,8 @@ contains @assertEqual([-3.0_rk, 3.0_rk],[grid%get_x(0, 0), grid%get_x(6, 0)]) @assertEqual(2.0_rk, grid%get_cell_volumes(1, 1)) - @assertEqual(2.0_rk, grid%get_cell_volumes(grid%cell_lbounds_halo(1), grid%cell_lbounds_halo(2))) - @assertEqual(2.0_rk, grid%get_cell_volumes(grid%cell_ubounds_halo(1), grid%cell_ubounds_halo(2))) + @assertEqual(2.0_rk, grid%get_cell_volumes(grid%lbounds_halo(1), grid%lbounds_halo(2))) + @assertEqual(2.0_rk, grid%get_cell_volumes(grid%ubounds_halo(1), grid%ubounds_halo(2))) i = 1; j = 1 ! Corners