From 2f5271e3e16ae2cdddd9970ebeaad1307f9698e2 Mon Sep 17 00:00:00 2001 From: diener Date: Fri, 16 May 2003 13:19:04 +0000 Subject: Include files for a new experimental upwinding scheme based on the characteristics of the evolution equation. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@105 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/include/metric.h | 19 ++++ src/include/upwind_characteristic_second2.h | 146 ++++++++++++++++++++++++++++ 2 files changed, 165 insertions(+) create mode 100644 src/include/metric.h create mode 100644 src/include/upwind_characteristic_second2.h diff --git a/src/include/metric.h b/src/include/metric.h new file mode 100644 index 0000000..f3a22ff --- /dev/null +++ b/src/include/metric.h @@ -0,0 +1,19 @@ +! Calculation of inverse three metric, lapse squared and shift +! $Header$ + +if ( eh_mask(i,j,k) .ge. 0 ) then + alp2 = alp(i,j,k)**2 + + tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2 + tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k) + tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k) + + idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + + g3xx = tmp1 * idetg + g3xy = tmp2 * idetg + g3xz = tmp3 * idetg + g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg + g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg +end if 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 -- cgit v1.2.3