From 22f71ac94e898afa4785f3703e7b7c66c8ad748d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 2 Apr 2014 15:25:08 -0300 Subject: [PATCH 01/41] DRIVER: Set the maximum number of iterations to a huge number. Signed-off-by: Grzegorz Kowal --- src/driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/driver.F90 b/src/driver.F90 index a37a866..4f17446 100644 --- a/src/driver.F90 +++ b/src/driver.F90 @@ -81,7 +81,7 @@ program amun ! integer, dimension(3) :: div = 1 logical, dimension(3) :: per = .true. - integer :: nmax = 0, ndat = 1 + integer :: nmax = huge(1), ndat = 1 real :: tmax = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01 real(kind=8) :: dtnext = 0.0d+00 From f9529def8c98a42d58a5f146c3cd70d0e7c9560f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 2 Apr 2014 16:17:30 -0300 Subject: [PATCH 02/41] PROBLEMS: Implement implosion test problem. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 225 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 225 insertions(+) diff --git a/src/problems.F90 b/src/problems.F90 index 91f925c..ef33ae7 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -131,6 +131,9 @@ module problems case("blast") setup_problem => setup_problem_blast + case("implosion") + setup_problem => setup_problem_implosion + case("kh", "kelvinhelmholtz", "kelvin-helmholtz") setup_problem => setup_problem_kelvin_helmholtz @@ -507,6 +510,228 @@ module problems ! !=============================================================================== ! +! subroutine SETUP_PROBLEM_IMPLOSION: +! ---------------------------------- +! +! Subroutine sets the initial conditions for the implosion problem. +! +! Arguments: +! +! pdata - pointer to the datablock structure of the currently initialized +! block; +! +!=============================================================================== +! + subroutine setup_problem_implosion(pdata) + +! include external procedures and variables +! + use blocks , only : block_data, ndims + use constants , only : d2r + use coordinates, only : im, jm, km + use coordinates, only : ax, ay, az, adx, ady, adz + use equations , only : prim2cons + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use parameters , only : get_parameter_real + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + type(block_data), pointer, intent(inout) :: pdata + +! default parameter values +! + real(kind=8), save :: sline = 1.50d-01 + real(kind=8), save :: adens = 1.00d+00 + real(kind=8), save :: apres = 1.00d+00 + real(kind=8), save :: drat = 1.25d-01 + real(kind=8), save :: prat = 1.40d-01 + real(kind=8), save :: buni = 1.00d+00 + real(kind=8), save :: bgui = 0.00d+00 + real(kind=8), save :: angle = 0.00d+00 + +! local saved parameters +! + logical , save :: first = .true. + real(kind=8), save :: odens = 1.25d-01 + real(kind=8), save :: opres = 1.40d-01 + +! local variables +! + integer :: i, j, k + real(kind=8) :: rl, ru, dx, dy, dz, dxh, dyh, dzh, ds, dl, dr + real(kind=8) :: sn, cs + +! local arrays +! + real(kind=8), dimension(nv,im) :: q, u + real(kind=8), dimension(im) :: x, xl, xu + real(kind=8), dimension(jm) :: y, yl, yu + real(kind=8), dimension(km) :: z, zl, zu +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the problem setup +! + call start_timer(imu) +#endif /* PROFILE */ + +! prepare problem constants during the first subroutine call +! + if (first) then + +! get problem parameters +! + call get_parameter_real("shock_line" , sline ) + call get_parameter_real("ambient_density" , adens ) + call get_parameter_real("ambient_pressure", apres ) + call get_parameter_real("density_ratio" , drat ) + call get_parameter_real("pressure_ratio" , prat ) + call get_parameter_real("buni" , buni ) + call get_parameter_real("bgui" , bgui ) + call get_parameter_real("angle" , angle ) + +! calculate parameters +! + odens = drat * adens + opres = prat * apres + +! reset the first execution flag +! + first = .false. + + end if ! first call + +! prepare block coordinates +! + x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) + y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) +#if NDIMS == 3 + z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km) +#endif /* NDIMS == 3 */ + +! calculate mesh intervals and areas +! + dx = adx(pdata%meta%level) + dy = ady(pdata%meta%level) +#if NDIMS == 3 + dz = adz(pdata%meta%level) +#endif /* NDIMS == 3 */ + dxh = 0.5d+00 * dx + dyh = 0.5d+00 * dy +#if NDIMS == 3 + dzh = 0.5d+00 * dz +#endif /* NDIMS == 3 */ + +! calculate edge coordinates +! + xl(:) = x(:) - dxh + xu(:) = x(:) + dxh + yl(:) = y(:) - dyh + yu(:) = y(:) + dyh +#if NDIMS == 3 + zl(:) = z(:) - dzh + zu(:) = z(:) + dzh +#endif /* NDIMS == 3 */ + +! reset velocity components +! + q(ivx,:) = 0.0d+00 + q(ivy,:) = 0.0d+00 + q(ivz,:) = 0.0d+00 + +! if magnetic field is present, set it to be uniform with the desired strength +! and orientation +! + if (ibx > 0) then + +! calculate the orientation angles +! + sn = sin(d2r * angle) + cs = sqrt(1.0d+00 - sn * sn) + +! set magnetic field components +! + q(ibx,:) = buni * cs + q(iby,:) = buni * sn + q(ibz,:) = bgui + q(ibp,:) = 0.0d+00 + + end if + +! iterate over all positions +! + do k = 1, km + do j = 1, jm + do i = 1, im + +! calculate the distance from the origin +! +#if NDIMS == 3 + rl = xl(i) + yl(j) + zl(k) + ru = xu(i) + yu(j) + zu(k) +#else /* NDIMS == 3 */ + rl = xl(i) + yl(j) + ru = xu(i) + yu(j) +#endif /* NDIMS == 3 */ + +! initialize density and pressure +! + if (ru <= sline) then + q(idn,i) = odens + if (ipr > 0) q(ipr,i) = opres + else if (rl >= sline) then + q(idn,i) = adens + if (ipr > 0) q(ipr,i) = apres + else + ds = (sline - rl) / dx + if (ds <= 1.0d+00) then + dl = 5.0d-01 * ds**ndims + dr = 1.0d+00 - dl + else + ds = (ru - sline) / dx + dr = 5.0d-01 * ds**ndims + dl = 1.0d+00 - dr + end if + + q(idn,i) = adens * dl + odens * dr + if (ipr > 0) q(ipr,i) = apres * dl + opres * dr + end if + + end do ! i = 1, im + +! convert the primitive variables to conservative ones +! + call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im)) + +! copy the conserved variables to the current block +! + pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im) + +! copy the primitive variables to the current block +! + pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im) + + end do ! j = 1, jm + end do ! k = 1, km + +#ifdef PROFILE +! stop accounting time for the problems setup +! + call stop_timer(imu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine setup_problem_implosion +! +!=============================================================================== +! ! subroutine SETUP_PROBLEM_KELVIN_HELMHOLTZ: ! ----------------------------------------- ! From 7255561d2489533f5a7b075f94b37ae14b950987 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 08:45:44 -0300 Subject: [PATCH 03/41] EQUATIONS: Conserved to primitive variables only in interioris. Conversion of the conserved to primitive variables can be expensive for some sets of equations, such as the relativistic MHD. Therefore, we will convert only the interior parts of blocks and update the ghost zones later in the boundary update. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index bde35e8..71cd7e8 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -605,7 +605,8 @@ module equations ! include external procedures and variables ! - use coordinates, only : im, jm, km + use coordinates, only : im, jm, km, in, jn, kn + use coordinates, only : ib, jb, kb, ie, je, ke ! local variables are not implicit by default ! @@ -618,29 +619,21 @@ module equations ! temporary variables ! - integer :: j, k - -! temporary array to store conserved variable vector -! - real(kind=8), dimension(nv,im) :: u + integer :: j, k ! !------------------------------------------------------------------------------- ! ! update primitive variables ! - do k = 1, km - do j = 1, jm - -! copy variables to temporary array of conserved variables -! - u(1:nv,1:im) = uu(1:nv,1:im,j,k) + do k = kb, ke + do j = jb, je ! convert conserved variables to primitive ones ! - call cons2prim(im, u(1:nv,1:im), qq(1:nv,1:im,j,k)) + call cons2prim(in, uu(1:nv,ib:ie,j,k), qq(1:nv,ib:ie,j,k)) - end do ! j = 1, jm - end do ! k = 1, km + end do ! j = jb, je + end do ! k = kb, ke !------------------------------------------------------------------------------- ! From f87451167316125038f601d61753c6e33390ea41 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 08:51:04 -0300 Subject: [PATCH 04/41] EVOLUTION: Update primitive variables before boundaries. The new order is to first, convert the updated conserved variables to their primitive representation in the block interiors, update boundaries of the primitive variables, and finally, convert the primitive variables at the ghost cells to conserved ones. The change updates the primitive variables in the block interiors before the boundary update. Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index 9091111..eadb831 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -249,14 +249,14 @@ module evolution ! call update_mesh() -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + ! set all meta blocks to be updated ! call set_blocks_update(.true.) @@ -467,14 +467,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + !------------------------------------------------------------------------------- ! end subroutine evolve_euler @@ -552,14 +552,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + ! update fluxes for the second step of the RK2 integration ! call update_fluxes() @@ -598,14 +598,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 From 125e267f9710da1d0999e6df26c6c86243d03747 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 09:35:15 -0300 Subject: [PATCH 05/41] BOUNDARIES: Add subroutine to update conservative variables. This subroutine is called after updating the boundaries of primitive variables. It converts those primitive variables to their conservative representation in all ghost cells. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 108 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 5b895bd..a2953e9 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -326,6 +326,10 @@ module boundaries ! call update_corners() +! convert updated primitive variables to conservative ones in all ghost cells +! + call update_ghost_cells() + #ifdef PROFILE ! stop accounting time for variable boundary update ! @@ -1021,6 +1025,110 @@ module boundaries ! !=============================================================================== ! +! subroutine UPDATE_GHOST_CELLS: +! ----------------------------- +! +! Subroutine updates conservative variables in all ghost cells from +! already updated primitive variables. +! +! +!=============================================================================== +! + subroutine update_ghost_cells() + +! include external variables +! + use blocks , only : block_data, list_data + use coordinates , only : im , jm , km , in , jn , kn + use coordinates , only : ib , jb , kb , ie , je , ke + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : prim2cons + +! local variables are not implicit by default +! + implicit none + +! local variables +! + integer :: i, j, k + +! local pointers +! + type(block_data), pointer :: pdata +! +!------------------------------------------------------------------------------- +! +! assign the pointer to the first block on the list +! + pdata => list_data + +! scan all data blocks until the last is reached +! + do while(associated(pdata)) + +! update the X and Y boundary ghost cells +! + do k = 1, km + +! update lower layers of the Y boundary +! + do j = 1, jbl + call prim2cons(im, pdata%q(1:nv,1:im,j,k), pdata%u(1:nv,1:im,j,k)) + end do ! j = 1, jbl + +! update upper layers of the Y boundary +! + do j = jeu, jm + call prim2cons(im, pdata%q(1:nv,1:im,j,k), pdata%u(1:nv,1:im,j,k)) + end do ! j = jeu, jm + +! update remaining left layers of the X boundary +! + do i = 1, ibl + call prim2cons(jn, pdata%q(1:nv,i,jb:je,k), pdata%u(1:nv,i,jb:je,k)) + end do ! i = 1, ibl + +! update remaining right layers of the X boundary +! + do i = ieu, im + call prim2cons(jn, pdata%q(1:nv,i,jb:je,k), pdata%u(1:nv,i,jb:je,k)) + end do ! i = 1, ibl + + end do ! k = 1, km + +#if NDIMS == 3 +! update the Z boundary ghost cells +! + do j = jb, je + +! update the remaining front layers of the Z boundary +! + do k = 1, kbl + call prim2cons(in, pdata%q(1:nv,ib:ie,j,k), pdata%u(1:nv,ib:ie,j,k)) + end do ! k = 1, kbl + +! update the remaining back layers of the Z boundary +! + do k = keu, km + call prim2cons(in, pdata%q(1:nv,ib:ie,j,k), pdata%u(1:nv,ib:ie,j,k)) + end do ! k = keu, km + + end do ! j = jb, je +#endif /* NDIMS == 3 */ + +! assign the pointer to the next block on the list +! + pdata => pdata%next + + end do ! data blocks + +!------------------------------------------------------------------------------- +! + end subroutine update_ghost_cells +! +!=============================================================================== +! ! subroutine SPECIFIC_BOUNDARIES: ! ------------------------------ ! From 722d623e0bb9b7679a324ab2b76af74a2408e873 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 09:40:33 -0300 Subject: [PATCH 06/41] BOUNDARIES: Convert update_corners() to update primitive variables. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index a2953e9..9ab6b2f 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -947,7 +947,7 @@ module boundaries ! local variables ! - integer :: n, i, j, k + integer :: i, j, k, p ! local pointers ! @@ -965,50 +965,50 @@ module boundaries ! iterate over all variables ! - do n = 1, nv + do p = 1, nv ! edges ! #if NDIMS == 3 do i = 1, im - pdata%u(n,i, 1:nh, 1:nh) = pdata%u(n,i,jbl,kbl) - pdata%u(n,i,jt:jm, 1:nh) = pdata%u(n,i,jeu,kbl) - pdata%u(n,i, 1:nh,kt:km) = pdata%u(n,i,jbl,keu) - pdata%u(n,i,jt:jm,kt:km) = pdata%u(n,i,jeu,keu) + pdata%q(p,i, 1:nh, 1:nh) = pdata%q(p,i,jbl,kbl) + pdata%q(p,i,jt:jm, 1:nh) = pdata%q(p,i,jeu,kbl) + pdata%q(p,i, 1:nh,kt:km) = pdata%q(p,i,jbl,keu) + pdata%q(p,i,jt:jm,kt:km) = pdata%q(p,i,jeu,keu) end do do j = 1, jm - pdata%u(n, 1:nh,j, 1:nh) = pdata%u(n,ibl,j,kbl) - pdata%u(n,it:im,j, 1:nh) = pdata%u(n,ieu,j,kbl) - pdata%u(n, 1:nh,j,kt:km) = pdata%u(n,ibl,j,keu) - pdata%u(n,it:im,j,kt:km) = pdata%u(n,ieu,j,keu) + pdata%q(p, 1:nh,j, 1:nh) = pdata%q(p,ibl,j,kbl) + pdata%q(p,it:im,j, 1:nh) = pdata%q(p,ieu,j,kbl) + pdata%q(p, 1:nh,j,kt:km) = pdata%q(p,ibl,j,keu) + pdata%q(p,it:im,j,kt:km) = pdata%q(p,ieu,j,keu) end do #endif /* == 3 */ do k = 1, km - pdata%u(n, 1:nh, 1:nh,k) = pdata%u(n,ibl,jbl,k) - pdata%u(n,it:im, 1:nh,k) = pdata%u(n,ieu,jbl,k) - pdata%u(n, 1:nh,jt:jm,k) = pdata%u(n,ibl,jeu,k) - pdata%u(n,it:im,jt:jm,k) = pdata%u(n,ieu,jeu,k) + pdata%q(p, 1:nh, 1:nh,k) = pdata%q(p,ibl,jbl,k) + pdata%q(p,it:im, 1:nh,k) = pdata%q(p,ieu,jbl,k) + pdata%q(p, 1:nh,jt:jm,k) = pdata%q(p,ibl,jeu,k) + pdata%q(p,it:im,jt:jm,k) = pdata%q(p,ieu,jeu,k) end do ! corners ! #if NDIMS == 3 - pdata%u(n, 1:nh, 1:nh, 1:nh) = pdata%u(n,ibl,jbl,kbl) - pdata%u(n,it:im, 1:nh, 1:nh) = pdata%u(n,ieu,jbl,kbl) - pdata%u(n, 1:nh,jt:jm, 1:nh) = pdata%u(n,ibl,jeu,kbl) - pdata%u(n,it:im,jt:jm, 1:nh) = pdata%u(n,ieu,jeu,kbl) - pdata%u(n, 1:nh, 1:nh,kt:km) = pdata%u(n,ibl,jbl,keu) - pdata%u(n,it:im, 1:nh,kt:km) = pdata%u(n,ieu,jbl,keu) - pdata%u(n, 1:nh,jt:jm,kt:km) = pdata%u(n,ibl,jeu,keu) - pdata%u(n,it:im,jt:jm,kt:km) = pdata%u(n,ieu,jeu,keu) + pdata%q(p, 1:nh, 1:nh, 1:nh) = pdata%q(p,ibl,jbl,kbl) + pdata%q(p,it:im, 1:nh, 1:nh) = pdata%q(p,ieu,jbl,kbl) + pdata%q(p, 1:nh,jt:jm, 1:nh) = pdata%q(p,ibl,jeu,kbl) + pdata%q(p,it:im,jt:jm, 1:nh) = pdata%q(p,ieu,jeu,kbl) + pdata%q(p, 1:nh, 1:nh,kt:km) = pdata%q(p,ibl,jbl,keu) + pdata%q(p,it:im, 1:nh,kt:km) = pdata%q(p,ieu,jbl,keu) + pdata%q(p, 1:nh,jt:jm,kt:km) = pdata%q(p,ibl,jeu,keu) + pdata%q(p,it:im,jt:jm,kt:km) = pdata%q(p,ieu,jeu,keu) #endif /* == 3 */ end do From 10072f5c04db43b6d3e5a9c8ccc8cee77faf3e2d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 09:45:59 -0300 Subject: [PATCH 07/41] BOUNDARIES: Copy primitive variables instead of conservative. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 9ab6b2f..409e729 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -1390,27 +1390,27 @@ module boundaries case(1) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,iel:ie,:,:), idir, iside) + , pneigh%data%q(:,iel:ie,:,:), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,ib:ibu,:,:), idir, iside) + , pneigh%data%q(:,ib:ibu,:,:), idir, iside) end if case(2) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,:,jel:je,:), idir, iside) + , pneigh%data%q(:,:,jel:je,:), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,:,jb:jbu,:), idir, iside) + , pneigh%data%q(:,:,jb:jbu,:), idir, iside) end if #if NDIMS == 3 case(3) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,:,:,kel:ke), idir, iside) + , pneigh%data%q(:,:,:,kel:ke), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,:,:,kb:kbu), idir, iside) + , pneigh%data%q(:,:,:,kb:kbu), idir, iside) end if #endif /* NDIMS == 3 */ end select @@ -1539,9 +1539,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,iel:ie,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,iel:ie,:,:) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,ib:ibu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,ib:ibu,:,:) end if ! associate the pointer with the next block @@ -1567,9 +1567,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jel:je,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jel:je,:) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jb:jbu,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jb:jbu,:) end if ! associate the pointer with the next block @@ -1596,9 +1596,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kel:ke) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kel:ke) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kb:kbu) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kb:kbu) end if ! associate the pointer with the next block @@ -2984,13 +2984,13 @@ module boundaries ! Arguments: ! ! pdata - the pointer to modified data block; -! u - the variable array from which boundaries are updated; +! q - the variable array from which boundaries are updated; ! idir - the direction to be processed; ! iside - the side to be processed; ! !=============================================================================== ! - subroutine boundary_copy(pdata, u, idir, iside) + subroutine boundary_copy(pdata, q, idir, iside) ! import external procedures and variables ! @@ -3005,7 +3005,7 @@ module boundaries ! subroutine arguments ! type(block_data), pointer , intent(inout) :: pdata - real , dimension(:,:,:,:), intent(in) :: u + real , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside ! !------------------------------------------------------------------------------- @@ -3017,26 +3017,26 @@ module boundaries case(1) if (iside == 1) then - pdata%u(1:nv, 1:ibl,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km) + pdata%q(1:nv, 1:ibl,1:jm,1:km) = q(1:nv,1:ng,1:jm,1:km) else - pdata%u(1:nv,ieu:im ,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km) + pdata%q(1:nv,ieu:im ,1:jm,1:km) = q(1:nv,1:ng,1:jm,1:km) end if case(2) if (iside == 1) then - pdata%u(1:nv,1:im, 1:jbl,1:km) = u(1:nv,1:im,1:ng,1:km) + pdata%q(1:nv,1:im, 1:jbl,1:km) = q(1:nv,1:im,1:ng,1:km) else - pdata%u(1:nv,1:im,jeu:jm ,1:km) = u(1:nv,1:im,1:ng,1:km) + pdata%q(1:nv,1:im,jeu:jm ,1:km) = q(1:nv,1:im,1:ng,1:km) end if #if NDIMS == 3 case(3) if (iside == 1) then - pdata%u(1:nv,1:im,1:jm, 1:kbl) = u(1:nv,1:im,1:jm,1:ng) + pdata%q(1:nv,1:im,1:jm, 1:kbl) = q(1:nv,1:im,1:jm,1:ng) else - pdata%u(1:nv,1:im,1:jm,keu:km ) = u(1:nv,1:im,1:jm,1:ng) + pdata%q(1:nv,1:im,1:jm,keu:km ) = q(1:nv,1:im,1:jm,1:ng) end if #endif /* NDIMS == 3 */ From 0876badc169c5542b5feb4f59421686cc07f3e94 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 09:49:21 -0300 Subject: [PATCH 08/41] BOUNDARIES: Restrict primitive variables instead of conservative. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 409e729..6d0f76b 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -1914,7 +1914,7 @@ module boundaries ! update boundaries of the current block ! call boundary_restrict(pdata & - , pneigh%data%u(:,il:iu,jl:ju,kl:ku) & + , pneigh%data%q(:,il:iu,jl:ju,kl:ku) & , idir, iside, iface) #ifdef MPI @@ -2046,7 +2046,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,il:iu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,il:iu,:,:) ! associate the pointer with the next block ! @@ -2081,7 +2081,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jl:ju,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jl:ju,:) ! associate the pointer with the next block ! @@ -2117,7 +2117,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kl:ku) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kl:ku) ! associate the pointer with the next block ! @@ -3058,12 +3058,12 @@ module boundaries ! Arguments: ! ! pdata - the input data block; -! u - the conserved array; +! q - the variable array from which boundaries are updated; ! idir, iside, iface - the positions of the neighbor block; ! !=============================================================================== ! - subroutine boundary_restrict(pdata, u, idir, iside, iface) + subroutine boundary_restrict(pdata, q, idir, iside, iface) ! import external procedures and variables ! @@ -3080,7 +3080,7 @@ module boundaries ! subroutine arguments ! type(block_data) , pointer, intent(inout) :: pdata - real(kind=8) , dimension(:,:,:,:), intent(in) :: u + real(kind=8) , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside, iface ! local variables @@ -3220,22 +3220,22 @@ module boundaries ! update boundaries of the conserved variables ! #if NDIMS == 2 - pdata%u(:,is:it,js:jt, 1 ) = & - 2.50d-01 * ((u(1:nv,il:iu:2,jl:ju:2, 1 ) & - + u(1:nv,ip:iu:2,jp:ju:2, 1 )) & - + (u(1:nv,il:iu:2,jp:ju:2, 1 ) & - + u(1:nv,ip:iu:2,jl:ju:2, 1 ))) + pdata%q(:,is:it,js:jt, 1 ) = & + 2.50d-01 * ((q(1:nv,il:iu:2,jl:ju:2, 1 ) & + + q(1:nv,ip:iu:2,jp:ju:2, 1 )) & + + (q(1:nv,il:iu:2,jp:ju:2, 1 ) & + + q(1:nv,ip:iu:2,jl:ju:2, 1 ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 - pdata%u(:,is:it,js:jt,ks:kt) = & - 1.25d-01 * (((u(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & - + u(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & - + (u(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & - + u(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & - + ((u(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & - + u(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & - + (u(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & - + u(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) + pdata%q(:,is:it,js:jt,ks:kt) = & + 1.25d-01 * (((q(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + + q(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + + (q(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + + q(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + + ((q(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + + q(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + + (q(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + + q(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- From 842f6bdd51f207663e37917dc3ba01c5efba532a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 09:54:12 -0300 Subject: [PATCH 09/41] BOUNDARIES: Prolong primitive variables instead of conservative. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 62 +++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 6d0f76b..45e214b 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -2423,7 +2423,7 @@ module boundaries ! update boundaries of the current block from its neighbor ! call boundary_prolong(pdata & - , pneigh%data%u(:,il:iu,jl:ju,kl:ku) & + , pneigh%data%q(:,il:iu,jl:ju,kl:ku) & , idir, iside, nface) #ifdef MPI @@ -2557,7 +2557,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,il:iu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,il:iu,:,:) ! associate the pointer with the next block ! @@ -2592,7 +2592,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jl:ju,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jl:ju,:) ! associate the pointer with the next block ! @@ -2628,7 +2628,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kl:ku) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kl:ku) ! associate the pointer with the next block ! @@ -3254,12 +3254,12 @@ module boundaries ! Arguments: ! ! pdata - the input data block; -! u - the conserved array; +! q - the variable array from which boundaries are updated; ! idir, iside, iface - the positions of the neighbor block; ! !=============================================================================== ! - subroutine boundary_prolong(pdata, u, idir, iside, iface) + subroutine boundary_prolong(pdata, q, idir, iside, iface) ! import external procedures and variables ! @@ -3277,16 +3277,16 @@ module boundaries ! subroutine arguments ! type(block_data), pointer , intent(inout) :: pdata - real , dimension(:,:,:,:), intent(in) :: u + real , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside, iface ! local variables ! - integer :: i, j, k, q + integer :: i, j, k, p integer :: ic, jc, kc, ip, jp, kp integer :: il, jl, kl, iu, ju, ku integer :: is, js, ks, it, jt, kt - real :: dul, dur, dux, duy, duz + real :: dql, dqr, dqx, dqy, dqz ! !------------------------------------------------------------------------------- ! @@ -3403,37 +3403,37 @@ module boundaries ! iterate over all variables ! - do q = 1, nv + do p = 1, nv - dul = u(q,i ,j,k) - u(q,i-1,j,k) - dur = u(q,i+1,j,k) - u(q,i ,j,k) - dux = limiter(0.25d+00, dul, dur) + dql = q(p,i ,j,k) - q(p,i-1,j,k) + dqr = q(p,i+1,j,k) - q(p,i ,j,k) + dqx = limiter(0.25d+00, dql, dqr) - dul = u(q,i,j ,k) - u(q,i,j-1,k) - dur = u(q,i,j+1,k) - u(q,i,j ,k) - duy = limiter(0.25d+00, dul, dur) + dql = q(p,i,j ,k) - q(p,i,j-1,k) + dqr = q(p,i,j+1,k) - q(p,i,j ,k) + dqy = limiter(0.25d+00, dql, dqr) #if NDIMS == 3 - dul = u(q,i,j,k ) - u(q,i,j,k-1) - dur = u(q,i,j,k+1) - u(q,i,j,k ) - duz = limiter(0.25d+00, dul, dur) + dql = q(p,i,j,k ) - q(p,i,j,k-1) + dqr = q(p,i,j,k+1) - q(p,i,j,k ) + dqz = limiter(0.25d+00, dql, dqr) #endif /* NDIMS == 3 */ #if NDIMS == 2 - pdata%u(q,it,jt,kt) = u(q,i,j,k) - (dux + duy) - pdata%u(q,ip,jt,kt) = u(q,i,j,k) + (dux - duy) - pdata%u(q,it,jp,kt) = u(q,i,j,k) + (duy - dux) - pdata%u(q,ip,jp,kt) = u(q,i,j,k) + (dux + duy) + pdata%q(p,it,jt,kt) = q(p,i,j,k) - (dqx + dqy) + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + (dqx - dqy) + pdata%q(p,it,jp,kt) = q(p,i,j,k) + (dqy - dqx) + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + (dqx + dqy) #endif /* NDIMS == 2 */ #if NDIMS == 3 - pdata%u(q,it,jt,kt) = u(q,i,j,k) - (dux + duy + duz) - pdata%u(q,ip,jt,kt) = u(q,i,j,k) + (dux - duy - duz) - pdata%u(q,it,jp,kt) = u(q,i,j,k) - (dux - duy + duz) - pdata%u(q,ip,jp,kt) = u(q,i,j,k) + (dux + duy - duz) - pdata%u(q,it,jt,kp) = u(q,i,j,k) - (dux + duy - duz) - pdata%u(q,ip,jt,kp) = u(q,i,j,k) + (dux - duy + duz) - pdata%u(q,it,jp,kp) = u(q,i,j,k) - (dux - duy - duz) - pdata%u(q,ip,jp,kp) = u(q,i,j,k) + (dux + duy + duz) + pdata%q(p,it,jt,kt) = q(p,i,j,k) - (dqx + dqy + dqz) + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + (dqx - dqy - dqz) + pdata%q(p,it,jp,kt) = q(p,i,j,k) - (dqx - dqy + dqz) + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + (dqx + dqy - dqz) + pdata%q(p,it,jt,kp) = q(p,i,j,k) - (dqx + dqy - dqz) + pdata%q(p,ip,jt,kp) = q(p,i,j,k) + (dqx - dqy + dqz) + pdata%q(p,it,jp,kp) = q(p,i,j,k) - (dqx - dqy - dqz) + pdata%q(p,ip,jp,kp) = q(p,i,j,k) + (dqx + dqy + dqz) #endif /* NDIMS == 3 */ end do ! q - 1, nv From b528322c63de969a0162cfac9a07d52ccf40244f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 10:05:39 -0300 Subject: [PATCH 10/41] BOUNDARIES: Apply specific boundaries to primitive variables. Convert the specific boundaries update to use primitive variables instead of conservative ones. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 43 ++++++++++++++++++++++--------------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 45e214b..d8c9563 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -2772,10 +2772,11 @@ module boundaries ! import external procedures and variables ! use blocks , only : block_data - use coordinates , only : ng, im, jm, km, ib, ibl, ie, ieu, jb & - , jbl, je, jeu, kb, kbl, ke, keu + use coordinates , only : im , jm , km , ng + use coordinates , only : ib , jb , kb , ie , je , ke + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu use equations , only : nv - use equations , only : idn, imx, imy, imz, ibx, iby, ibz, ibp + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp use error , only : print_warning ! local variables are not implicit by default @@ -2816,15 +2817,15 @@ module boundaries it = ib - i is = ibl + i - pdata%u( :,it,:,:) = pdata%u( :,is,:,:) - pdata%u(imx,it,:,:) = - pdata%u(imx,is,:,:) + pdata%q( :,it,:,:) = pdata%q( :,is,:,:) + pdata%q(ivx,it,:,:) = - pdata%q(ivx,is,:,:) end do case default ! "open" as default boundary conditions do i = 1, ng - pdata%u( :,i,:,:) = pdata%u(:,ib,:,:) + pdata%q( :,i,:,:) = pdata%q(:,ib,:,:) end do end select @@ -2843,14 +2844,14 @@ module boundaries it = ie + i is = ieu - i - pdata%u( :,it,:,:) = pdata%u( :,is,:,:) - pdata%u(imx,it,:,:) = - pdata%u(imx,is,:,:) + pdata%q( :,it,:,:) = pdata%q( :,is,:,:) + pdata%q(ivx,it,:,:) = - pdata%q(ivx,is,:,:) end do case default ! "open" as default boundary conditions do i = ieu, im - pdata%u( :,i,:,:) = pdata%u(:,ie,:,:) + pdata%q( :,i ,:,:) = pdata%q( :,ie,:,:) end do end select @@ -2869,14 +2870,14 @@ module boundaries jt = jb - j js = jbl + j - pdata%u( :,:,jt,:) = pdata%u( :,:,js,:) - pdata%u(imy,:,jt,:) = - pdata%u(imy,:,js,:) + pdata%q( :,:,jt,:) = pdata%q( :,:,js,:) + pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do case default ! "open" as default boundary conditions do j = 1, ng - pdata%u( :,:,j,:) = pdata%u(:,:,jb,:) + pdata%q( :,:,j ,:) = pdata%q( :,:,jb,:) end do end select @@ -2895,14 +2896,14 @@ module boundaries jt = je + j js = jeu - j - pdata%u( :,:,jt,:) = pdata%u( :,:,js,:) - pdata%u(imy,:,jt,:) = - pdata%u(imy,:,js,:) + pdata%q( :,:,jt,:) = pdata%q( :,:,js,:) + pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do case default ! "open" as default boundary conditions do j = jeu, jm - pdata%u( :,:,j,:) = pdata%u(:,:,je,:) + pdata%q( :,:,j ,:) = pdata%q( :,:,je,:) end do end select @@ -2922,14 +2923,14 @@ module boundaries kt = kb - k ks = kbl + k - pdata%u( :,:,:,kt) = pdata%u( :,:,:,ks) - pdata%u(imz,:,:,kt) = - pdata%u(imz,:,:,ks) + pdata%q( :,:,:,kt) = pdata%q( :,:,:,ks) + pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do case default ! "open" as default boundary conditions do k = 1, ng - pdata%u( :,:,:,k) = pdata%u(:,:,:,kb) + pdata%q( :,:,:,k ) = pdata%q( :,:,:,kb) end do end select @@ -2948,14 +2949,14 @@ module boundaries kt = ke + k ks = keu - k - pdata%u( :,:,:,kt) = pdata%u( :,:,:,ks) - pdata%u(imz,:,:,kt) = - pdata%u(imz,:,:,ks) + pdata%q( :,:,:,kt) = pdata%q( :,:,:,ks) + pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do case default ! "open" as default boundary conditions do k = keu, km - pdata%u( :,:,:,k) = pdata%u(:,:,:,ke) + pdata%q( :,:,:,k ) = pdata%q( :,:,:,ke) end do end select From 9b7f0232cdfdf88f31c1dd06bc9c3752d889f110 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 11:01:20 -0300 Subject: [PATCH 11/41] PROBLEMS: Make implosion problem symmetric about axes. Since we use the reflective boundary conditions for this problem, the initial profile should respect that and also apply reflection at the boundaries. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/problems.F90 b/src/problems.F90 index ef33ae7..3d7f216 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -630,13 +630,13 @@ module problems ! calculate edge coordinates ! - xl(:) = x(:) - dxh - xu(:) = x(:) + dxh - yl(:) = y(:) - dyh - yu(:) = y(:) + dyh + xl(:) = abs(x(:)) - dxh + xu(:) = abs(x(:)) + dxh + yl(:) = abs(y(:)) - dyh + yu(:) = abs(y(:)) + dyh #if NDIMS == 3 - zl(:) = z(:) - dzh - zu(:) = z(:) + dzh + zl(:) = abs(z(:)) - dzh + zu(:) = abs(z(:)) + dzh #endif /* NDIMS == 3 */ ! reset velocity components From f85db888ffda309a99d91f0618b01c6cb65e66c3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 9 Apr 2014 13:48:05 -0300 Subject: [PATCH 12/41] BOUNDARIES: Replace string boundary types with integer ones. Now, the boundary type is stored in an integer array instead of string variables. This should reduce the comparison time in boundary_specific(). Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 270 ++++++++++++++++++++++++++++++--------------- src/domains.F90 | 9 +- 2 files changed, 182 insertions(+), 97 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index d8c9563..5a41eb8 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -42,17 +42,18 @@ module boundaries #ifdef PROFILE ! timer indices ! - integer , save :: imi, imv, imf, ims, imc, imr, imp + integer , save :: imi, imv, imf, ims, imc, imr, imp #endif /* PROFILE */ -! module parameters for the boundary update order and boundary type +! parameters corresponding to the boundary type ! - character(len = 32), save :: xlbndry = "periodic" - character(len = 32), save :: xubndry = "periodic" - character(len = 32), save :: ylbndry = "periodic" - character(len = 32), save :: yubndry = "periodic" - character(len = 32), save :: zlbndry = "periodic" - character(len = 32), save :: zubndry = "periodic" + integer, parameter :: bnd_periodic = 0 + integer, parameter :: bnd_open = 1 + integer, parameter :: bnd_reflective = 2 + +! variable to store boundary type flags +! + integer, dimension(3,2), save :: bnd_type = bnd_periodic ! by default everything is private ! @@ -62,7 +63,7 @@ module boundaries ! public :: initialize_boundaries, finalize_boundaries public :: boundary_variables, boundary_fluxes - public :: xlbndry, ylbndry, zlbndry, xubndry, yubndry, zubndry + public :: bnd_type, bnd_periodic !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -105,6 +106,15 @@ module boundaries ! logical, intent(in) :: verbose integer, intent(inout) :: iret + +! module parameters for the boundary update order and boundary type +! + character(len = 32) :: xlbndry = "periodic" + character(len = 32) :: xubndry = "periodic" + character(len = 32) :: ylbndry = "periodic" + character(len = 32) :: yubndry = "periodic" + character(len = 32) :: zlbndry = "periodic" + character(len = 32) :: zubndry = "periodic" ! !------------------------------------------------------------------------------- ! @@ -133,6 +143,62 @@ module boundaries call get_parameter_string ("zlbndry" , zlbndry) call get_parameter_string ("zubndry" , zubndry) +! fill the boundary type flags +! + select case(xlbndry) + case("open") + bnd_type(1,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(1,1) = bnd_reflective + case default + bnd_type(1,1) = bnd_periodic + end select + + select case(xubndry) + case("open") + bnd_type(1,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(1,2) = bnd_reflective + case default + bnd_type(1,2) = bnd_periodic + end select + + select case(ylbndry) + case("open") + bnd_type(2,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(2,1) = bnd_reflective + case default + bnd_type(2,1) = bnd_periodic + end select + + select case(yubndry) + case("open") + bnd_type(2,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(2,2) = bnd_reflective + case default + bnd_type(2,2) = bnd_periodic + end select + + select case(zlbndry) + case("open") + bnd_type(3,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(3,1) = bnd_reflective + case default + bnd_type(3,1) = bnd_periodic + end select + + select case(zubndry) + case("open") + bnd_type(3,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(3,2) = bnd_reflective + case default + bnd_type(3,2) = bnd_periodic + end select + ! print information about the boundary conditions ! if (verbose) then @@ -148,52 +214,6 @@ module boundaries end if -#ifdef MPI -! change the internal boundaries to "exchange" type for the MPI update -! - if (pdims(1) > 1) then - if (periodic(1)) then - xlbndry = "exchange" - xubndry = "exchange" - else - if (pcoords(1) > 0 ) then - xlbndry = "exchange" - end if - if (pcoords(1) < pdims(1)-1) then - xubndry = "exchange" - end if - end if - end if - - if (pdims(2) > 1) then - if (periodic(2)) then - ylbndry = "exchange" - yubndry = "exchange" - else - if (pcoords(2) > 0 ) then - ylbndry = "exchange" - end if - if (pcoords(2) < pdims(2)-1) then - yubndry = "exchange" - end if - end if - end if - - if (pdims(3) > 1) then - if (periodic(3)) then - zlbndry = "exchange" - zubndry = "exchange" - else - if (pcoords(3) > 0 ) then - zlbndry = "exchange" - end if - if (pcoords(3) < pdims(3)-1) then - zubndry = "exchange" - end if - end if - end if -#endif /* MPI */ - #ifdef PROFILE ! stop accounting time for module initialization/finalization ! @@ -2777,7 +2797,7 @@ module boundaries use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp - use error , only : print_warning + use error , only : print_error, print_warning ! local variables are not implicit by default ! @@ -2808,9 +2828,19 @@ module boundaries ! apply selected boundary condition ! - select case(xlbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do i = 1, ng + pdata%q( :,i,:,:) = pdata%q(:,ib,:,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do i = 1, ng @@ -2822,11 +2852,12 @@ module boundaries end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do i = 1, ng - pdata%q( :,i,:,:) = pdata%q(:,ib,:,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left X boundary type!") end select @@ -2836,9 +2867,19 @@ module boundaries ! apply selected boundary condition ! - select case(xubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do i = ieu, im + pdata%q( :,i ,:,:) = pdata%q( :,ie,:,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do i = 1, ng it = ie + i @@ -2848,11 +2889,12 @@ module boundaries pdata%q(ivx,it,:,:) = - pdata%q(ivx,is,:,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do i = ieu, im - pdata%q( :,i ,:,:) = pdata%q( :,ie,:,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right X boundary type!") end select @@ -2862,9 +2904,19 @@ module boundaries ! apply selected boundary condition ! - select case(ylbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do j = 1, ng + pdata%q( :,:,j ,:) = pdata%q( :,:,jb,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do j = 1, ng jt = jb - j @@ -2874,11 +2926,12 @@ module boundaries pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do j = 1, ng - pdata%q( :,:,j ,:) = pdata%q( :,:,jb,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left Y boundary type!") end select @@ -2888,9 +2941,19 @@ module boundaries ! apply selected boundary condition ! - select case(yubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do j = jeu, jm + pdata%q( :,:,j ,:) = pdata%q( :,:,je,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do j = 1, ng jt = je + j @@ -2900,11 +2963,12 @@ module boundaries pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do j = jeu, jm - pdata%q( :,:,j ,:) = pdata%q( :,:,je,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right Y boundary type!") end select @@ -2915,9 +2979,19 @@ module boundaries ! apply selected boundary condition ! - select case(zlbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do k = 1, ng + pdata%q( :,:,:,k ) = pdata%q( :,:,:,kb) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do k = 1, ng kt = kb - k @@ -2927,11 +3001,12 @@ module boundaries pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do k = 1, ng - pdata%q( :,:,:,k ) = pdata%q( :,:,:,kb) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left Z boundary type!") end select @@ -2941,9 +3016,19 @@ module boundaries ! apply selected boundary condition ! - select case(zubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do k = keu, km + pdata%q( :,:,:,k ) = pdata%q( :,:,:,ke) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do k = 1, ng kt = ke + k @@ -2953,11 +3038,12 @@ module boundaries pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do k = keu, km - pdata%q( :,:,:,k ) = pdata%q( :,:,:,ke) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right Z boundary type!") end select #endif /* NDIMS == 3 */ diff --git a/src/domains.F90 b/src/domains.F90 index f32de85..8c3eda3 100644 --- a/src/domains.F90 +++ b/src/domains.F90 @@ -143,8 +143,7 @@ module domains use blocks , only : metablock_set_configuration use blocks , only : metablock_set_coordinates, metablock_set_bounds use blocks , only : nsides, nfaces - use boundaries , only : xlbndry, ylbndry, zlbndry - use boundaries , only : xubndry, yubndry, zubndry + use boundaries , only : bnd_type, bnd_periodic use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : ir, jr, kr @@ -349,7 +348,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (xlbndry == 'periodic' .and. xubndry == 'periodic') then + if (bnd_type(1,1) == bnd_periodic .and. bnd_type(1,2) == bnd_periodic) then do k = 1, kr do j = 1, jr @@ -395,7 +394,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (ylbndry == 'periodic' .and. yubndry == 'periodic') then + if (bnd_type(2,1) == bnd_periodic .and. bnd_type(2,2) == bnd_periodic) then do k = 1, kr do i = 1, ir @@ -442,7 +441,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (zlbndry == 'periodic' .and. zubndry == 'periodic') then + if (bnd_type(3,1) == bnd_periodic .and. bnd_type(3,2) == bnd_periodic) then do j = 1, jr do i = 1, ir From a4805761db74fac0c908e27197bfb55571a8bbda Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 17 Apr 2014 11:31:33 -0300 Subject: [PATCH 13/41] BOUNDARIES: Interpolation symmetry in boundary_prolong(). Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 5a41eb8..ee922d8 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -3509,7 +3509,7 @@ module boundaries #if NDIMS == 2 pdata%q(p,it,jt,kt) = q(p,i,j,k) - (dqx + dqy) pdata%q(p,ip,jt,kt) = q(p,i,j,k) + (dqx - dqy) - pdata%q(p,it,jp,kt) = q(p,i,j,k) + (dqy - dqx) + pdata%q(p,it,jp,kt) = q(p,i,j,k) - (dqx - dqy) pdata%q(p,ip,jp,kt) = q(p,i,j,k) + (dqx + dqy) #endif /* NDIMS == 2 */ #if NDIMS == 3 From b226959097629165062625d835c9ab82a57891d3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 17 Apr 2014 11:35:15 -0300 Subject: [PATCH 14/41] MESH: Interpolation symmetry in prolong_block(). Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 236e14d..e247a76 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1588,19 +1588,19 @@ module mesh #if NDIMS == 2 u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy) u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy) - u(p,ic,jp,kc) = pdata%u(p,i,j,k) + (duy - dux) + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - (dux - duy) u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy) #endif /* NDIMS == 2 */ #if NDIMS == 3 - u(p,ic,jc,kc) = pdata%u(p,i,j,k) - dux - duy - duz - u(p,ip,jc,kc) = pdata%u(p,i,j,k) + dux - duy - duz - u(p,ic,jp,kc) = pdata%u(p,i,j,k) - dux + duy - duz - u(p,ip,jp,kc) = pdata%u(p,i,j,k) + dux + duy - duz - u(p,ic,jc,kp) = pdata%u(p,i,j,k) - dux - duy + duz - u(p,ip,jc,kp) = pdata%u(p,i,j,k) + dux - duy + duz - u(p,ic,jp,kp) = pdata%u(p,i,j,k) - dux + duy + duz - u(p,ip,jp,kp) = pdata%u(p,i,j,k) + dux + duy + duz + u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy + duz) + u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy - duz) + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - (dux - duy + duz) + u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy - duz) + u(p,ic,jc,kp) = pdata%u(p,i,j,k) - (dux + duy - duz) + u(p,ip,jc,kp) = pdata%u(p,i,j,k) + (dux - duy + duz) + u(p,ic,jp,kp) = pdata%u(p,i,j,k) - (dux - duy - duz) + u(p,ip,jp,kp) = pdata%u(p,i,j,k) + (dux + duy + duz) #endif /* NDIMS == 3 */ end do end do From d3a5055cf5ac9ba51452cf3b719f3e6120f05b4a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 17 Apr 2014 11:44:47 -0300 Subject: [PATCH 15/41] BOUNDARIES: Micro optimization in boundary_prolong(). Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index ee922d8..f59ac04 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -3373,7 +3373,7 @@ module boundaries integer :: ic, jc, kc, ip, jp, kp integer :: il, jl, kl, iu, ju, ku integer :: is, js, ks, it, jt, kt - real :: dql, dqr, dqx, dqy, dqz + real :: dql, dqr, dqx, dqy, dqz, dq1, dq2, dq3, dq4 ! !------------------------------------------------------------------------------- ! @@ -3507,20 +3507,26 @@ module boundaries #endif /* NDIMS == 3 */ #if NDIMS == 2 - pdata%q(p,it,jt,kt) = q(p,i,j,k) - (dqx + dqy) - pdata%q(p,ip,jt,kt) = q(p,i,j,k) + (dqx - dqy) - pdata%q(p,it,jp,kt) = q(p,i,j,k) - (dqx - dqy) - pdata%q(p,ip,jp,kt) = q(p,i,j,k) + (dqx + dqy) + dq1 = dqx + dqy + dq2 = dqx - dqy + pdata%q(p,it,jt,kt) = q(p,i,j,k) - dq1 + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + dq2 + pdata%q(p,it,jp,kt) = q(p,i,j,k) - dq2 + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 - pdata%q(p,it,jt,kt) = q(p,i,j,k) - (dqx + dqy + dqz) - pdata%q(p,ip,jt,kt) = q(p,i,j,k) + (dqx - dqy - dqz) - pdata%q(p,it,jp,kt) = q(p,i,j,k) - (dqx - dqy + dqz) - pdata%q(p,ip,jp,kt) = q(p,i,j,k) + (dqx + dqy - dqz) - pdata%q(p,it,jt,kp) = q(p,i,j,k) - (dqx + dqy - dqz) - pdata%q(p,ip,jt,kp) = q(p,i,j,k) + (dqx - dqy + dqz) - pdata%q(p,it,jp,kp) = q(p,i,j,k) - (dqx - dqy - dqz) - pdata%q(p,ip,jp,kp) = q(p,i,j,k) + (dqx + dqy + dqz) + dq1 = dqx + dqy + dqz + dq2 = dqx - dqy - dqz + dq3 = dqx - dqy + dqz + dq4 = dqx + dqy - dqz + pdata%q(p,it,jt,kt) = q(p,i,j,k) - dq1 + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + dq2 + pdata%q(p,it,jp,kt) = q(p,i,j,k) - dq3 + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + dq4 + pdata%q(p,it,jt,kp) = q(p,i,j,k) - dq4 + pdata%q(p,ip,jt,kp) = q(p,i,j,k) + dq3 + pdata%q(p,it,jp,kp) = q(p,i,j,k) - dq2 + pdata%q(p,ip,jp,kp) = q(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ end do ! q - 1, nv From 0a81a0124247dd168c4e99cec8ce7c8101552ec3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 17 Apr 2014 11:52:06 -0300 Subject: [PATCH 16/41] MESH: Micro optimization in prolong_block(). Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index e247a76..0625515 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1502,7 +1502,7 @@ module mesh integer :: i, j, k, q, p integer :: il, iu, jl, ju, kl, ku integer :: ic, jc, kc, ip, jp, kp - real :: dul, dur, dux, duy, duz + real :: dul, dur, dux, duy, duz, du1, du2, du3, du4 ! local pointers ! @@ -1586,21 +1586,27 @@ module mesh #endif /* NDIMS == 3 */ #if NDIMS == 2 - u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy) - u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy) - u(p,ic,jp,kc) = pdata%u(p,i,j,k) - (dux - duy) - u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy) + du1 = dux + duy + du2 = dux - duy + u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 + u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2 + u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 2 */ #if NDIMS == 3 - u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy + duz) - u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy - duz) - u(p,ic,jp,kc) = pdata%u(p,i,j,k) - (dux - duy + duz) - u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy - duz) - u(p,ic,jc,kp) = pdata%u(p,i,j,k) - (dux + duy - duz) - u(p,ip,jc,kp) = pdata%u(p,i,j,k) + (dux - duy + duz) - u(p,ic,jp,kp) = pdata%u(p,i,j,k) - (dux - duy - duz) - u(p,ip,jp,kp) = pdata%u(p,i,j,k) + (dux + duy + duz) + du1 = dux + duy + duz + du2 = dux - duy - duz + du3 = dux - duy + duz + du4 = dux + duy - duz + u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 + u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3 + u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4 + u(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4 + u(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3 + u(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2 + u(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 3 */ end do end do From f012eefe9f1b740f5f14b5c57e1188b07f7645ea Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 19 Apr 2014 13:50:44 -0300 Subject: [PATCH 17/41] IO: Slightly rewrite write_restart_snapshot(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index cecf5a3..c602af2 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -58,7 +58,7 @@ module io ! irest - the local counter for the restart snapshots; ! isnap - the local counter for the regular snapshots; ! hrest - the execution time interval for restart snapshot storing -! (in hours); +! (in hours); the minimum allowed value is 3 minutes; ! hsnap - the problem time interval for regular snapshot storing; ! character(len=255), save :: respath = "./" @@ -66,7 +66,7 @@ module io integer , save :: nrest = -1 integer(kind=4) , save :: irest = 1 integer(kind=4) , save :: isnap = -1 - real , save :: hrest = 6.0d+00 + real , save :: hrest = 6.0e+00 real , save :: hsnap = 1.0d+00 ! flags to determine the way of data writing @@ -302,14 +302,14 @@ module io ! subroutine WRITE_RESTART_SNAPSHOTS: ! ---------------------------------- ! -! Subroutine stores current snapshot in restart files. -! This is a wrapper calling specific format subroutine. +! Subroutine stores current restart snapshot files. This is a wrapper +! calling specific format subroutine. ! ! Arguments: ! -! thrs - the interval time in hours for storing the restart snapshots; -! nrun - the restart snapshot number; -! iret - the return flag to inform if subroutine succeeded or failed; +! thrs - the current execution time in hours; +! nrun - the run number; +! iret - the return flag; ! !=============================================================================== ! @@ -327,28 +327,25 @@ module io ! !------------------------------------------------------------------------------- ! +! check if the condition for storing the restart snapshot has been met +! + if (hrest < 5.0e-02 .or. thrs < irest * hrest) return + #ifdef PROFILE ! start accounting time for the data writing ! call start_timer(iow) #endif /* PROFILE */ -! store restart snapshots at constant interval or when the job is done only -! if hrest is a positive number -! - if (hrest > 0.0d+00 .and. thrs >= irest * hrest) then - #ifdef HDF5 ! store restart file ! - call write_restart_snapshot_h5(nrun, iret) + call write_restart_snapshot_h5(nrun, iret) #endif /* HDF5 */ ! increase the restart snapshot counter ! - irest = irest + 1 - - end if + irest = irest + 1 #ifdef PROFILE ! stop accounting time for the data writing From 67a9beaf65d88fe34fd59005565b1e887989d197 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 11:49:55 -0300 Subject: [PATCH 18/41] IO: Make the snapshot number continuously increasing. If we restart a job with a different snapshot interval, the snapshot number will continue from a different value corresponding to the simulation time. The snapshot number shouldn't be related to the current simulation time. It should be always increasing. The simulation time is stored as an attribute, anyway. This patch makes the snapshot number always increasing continuously, without any gaps for restarted jobs. Signed-off-by: Grzegorz Kowal --- src/io.F90 | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index c602af2..09576f7 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -57,17 +57,22 @@ module io ! nrest - for job restarting, this is the number of restart snapshot; ! irest - the local counter for the restart snapshots; ! isnap - the local counter for the regular snapshots; +! ishift - the shift of the snapshot counter for restarting job with +! different snapshot interval; ! hrest - the execution time interval for restart snapshot storing ! (in hours); the minimum allowed value is 3 minutes; ! hsnap - the problem time interval for regular snapshot storing; +! tsnap - the next snapshot time; ! character(len=255), save :: respath = "./" character , save :: ftype = "p" integer , save :: nrest = -1 integer(kind=4) , save :: irest = 1 - integer(kind=4) , save :: isnap = -1 + integer(kind=4) , save :: isnap = 0 + integer(kind=4) , save :: ishift = 0 real , save :: hrest = 6.0e+00 - real , save :: hsnap = 1.0d+00 + real , save :: hsnap = 1.0e+00 + real , save :: tsnap = 0.0e+00 ! flags to determine the way of data writing ! @@ -283,9 +288,10 @@ module io call read_restart_snapshot_h5() #endif /* HDF5 */ -! calculate the next snapshot number +! calculate the shift of the snapshot counter, and the next snapshot time ! - isnap = int(time / hsnap) + ishift = int(time / hsnap) - isnap + 1 + tsnap = (ishift + isnap) * hsnap #ifdef PROFILE ! stop accounting time for the data reading @@ -327,7 +333,7 @@ module io ! !------------------------------------------------------------------------------- ! -! check if the condition for storing the restart snapshot has been met +! check if conditions for storing the restart snapshot have been met ! if (hrest < 5.0e-02 .or. thrs < irest * hrest) return @@ -381,9 +387,9 @@ module io ! !------------------------------------------------------------------------------- ! -! exit the subroutine, if the time of the next snapshot is not reached +! check if conditions for storing the regular snapshot have been met ! - if (hsnap <= 0.0d+00 .or. isnap >= (int(time / hsnap))) return + if (hsnap <= 0.0e+00 .or. time < tsnap) return #ifdef PROFILE ! start accounting time for the data writing @@ -391,16 +397,17 @@ module io call start_timer(iow) #endif /* PROFILE */ -! increase the file counter -! - isnap = isnap + 1 - #ifdef HDF5 ! store variable snapshot file ! call write_snapshot_h5() #endif /* HDF5 */ +! increase the snapshot counter and calculate the next snapshot time +! + isnap = isnap + 1 + tsnap = (ishift + isnap) * hsnap + #ifdef PROFILE ! stop accounting time for the data writing ! @@ -430,7 +437,7 @@ module io !------------------------------------------------------------------------------- ! if (hsnap > 0.0d+00) then - next_tout = (isnap + 1) * hsnap + next_tout = tsnap else next_tout = huge(hsnap) end if @@ -1092,6 +1099,7 @@ module io call write_attribute_integer_h5(gid, 'nproc' , nproc) call write_attribute_integer_h5(gid, 'nseeds' , nseeds) call write_attribute_integer_h5(gid, 'step' , step ) + call write_attribute_integer_h5(gid, 'isnap' , isnap ) ! store the real attributes ! @@ -1317,6 +1325,8 @@ module io end if case('step') call read_attribute_integer_h5(aid, aname, step) + case('isnap') + call read_attribute_integer_h5(aid, aname, isnap) case('time') call read_attribute_double_h5(aid, aname, time) case('dt') From 3132fb8da8e605612387c3fc687f151b40a04c8d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 11:56:48 -0300 Subject: [PATCH 19/41] IO: Correct the description of function next_tout(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io.F90 b/src/io.F90 index 09576f7..3b510ab 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -423,7 +423,7 @@ module io ! function NEXT_TOUT: ! ------------------ ! -! Function returns time of the next data snapshot. +! Function returns the next data snapshot time. ! ! !=============================================================================== From b746cd1c7ea2d8bdccfd2a9c0eee4ad202791dcf Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 12:58:12 -0300 Subject: [PATCH 20/41] IO: Rewrite write_attributes_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 195 +++++++++++++++++++++++++++-------------------------- 1 file changed, 98 insertions(+), 97 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 3b510ab..381b3f0 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1030,46 +1030,49 @@ module io ! !=============================================================================== ! -! write_attributes_h5: subroutine writes attributes in the HDF5 format -! connected to the provided identificator +! subroutine WRITE_ATTRIBUTES_H5: +! ------------------------------ ! -! info: this subroutine stores only the global attributes +! Subroutine stores global attributes in the HDF5 file provided by an +! identificator. ! -! arguments: -! fid - the HDF5 file identificator +! Arguments: +! +! fid - the HDF5 file identificator; ! !=============================================================================== ! subroutine write_attributes_h5(fid) -! references to other modules +! import external procedures and variables ! - use blocks , only : get_mblocks, get_dblocks, get_nleafs - use blocks , only : get_last_id - use coordinates, only : nn, ng, in, jn, kn, minlev, maxlev, toplev, ir, jr, kr - use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax - use error , only : print_error - use evolution, only : step, time, dt, dtn - use hdf5 , only : hid_t - use hdf5 , only : h5gcreate_f, h5gclose_f - use mpitools , only : nprocs, nproc - use random , only : nseeds, get_seeds + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use blocks , only : get_last_id + use coordinates , only : minlev, maxlev, toplev + use coordinates , only : nn, ng, in, jn, kn, ir, jr, kr + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use error , only : print_error + use evolution , only : step, time, dt, dtn + use hdf5 , only : hid_t + use hdf5 , only : h5gcreate_f, h5gclose_f + use mpitools , only : nprocs, nproc + use random , only : nseeds, get_seeds -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! - integer(hid_t) :: gid - integer(kind=4) :: dm(3) - integer :: err + integer(hid_t) :: gid + integer :: err + integer(kind=4), dimension(3) :: dm -! allocatable arrays +! local allocatable arrays ! integer(kind=4), dimension(:), allocatable :: seeds ! @@ -1081,81 +1084,7 @@ module io ! check if the group has been created successfuly ! - if (err .ge. 0) then - -! store the integer attributes -! - call write_attribute_integer_h5(gid, 'ndims' , NDIMS) - call write_attribute_integer_h5(gid, 'last_id', get_last_id()) - call write_attribute_integer_h5(gid, 'mblocks', get_mblocks()) - call write_attribute_integer_h5(gid, 'dblocks', get_dblocks()) - call write_attribute_integer_h5(gid, 'nleafs' , get_nleafs()) - call write_attribute_integer_h5(gid, 'ncells' , nn) - call write_attribute_integer_h5(gid, 'nghost' , ng) - call write_attribute_integer_h5(gid, 'minlev' , minlev) - call write_attribute_integer_h5(gid, 'maxlev' , maxlev) - call write_attribute_integer_h5(gid, 'toplev' , toplev) - call write_attribute_integer_h5(gid, 'nprocs' , nprocs) - call write_attribute_integer_h5(gid, 'nproc' , nproc) - call write_attribute_integer_h5(gid, 'nseeds' , nseeds) - call write_attribute_integer_h5(gid, 'step' , step ) - call write_attribute_integer_h5(gid, 'isnap' , isnap ) - -! store the real attributes -! - call write_attribute_double_h5(gid, 'xmin', xmin) - call write_attribute_double_h5(gid, 'xmax', xmax) - call write_attribute_double_h5(gid, 'ymin', ymin) - call write_attribute_double_h5(gid, 'ymax', ymax) - call write_attribute_double_h5(gid, 'zmin', zmin) - call write_attribute_double_h5(gid, 'zmax', zmax) - call write_attribute_double_h5(gid, 'time', time) - call write_attribute_double_h5(gid, 'dt' , dt ) - call write_attribute_double_h5(gid, 'dtn' , dtn ) - -! store the vector attributes -! - dm(:) = (/ in, jn, kn /) - call write_attribute_vector_integer_h5(gid, 'dims' , 3, dm(:)) - call write_attribute_vector_integer_h5(gid, 'rdims', 3, (/ ir, jr, kr /)) - -! store random number generator seed values -! - if (nseeds .gt. 0) then - -! allocate space for seeds -! - allocate(seeds(nseeds)) - -! get the seed values -! - call get_seeds(seeds) - -! store them in the current group -! - call write_attribute_vector_integer_h5(gid, 'seeds', nseeds, seeds(:)) - -! deallocate seed array -! - deallocate(seeds) - - end if ! nseeds > 0 - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::write_attributes_h5", "Cannot close the group!") - - end if - - else + if (err < 0) then ! print error about the problem with creating the group ! @@ -1163,6 +1092,78 @@ module io end if +! store the integer attributes +! + call write_attribute_integer_h5(gid, 'ndims' , NDIMS ) + call write_attribute_integer_h5(gid, 'last_id', get_last_id()) + call write_attribute_integer_h5(gid, 'mblocks', get_mblocks()) + call write_attribute_integer_h5(gid, 'dblocks', get_dblocks()) + call write_attribute_integer_h5(gid, 'nleafs' , get_nleafs() ) + call write_attribute_integer_h5(gid, 'ncells' , nn ) + call write_attribute_integer_h5(gid, 'nghost' , ng ) + call write_attribute_integer_h5(gid, 'minlev' , minlev ) + call write_attribute_integer_h5(gid, 'maxlev' , maxlev ) + call write_attribute_integer_h5(gid, 'toplev' , toplev ) + call write_attribute_integer_h5(gid, 'nprocs' , nprocs ) + call write_attribute_integer_h5(gid, 'nproc' , nproc ) + call write_attribute_integer_h5(gid, 'nseeds' , nseeds ) + call write_attribute_integer_h5(gid, 'step' , step ) + call write_attribute_integer_h5(gid, 'isnap' , isnap ) + +! store the real attributes +! + call write_attribute_double_h5(gid, 'xmin', xmin) + call write_attribute_double_h5(gid, 'xmax', xmax) + call write_attribute_double_h5(gid, 'ymin', ymin) + call write_attribute_double_h5(gid, 'ymax', ymax) + call write_attribute_double_h5(gid, 'zmin', zmin) + call write_attribute_double_h5(gid, 'zmax', zmax) + call write_attribute_double_h5(gid, 'time', time) + call write_attribute_double_h5(gid, 'dt' , dt ) + call write_attribute_double_h5(gid, 'dtn' , dtn ) + +! store the vector attributes +! + dm(:) = (/ in, jn, kn /) + call write_attribute_vector_integer_h5(gid, 'dims' , 3, dm(:)) + call write_attribute_vector_integer_h5(gid, 'rdims', 3, (/ ir, jr, kr /)) + +! store random number generator seed values +! + if (nseeds > 0) then + +! allocate space for seeds +! + allocate(seeds(nseeds)) + +! get the seed values +! + call get_seeds(seeds) + +! store them in the current group +! + call write_attribute_vector_integer_h5(gid, 'seeds', nseeds, seeds(:)) + +! deallocate seed array +! + deallocate(seeds) + + end if ! nseeds > 0 + +! close the group +! + call h5gclose_f(gid, err) + +! check if the group has been closed successfuly +! + if (err < 0) then + +! print error about the problem with closing the group +! + call print_error("io::write_attributes_h5", "Cannot close the group!") + + end if + !------------------------------------------------------------------------------- ! end subroutine write_attributes_h5 From 29c842155fc790da4986b9d5fa271259614b9194 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 13:13:44 -0300 Subject: [PATCH 21/41] IO: Quit write_attributes_h5() if cannot create group. Signed-off-by: Grzegorz Kowal --- src/io.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/io.F90 b/src/io.F90 index 381b3f0..1c4c21f 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1090,6 +1090,10 @@ module io ! call print_error("io::write_attributes_h5", "Cannot create the group!") +! return from the subroutine +! + return + end if ! store the integer attributes From 740fedede929ee2d1ae4f113149560d01c981d9f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 15:46:19 -0300 Subject: [PATCH 22/41] IO: Rewrite read_attribute_integer_h5(). The attribute, provided by the group identifier to which it is linked and its name, is opened inside the subroutine read_attribute_integer_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 177 ++++++++++++++++++++++++++++++++--------------------- 1 file changed, 108 insertions(+), 69 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 1c4c21f..de53994 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1268,7 +1268,7 @@ module io ! select case(trim(aname)) case('ndims') - call read_attribute_integer_h5(aid, aname, lndims) + call read_attribute_integer_h5(gid, aname, lndims) ! check if the restart file and compiled program have the same number of ! dimensions @@ -1278,7 +1278,7 @@ module io , "File and program dimensions are incompatible!") end if case('maxlev') - call read_attribute_integer_h5(aid, aname, lmaxlev) + call read_attribute_integer_h5(gid, aname, lmaxlev) if (lmaxlev .gt. toplev) then ! subtitute the new value of toplev @@ -1301,17 +1301,17 @@ module io ucor = 2**(maxlev - lmaxlev) end if case('nprocs') - call read_attribute_integer_h5(aid, aname, nfiles) + call read_attribute_integer_h5(gid, aname, nfiles) case('nproc') - call read_attribute_integer_h5(aid, aname, lnproc) + call read_attribute_integer_h5(gid, aname, lnproc) case('last_id') - call read_attribute_integer_h5(aid, aname, llast_id) + call read_attribute_integer_h5(gid, aname, llast_id) case('mblocks') - call read_attribute_integer_h5(aid, aname, lmblocks) + call read_attribute_integer_h5(gid, aname, lmblocks) case('nleafs') - call read_attribute_integer_h5(aid, aname, lnleafs) + call read_attribute_integer_h5(gid, aname, lnleafs) case('ncells') - call read_attribute_integer_h5(aid, aname, lncells) + call read_attribute_integer_h5(gid, aname, lncells) ! check if the block dimensions are compatible ! @@ -1320,7 +1320,7 @@ module io , "File and program block dimensions are incompatible!") end if case('nghost') - call read_attribute_integer_h5(aid, aname, lnghost) + call read_attribute_integer_h5(gid, aname, lnghost) ! check if the ghost layers are compatible ! @@ -1329,9 +1329,9 @@ module io , "File and program block ghost layers are incompatible!") end if case('step') - call read_attribute_integer_h5(aid, aname, step) + call read_attribute_integer_h5(gid, aname, step) case('isnap') - call read_attribute_integer_h5(aid, aname, isnap) + call read_attribute_integer_h5(gid, aname, isnap) case('time') call read_attribute_double_h5(aid, aname, time) case('dt') @@ -1351,7 +1351,7 @@ module io case('zmax') call read_attribute_double_h5(aid, aname, zmax) case('nseeds') - call read_attribute_integer_h5(aid, aname, lnseeds) + call read_attribute_integer_h5(gid, aname, lnseeds) ! ! check if the numbers of seeds are compatible ! ! @@ -1547,7 +1547,7 @@ module io ! obtain the number of data blocks ! - call read_attribute_integer_h5(aid, aname, ldblocks) + call read_attribute_integer_h5(gid, aname, ldblocks) case('dims') @@ -1559,7 +1559,7 @@ module io ! obtain the number of data blocks ! - call read_attribute_integer_h5(aid, aname, lnghost) + call read_attribute_integer_h5(gid, aname, lnghost) case default end select @@ -1634,6 +1634,100 @@ module io end subroutine read_datablock_dims_h5 ! !=============================================================================== +!! +!!--- ATTRIBUTE SUBROUTINES -------------------------------------------------- +!! +!=============================================================================== +! +!=============================================================================== +! +! subroutine READ_ATTRIBUTES_INTEGER_H5: +! ------------------------------------- +! +! Subroutine reads a value of the integer value provided by the group +! identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_attribute_integer_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5aexists_by_name_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! input variables +! + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + integer , intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_integer_h5" & + , "Attribute does not exist :" // trim(aname)) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_INTEGER, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot read attribute :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot close attribute :" // trim(aname)) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_attribute_integer_h5 +! +!=============================================================================== ! ! write_attribute_integer_h5: subroutine writes an attribute with the integer ! value in a group given by its indentificator @@ -1757,61 +1851,6 @@ module io ! !=============================================================================== ! -! read_attribute_integer_h5: subroutine read an integer value from an attribute -! given by the identificator -! -! arguments: -! aid - the HDF5 attribute identificator -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_integer_h5(aid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - integer , intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_INTEGER, value, am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute -! - call print_error("io::read_attribute_integer_h5" & - , "Cannot read attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine read_attribute_integer_h5 -! -!=============================================================================== -! ! write_attribute_vector_integer_h5: subroutine writes an vector intiger ! attribute in a group given by its indentificator ! From 624ee8be66caf5914827c0700b59a5dec9bf8d28 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 15:49:12 -0300 Subject: [PATCH 23/41] IO: Spell corrections in read_attribute_integer_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io.F90 b/src/io.F90 index de53994..e132701 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1644,7 +1644,7 @@ module io ! subroutine READ_ATTRIBUTES_INTEGER_H5: ! ------------------------------------- ! -! Subroutine reads a value of the integer value provided by the group +! Subroutine reads a value of the integer attribute provided by the group ! identifier to which it is linked and its name. ! ! Arguments: From 6e6e4141e06bf92e263c0e1abd48be406f3eb323 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 16:06:42 -0300 Subject: [PATCH 24/41] IO: Rewrite write_attribute_integer_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 134 ++++++++++++++++++++--------------------------------- 1 file changed, 51 insertions(+), 83 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index e132701..8cc87fb 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1641,8 +1641,8 @@ module io ! !=============================================================================== ! -! subroutine READ_ATTRIBUTES_INTEGER_H5: -! ------------------------------------- +! subroutine READ_ATTRIBUTE_INTEGER_H5: +! ------------------------------------ ! ! Subroutine reads a value of the integer attribute provided by the group ! identifier to which it is linked and its name. @@ -1668,7 +1668,7 @@ module io ! implicit none -! input variables +! attribute arguments ! integer(hid_t) , intent(in) :: gid character(len=*), intent(in) :: aname @@ -1729,120 +1729,88 @@ module io ! !=============================================================================== ! -! write_attribute_integer_h5: subroutine writes an attribute with the integer -! value in a group given by its indentificator +! subroutine WRITE_ATTRIBUTE_INTEGER_H5: +! ------------------------------------- ! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! value - the attribute value +! Subroutine stores a value of the integer attribute in the group provided +! by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute value; ! !=============================================================================== ! - subroutine write_attribute_integer_h5(gid, name, value) + subroutine write_attribute_integer_h5(gid, aname, avalue) -! references to other modules +! import external procedures and variables ! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5screate_simple_f, h5sclose_f & - , h5acreate_f, h5awrite_f, h5aclose_f + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! attribute arguments ! integer(hid_t) , intent(in) :: gid - character(len=*), intent(in) :: name - integer , intent(in) :: value + character(len=*), intent(in) :: aname + integer , intent(in) :: avalue ! local variables ! integer(hid_t) :: sid, aid integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err + integer :: ierr ! !------------------------------------------------------------------------------- ! -! create space for the attribute +! create space for the attribute value ! - call h5screate_simple_f(1, am, sid, err) + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if -! check if the space has been created successfuly +! create the attribute in the given group ! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_INTEGER, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then + call h5acreate_f(gid, aname, H5T_NATIVE_INTEGER, sid, aid, ierr) + if (ierr == 0) then ! write the attribute data ! - call h5awrite_f(aid, H5T_NATIVE_INTEGER, value, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if + call h5awrite_f(aid, H5T_NATIVE_INTEGER, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if ! close the attribute ! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot close space for attribute :" // trim(name)) - + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot close attribute :" // trim(aname)) end if else + call print_error("io::write_attribute_integer_h5" & + , "Cannot create attribute :" // trim(aname)) + end if -! print error about the problem with creating the space for the attribute +! release the space ! - call print_error("io::write_attribute_integer_h5" & - , "Cannot create space for attribute :" // trim(name)) - + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot close space for attribute :" // trim(aname)) end if !------------------------------------------------------------------------------- From 3cd957568542e2d40c2be255064dd775384f4a14 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 16:13:36 -0300 Subject: [PATCH 25/41] IO: Rewrite subroutine read_attribute_double_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 160 ++++++++++++++++++++++++++++++++--------------------- 1 file changed, 97 insertions(+), 63 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 8cc87fb..21c390a 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1333,23 +1333,23 @@ module io case('isnap') call read_attribute_integer_h5(gid, aname, isnap) case('time') - call read_attribute_double_h5(aid, aname, time) + call read_attribute_double_h5(gid, aname, time) case('dt') - call read_attribute_double_h5(aid, aname, dt) + call read_attribute_double_h5(gid, aname, dt) case('dtn') - call read_attribute_double_h5(aid, aname, dtn) + call read_attribute_double_h5(gid, aname, dtn) case('xmin') - call read_attribute_double_h5(aid, aname, xmin) + call read_attribute_double_h5(gid, aname, xmin) case('xmax') - call read_attribute_double_h5(aid, aname, xmax) + call read_attribute_double_h5(gid, aname, xmax) case('ymin') - call read_attribute_double_h5(aid, aname, ymin) + call read_attribute_double_h5(gid, aname, ymin) case('ymax') - call read_attribute_double_h5(aid, aname, ymax) + call read_attribute_double_h5(gid, aname, ymax) case('zmin') - call read_attribute_double_h5(aid, aname, zmin) + call read_attribute_double_h5(gid, aname, zmin) case('zmax') - call read_attribute_double_h5(aid, aname, zmax) + call read_attribute_double_h5(gid, aname, zmax) case('nseeds') call read_attribute_integer_h5(gid, aname, lnseeds) @@ -1729,6 +1729,94 @@ module io ! !=============================================================================== ! +! subroutine READ_ATTRIBUTE_DOUBLE_H5: +! ----------------------------------- +! +! Subroutine reads a value of the double precision attribute provided by +! the group identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_attribute_double_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE + use hdf5 , only : h5aexists_by_name_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + real(kind=8) , intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_double_h5" & + , "Attribute does not exist :" // trim(aname)) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_DOUBLE, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot read attribute :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot close attribute :" // trim(aname)) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_attribute_double_h5 +! +!=============================================================================== +! ! subroutine WRITE_ATTRIBUTE_INTEGER_H5: ! ------------------------------------- ! @@ -2144,60 +2232,6 @@ module io ! !=============================================================================== ! -! read_attribute_double_h5: subroutine reads a double precision attribute -! -! arguments: -! aid - the HDF5 attribute identificator -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_double_h5(aid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - real(kind=8) , intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_DOUBLE, value, am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute -! - call print_error("io::read_attribute_double_h5" & - , "Cannot read attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine read_attribute_double_h5 -! -!=============================================================================== -! ! write_metablocks_h5: subroutine writes metablocks in the HDF5 format connected ! to the provided identificator ! From 07ca67549715aae7e1b5647546493c5695d443f7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 16:18:11 -0300 Subject: [PATCH 26/41] IO: Rewrite subroutine write_attribute_double_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 212 +++++++++++++++++++++++------------------------------ 1 file changed, 90 insertions(+), 122 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 21c390a..78034ae 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1907,6 +1907,96 @@ module io ! !=============================================================================== ! +! subroutine WRITE_ATTRIBUTE_DOUBLE_H5: +! ------------------------------------ +! +! Subroutine stores a value of the double precision attribute in the group +! provided by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine write_attribute_double_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + real(kind=8) , intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! create space for the attribute value +! + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_NATIVE_DOUBLE, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_NATIVE_DOUBLE, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot close attribute :" // trim(aname)) + end if + + else + call print_error("io::write_attribute_double_h5" & + , "Cannot create attribute :" // trim(aname)) + end if + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot close space for attribute :" // trim(aname)) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_double_h5 +! +!=============================================================================== +! ! write_attribute_vector_integer_h5: subroutine writes an vector intiger ! attribute in a group given by its indentificator ! @@ -2110,128 +2200,6 @@ module io ! !=============================================================================== ! -! write_attribute_double_h5: subroutine writes a double precision attribute in -! a group given by its indentificator -! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! value - the attribute value -! -!=============================================================================== -! - subroutine write_attribute_double_h5(gid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE - use hdf5 , only : h5screate_simple_f, h5sclose_f, h5acreate_f, h5awrite_f & - , h5aclose_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: gid - character(len=*), intent(in) :: name - real(kind=8) , intent(in) :: value - -! local variables -! - integer(hid_t) :: sid, aid - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! create space for the attribute -! - call h5screate_simple_f(1, am, sid, err) - -! check if the space has been created successfuly -! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_DOUBLE, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then - -! write the attribute data -! - call h5awrite_f(aid, H5T_NATIVE_DOUBLE, value, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_double_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if - -! close the attribute -! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_double_h5" & - , "Cannot close space for attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the space for the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot create space for attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine write_attribute_double_h5 -! -!=============================================================================== -! ! write_metablocks_h5: subroutine writes metablocks in the HDF5 format connected ! to the provided identificator ! From a169a55085f6de6cc3f0164578cb99ad45ac97cc Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 16:55:14 -0300 Subject: [PATCH 27/41] IO: Rewrite subroutine read_attribute_vector_integer_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 198 ++++++++++++++++++++++++++++++++--------------------- 1 file changed, 120 insertions(+), 78 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 78034ae..3bd712a 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1371,8 +1371,7 @@ module io ! store them in the current group ! - call read_attribute_vector_integer_h5(aid, aname, nseeds & - , seeds(:)) + call read_attribute_vector_integer_h5(gid, aname, seeds(:)) ! set the seed values ! @@ -1553,7 +1552,7 @@ module io ! obtain the block dimensions ! - call read_attribute_vector_integer_h5(aid, aname, 3, lm(:)) + call read_attribute_vector_integer_h5(gid, aname, lm(:)) case('nghost') @@ -1817,6 +1816,124 @@ module io ! !=============================================================================== ! +! subroutine READ_ATTRIBUTE_VECTOR_INTEGER_H5: +! ------------------------------------------- +! +! Subroutine reads a vector of the integer attribute provided by the group +! identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_attribute_vector_integer_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5aexists_by_name_f, h5aget_space_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + integer , dimension(:), intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid, sid + integer(hsize_t), dimension(1) :: am, bm + integer(hsize_t) :: alen + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Attribute does not exist :" // trim(aname)) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if + +! get the attribute space +! + call h5aget_space_f(aid, sid, ierr) + if (ierr == 0) then + call h5sget_simple_extent_dims_f(sid, am, bm, ierr) + if (ierr /= 1) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot get attribute dimensions :" // trim(aname)) + end if + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot close the attribute space :" // trim(aname)) + end if + else + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot get the attribute space :" // trim(aname)) + return + end if + +! check if the output array is large enough +! + if (am(1) > size(avalue)) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Attribute too large for output argument :" // trim(aname)) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_INTEGER, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot read attribute :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot close attribute :" // trim(aname)) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_attribute_vector_integer_h5 +! +!=============================================================================== +! ! subroutine WRITE_ATTRIBUTE_INTEGER_H5: ! ------------------------------------- ! @@ -2125,81 +2242,6 @@ module io ! !=============================================================================== ! -! read_attribute_vector_integer_h5: subroutine reads a vector of integer values -! from an attribute given by the identificator -! -! arguments: -! aid - the HDF5 attribute identificator -! name - the attribute name -! length - the vector length -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_vector_integer_h5(aid, name, length, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - integer , intent(in) :: length - integer, dimension(:), intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am - integer :: err -! -!------------------------------------------------------------------------------- -! -! check if the length is larger than 0 -! - if (length .gt. 0) then - -! prepare dimension array -! - am(1) = length - -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_INTEGER, value(:), am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute -! - call print_error("io::read_attribute_vector_integer_h5" & - , "Cannot read attribute :" // trim(name)) - - end if - - else ! length > 0 - -! print error about the wrong vector size -! - call print_error("io::read_attribute_vector_integer_h5" & - , "Wrong length of vector") - - end if ! length > 0 - -!------------------------------------------------------------------------------- -! - end subroutine read_attribute_vector_integer_h5 -! -!=============================================================================== -! ! write_metablocks_h5: subroutine writes metablocks in the HDF5 format connected ! to the provided identificator ! From 140daaea04e3e443e92c264d4efa89cef543f517 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 21 Apr 2014 17:03:26 -0300 Subject: [PATCH 28/41] IO: Rewrite subroutine write_attribute_vector_integer_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 142 ++++++++++++++++++++--------------------------------- 1 file changed, 54 insertions(+), 88 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 3bd712a..17abc7f 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1129,8 +1129,8 @@ module io ! store the vector attributes ! dm(:) = (/ in, jn, kn /) - call write_attribute_vector_integer_h5(gid, 'dims' , 3, dm(:)) - call write_attribute_vector_integer_h5(gid, 'rdims', 3, (/ ir, jr, kr /)) + call write_attribute_vector_integer_h5(gid, 'dims' , dm(:)) + call write_attribute_vector_integer_h5(gid, 'rdims', (/ ir, jr, kr /)) ! store random number generator seed values ! @@ -1146,7 +1146,7 @@ module io ! store them in the current group ! - call write_attribute_vector_integer_h5(gid, 'seeds', nseeds, seeds(:)) + call write_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) ! deallocate seed array ! @@ -2114,126 +2114,92 @@ module io ! !=============================================================================== ! -! write_attribute_vector_integer_h5: subroutine writes an vector intiger -! attribute in a group given by its indentificator +! subroutine WRITE_ATTRIBUTE_VECTOR_INTEGER_H5: +! -------------------------------------------- ! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! length - the vector length -! data - the attribute value +! Subroutine stores a vector of the integer attribute in the group provided +! by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute values; ! !=============================================================================== ! - subroutine write_attribute_vector_integer_h5(gid, name, length, data) + subroutine write_attribute_vector_integer_h5(gid, aname, avalue) -! references to other modules +! import external procedures and variables ! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5screate_simple_f, h5sclose_f & - , h5acreate_f, h5awrite_f, h5aclose_f + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! attribute arguments ! integer(hid_t) , intent(in) :: gid - character(len=*) , intent(in) :: name - integer , intent(in) :: length - integer(kind=4) , dimension(:), intent(in) :: data + character(len=*) , intent(in) :: aname + integer , dimension(:), intent(in) :: avalue ! local variables ! integer(hid_t) :: sid, aid - integer(hsize_t), dimension(1) :: am - integer :: err + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr ! !------------------------------------------------------------------------------- ! -! prepare the space dimensions +! set the proper attribute length ! - am(1) = length + am(1) = size(avalue) -! create space for the attribute +! create space for the attribute value ! - call h5screate_simple_f(1, am, sid, err) + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if -! check if the space has been created successfuly +! create the attribute in the given group ! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_INTEGER, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then + call h5acreate_f(gid, aname, H5T_NATIVE_INTEGER, sid, aid, ierr) + if (ierr == 0) then ! write the attribute data ! - call h5awrite_f(aid, H5T_NATIVE_INTEGER, data, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if + call h5awrite_f(aid, H5T_NATIVE_INTEGER, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if ! close the attribute ! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot close space for attribute :" // trim(name)) - + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot close attribute :" // trim(aname)) end if else + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot create attribute :" // trim(aname)) + end if -! print error about the problem with creating the space for the attribute +! release the space ! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot create space for attribute :" // trim(name)) - + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot close space for attribute :" // trim(aname)) end if !------------------------------------------------------------------------------- From 5180d03f2f4a0ee02ca766427c8f474b63b17568 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 22 Apr 2014 09:35:28 -0300 Subject: [PATCH 29/41] IO, MPITOOLS: Spell checking. Signed-off-by: Grzegorz Kowal --- src/io.F90 | 66 ++++++++++++++++++++++++------------------------ src/mpitools.F90 | 2 +- 2 files changed, 34 insertions(+), 34 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 17abc7f..bd2fd85 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -369,7 +369,7 @@ module io ! ------------------------- ! ! Subroutine stores block data in snapshots. Block variables are grouped -! todether and stored in big 4D arrays separately. This is a wrapper for +! together and stored in big 4D arrays separately. This is a wrapper for ! specific format storing. ! ! @@ -544,7 +544,7 @@ module io return end if -! opent the HDF5 file +! open the HDF5 file ! call h5fopen_f(fl, H5F_ACC_RDONLY_F, fid, err) @@ -1034,11 +1034,11 @@ module io ! ------------------------------ ! ! Subroutine stores global attributes in the HDF5 file provided by an -! identificator. +! identifier. ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -1175,12 +1175,12 @@ module io !=============================================================================== ! ! read_attributes_h5: subroutine restores attributes from an HDF5 file linked -! to the HDF5 file identificator +! to the HDF5 file identifier ! ! info: this subroutine restores only the global attributes ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! !=============================================================================== ! @@ -1457,10 +1457,10 @@ module io ! ! read_datablock_dims_h5: subroutine reads the data block dimensions from the ! attributes group of the file given by the file -! identificator +! identifier ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! dm - the data block dimensions ! !=============================================================================== @@ -2209,12 +2209,12 @@ module io !=============================================================================== ! ! write_metablocks_h5: subroutine writes metablocks in the HDF5 format connected -! to the provided identificator +! to the provided identifier ! ! info: this subroutine stores only the metablocks ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! !=============================================================================== ! @@ -2440,7 +2440,7 @@ module io ! info: this subroutine restores metablocks only ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! !=============================================================================== ! @@ -2684,7 +2684,7 @@ module io ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -2754,7 +2754,7 @@ module io dm(4) = jm dm(5) = km -! allocate arrays to store associated meta block identificators, conserved and +! allocate arrays to store associated meta block identifiers, conserved and ! primitive variables ! allocate(id(am(1))) @@ -2870,7 +2870,7 @@ module io ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3006,7 +3006,7 @@ module io ! info: this subroutine stores coordinates ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3168,11 +3168,11 @@ module io ! ! Subroutine groups each primitive variable from all data blocks and writes ! it as an array in the HDF5 dataset connected to the input HDF file -! identificator. +! identifier. ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3340,11 +3340,11 @@ module io ! ! Subroutine groups each conservative variable from all data blocks and writes ! it as an array in the HDF5 dataset connected to the input HDF file -! identificator. +! identifier. ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3510,7 +3510,7 @@ module io ! write_vector_integer_h5: subroutine stores a 1D integer vector in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! length - the vector length ! value - the data @@ -3637,7 +3637,7 @@ module io ! read_vector_integer_h5: subroutine reads a 1D integer vector ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! length - the vector length ! value - the data @@ -3726,7 +3726,7 @@ module io ! write_array2_integer_h5: subroutine stores a 2D integer array in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -3936,7 +3936,7 @@ module io ! read_array2_integer_h5: subroutine reads a 2D integer array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4025,7 +4025,7 @@ module io ! write_array4_integer_h5: subroutine stores a 4D integer array in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4235,7 +4235,7 @@ module io ! read_array4_integer_h5: subroutine reads a 4D integer array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4325,7 +4325,7 @@ module io ! a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! !=============================================================================== ! @@ -4449,7 +4449,7 @@ module io ! read_vector_double_h5: subroutine reads a 1D double precision vector ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the vector dimensions ! value - the data @@ -4538,7 +4538,7 @@ module io ! write_array4_float_h5: subroutine stores a 4D single precision array ! ! arguments: -! gid - the HDF5 group identificator where the dataset should be located +! gid - the HDF5 group identifier where the dataset should be located ! name - the string name representing the dataset ! dm - the dataset dimensions ! value - the dataset values @@ -4692,7 +4692,7 @@ module io ! write_array3_double_h5: subroutine stores a 3D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4902,7 +4902,7 @@ module io ! write_array4_double_h5: subroutine stores a 4D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5112,7 +5112,7 @@ module io ! write_array5_double_h5: subroutine stores a 5D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5322,7 +5322,7 @@ module io ! read_array5_double_h5: subroutine reads a 5D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5411,7 +5411,7 @@ module io ! write_array6_double_h5: subroutine stores a 6D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data diff --git a/src/mpitools.F90 b/src/mpitools.F90 index 0d137bc..93e74f5 100644 --- a/src/mpitools.F90 +++ b/src/mpitools.F90 @@ -141,7 +141,7 @@ module mpitools stop end if -! obtain the current process identificator +! obtain the current process identifier ! call mpi_comm_rank(mpi_comm_world, nproc , iret) From 1e6b6f71e20b93546253b80c85e51e1bb9bea401 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 23 Apr 2014 14:39:01 -0300 Subject: [PATCH 30/41] IO: Rewrite subroutine read_attributes_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 347 ++++++++++++++++++++--------------------------------- 1 file changed, 133 insertions(+), 214 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index bd2fd85..7e6c71a 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1174,57 +1174,55 @@ module io ! !=============================================================================== ! -! read_attributes_h5: subroutine restores attributes from an HDF5 file linked -! to the HDF5 file identifier +! subroutine READ_ATTRIBUTES_H5: +! ----------------------------- ! -! info: this subroutine restores only the global attributes +! Subroutine restores global attributes from an HDF5 file provided by its +! identifier. ! -! arguments: -! fid - the HDF5 file identifier +! Arguments: +! +! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine read_attributes_h5(fid) -! references to other modules +! import external procedures and variables ! - use blocks , only : block_meta, block_data - use blocks , only : append_metablock - use blocks , only : set_last_id, get_last_id, get_mblocks, get_dblocks & - , get_nleafs - use coordinates, only : nn, ng, in, jn, kn, maxlev, toplev, ir, jr, kr - use coordinates, only : initialize_coordinates, finalize_coordinates - use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax - use error , only : print_error - use evolution, only : step, time, dt, dtn - use hdf5 , only : hid_t, hsize_t - use hdf5 , only : h5gopen_f, h5gclose_f, h5aget_num_attrs_f & - , h5aopen_idx_f, h5aclose_f, h5aget_name_f - use mpitools , only : nprocs, nproc - use random , only : nseeds, set_seeds + use blocks , only : block_meta + use blocks , only : append_metablock + use blocks , only : set_last_id, get_last_id + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use coordinates , only : nn, ng, in, jn, kn, ir, jr, kr + use coordinates , only : maxlev, toplev + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use coordinates , only : initialize_coordinates, finalize_coordinates + use error , only : print_error + use evolution , only : step, time, dt, dtn + use hdf5 , only : hid_t, hsize_t + use hdf5 , only : h5gopen_f, h5gclose_f + use mpitools , only : nprocs, nproc + use random , only : nseeds, set_seeds -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! - character(len=16) :: aname - integer(hid_t) :: gid, aid - integer(hsize_t) :: alen = 16 - integer(kind=4) :: dm(3) - integer :: err, i, l - integer :: nattrs, lndims, llast_id, lmblocks, lnleafs & - , lncells, lnghost, lnseeds, lmaxlev, lnproc + integer(hid_t) :: gid + integer :: ierr, l + integer :: lndims, lmaxlev, lmblocks, lnleafs, llast_id + integer :: lncells, lnghost, lnproc, lnseeds ! local pointers ! type(block_meta), pointer :: pmeta - type(block_data), pointer :: pdata ! allocatable arrays ! @@ -1234,221 +1232,142 @@ module io ! ! open the global attributes group ! - call h5gopen_f(fid, 'attributes', gid, err) + call h5gopen_f(fid, 'attributes', gid, ierr) ! check if the group has been opened successfuly ! - if (err .ge. 0) then + if (ierr < 0) then + call print_error("io::read_attributes_h5", "Cannot open the group!") + return + end if -! read the number of global attributes +! restore integer attributes ! - call h5aget_num_attrs_f(gid, nattrs, err) + call read_attribute_integer_h5(gid, 'ndims' , lndims ) + call read_attribute_integer_h5(gid, 'maxlev' , lmaxlev ) + call read_attribute_integer_h5(gid, 'nprocs' , nfiles ) + call read_attribute_integer_h5(gid, 'nproc' , lnproc ) + call read_attribute_integer_h5(gid, 'mblocks', lmblocks) + call read_attribute_integer_h5(gid, 'nleafs' , lnleafs ) + call read_attribute_integer_h5(gid, 'last_id', llast_id) + call read_attribute_integer_h5(gid, 'ncells' , lncells ) + call read_attribute_integer_h5(gid, 'nghost' , lnghost ) + call read_attribute_integer_h5(gid, 'nseeds' , lnseeds ) + call read_attribute_integer_h5(gid, 'step' , step ) + call read_attribute_integer_h5(gid, 'isnap' , isnap ) -! check if the number of attributes has been read successfuly +! restore double precision attributes ! - if (err .ge. 0) then + call read_attribute_double_h5(gid, 'xmin', xmin) + call read_attribute_double_h5(gid, 'xmax', xmax) + call read_attribute_double_h5(gid, 'ymin', ymin) + call read_attribute_double_h5(gid, 'ymax', ymax) + call read_attribute_double_h5(gid, 'zmin', zmin) + call read_attribute_double_h5(gid, 'zmax', zmax) + call read_attribute_double_h5(gid, 'time', time) + call read_attribute_double_h5(gid, 'dt' , dt ) + call read_attribute_double_h5(gid, 'dtn' , dtn ) -! iterate over all attributes +! check the number of dimensions ! - do i = 0, nattrs - 1 + if (lndims /= NDIMS) then + call print_error("io::read_attributes_h5" & + , "The number of dimensions does not match!") + return + end if -! open the current attribute +! check the block dimensions ! - call h5aopen_idx_f(gid, i, aid, err) + if (lncells /= nn) then + call print_error("io::read_attributes_h5" & + , "The block dimensions do not match!") + end if -! check if the attribute has been opened successfuly +! check the number of ghost layers ! - if (err .ge. 0) then + if (lnghost /= ng) then + call print_error("io::read_attributes_h5" & + , "The number of ghost layers does not match!") + end if -! obtain the attribute name +! prepare coordinates and rescaling factors if the maximum level has changed ! - call h5aget_name_f(aid, alen, aname, err) - -! depending on the attribute name use proper subroutine to read its value -! - select case(trim(aname)) - case('ndims') - call read_attribute_integer_h5(gid, aname, lndims) - -! check if the restart file and compiled program have the same number of -! dimensions -! - if (lndims .ne. NDIMS) then - call print_error("io::read_attributes_h5" & - , "File and program dimensions are incompatible!") - end if - case('maxlev') - call read_attribute_integer_h5(gid, aname, lmaxlev) - if (lmaxlev .gt. toplev) then + if (lmaxlev > toplev) then ! subtitute the new value of toplev ! - toplev = lmaxlev + toplev = lmaxlev ! regenerate coordinates ! - call finalize_coordinates(err) - call initialize_coordinates(.false., err) + call finalize_coordinates(ierr) + call initialize_coordinates(.false., ierr) ! calculate a factor to rescale the block coordinates ! - dcor = 2**(toplev - maxlev) - - else - -! calculate a factor to rescale the block coordinates -! - ucor = 2**(maxlev - lmaxlev) - end if - case('nprocs') - call read_attribute_integer_h5(gid, aname, nfiles) - case('nproc') - call read_attribute_integer_h5(gid, aname, lnproc) - case('last_id') - call read_attribute_integer_h5(gid, aname, llast_id) - case('mblocks') - call read_attribute_integer_h5(gid, aname, lmblocks) - case('nleafs') - call read_attribute_integer_h5(gid, aname, lnleafs) - case('ncells') - call read_attribute_integer_h5(gid, aname, lncells) - -! check if the block dimensions are compatible -! - if (lncells .ne. nn) then - call print_error("io::read_attributes_h5" & - , "File and program block dimensions are incompatible!") - end if - case('nghost') - call read_attribute_integer_h5(gid, aname, lnghost) - -! check if the ghost layers are compatible -! - if (lnghost .ne. ng) then - call print_error("io::read_attributes_h5" & - , "File and program block ghost layers are incompatible!") - end if - case('step') - call read_attribute_integer_h5(gid, aname, step) - case('isnap') - call read_attribute_integer_h5(gid, aname, isnap) - case('time') - call read_attribute_double_h5(gid, aname, time) - case('dt') - call read_attribute_double_h5(gid, aname, dt) - case('dtn') - call read_attribute_double_h5(gid, aname, dtn) - case('xmin') - call read_attribute_double_h5(gid, aname, xmin) - case('xmax') - call read_attribute_double_h5(gid, aname, xmax) - case('ymin') - call read_attribute_double_h5(gid, aname, ymin) - case('ymax') - call read_attribute_double_h5(gid, aname, ymax) - case('zmin') - call read_attribute_double_h5(gid, aname, zmin) - case('zmax') - call read_attribute_double_h5(gid, aname, zmax) - case('nseeds') - call read_attribute_integer_h5(gid, aname, lnseeds) - -! ! check if the numbers of seeds are compatible -! ! -! if (lnseeds .ne. nseeds) then -! call print_error("io::read_attributes_h5" & -! , "The number of seeds from file and program are incompatible!") -! end if - case('seeds') - -! check if the numbers of seeds are compatible -! - if (lnseeds .eq. nseeds) then - -! allocate space for seeds -! - allocate(seeds(nseeds)) - -! store them in the current group -! - call read_attribute_vector_integer_h5(gid, aname, seeds(:)) - -! set the seed values -! - call set_seeds(nseeds, seeds(:)) - -! deallocate seed array -! - deallocate(seeds) - - end if - case default - end select - -! close the current attribute -! - call h5aclose_f(aid, err) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_attributes_h5", & - "Cannot open the current attribute!") - - end if - - end do - -! allocate all metablocks -! - do l = 1, lmblocks - call append_metablock(pmeta) - end do - -! check if the number of created metablocks is equal to lbmcloks -! - if (lmblocks .ne. get_mblocks()) then - call print_error("io::read_attributes_h5" & - , "Number of metablocks doesn't match!") - end if - -! allocate an array of pointers with the size llast_id -! - allocate(block_array(llast_id)) - call set_last_id(llast_id) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_attributes_h5", & - "Cannot get the number of global attributes!") - - end if - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_attributes_h5", "Cannot close the group!") - - end if + dcor = 2**(toplev - maxlev) else -! print error about the problem with creating the group +! calculate a factor to rescale the block coordinates ! - call print_error("io::read_attributes_h5", "Cannot open the group!") + ucor = 2**(maxlev - lmaxlev) end if +! check if the numbers of seeds are compatible +! + if (lnseeds == nseeds) then + +! allocate space for seeds +! + allocate(seeds(nseeds)) + +! store them in the current group +! + call read_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) + +! set the seed values +! + call set_seeds(nseeds, seeds(:)) + +! deallocate seed array +! + deallocate(seeds) + + end if + +! allocate all metablocks +! + do l = 1, lmblocks + call append_metablock(pmeta) + end do + +! check if the number of created metablocks is equal to lbmcloks +! + if (lmblocks /= get_mblocks()) then + call print_error("io::read_attributes_h5" & + , "Number of metablocks does not match!") + end if + +! allocate an array of pointers with the size llast_id +! + allocate(block_array(llast_id)) + +! set the last_id +! + call set_last_id(llast_id) + +! close the group +! + call h5gclose_f(gid, ierr) + +! check if the group has been closed successfuly +! + if (ierr /= 0) then + call print_error("io::read_attributes_h5", "Cannot close the group!") + end if + !------------------------------------------------------------------------------- ! end subroutine read_attributes_h5 From 64d360ad3d0da0e63e0af5483aa18523bc968110 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 23 Apr 2014 14:56:23 -0300 Subject: [PATCH 31/41] RANDOM: Allow for seed expanding or shrinking in set_seeds(). Signed-off-by: Grzegorz Kowal --- src/random.F90 | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/random.F90 b/src/random.F90 index 125112f..cf95008 100644 --- a/src/random.F90 +++ b/src/random.F90 @@ -214,6 +214,11 @@ module random ! integer , intent(in) :: np integer(kind=4), dimension(0:np-1), intent(in) :: seed + +! local variables +! + integer :: i, l + real :: r ! !------------------------------------------------------------------------------- ! @@ -225,10 +230,34 @@ module random ! set the seeds only if the input array and seeds have the same sizes ! - if (np .eq. nseeds) then + if (np == nseeds) then seeds(0:lseed) = seed(0:lseed) + else + +! if the input array and seeds have different sizes, expand or shrink seeds +! + select case(gentype) + case('random') + l = min(lseed, np - 1) + seeds(0:l) = seed(0:l) + if (l < lseed) then + do i = l + 1, lseed + call random_number(r) + seeds(i) = 123456789 * r + end do + end if + case default + l = nseeds / 2 + do i = 0, l - 1 + seeds(i) = seed(0) + end do + do i = l, lseed + seeds(i) = seed(np-1) + end do + end select + end if #ifdef PROFILE From 14727d3afbd79454909a15e960194085b6eaea89 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 23 Apr 2014 15:02:56 -0300 Subject: [PATCH 32/41] IO: Restore seeds for all cases in read_attributes_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 7e6c71a..9478b17 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1315,27 +1315,21 @@ module io end if -! check if the numbers of seeds are compatible -! - if (lnseeds == nseeds) then - ! allocate space for seeds ! - allocate(seeds(nseeds)) + allocate(seeds(lnseeds)) ! store them in the current group ! - call read_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) + call read_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) ! set the seed values ! - call set_seeds(nseeds, seeds(:)) + call set_seeds(lnseeds, seeds(:)) ! deallocate seed array ! - deallocate(seeds) - - end if + deallocate(seeds) ! allocate all metablocks ! From 962c5fce94cdbf5771729b0462371e2163d496d4 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 23 Apr 2014 15:08:58 -0300 Subject: [PATCH 33/41] MESH: Remove calling check_mesh() if DEBUG=Y. Subroutine check_mesh() has been removed. Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 0625515..218d879 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -819,12 +819,6 @@ module mesh call start_timer(imu) #endif /* PROFILE */ -#ifdef DEBUG -! check the mesh when debugging -! - call check_mesh('before update_mesh') -#endif /* DEBUG */ - !! DETERMINE THE REFINEMENT OF ALL DATA BLOCKS !! ! set the pointer to the first block on the data block list @@ -1263,12 +1257,6 @@ module mesh call redistribute_blocks() #endif /* MPI */ -#ifdef DEBUG -! check mesh -! - call check_mesh('after update_mesh') -#endif /* DEBUG */ - #ifdef PROFILE ! stop accounting time for the adaptive mesh refinement update ! From 75eb97125e479f0ceba773b9de11e4d229374c01 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 23 Apr 2014 17:33:44 -0300 Subject: [PATCH 34/41] BLOCKS, MESH, IO: Fix restart with less processes. The job restart with smaller number of processors than the number of files didn't work well, since all metablocks of all remaining files where set to the last process before actually reading the corresponding data blocks. This resulted in a crash in redistribute_blocks() due to pointer not associated. The fix lets the meta blocks to be restored with the original process numbers. While reading files with the indices larger than the index of the last process, all meta blocks which corresponding data blocks are stored in the read file, are set to the process number of the reading process. After reading data blocks from each file, the blocks are redistributed. However, only data blocks belonging to the active processes are redistributed. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 53 +++++++++++++++++++++++++++++++++++++++++++++ src/io.F90 | 7 +++++- src/mesh.F90 | 58 ++++++++++++++++++++++++++++---------------------- 3 files changed, 91 insertions(+), 27 deletions(-) diff --git a/src/blocks.F90 b/src/blocks.F90 index 0a47ec4..5b77e6f 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -271,6 +271,7 @@ module blocks public :: refine_block, derefine_block public :: set_last_id, get_last_id, get_mblocks, get_dblocks, get_nleafs public :: set_blocks_update + public :: change_blocks_process public :: set_neighbors_refine public :: metablock_set_id, metablock_set_process, metablock_set_level public :: metablock_set_configuration, metablock_set_refinement @@ -2152,6 +2153,58 @@ module blocks ! !=============================================================================== ! +! subroutine CHANGE_BLOCKS_PROCESS: +! -------------------------------- +! +! Subroutine switches meta blocks which belong to old process to the new one. +! +! Arguments: +! +! npold - the old process number; +! npnew - the new process number; +! +!=============================================================================== +! + subroutine change_blocks_process(npold, npnew) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer, intent(in) :: npold, npnew + +! local pointers +! + type(block_meta), pointer :: pmeta +! +!------------------------------------------------------------------------------- +! +! associate the pointer with the first block on the meta block list +! + pmeta => list_meta + +! iterate over all blocks in the list +! + do while(associated(pmeta)) + +! if the meta block belongs to process npold, switch it to process npnew +! + if (pmeta%process == npold) pmeta%process = npnew + +! associate the pointer with the next block on the list +! + pmeta => pmeta%next + + end do ! meta blocks + +!------------------------------------------------------------------------------- +! + end subroutine change_blocks_process +! +!=============================================================================== +! ! subroutine METABLOCK_SET_ID: ! --------------------------- ! diff --git a/src/io.F90 b/src/io.F90 index 9478b17..df8d9ed 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -469,6 +469,7 @@ module io ! import external procedures and variables ! + use blocks , only : change_blocks_process use error , only : print_error use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F @@ -669,6 +670,10 @@ module io ! do lfile = nprocs, nfiles - 1 +! switch meta blocks from the read file to belong to the reading process +! + call change_blocks_process(lfile, nprocs - 1) + ! read the remaining files by the last process only ! if (nproc == nprocs - 1) then @@ -2497,7 +2502,7 @@ module io block_array(id(l))%ptr => pmeta call metablock_set_id (pmeta, id (l)) - call metablock_set_process (pmeta, min(lcpu, cpu(l))) + call metablock_set_process (pmeta, cpu(l)) call metablock_set_refinement (pmeta, ref(l)) call metablock_set_configuration(pmeta, cfg(l)) call metablock_set_level (pmeta, lev(l)) diff --git a/src/mesh.F90 b/src/mesh.F90 index 218d879..d640bf9 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1347,92 +1347,98 @@ module mesh ! do while (associated(pmeta)) +! consider only meta blocks which belong to active processes +! + if (pmeta%process < nprocs) then + ! check if the block belongs to another process ! - if (pmeta%process /= np) then + if (pmeta%process /= np) then ! check if the block is the leaf ! - if (pmeta%leaf) then + if (pmeta%leaf) then ! generate a tag for communication ! - itag = pmeta%process * nprocs + np + nprocs + 1 + itag = pmeta%process * nprocs + np + nprocs + 1 ! sends the block to the right process ! - if (nproc == pmeta%process) then + if (nproc == pmeta%process) then ! copy data to buffer ! - rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:) - rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:) + rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:) + rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:) ! send data ! - call send_real_array(size(rbuf), np, itag, rbuf, iret) + call send_real_array(size(rbuf), np, itag, rbuf, iret) ! remove data block from the current process ! - call remove_datablock(pmeta%data) + call remove_datablock(pmeta%data) ! send data block ! - end if ! nproc == pmeta%process + end if ! nproc == pmeta%process ! receive the block from another process ! - if (nproc == np) then + if (nproc == np) then ! allocate a new data block and link it with the current meta block ! - call append_datablock(pdata) - call link_blocks(pmeta, pdata) + call append_datablock(pdata) + call link_blocks(pmeta, pdata) ! receive the data ! - call receive_real_array(size(rbuf), pmeta%process, itag, rbuf, iret) + call receive_real_array(size(rbuf), pmeta%process, itag, rbuf, iret) ! coppy the buffer to data block ! - pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:) - pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:) + pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:) + pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:) - end if ! nproc == n + end if ! nproc == n - end if ! leaf + end if ! leaf ! set new processor number ! - pmeta%process = np + pmeta%process = np - end if ! pmeta%process /= np + end if ! pmeta%process /= np ! increase the number of blocks on the current process; if it exceeds the ! allowed number reset the counter and increase the processor number ! - if (pmeta%leaf) then + if (pmeta%leaf) then ! increase the number of leafs for the current process ! - nl = nl + 1 + nl = nl + 1 ! if the number of leafs for the current process exceeds the number of assigned ! blocks, reset the counter and increase the process number ! - if (nl >= lb(np)) then + if (nl >= lb(np)) then ! reset the leaf counter for the current process ! - nl = 0 + nl = 0 ! increase the process number ! - np = min(nprocs - 1, np + 1) + np = min(nprocs - 1, np + 1) - end if ! l >= lb(n) + end if ! l >= lb(n) - end if ! leaf + end if ! leaf + + end if ! pmeta%process < nprocs ! assign the pointer to the next meta block ! From aafbca8203e4fd97a4f63ee913a4f4aad8573961 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 24 Apr 2014 13:11:25 -0300 Subject: [PATCH 35/41] IO: Reduce local variables in write_attributes_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index df8d9ed..7decd7a 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1075,7 +1075,6 @@ module io ! integer(hid_t) :: gid integer :: err - integer(kind=4), dimension(3) :: dm ! local allocatable arrays ! @@ -1133,8 +1132,7 @@ module io ! store the vector attributes ! - dm(:) = (/ in, jn, kn /) - call write_attribute_vector_integer_h5(gid, 'dims' , dm(:)) + call write_attribute_vector_integer_h5(gid, 'dims' , (/ in, jn, kn /)) call write_attribute_vector_integer_h5(gid, 'rdims', (/ ir, jr, kr /)) ! store random number generator seed values From e8c2fac09f63c5005f245981f411cbe75d50d64d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 24 Apr 2014 15:38:03 -0300 Subject: [PATCH 36/41] IO: Rewrite subroutine read_datablocks_h5(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 290 +++++++++++------------------------------------------ 1 file changed, 59 insertions(+), 231 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 7decd7a..d9948c3 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -1370,185 +1370,6 @@ module io end subroutine read_attributes_h5 ! !=============================================================================== -! -! read_datablock_dims_h5: subroutine reads the data block dimensions from the -! attributes group of the file given by the file -! identifier -! -! arguments: -! fid - the HDF5 file identifier -! dm - the data block dimensions -! -!=============================================================================== -! - subroutine read_datablock_dims_h5(fid, dm) - -! references to other modules -! - use error , only : print_error - use hdf5 , only : hid_t, hsize_t - use hdf5 , only : h5gopen_f, h5gclose_f, h5aget_num_attrs_f & - , h5aopen_idx_f, h5aclose_f, h5aget_name_f - use equations, only : nv - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: fid - integer(hsize_t), dimension(5), intent(out) :: dm - -! local variables -! - character(len=16) :: aname - integer(hid_t) :: gid, aid - integer(hsize_t) :: alen = 16 - integer :: err, i - integer :: nattrs, ldblocks, lnghost - -! local arrays -! - integer(kind=4), dimension(3) :: lm -! -!------------------------------------------------------------------------------- -! -! initiate the output vector -! - dm(:) = 0 - dm(2) = nv - -! open the global attributes group -! - call h5gopen_f(fid, 'attributes', gid, err) - -! check if the group has been opened successfuly -! - if (err .ge. 0) then - -! read the number of global attributes -! - call h5aget_num_attrs_f(gid, nattrs, err) - -! check if the number of attributes has been read successfuly -! - if (err .ge. 0) then - -! iterate over all attributes -! - do i = 0, nattrs - 1 - -! open the current attribute -! - call h5aopen_idx_f(gid, i, aid, err) - -! check if the attribute has been opened successfuly -! - if (err .ge. 0) then - -! obtain the attribute name -! - call h5aget_name_f(aid, alen, aname, err) - -! check if the attribute name has been read successfuly -! - if (err .ge. 0) then - -! depending on the attribute name use proper subroutine to read its value -! - select case(trim(aname)) - case('dblocks') - -! obtain the number of data blocks -! - call read_attribute_integer_h5(gid, aname, ldblocks) - - case('dims') - -! obtain the block dimensions -! - call read_attribute_vector_integer_h5(gid, aname, lm(:)) - - case('nghost') - -! obtain the number of data blocks -! - call read_attribute_integer_h5(gid, aname, lnghost) - - case default - end select - - else - -! print error about the problem with reading the attribute name -! - call print_error("io::read_datablock_dims_h5", & - "Cannot read the current attribute name!") - - end if - -! close the current attribute -! - call h5aclose_f(aid, err) - - else - -! print error about the problem with opening the current attribute -! - call print_error("io::read_datablock_dims_h5", & - "Cannot open the current attribute!") - - end if - - end do ! i = 0, nattrs - 1 - -! prepare the output array -! - dm(1) = ldblocks - do i = 1, 3 - if (lm(i) .gt. 1) lm(i) = lm(i) + 2 * lnghost - end do - dm(3:5) = lm(1:3) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_datablock_dims_h5", & - "Cannot get the number of global attributes!") - - end if - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_datablock_dims_h5" & - , "Cannot close the attributes group!") - - end if - - else - -! print error about the problem with creating the group -! - call print_error("io::read_datablock_dims_h5" & - , "Cannot open the attributes group!") - - end if - -!------------------------------------------------------------------------------- -! - end subroutine read_datablock_dims_h5 -! -!=============================================================================== !! !!--- ATTRIBUTE SUBROUTINES -------------------------------------------------- !! @@ -2797,6 +2618,7 @@ module io use blocks , only : block_meta, block_data, list_data use blocks , only : append_datablock, link_blocks use coordinates , only : im, jm, km + use equations , only : nv use error , only : print_error use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gopen_f, h5gclose_f @@ -2817,7 +2639,7 @@ module io ! integer(hid_t) :: gid integer(kind=4) :: l - integer :: err + integer :: dblocks, ierr ! local arrays ! @@ -2830,84 +2652,90 @@ module io ! !------------------------------------------------------------------------------- ! -! get datablock array dimensions +! read the number of data blocks ! - call read_datablock_dims_h5(fid, dm(:)) - -! open the group 'datablocks' -! - call h5gopen_f(fid, 'datablocks', gid, err) - -! check if the datablock group has been opened successfuly -! - if (err >= 0) then + call h5gopen_f(fid, 'attributes', gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot open the attribute group!") + return + end if + call read_attribute_integer_h5(gid, 'dblocks', dblocks) + call h5gclose_f(gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot close the attribute group!") + return + end if ! restore data blocks only if there are any ! - if (dm(1) > 0) then + if (dblocks > 0) then + +! fill out dimensions dm(:) +! + dm(1) = dblocks + dm(2) = nv + dm(3) = im + dm(4) = jm + dm(5) = km ! allocate arrays to read data ! - allocate(id(dm(1))) - allocate(uv(dm(1),dm(2),dm(3),dm(4),dm(5))) - allocate(qv(dm(1),dm(2),dm(3),dm(4),dm(5))) + allocate(id(dm(1))) + allocate(uv(dm(1),dm(2),dm(3),dm(4),dm(5))) + allocate(qv(dm(1),dm(2),dm(3),dm(4),dm(5))) + +! open the group 'datablocks' +! + call h5gopen_f(fid, 'datablocks', gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot open the data block group!") + return + end if ! read array data from the HDF5 file ! - call read_vector_integer_h5(gid, 'meta', dm(1:1), id(:) ) - call read_array5_double_h5 (gid, 'uvar', dm(1:5), uv(:,:,:,:,:)) - call read_array5_double_h5 (gid, 'qvar', dm(1:5), qv(:,:,:,:,:)) + call read_vector_integer_h5(gid, 'meta', dm(1:1), id(:) ) + call read_array5_double_h5 (gid, 'uvar', dm(1:5), uv(:,:,:,:,:)) + call read_array5_double_h5 (gid, 'qvar', dm(1:5), qv(:,:,:,:,:)) + +! close the data block group +! + call h5gclose_f(gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot close the data block group!") + return + end if ! iterate over data blocks, allocate them and fill out their fields ! - do l = 1, dm(1) + do l = 1, dm(1) ! allocate and append to the end of the list a new datablock ! - call append_datablock(pdata) + call append_datablock(pdata) ! associate the corresponding meta block with the current data block ! - call link_blocks(block_array(id(l))%ptr, pdata) + call link_blocks(block_array(id(l))%ptr, pdata) ! fill out the array of conservative and primitive variables ! - pdata%u(:,:,:,:) = uv(l,:,:,:,:) - pdata%q(:,:,:,:) = qv(l,:,:,:,:) + pdata%u(:,:,:,:) = uv(l,:,:,:,:) + pdata%q(:,:,:,:) = qv(l,:,:,:,:) - end do ! l = 1, dm(1) + end do ! l = 1, dm(1) ! deallocate allocatable arrays ! - if (allocated(id)) deallocate(id) - if (allocated(uv)) deallocate(uv) - if (allocated(qv)) deallocate(qv) + if (allocated(id)) deallocate(id) + if (allocated(uv)) deallocate(uv) + if (allocated(qv)) deallocate(qv) - end if ! dblocks > 0 - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err > 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_datablocks_h5" & - , "Cannot close data block group!") - - end if - - else - -! print error about the problem with opening the group -! - call print_error("io::read_datablocks_h5" & - , "Cannot open data block group!") - - end if + end if ! dblocks > 0 !------------------------------------------------------------------------------- ! From 9edf7111b5e1b78bc82bd4918c8279767f74594f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 25 Apr 2014 12:45:00 -0300 Subject: [PATCH 37/41] BLOCKS: Iterate over corners from all edges in 2D iterate_over_neighbors(). In 2D case, if we update corners using only X-edge neighbors, we get some rare problems with symmetry. If we additionally use Y-edge neighbors, to update corners, those problems disappear. This should be additionally tested for both, 2D and 3D cases. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/blocks.F90 b/src/blocks.F90 index 5b77e6f..be3b983 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -3521,32 +3521,52 @@ module blocks end if end if -! (0,0)-(½,0) edge +! (0,0)-(½,0) edge and (0,0) corner ! pneigh => pmeta%neigh(2,1,1)%ptr if (associated(pneigh)) then call pprocedure(pmeta, pneigh) + + pneigh => pneigh%neigh(1,1,2)%ptr + if (associated(pneigh)) then + call pprocedure(pmeta, pneigh) + end if end if -! (½,0)-(1,0) edge +! (½,0)-(1,0) edge and (1,0) corner ! pneigh => pmeta%neigh(2,1,2)%ptr if (associated(pneigh)) then call pprocedure(pmeta, pneigh) + + pneigh => pneigh%neigh(1,2,2)%ptr + if (associated(pneigh)) then + call pprocedure(pmeta, pneigh) + end if end if -! (0,1)-(½,1) edge +! (0,1)-(½,1) edge and (0,1) corner ! pneigh => pmeta%neigh(2,2,1)%ptr if (associated(pneigh)) then call pprocedure(pmeta, pneigh) + + pneigh => pneigh%neigh(1,1,1)%ptr + if (associated(pneigh)) then + call pprocedure(pmeta, pneigh) + end if end if -! (½,1)-(1,1) edge +! (½,1)-(1,1) edge and (1,1) corner ! pneigh => pmeta%neigh(2,2,2)%ptr if (associated(pneigh)) then call pprocedure(pmeta, pneigh) + + pneigh => pneigh%neigh(1,2,1)%ptr + if (associated(pneigh)) then + call pprocedure(pmeta, pneigh) + end if end if #endif /* NDIMS == 3 */ From 91a4983660c3b42d56af3d9e6b2840310ae7c829 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 26 Apr 2014 13:36:29 -0300 Subject: [PATCH 38/41] BOUNDARIES: Fix asymmetries in flux boundary update. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index f59ac04..0a3f837 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -3628,8 +3628,8 @@ module boundaries k1 = 2 * (k - kl) + kb k2 = k1 + 1 - pdata%f(idir,:,it,j,k) = 2.5d-01 * (f(:,j1,k1) + f(:,j2,k1) & - + f(:,j1,k2) + f(:,j2,k2)) + pdata%f(idir,:,it,j,k) = 2.5d-01 * ((f(:,j1,k1) + f(:,j2,k2)) & + + (f(:,j2,k1) + f(:,j1,k2))) end do #endif /* NDIMS == 3 */ end do @@ -3676,8 +3676,8 @@ module boundaries k1 = 2 * (k - kl) + kb k2 = k1 + 1 - pdata%f(idir,:,i,jt,k) = 2.5d-01 * (f(:,i1,k1) + f(:,i2,k1) & - + f(:,i1,k2) + f(:,i2,k2)) + pdata%f(idir,:,i,jt,k) = 2.5d-01 * ((f(:,i1,k1) + f(:,i2,k2)) & + + (f(:,i2,k1) + f(:,i1,k2))) end do #endif /* NDIMS == 3 */ end do @@ -3717,8 +3717,8 @@ module boundaries j1 = 2 * (j - jl) + jb j2 = j1 + 1 - pdata%f(idir,:,i,j,kt) = 2.5d-01 * (f(:,i1,j1) + f(:,i2,j1) & - + f(:,i1,j2) + f(:,i2,j2)) + pdata%f(idir,:,i,j,kt) = 2.5d-01 * ((f(:,i1,j1) + f(:,i2,j2)) & + + (f(:,i2,j1) + f(:,i1,j2))) end do end do #endif /* NDIMS == 3 */ From e06a7edfdd0c77f176bc59f2def9ad473b1e9837 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 27 Apr 2014 21:00:08 -0300 Subject: [PATCH 39/41] BLOCKS: Check neighbor consistency if DEBUG=Y. Add subroutines to iterate over all block neighbors and check if they reference to correct meta blocks. The neighbor consistency check is done after initial mesh generation and each time the mesh has been updated. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 207 +++++++++++++++++++++++++++++++++++++++++++++++++ src/mesh.F90 | 18 +++++ 2 files changed, 225 insertions(+) diff --git a/src/blocks.F90 b/src/blocks.F90 index be3b983..c6533d0 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -277,6 +277,9 @@ module blocks public :: metablock_set_configuration, metablock_set_refinement public :: metablock_set_position, metablock_set_coordinates public :: metablock_set_bounds, metablock_set_leaf, metablock_unset_leaf +#ifdef DEBUG + public :: check_neighbors +#endif /* DEBUG */ !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -2711,6 +2714,52 @@ module blocks !------------------------------------------------------------------------------- ! end subroutine set_neighbors_refine +#ifdef DEBUG +! +!=============================================================================== +! +! subroutine CHECK_NEIGHBORS: +! -------------------------- +! +! Subroutine iterates over all blocks and checks if the pointers to their +! neighbors are consistent. +! +!=============================================================================== +! + subroutine check_neighbors() + +! local variables are not implicit by default +! + implicit none + +! local pointers +! + type(block_meta), pointer :: pmeta +! +!------------------------------------------------------------------------------- +! +! associate the pointer with the first block on the meta block list +! + pmeta => list_meta + +! iterate over all blocks in the list +! + do while(associated(pmeta)) + +! check the block neighbors +! + call check_block_neighbors(pmeta) + +! associate the pointer with the next block on the list +! + pmeta => pmeta%next + + end do ! meta blocks + +!------------------------------------------------------------------------------- +! + end subroutine check_neighbors +#endif /* DEBUG */ ! !=============================================================================== !! @@ -3573,6 +3622,164 @@ module blocks !------------------------------------------------------------------------------- ! end subroutine iterate_over_neighbors +#ifdef DEBUG +! +!=============================================================================== +! +! subroutine CHECK_BLOCK_NEIGHBORS: +! -------------------------------- +! +! Subroutine iterates over all neighbor blocks and checks if their pointers +! are consistent, i.e. if their corresponding neighbor pointers point to +! the current block. +! +! Arguments: +! +! pmeta - a pointer to the meta block; +! +!=============================================================================== +! + subroutine check_block_neighbors(pmeta) + +! import external procedures and variables +! + use error , only : print_warning + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_meta), pointer, intent(in) :: pmeta + +! local pointers +! + type(block_meta), pointer :: pneigh, pnneigh + +! local variables +! + integer :: i, j, k, l, m +! +!------------------------------------------------------------------------------- +! +! return if the block is not leaf +! + if (.not. pmeta%leaf) return + +! iterate over all face neighbors +! + do i = 1, ndims + do j = 1, nsides + m = 3 - j + do k = 1, nfaces + +! assign pointer with the neighbor +! + pneigh => pmeta%neigh(i,j,k)%ptr + +! check if the neighbor is associated +! + if (associated(pneigh)) then + +! check neighbors on the same levels +! + if (pmeta%level == pneigh%level) then + +! assign pointer to the neighbor of the neighbor pointing to the current meta +! block +! + pnneigh => pneigh%neigh(i,m,k)%ptr + +! check if it is associated +! + if (associated(pnneigh)) then + +! check if the pointer of the neighbor points to the current meta block +! + if (pmeta%id /= pnneigh%id) then + +! print warning, since the blocks differ +! + call print_warning("blocks::check_block_neighbors" & + , "Inconsistent same level neighbors!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k + + end if + + else ! pnneigh associated + +! print warning, since the pointer should be associated +! + call print_warning("blocks::check_block_neighbors" & + , "Same level neighbor not associated!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k + + end if ! pnneigh associated + + end if ! the same levels + +! check neighbors on the level higher than the meta block's level; it also +! covers the other way around, since we iterate over all neighbor faces +! + if (pmeta%level < pneigh%level) then + +! iterate over whole face of the corresponding neighbor +! + do l = 1, nfaces + +! assign pointer to the corresponding neighbor of the neighbor +! + pnneigh => pneigh%neigh(i,m,l)%ptr + +! check if it is associated +! + if (associated(pnneigh)) then + +! check if the pointer of the neighbor points to the current meta block +! + if (pmeta%id /= pnneigh%id) then + +! print warning, since the blocks differ +! + call print_warning("blocks::check_block_neighbors" & + , "Inconsistent different level neighbors!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k, l + + end if + + else ! pnneigh associated + +! print warning, since the pointer should be associated +! + call print_warning("blocks::check_block_neighbors" & + , "Different level neighbor not associated!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k, l + + end if ! pnneigh associated + + end do ! l = 1, nfaces + + end if ! pmeta's level < pneigh's level + + end if ! pneigh associated + + end do ! nfaces + end do ! nsides + end do ! ndims + +!------------------------------------------------------------------------------- +! + end subroutine check_block_neighbors +#endif /* DEBUG */ !=============================================================================== ! diff --git a/src/mesh.F90 b/src/mesh.F90 index d640bf9..baa8ac2 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -462,6 +462,9 @@ module mesh use blocks , only : link_blocks, unlink_blocks, refine_block use blocks , only : get_mblocks, get_nleafs use blocks , only : set_neighbors_refine +#ifdef DEBUG + use blocks , only : check_neighbors +#endif /* DEBUG */ use coordinates , only : minlev, maxlev use domains , only : setup_domain use error , only : print_error @@ -735,6 +738,12 @@ module mesh end do ! pmeta +#ifdef DEBUG +! check if neighbors are consistent after mesh generation +! + call check_neighbors() +#endif /* DEBUG */ + #ifdef PROFILE ! stop accounting time for the initial mesh generation ! @@ -768,6 +777,9 @@ module mesh use blocks , only : refine_block, derefine_block use blocks , only : append_datablock, remove_datablock, link_blocks use blocks , only : set_neighbors_refine +#ifdef DEBUG + use blocks , only : check_neighbors +#endif /* DEBUG */ use coordinates , only : minlev, maxlev, toplev, im, jm, km use equations , only : nv use error , only : print_error @@ -1257,6 +1269,12 @@ module mesh call redistribute_blocks() #endif /* MPI */ +#ifdef DEBUG +! check if neighbors are consistent after mesh refinement +! + call check_neighbors() +#endif /* DEBUG */ + #ifdef PROFILE ! stop accounting time for the adaptive mesh refinement update ! From c1d1d54b425210bc8cd8566435c584e9420d5700 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 28 Apr 2014 10:43:37 -0300 Subject: [PATCH 40/41] BLOCKS: Rewrite and micro optimize iterate_over_neighbors(). We prepare a set of indices for all faces in the block. Those indices point to edges and corners for each face, and are set one during the first execution. Then, we simply iterate over faces and call iterate_over_face() in order to apply subroutine pprocedure to edges and corners connected to a given face. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 755 ++++++++++++++++++++++++------------------------- 1 file changed, 366 insertions(+), 389 deletions(-) diff --git a/src/blocks.F90 b/src/blocks.F90 index c6533d0..b742b5b 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -54,11 +54,13 @@ module blocks ! nsides - the number of sides along each direction (2); ! nfaces - the number of faces at each side (2 for 2D, 4 for 3D); ! nchildren - the number of child blocks for each block (4 for 2D, 8 for 3D); +! mfaces - the number of faces in block (8 for 2D, 24 for 3D); ! integer(kind=4), parameter :: ndims = NDIMS integer(kind=4), parameter :: nsides = 2 integer(kind=4), parameter :: nfaces = 2**(ndims - 1) integer(kind=4), parameter :: nchildren = 2**ndims + integer(kind=4), parameter :: mfaces = nsides * nfaces * ndims ! MODULE VARIABLES: ! ================ @@ -3166,462 +3168,437 @@ module blocks type(block_meta) , pointer, intent(inout) :: pmeta procedure(reset_neighbors_update), pointer, intent(in) :: pprocedure -! local pointers +! local saved variables ! - type(block_meta), pointer :: pneigh + logical, save :: first = .true. + +! local saved arrays +! + integer, dimension(mfaces,3) , save :: fidx + integer, dimension(mfaces,2,3), save :: eidx + integer, dimension(mfaces,2,3), save :: cidx + +! local variables +! + integer :: l ! !------------------------------------------------------------------------------- ! -#if NDIMS == 3 -!! 3D CASE: WALK AROUND THE CORNERS -!! +! prepare indices +! + if (first) then + +! inicialize indices +! + fidx(: ,:) = 0 + eidx(:,:,:) = 0 + cidx(:,:,:) = 0 + +! prepare indices to get proper face, edge and corner neighbors +! +#if NDIMS == 2 ! around (0,0,0) corner ! -! (1,1,1) - X = (1,2,1) - Y = (2,1,2) +! [1,1,1]:X:(1,2,1) <- Y = [2,1,2] ! - pneigh => pmeta%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 1, :) = (/ 1, 1, 1 /) + eidx( 1,1,:) = (/ 2, 1, 2 /) - pneigh => pneigh%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,1) - Z = (3,2,1) - X = (1,1,3) +! [2,1,1]:Y:(2,2,1) <- X = [1,1,2] ! - pneigh => pmeta%neigh(3,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,1,1) - Y = (2,2,1) - Z = (3,1,3) - X = [1,1,4] -! - pneigh => pmeta%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 2, :) = (/ 2, 1, 1 /) + eidx( 2,1,:) = (/ 1, 1, 2 /) ! around (0,1,0) corner ! -! (2,2,1) + Y = (2,1,1) - X = (1,1,1) +! [1,1,2]:X:(1,2,2) -> Y = [2,2,2] ! - pneigh => pmeta%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 3, :) = (/ 1, 1, 2 /) + eidx( 3,1,:) = (/ 2, 2, 2 /) - pneigh => pneigh%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,3) - Z = (3,2,3) + Y = (2,2,3) +! [2,2,1]:Y:(2,1,2) <- X = [1,1,1] ! - pneigh => pmeta%neigh(3,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,1,2) - X = (1,2,2) - Z = (3,1,4) + Y = [2,2,4] -! - pneigh => pmeta%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 4, :) = (/ 2, 2, 1 /) + eidx( 4,1,:) = (/ 1, 1, 1 /) ! around (1,1,0) corner ! -! (1,2,2) + X = (1,1,2) + Y = (2,2,1) +! [1,2,2]:X:(1,1,2) -> Y = [2,2,1] ! - pneigh => pmeta%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 5, :) = (/ 1, 2, 2 /) + eidx( 5,1,:) = (/ 2, 2, 1 /) - pneigh => pneigh%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,4) - Z = (3,2,4) + X = (1,2,4) +! [2,2,2]:Y:(2,1,2) -> X = [1,2,1] ! - pneigh => pmeta%neigh(3,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,2,2) + Y = (2,1,2) - Z = (3,1,2) + X = [1,2,3] -! - pneigh => pmeta%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 6, :) = (/ 2, 2, 2 /) + eidx( 6,1,:) = (/ 1, 2, 1 /) ! around (1,0,0) corner ! -! (2,1,2) - Y = (2,2,2) + X = (1,2,2) +! [1,2,1]:X:(1,1,1) <- Y = [2,1,1] ! - pneigh => pmeta%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 7, :) = (/ 1, 2, 1 /) + eidx( 7,1,:) = (/ 2, 1, 1 /) - pneigh => pneigh%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,2) - Z = (3,2,2) - Y = (2,1,4) +! [2,1,2]:Y:(2,2,2) -> X = [1,2,2] ! - pneigh => pmeta%neigh(3,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,2,1) + X = (1,1,1) - Z = (3,1,1) - Y = [2,1,3] + fidx( 8, :) = (/ 2, 1, 2 /) + eidx( 8,1,:) = (/ 1, 2, 2 /) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 +! around (0,0,0) corner ! - pneigh => pmeta%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) +! [1,1,1]:X:(1,2,1) <- Y = [2,1,2] <- Z = [3,1,4] +! [1,1,1]:X:(1,2,1) <- Z = [3,1,2] <- Y = [2,1,4] +! + fidx( 1, :) = (/ 1, 1, 1 /) + eidx( 1,1,:) = (/ 2, 1, 2 /) + eidx( 1,2,:) = (/ 3, 1, 2 /) + cidx( 1,1,:) = (/ 3, 1, 4 /) + cidx( 1,2,:) = (/ 2, 1, 4 /) - pneigh => pneigh%neigh(3,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) +! [2,1,1]:Y:(2,2,1) <- Z = [3,1,3] <- X = [1,1,4] +! [2,1,1]:Y:(2,2,1) <- X = [1,1,2] <- Z = [3,1,4] +! + fidx( 2, :) = (/ 2, 1, 1 /) + eidx( 2,1,:) = (/ 3, 1, 3 /) + eidx( 2,2,:) = (/ 1, 1, 2 /) + cidx( 2,1,:) = (/ 1, 1, 4 /) + cidx( 2,2,:) = (/ 3, 1, 4 /) - pneigh => pneigh%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if +! [3,1,1]:Z:(3,2,1) <- X = [1,1,3] <- Y = [2,1,4] +! [3,1,1]:Z:(3,2,1) <- Y = [2,1,3] <- X = [1,1,4] +! + fidx( 3, :) = (/ 3, 1, 1 /) + eidx( 3,1,:) = (/ 1, 1, 3 /) + eidx( 3,2,:) = (/ 2, 1, 3 /) + cidx( 3,1,:) = (/ 2, 1, 4 /) + cidx( 3,2,:) = (/ 1, 1, 4 /) + +! around (0,1,0) corner +! +! [1,1,2]:X:(1,2,2) -> Y = [2,2,2] <- Z = [3,1,2] +! [1,1,2]:X:(1,2,2) <- Z = [3,1,4] -> Y = [2,2,4] +! + fidx( 4, :) = (/ 1, 1, 2 /) + eidx( 4,1,:) = (/ 2, 2, 2 /) + eidx( 4,2,:) = (/ 3, 1, 4 /) + cidx( 4,1,:) = (/ 3, 1, 2 /) + cidx( 4,2,:) = (/ 2, 2, 4 /) + +! [2,2,1]:Y:(2,1,2) <- Z = [3,1,1] <- X = [1,1,3] +! [2,2,1]:Y:(2,1,2) <- X = [1,1,1] <- Z = [3,1,2] +! + fidx( 5, :) = (/ 2, 2, 1 /) + eidx( 5,1,:) = (/ 3, 1, 1 /) + eidx( 5,2,:) = (/ 1, 1, 1 /) + cidx( 5,1,:) = (/ 1, 1, 3 /) + cidx( 5,2,:) = (/ 3, 1, 2 /) + +! [3,1,3]:Z:(3,2,3) <- X = [1,1,4] -> Y = [2,2,4] +! [3,1,3]:Z:(3,2,3) -> Y = [2,2,3] <- X = [1,1,3] +! + fidx( 6, :) = (/ 3, 1, 3 /) + eidx( 6,1,:) = (/ 1, 1, 4 /) + eidx( 6,2,:) = (/ 2, 2, 3 /) + cidx( 6,1,:) = (/ 2, 2, 4 /) + cidx( 6,2,:) = (/ 1, 1, 3 /) + +! around (1,1,0) corner +! +! [1,2,2]:X:(1,1,2) -> Y = [2,2,1] <- Z = [3,1,1] +! [1,2,2]:X:(1,1,2) <- Z = [3,1,3] -> Y = [2,2,3] +! + fidx( 7, :) = (/ 1, 2, 2 /) + eidx( 7,1,:) = (/ 2, 2, 1 /) + eidx( 7,2,:) = (/ 3, 1, 3 /) + cidx( 7,1,:) = (/ 3, 1, 1 /) + cidx( 7,2,:) = (/ 2, 2, 3 /) + +! [2,2,2]:Y:(2,1,2) <- Z = [3,1,2] -> X = [1,2,3] +! [2,2,2]:Y:(2,1,2) -> X = [1,2,1] <- Z = [3,1,1] +! + fidx( 8, :) = (/ 2, 2, 2 /) + eidx( 8,1,:) = (/ 3, 1, 2 /) + eidx( 8,2,:) = (/ 1, 2, 1 /) + cidx( 8,1,:) = (/ 1, 2, 3 /) + cidx( 8,2,:) = (/ 3, 1, 1 /) + +! [3,1,4]:Z:(3,2,4) -> X = [1,2,4] -> Y = [2,2,3] +! [3,1,4]:Z:(3,2,4) -> Y = [2,2,4] -> X = [1,2,3] +! + fidx( 9, :) = (/ 3, 1, 4 /) + eidx( 9,1,:) = (/ 1, 2, 4 /) + eidx( 9,2,:) = (/ 2, 2, 4 /) + cidx( 9,1,:) = (/ 2, 2, 3 /) + cidx( 9,2,:) = (/ 1, 2, 3 /) + +! around (1,0,0) corner +! +! [1,2,1]:X:(1,1,1) <- Y = [2,1,1] <- Z = [3,1,3] +! [1,2,1]:X:(1,1,1) <- Z = [3,1,1] <- Y = [2,1,3] +! + fidx(10, :) = (/ 1, 2, 1 /) + eidx(10,1,:) = (/ 2, 1, 1 /) + eidx(10,2,:) = (/ 3, 1, 1 /) + cidx(10,1,:) = (/ 3, 1, 3 /) + cidx(10,2,:) = (/ 2, 1, 3 /) + +! [2,1,2]:Y:(2,2,2) <- Z = [3,1,4] -> X = [1,2,4] +! [2,1,2]:Y:(2,2,2) -> X = [1,2,2] <- Z = [3,1,3] +! + fidx(11, :) = (/ 2, 1, 2 /) + eidx(11,1,:) = (/ 3, 1, 4 /) + eidx(11,2,:) = (/ 1, 2, 2 /) + cidx(11,1,:) = (/ 1, 2, 4 /) + cidx(11,2,:) = (/ 3, 1, 3 /) + +! [3,1,2]:Z:(3,2,2) -> X = [1,2,3] <- Y = [2,1,3] +! [3,1,2]:Z:(3,2,2) <- Y = [2,1,4] -> X = [1,2,4] +! + fidx(12, :) = (/ 3, 1, 2 /) + eidx(12,1,:) = (/ 1, 2, 3 /) + eidx(12,2,:) = (/ 2, 1, 4 /) + cidx(12,1,:) = (/ 2, 1, 3 /) + cidx(12,2,:) = (/ 1, 2, 4 /) ! around (0,0,1) corner ! -! (1,1,3) - X = (1,2,3) + Z = (3,2,2) +! [1,1,3]:X:(1,2,3) <- Y = [2,1,4] -> Z = [3,2,4] +! [1,1,3]:X:(1,2,3) -> Z = [3,2,2] <- Y = [2,1,2] ! - pneigh => pmeta%neigh(1,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(13, :) = (/ 1, 1, 3 /) + eidx(13,1,:) = (/ 2, 1, 4 /) + eidx(13,2,:) = (/ 3, 2, 2 /) + cidx(13,1,:) = (/ 3, 2, 4 /) + cidx(13,2,:) = (/ 2, 1, 2 /) - pneigh => pneigh%neigh(3,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,1,3) - Y = (2,2,3) - X = (1,1,4) +! [2,1,3]:Y:(2,2,3) -> Z = [3,2,3] <- X = [1,1,2] +! [2,1,3]:Y:(2,2,3) <- X = [1,1,4] -> Z = [3,2,1] ! - pneigh => pmeta%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(14, :) = (/ 2, 1, 3 /) + eidx(14,1,:) = (/ 3, 2, 3 /) + eidx(14,2,:) = (/ 1, 1, 4 /) + cidx(14,1,:) = (/ 1, 1, 2 /) + cidx(14,2,:) = (/ 3, 2, 1 /) - pneigh => pneigh%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,1) + Z = (3,1,1) - Y = (2,1,1) - X = [1,1,2] +! [3,2,1]:Z:(3,1,1) <- X = [1,1,1] <- Y = [2,1,2] +! [3,2,1]:Z:(3,1,1) <- Y = [2,1,1] <- X = [1,1,2] ! - pneigh => pmeta%neigh(3,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(15, :) = (/ 3, 2, 1 /) + eidx(15,1,:) = (/ 1, 1, 1 /) + eidx(15,2,:) = (/ 2, 1, 1 /) + cidx(15,1,:) = (/ 2, 1, 2 /) + cidx(15,2,:) = (/ 1, 1, 2 /) ! around (0,1,1) corner ! -! (2,2,3) + Y = (2,1,3) + Z = (3,2,1) +! [1,1,4]:X:(1,2,4) -> Y = [2,2,4] -> Z = [3,2,2] +! [1,1,4]:X:(1,2,4) -> Z = [3,2,4] -> Y = [2,2,2] ! - pneigh => pmeta%neigh(2,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(16, :) = (/ 1, 1, 4 /) + eidx(16,1,:) = (/ 2, 2, 4 /) + eidx(16,2,:) = (/ 3, 2, 4 /) + cidx(16,1,:) = (/ 3, 2, 2 /) + cidx(16,2,:) = (/ 2, 2, 2 /) - pneigh => pneigh%neigh(3,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,1,4) - X = (1,2,4) + Y = (2,2,4) +! [2,2,3]:Y:(2,1,3) -> Z = [3,2,1] <- X = [1,1,1] +! [2,2,3]:Y:(2,1,3) <- X = [1,1,3] -> Z = [3,2,2] ! - pneigh => pmeta%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(17, :) = (/ 2, 2, 3 /) + eidx(17,1,:) = (/ 3, 2, 1 /) + eidx(17,2,:) = (/ 1, 1, 3 /) + cidx(17,1,:) = (/ 1, 1, 1 /) + cidx(17,2,:) = (/ 3, 2, 2 /) - pneigh => pneigh%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,3) + Z = (3,1,3) - X = (1,1,2) + Z = [2,2,2] +! [3,2,3]:Z:(3,1,3) <- X = [1,1,2] -> Y = [2,2,2] +! [3,2,3]:Z:(3,1,3) -> Y = [2,2,1] <- X = [1,1,1] ! - pneigh => pmeta%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(18, :) = (/ 3, 2, 3 /) + eidx(18,1,:) = (/ 1, 1, 2 /) + eidx(18,2,:) = (/ 2, 2, 1 /) + cidx(18,1,:) = (/ 2, 2, 2 /) + cidx(18,2,:) = (/ 1, 1, 1 /) ! around (1,1,1) corner ! -! (1,2,4) + X = (1,1,4) + Z = (3,2,3) +! [1,2,4]:X:(1,1,4) -> Y = [2,2,3] -> Z = [3,2,1] +! [1,2,4]:X:(1,1,4) -> Z = [3,2,3] -> Y = [2,2,1] ! - pneigh => pmeta%neigh(1,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(19, :) = (/ 1, 2, 4 /) + eidx(19,1,:) = (/ 2, 2, 3 /) + eidx(19,2,:) = (/ 3, 2, 3 /) + cidx(19,1,:) = (/ 3, 2, 1 /) + cidx(19,2,:) = (/ 2, 2, 1 /) - pneigh => pneigh%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,2,4) + Y = (2,1,4) + X = (1,2,3) +! [2,2,4]:Y:(2,1,4) -> Z = [3,2,2] -> X = [1,2,1] +! [2,2,4]:Y:(2,1,4) -> X = [1,2,3] -> Z = [3,2,1] ! - pneigh => pmeta%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(20, :) = (/ 2, 2, 4 /) + eidx(20,1,:) = (/ 3, 2, 2 /) + eidx(20,2,:) = (/ 1, 2, 3 /) + cidx(20,1,:) = (/ 1, 2, 1 /) + cidx(20,2,:) = (/ 3, 2, 1 /) - pneigh => pneigh%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,4) + Z = (3,1,4) + Y = (2,2,2) + X = [1,2,1] +! [3,2,4]:Z:(3,1,4) -> X = [1,2,2] -> Y = [2,2,1] +! [3,2,4]:Z:(3,1,4) -> Y = [2,2,2] -> X = [1,2,1] ! - pneigh => pmeta%neigh(3,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(21, :) = (/ 3, 2, 4 /) + eidx(21,1,:) = (/ 1, 2, 2 /) + eidx(21,2,:) = (/ 2, 2, 2 /) + cidx(21,1,:) = (/ 2, 2, 1 /) + cidx(21,2,:) = (/ 1, 2, 1 /) ! around (1,0,1) corner ! -! (2,1,4) - Y = (2,2,4) + Z = (3,2,4) +! [1,2,3]:X:(1,1,3) <- Y = [2,1,3] -> Z = [3,2,3] +! [1,2,3]:X:(1,1,3) -> Z = [3,2,1] <- Y = [2,1,1] ! - pneigh => pmeta%neigh(2,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(22, :) = (/ 1, 2, 3 /) + eidx(22,1,:) = (/ 2, 1, 3 /) + eidx(22,2,:) = (/ 3, 2, 1 /) + cidx(22,1,:) = (/ 3, 2, 3 /) + cidx(22,2,:) = (/ 2, 1, 1 /) - pneigh => pneigh%neigh(3,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,2) + Z = (3,1,2) + X = (1,2,1) +! [2,1,4]:Y:(2,2,4) -> Z = [3,2,4] -> X = [1,2,2] +! [2,1,4]:Y:(2,2,4) -> X = [1,2,4] -> Z = [3,2,3] ! - pneigh => pmeta%neigh(3,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(23, :) = (/ 2, 1, 4 /) + eidx(23,1,:) = (/ 3, 2, 4 /) + eidx(23,2,:) = (/ 1, 2, 4 /) + cidx(23,1,:) = (/ 1, 2, 2 /) + cidx(23,2,:) = (/ 3, 2, 3 /) - pneigh => pneigh%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,2,3) + X = (1,1,3) - Y = (2,1,3) + Z = [3,2,3] +! [3,2,2]:Z:(3,1,2) -> X = [1,2,1] <- Y = [2,1,3] +! [3,2,2]:Z:(3,1,2) <- Y = [2,1,2] -> X = [1,2,4] ! - pneigh => pmeta%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if - -#else /* NDIMS == 3 */ -!! 2D CASE -!! -! (0,0)-(0,½) edge and (0,0) corner -! - pneigh => pmeta%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (0,½)-(0,1) edge and (0,1) corner -! - pneigh => pmeta%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,0)-(1,½) edge and (1,0) corner -! - pneigh => pmeta%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,½)-(1,1) edge and (1,1) corner -! - pneigh => pmeta%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (0,0)-(½,0) edge and (0,0) corner -! - pneigh => pmeta%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (½,0)-(1,0) edge and (1,0) corner -! - pneigh => pmeta%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (0,1)-(½,1) edge and (0,1) corner -! - pneigh => pmeta%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (½,1)-(1,1) edge and (1,1) corner -! - pneigh => pmeta%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if + fidx(24, :) = (/ 3, 2, 2 /) + eidx(24,1,:) = (/ 1, 2, 1 /) + eidx(24,2,:) = (/ 2, 1, 2 /) + cidx(24,1,:) = (/ 2, 1, 3 /) + cidx(24,2,:) = (/ 1, 2, 3 /) #endif /* NDIMS == 3 */ +! reset the first time execution flag +! + first = .false. + + end if + +! iterate over all block faces (or edges in the 2D case) +! + do l = 1, mfaces + call iterate_over_face(pmeta, pprocedure & + , fidx(l,:), eidx(l,:,:), cidx(l,:,:)) + end do + !------------------------------------------------------------------------------- ! end subroutine iterate_over_neighbors +! +!=============================================================================== +! +! subroutine ITERATE_OVER_FACE: +! ---------------------------- +! +! Subroutine iterates over all neighbors, edges and corners linked to +! the input meta block and executes a subroutine provided by the pointer. +! +! Arguments: +! +! pmeta - a pointer to the meta block which neighbors are iterated over; +! pproc - a pointer to the subroutine called with each pair (pmeta, pneigh); +! fidx - the index of face to process; +! eidx - the indices of faces connected with edges; +! cidx - the indices of faces connected with corners; +! +!=============================================================================== +! + subroutine iterate_over_face(pmeta, pprocedure, fidx, eidx, cidx) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_meta) , pointer, intent(inout) :: pmeta + procedure(reset_neighbors_update), pointer, intent(in) :: pprocedure + integer, dimension(3) , intent(in) :: fidx + integer, dimension(2,3) , intent(in) :: eidx + integer, dimension(2,3) , intent(in) :: cidx + +! local pointers +! + type(block_meta), pointer :: pneigh, pedge, pcorner +! +!------------------------------------------------------------------------------- +! +! associate a pointer with the neighbor +! + pneigh => pmeta%neigh(fidx(1),fidx(2),fidx(3))%ptr + +! check if the neighbors is associated +! + if (associated(pneigh)) then + +! call the procedure for the face neighbor +! + call pprocedure(pmeta, pneigh) + +! associate a pointer with the first edge +! + pedge => pneigh%neigh(eidx(1,1),eidx(1,2),eidx(1,3))%ptr + +! check if the edge pointer is associated +! + if (associated(pedge)) then + +! call the procedure for the edge neighbor +! + call pprocedure(pmeta, pedge) + +#if NDIMS == 3 +! associate a pointer with the first corner +! + pcorner => pedge%neigh(cidx(1,1),cidx(1,2),cidx(1,3))%ptr + +! call the procedure for the corner neighbor if it is associated +! + if (associated(pcorner)) call pprocedure(pmeta, pcorner) +#endif /* NDIMS == 3 */ + + end if ! pedge associated + +#if NDIMS == 3 +! associate a pointer with the second edge +! + pedge => pneigh%neigh(eidx(2,1),eidx(2,2),eidx(2,3))%ptr + +! check if the edge pointer is associated +! + if (associated(pedge)) then + +! call the procedure for the edge neighbor +! + call pprocedure(pmeta, pedge) + +! associate a pointer with the second corner +! + pcorner => pedge%neigh(cidx(2,1),cidx(2,2),cidx(2,3))%ptr + +! call the procedure for the corner neighbor if it is associated +! + if (associated(pcorner)) call pprocedure(pmeta, pcorner) + + end if ! pedge associated +#endif /* NDIMS == 3 */ + + end if ! pneigh associated + +!------------------------------------------------------------------------------- +! + end subroutine iterate_over_face #ifdef DEBUG ! !=============================================================================== From 471e8a690a98ca021ce914bd9ce206790b1e102b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 28 Apr 2014 10:50:41 -0300 Subject: [PATCH 41/41] BLOCKS: Add timer for neighbor consistency check. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/blocks.F90 b/src/blocks.F90 index b742b5b..7222a0b 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -45,6 +45,9 @@ module blocks ! timer indices ! integer, save :: imi, ima, imu, imp, imq, imr, imd +#ifdef DEBUG + integer, save :: imc +#endif /* DEBUG */ #endif /* PROFILE */ ! MODULE PARAMETERS: @@ -330,6 +333,9 @@ module blocks call set_timer('blocks:: data block deallocation', imq) call set_timer('blocks:: refine' , imr) call set_timer('blocks:: derefine' , imd) +#ifdef DEBUG + call set_timer('blocks:: check neighbors' , imc) +#endif /* DEBUG */ ! start accounting time for module initialization/finalization ! @@ -2740,6 +2746,12 @@ module blocks ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the neighbor consistency check +! + call start_timer(imc) +#endif /* PROFILE */ + ! associate the pointer with the first block on the meta block list ! pmeta => list_meta @@ -2758,6 +2770,12 @@ module blocks end do ! meta blocks +#ifdef PROFILE +! stop accounting time for the neighbor consistency check +! + call stop_timer(imc) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine check_neighbors