From f8fe8c48366da18ebee8611369bee8c40478bc3a Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 30 Sep 2003 09:08:40 +0000 Subject: Preparing for release. Added comments. Changed the way the inverse three metric is calculated, so I don't have to calculate it both for the characteristic upwinding and for the rhs evaluation and to make it easier to convert between conformal metric and physical metric. This requires a set of grid functions for the inverse 3 metric. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@144 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/include/metric.h | 17 ++- src/include/upwind_characteristic_second2.h | 14 +- src/include/upwind_first.h | 190 +++++++++++++--------------- 3 files changed, 108 insertions(+), 113 deletions(-) (limited to 'src') diff --git a/src/include/metric.h b/src/include/metric.h index d4890df..97e7a18 100644 --- a/src/include/metric.h +++ b/src/include/metric.h @@ -1,19 +1,18 @@ -! Calculation of inverse three metric, lapse squared and shift +! Calculation of inverse three metric. At this point there is no distinction +! between the case of a physical or static conformal metric. ! $Header$ if ( eh_mask(i,j,k,l) .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 + g3xx(i,j,k) = tmp1 * idetg + g3xy(i,j,k) = tmp2 * idetg + g3xz(i,j,k) = tmp3 * idetg + g3yy(i,j,k) = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + g3yz(i,j,k) = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg + g3zz(i,j,k) = ( 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 index 04b1d06..deec132 100644 --- a/src/include/upwind_characteristic_second2.h +++ b/src/include/upwind_characteristic_second2.h @@ -40,9 +40,17 @@ if ( eh_mask(i,j,k,l) .ge. 0 ) then ! Calculate the characteristic direction using these preliminary ! derivatives and multiply it with the timestep. - dfup(1) = g3xx * dfx(i,j,k,l) + g3xy * dfy(i,j,k,l) + g3xz * dfz(i,j,k,l) - dfup(2) = g3xy * dfx(i,j,k,l) + g3yy * dfy(i,j,k,l) + g3yz * dfz(i,j,k,l) - dfup(3) = g3xz * dfx(i,j,k,l) + g3yz * dfy(i,j,k,l) + g3zz * dfz(i,j,k,l) + dfup(1) = g3xx(i,j,k) * dfx(i,j,k,l) + & + g3xy(i,j,k) * dfy(i,j,k,l) + & + g3xz(i,j,k) * dfz(i,j,k,l) + dfup(2) = g3xy(i,j,k) * dfx(i,j,k,l) + & + g3yy(i,j,k) * dfy(i,j,k,l) + & + g3yz(i,j,k) * dfz(i,j,k,l) + dfup(3) = g3xz(i,j,k) * dfx(i,j,k,l) + & + g3yz(i,j,k) * dfy(i,j,k,l) + & + g3zz(i,j,k) * dfz(i,j,k,l) + + alp2 = alp(i,j,k)**2 cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k,l) + & dfup(2) * dfy(i,j,k,l) + & diff --git a/src/include/upwind_first.h b/src/include/upwind_first.h index 1f8c309..366da99 100644 --- a/src/include/upwind_first.h +++ b/src/include/upwind_first.h @@ -1,105 +1,93 @@ ! Calculate first order upwinded differences. ! $Header$ -idx = one / cctk_delta_space(1) -idy = one / cctk_delta_space(2) -idz = one / cctk_delta_space(3) +if ( eh_mask(i,j,k,l) .ge. 0 ) then + if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),1) ) ) then + al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) ) + ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) ) + else if ( btest(eh_mask(i,j,k,l),0) ) then + al = zero + ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) ) + else + al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) ) + ar = zero + end if + if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),3) ) ) then + bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) ) + br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) ) + else if ( btest(eh_mask(i,j,k,l),2) ) then + bl = zero + br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) ) + else + bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) ) + br = zero + end if + if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),5) ) ) then + cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) ) + cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) ) + else if ( btest(eh_mask(i,j,k,l),4) ) then + cl = zero + cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) ) + else + cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) ) + cr = zero + end if -do l = 1, eh_number_level_sets - do k = 1, nz - do j = 1, ny - do i = 1, nx - if ( eh_mask(i,j,k,l) .ge. 0 ) then - if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. & - ( .not. btest(eh_mask(i,j,k,l),1) ) ) then - al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) ) - ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) ) - else if ( btest(eh_mask(i,j,k,l),0) ) then - al = zero - ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) ) - else - al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) ) - ar = zero - end if - if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. & - ( .not. btest(eh_mask(i,j,k,l),3) ) ) then - bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) ) - br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) ) - else if ( btest(eh_mask(i,j,k,l),2) ) then - bl = zero - br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) ) - else - bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) ) - br = zero - end if - if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. & - ( .not. btest(eh_mask(i,j,k,l),5) ) ) then - cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) ) - cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) ) - else if ( btest(eh_mask(i,j,k,l),4) ) then - cl = zero - cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) ) - else - cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) ) - cr = zero - end if - - alminus = half*(al-abs(al)) - alplus = half*(al+abs(al)) - arminus = half*(ar-abs(ar)) - arplus = half*(ar+abs(ar)) - - blminus = half*(bl-abs(bl)) - blplus = half*(bl+abs(bl)) - brminus = half*(br-abs(br)) - brplus = half*(br+abs(br)) - - clminus = half*(cl-abs(cl)) - clplus = half*(cl+abs(cl)) - crminus = half*(cr-abs(cr)) - crplus = half*(cr+abs(cr)) - - if ( f(i,j,k,l) .gt. 0 ) then - if ( abs(alplus) .gt. abs(arminus) ) then - dfx(i,j,k,l) = alplus - else - dfx(i,j,k,l) = arminus - end if - if ( abs(blplus) .gt. abs(brminus) ) then - dfy(i,j,k,l) = blplus - else - dfy(i,j,k,l) = brminus - end if - if ( abs(clplus) .gt. abs(crminus) ) then - dfz(i,j,k,l) = clplus - else - dfz(i,j,k,l) = crminus - end if - endif - - if ( f(i,j,k,l) .lt. 0 ) then - if ( abs(alminus) .gt. abs(arplus) ) then - dfx(i,j,k,l) = alminus - else - dfx(i,j,k,l) = arplus - end if - if ( abs(blminus) .gt. abs(brplus) ) then - dfy(i,j,k,l) = blminus - else - dfy(i,j,k,l) = brplus - end if - if ( abs(clminus) .gt. abs(crplus) ) then - dfz(i,j,k,l) = clminus - else - dfz(i,j,k,l) = crplus - end if - end if - else - dfx(i,j,k,l) = zero - dfy(i,j,k,l) = zero - dfz(i,j,k,l) = zero - end if - end do - end do - end do -end do + alminus = half*(al-abs(al)) + alplus = half*(al+abs(al)) + arminus = half*(ar-abs(ar)) + arplus = half*(ar+abs(ar)) + + blminus = half*(bl-abs(bl)) + blplus = half*(bl+abs(bl)) + brminus = half*(br-abs(br)) + brplus = half*(br+abs(br)) + + clminus = half*(cl-abs(cl)) + clplus = half*(cl+abs(cl)) + crminus = half*(cr-abs(cr)) + crplus = half*(cr+abs(cr)) + + if ( f(i,j,k,l) .gt. 0 ) then + if ( abs(alplus) .gt. abs(arminus) ) then + dfx(i,j,k,l) = alplus + else + dfx(i,j,k,l) = arminus + end if + if ( abs(blplus) .gt. abs(brminus) ) then + dfy(i,j,k,l) = blplus + else + dfy(i,j,k,l) = brminus + end if + if ( abs(clplus) .gt. abs(crminus) ) then + dfz(i,j,k,l) = clplus + else + dfz(i,j,k,l) = crminus + end if + endif + + if ( f(i,j,k,l) .lt. 0 ) then + if ( abs(alminus) .gt. abs(arplus) ) then + dfx(i,j,k,l) = alminus + else + dfx(i,j,k,l) = arplus + end if + if ( abs(blminus) .gt. abs(brplus) ) then + dfy(i,j,k,l) = blminus + else + dfy(i,j,k,l) = brplus + end if + if ( abs(clminus) .gt. abs(crplus) ) then + dfz(i,j,k,l) = clminus + else + dfz(i,j,k,l) = crplus + end if + end if +else + dfx(i,j,k,l) = zero + dfy(i,j,k,l) = zero + dfz(i,j,k,l) = zero +end if -- cgit v1.2.3