diff options
Diffstat (limited to 'src/include/upwind_characteristic_second2.h')
-rw-r--r-- | src/include/upwind_characteristic_second2.h | 146 |
1 files changed, 146 insertions, 0 deletions
diff --git a/src/include/upwind_characteristic_second2.h b/src/include/upwind_characteristic_second2.h new file mode 100644 index 0000000..e851859 --- /dev/null +++ b/src/include/upwind_characteristic_second2.h @@ -0,0 +1,146 @@ +! Calculation of characteristic upwinded differences except at boundaries. +! $Header$ + +! If this is an active cell we have to compute its derivative. +if ( eh_mask(i,j,k) .ge. 0 ) then + + ! If it is not a boundary cell in the x-direction then calculate as a + ! starting point the centered derivative in the x-direction. + if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. & + ( .not. btest(eh_mask(i,j,k),1) ) ) then + dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) ) + +! Otherwise we will have to calculate the appropriate one sided derivative. + else if ( btest(eh_mask(i,j,k),0) ) then + dfx(i,j,k) = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) ) + else + dfx(i,j,k) = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) ) + end if + + ! Do the same for the y direction. + if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. & + ( .not. btest(eh_mask(i,j,k),3) ) ) then + dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) ) + else if ( btest(eh_mask(i,j,k),2) ) then + dfy(i,j,k) = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) ) + else + dfy(i,j,k) = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) ) + end if + + ! And for the z direction. + if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. & + ( .not. btest(eh_mask(i,j,k),5) ) ) then + dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) ) + else if ( btest(eh_mask(i,j,k),4) ) then + dfz(i,j,k) = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) ) + else + dfz(i,j,k) = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) ) + end if + +! Calculate the characteristic direction using these preliminary +! derivatives and multiply it with the timestep. + + dfup(1) = g3xx * dfx(i,j,k) + g3xy * dfy(i,j,k) + g3xz * dfz(i,j,k) + dfup(2) = g3xy * dfx(i,j,k) + g3yy * dfy(i,j,k) + g3yz * dfz(i,j,k) + dfup(3) = g3xz * dfx(i,j,k) + g3yz * dfy(i,j,k) + g3zz * dfz(i,j,k) + + cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k) + & + dfup(2) * dfy(i,j,k) + & + dfup(3) * dfz(i,j,k) ) ) + + cdx(1) = ( -betax(i,j,k) + alp2 * dfup(1) * cfactor ) * cctk_delta_time + cdx(2) = ( -betay(i,j,k) + alp2 * dfup(2) * cfactor ) * cctk_delta_time + cdx(3) = ( -betaz(i,j,k) + alp2 * dfup(3) * cfactor ) * cctk_delta_time + + if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. & + ( .not. btest(eh_mask(i,j,k),1) ) ) then + ! Look at the characteristic direction to figure out in which direction + ! the stencil should be chosen. + if ( cdx(1) .gt. zero ) then ! Choose the stencil to the left. + + ! If the neighbour to the left is also an active non-boundary cell its + ! safe to assume that i-2, i-1,i and i+1 have good values. + ! Otherwise the neigbour must be a boundary cell, so we can only use a + ! centered derivative, which we have already computed. + if ( eh_mask(i-1,j,k) .eq. 0 ) then + + ! Choose the stencil, that gives the smalles second derivative. + ! If it happens to be the centered, do not do anything since we have + ! already calculated the centered derivative. + if ( f(i-2,j,k) - two * f(i-1,j,k) + f(i,j,k) .le. & + f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) then + dfx(i,j,k) = idx * ( three * f(i,j,k) - & + four * f(i-1,j,k) + f(i-2,j,k) ) + end if + end if + + else if ( cdx(1) .lt. zero ) then ! Choose the stencil to the right. + + ! If the neighbour to the right is also an active non-boundary cell its + ! safe to assume that i-1, i,i+1 and i+2 have good values. + ! Otherwise the neigbour must be a boundary cell, so we can only use a + ! centered derivative, which we have already computed. + if ( eh_mask(i+1,j,k) .eq. 0 ) then + + ! Choose the stencil, that gives the smalles second derivative. + ! If it happens to be the centered, do not do anything since we have + ! already calculated the centered derivative. + if ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) .ge. & + f(i,j,k) - two * f(i+1,j,k) + f(i+2,j,k) ) then + dfx(i,j,k) = idx * ( -three * f(i,j,k) + & + four * f(i+1,j,k) - f(i+2,j,k) ) + end if + end if + end if + end if + + ! Do the same for the y direction. + if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. & + ( .not. btest(eh_mask(i,j,k),3) ) ) then + + if ( cdx(2) .gt. zero ) then ! Choose the stencil to the left. + if ( eh_mask(i,j-1,k) .eq. 0 ) then + if ( f(i,j-2,k) - two * f(i,j-1,k) + f(i,j,k) .le. & + f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) then + dfy(i,j,k) = idy * ( three * f(i,j,k) - & + four * f(i,j-1,k) + f(i,j-2,k) ) + end if + end if + else if ( cdx(2) .lt. zero ) then ! Choose the stencil to the left. + if ( eh_mask(i,j+1,k) .eq. 0 ) then + if ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) .ge. & + f(i,j,k) - two * f(i,j+1,k) + f(i,j+2,k) ) then + dfy(i,j,k) = idy * ( -three * f(i,j,k) + & + four * f(i,j+1,k) - f(i,j+2,k) ) + end if + end if + end if + end if + + ! And for the z direction. + if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. & + ( .not. btest(eh_mask(i,j,k),5) ) ) then + + if ( cdx(3) .gt. zero ) then + if ( eh_mask(i,j,k-1) .eq. 0 ) then + if ( f(i,j,k-2) - two * f(i,j,k-1) + f(i,j,k) .le. & + f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) then + dfz(i,j,k) = idz * ( three * f(i,j,k) - & + four * f(i,j,k-1) + f(i,j,k-2) ) + end if + end if + else if ( cdx(3) .lt. zero ) then + if ( eh_mask(i,j,k+1) .eq. 0 ) then + if ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) .ge. & + f(i,j,k) - two * f(i,j,k+1) + f(i,j,k+2) ) then + dfz(i,j,k) = idz * ( -three * f(i,j,k) + & + four * f(i,j,k+1) - f(i,j,k+2) ) + end if + end if + end if + end if +else + dfx(i,j,k) = zero + dfy(i,j,k) = zero + dfz(i,j,k) = zero +end if |