From a42e6a0a166e9ceb60a89db1763310c695623c9b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 5 Jun 2015 18:41:39 -0300 Subject: [PATCH 1/4] EQUATIONS: Fill the block borders in update_primitive_variables(). This should not be normally required, but it helps to avoid the unexpected appearance of NaNs and the code execution termination due to the floating point exceptions. This should also help with more advanced open boundary conditions (such as the conditions which keep the divergence or curl of vector fields zero). Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 0b59487..9d4d42e 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -1017,8 +1017,9 @@ module equations ! include external procedures and variables ! - use coordinates, only : im, jm, km, in, jn, kn - use coordinates, only : ib, jb, kb, ie, je, ke + 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 ! local variables are not implicit by default ! @@ -1079,6 +1080,29 @@ module equations end if +! fill out the borders +! + do i = ibl, 1, -1 + qq(1:nv, i,jb:je,kb:ke) = qq(1:nv,ib,jb:je,kb:ke) + end do + do i = ieu, im + qq(1:nv, i,jb:je,kb:ke) = qq(1:nv,ie,jb:je,kb:ke) + end do + do j = jbl, 1, -1 + qq(1:nv, 1:im, j,kb:ke) = qq(1:nv, 1:im,jb,kb:ke) + end do + do j = jeu, jm + qq(1:nv, 1:im, j,kb:ke) = qq(1:nv, 1:im,je,kb:ke) + end do +#if NDIMS == 3 + do k = kbl, 1, -1 + qq(1:nv, 1:im, 1:jm, k) = qq(1:nv, 1:im, 1:jm,kb) + end do + do k = keu, km + qq(1:nv, 1:im, 1:jm, k) = qq(1:nv, 1:im, 1:jm,ke) + end do +#endif /* NDIMS == 3 */ + !------------------------------------------------------------------------------- ! end subroutine update_primitive_variables From ff9483e95c575d124778fb921b6170d2d87147d3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 5 Jun 2015 18:59:31 -0300 Subject: [PATCH 2/4] MESH: Fill the block borders in restrict_block(). Since the borders are undefined, it is better to fill them with the copy of the last available boundary layers. Even though they will be later updated by the boundary update, it eliminates the possibilty of propagating NaNs inside the blocks. Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/mesh.F90 b/src/mesh.F90 index ed58cdf..e406a31 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1373,6 +1373,43 @@ module mesh end do ! p = 1, nchildren +! fill out the borders +! +#if NDIMS == 2 + do il = is - 1, 1, -1 + pparent%u(1:nv,il,js:jt, 1 ) = pparent%u(1:nv,is,js:jt, 1 ) + end do + do iu = it + 1, im + pparent%u(1:nv,iu,js:jt, 1 ) = pparent%u(1:nv,it,js:jt, 1 ) + end do + do jl = js - 1, 1, -1 + pparent%u(1:nv, 1:im,jl, 1 ) = pparent%u(1:nv, 1:im,js, 1 ) + end do + do ju = jt + 1, jm + pparent%u(1:nv, 1:im,ju, 1 ) = pparent%u(1:nv, 1:im,jt, 1 ) + end do +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + do il = is - 1, 1, -1 + pparent%u(1:nv,il,js:jt,ks:kt) = pparent%u(1:nv,is,js:jt,ks:kt) + end do + do iu = it + 1, im + pparent%u(1:nv,iu,js:jt,ks:kt) = pparent%u(1:nv,it,js:jt,ks:kt) + end do + do jl = js - 1, 1, -1 + pparent%u(1:nv, 1:im,jl,ks:kt) = pparent%u(1:nv, 1:im,js,ks:kt) + end do + do ju = jt + 1, jm + pparent%u(1:nv, 1:im,ju,ks:kt) = pparent%u(1:nv, 1:im,jt,ks:kt) + end do + do kl = ks - 1, 1, -1 + pparent%u(1:nv, 1:im, 1:jm,kl) = pparent%u(1:nv, 1:im, 1:jm,ks) + end do + do ku = kt + 1, km + pparent%u(1:nv, 1:im, 1:jm,ku) = pparent%u(1:nv, 1:im, 1:jm,kt) + end do +#endif /* NDIMS == 3 */ + #ifdef PROFILE ! stop accounting time for the block restriction ! From 95deee811397e1d3dc67a61f21b1e6ea6aaa67a3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 5 Jun 2015 19:07:30 -0300 Subject: [PATCH 3/4] SCHEMES: Initialize fluxes in update_flux() subroutines. Make the flux f(:,:,:,:) output argument only. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 72 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 48 insertions(+), 24 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index d6ed6ad..29ceb2c 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -592,10 +592,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -615,6 +615,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) @@ -1094,10 +1098,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -1117,6 +1121,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) @@ -1778,10 +1786,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -1801,6 +1809,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) @@ -3090,10 +3102,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -3113,6 +3125,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) @@ -4538,10 +4554,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -4561,6 +4577,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) @@ -5233,10 +5253,10 @@ module schemes ! input arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(out) :: f ! local variables ! @@ -5256,6 +5276,10 @@ module schemes call start_timer(imf) #endif /* PROFILE */ +! initialize fluxes +! + f(1:nv,1:im,1:jm,1:km) = 0.0d+00 + ! select the directional flux to compute ! select case(idir) From 4a5cec8168d4ac96b52ec50465d5728e4c073721 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 5 Jun 2015 19:09:46 -0300 Subject: [PATCH 4/4] EVOLUTION: No need to initialize fluxes in update_fluxes(). The fluxes are already initialized in update_flux(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index 8349adc..ff8146b 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -1848,10 +1848,6 @@ module evolution dx(2) = ady(pdata%meta%level) dx(3) = adz(pdata%meta%level) -! initialize fluxes -! - pdata%f(1:NDIMS,1:nv,1:im,1:jm,1:km) = 0.0d+00 - ! update fluxes for the current block ! do n = 1, NDIMS