diff options
Diffstat (limited to 'src/CalcK.F90')
| -rw-r--r-- | src/CalcK.F90 | 157 |
1 files changed, 157 insertions, 0 deletions
diff --git a/src/CalcK.F90 b/src/CalcK.F90 new file mode 100644 index 0000000..f77d500 --- /dev/null +++ b/src/CalcK.F90 @@ -0,0 +1,157 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine CalcK (CCTK_ARGUMENTS) + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_INT, parameter :: idummy=0 + integer, parameter :: ik=kind(idummy) + + CCTK_INT, parameter :: bndwidth = 1 + + CCTK_INT :: len_boundary, len_options + character(1000) :: str_boundary, str_options + integer :: options + + CCTK_REAL :: dt, dx(3) + + CCTK_REAL :: gama(3,3), gama_dot(3,3), dgama(3,3,3) + CCTK_REAL :: alfa + CCTK_REAL :: beta(3), dbeta(3,3) + CCTK_REAL :: kk(3,3) + + integer :: imin(3), imax(3) + integer :: i, j, k + integer :: a, b, c + integer :: ierr + + dt = CCTK_DELTA_TIME + dx(:) = CCTK_DELTA_SPACE(:) + + call CCTK_FortranString & + (len_boundary, int(extcurv_boundary,ik), str_boundary) + call CCTK_FortranString & + (len_options, int(extcurv_boundary_options,ik), str_options) + call Util_TableCreateFromString (options, str_options) + if (options<0) then + call CCTK_WARN (0, "Parameter ""extcurv_boundary_options"" has an illegal syntax") + end if + + imin(:) = 1+cctk_nghostzones(:) + imax(:) = cctk_lsh(:)-cctk_nghostzones(:) + where (cctk_bbox(1::2)/=0) imin(:) = 1+bndwidth + where (cctk_bbox(2::2)/=0) imax(:) = cctk_lsh(:)-bndwidth + + do k = imin(3), imax(3) + do j = imin(2), imax(2) + do i = imin(1), imax(1) + + gama(1,1) = gxx(i,j,k) + gama(1,2) = gxy(i,j,k) + gama(1,3) = gxz(i,j,k) + gama(2,2) = gyy(i,j,k) + gama(2,3) = gyz(i,j,k) + gama(3,3) = gzz(i,j,k) + gama(2,1) = gama(1,2) + gama(3,1) = gama(1,3) + gama(3,2) = gama(2,3) + +#if 0 + gama_dot(1,1) = (-3*gxx(i,j,k) + 4*gxx_prev(i,j,k) - gxx_prev2(i,j,k)) / (2*dt) + gama_dot(1,2) = (-3*gxy(i,j,k) + 4*gxy_prev(i,j,k) - gxy_prev2(i,j,k)) / (2*dt) + gama_dot(1,3) = (-3*gxz(i,j,k) + 4*gxz_prev(i,j,k) - gxz_prev2(i,j,k)) / (2*dt) + gama_dot(2,2) = (-3*gyy(i,j,k) + 4*gyy_prev(i,j,k) - gyy_prev2(i,j,k)) / (2*dt) + gama_dot(2,3) = (-3*gyz(i,j,k) + 4*gyz_prev(i,j,k) - gyz_prev2(i,j,k)) / (2*dt) + gama_dot(3,3) = (-3*gzz(i,j,k) + 4*gzz_prev(i,j,k) - gzz_prev2(i,j,k)) / (2*dt) +#else + gama_dot(1,1) = (gxx_next(i,j,k) - gxx_prev(i,j,k)) / (2*dt) + gama_dot(1,2) = (gxy_next(i,j,k) - gxy_prev(i,j,k)) / (2*dt) + gama_dot(1,3) = (gxz_next(i,j,k) - gxz_prev(i,j,k)) / (2*dt) + gama_dot(2,2) = (gyy_next(i,j,k) - gyy_prev(i,j,k)) / (2*dt) + gama_dot(2,3) = (gyz_next(i,j,k) - gyz_prev(i,j,k)) / (2*dt) + gama_dot(3,3) = (gzz_next(i,j,k) - gzz_prev(i,j,k)) / (2*dt) +#endif + gama_dot(2,1) = gama_dot(1,2) + gama_dot(3,1) = gama_dot(1,3) + gama_dot(3,2) = gama_dot(2,3) + + dgama(1,1,1) = (gxx(i+1,j,k) - gxx(i-1,j,k)) / (2*dx(1)) + dgama(1,2,1) = (gxy(i+1,j,k) - gxy(i-1,j,k)) / (2*dx(1)) + dgama(1,3,1) = (gxz(i+1,j,k) - gxz(i-1,j,k)) / (2*dx(1)) + dgama(2,2,1) = (gyy(i+1,j,k) - gyy(i-1,j,k)) / (2*dx(1)) + dgama(2,3,1) = (gyz(i+1,j,k) - gyz(i-1,j,k)) / (2*dx(1)) + dgama(3,3,1) = (gzz(i+1,j,k) - gzz(i-1,j,k)) / (2*dx(1)) + dgama(1,1,2) = (gxx(i,j+1,k) - gxx(i,j-1,k)) / (2*dx(2)) + dgama(1,2,2) = (gxy(i,j+1,k) - gxy(i,j-1,k)) / (2*dx(2)) + dgama(1,3,2) = (gxz(i,j+1,k) - gxz(i,j-1,k)) / (2*dx(2)) + dgama(2,2,2) = (gyy(i,j+1,k) - gyy(i,j-1,k)) / (2*dx(2)) + dgama(2,3,2) = (gyz(i,j+1,k) - gyz(i,j-1,k)) / (2*dx(2)) + dgama(3,3,2) = (gzz(i,j+1,k) - gzz(i,j-1,k)) / (2*dx(2)) + dgama(1,1,3) = (gxx(i,j,k+1) - gxx(i,j,k-1)) / (2*dx(3)) + dgama(1,2,3) = (gxy(i,j,k+1) - gxy(i,j,k-1)) / (2*dx(3)) + dgama(1,3,3) = (gxz(i,j,k+1) - gxz(i,j,k-1)) / (2*dx(3)) + dgama(2,2,3) = (gyy(i,j,k+1) - gyy(i,j,k-1)) / (2*dx(3)) + dgama(2,3,3) = (gyz(i,j,k+1) - gyz(i,j,k-1)) / (2*dx(3)) + dgama(3,3,3) = (gzz(i,j,k+1) - gzz(i,j,k-1)) / (2*dx(3)) + dgama(2,1,:) = dgama(1,2,:) + dgama(3,1,:) = dgama(1,3,:) + dgama(3,2,:) = dgama(2,3,:) + + alfa = alp(i,j,k) + + beta(1) = betax(i,j,k) + beta(2) = betay(i,j,k) + beta(3) = betaz(i,j,k) + + dbeta(1,1) = (betax(i+1,j,k) - betax(i-1,j,k)) / (2*dx(1)) + dbeta(2,1) = (betay(i+1,j,k) - betay(i-1,j,k)) / (2*dx(1)) + dbeta(3,1) = (betaz(i+1,j,k) - betaz(i-1,j,k)) / (2*dx(1)) + dbeta(1,2) = (betax(i,j+1,k) - betax(i,j-1,k)) / (2*dx(2)) + dbeta(2,2) = (betay(i,j+1,k) - betay(i,j-1,k)) / (2*dx(2)) + dbeta(3,2) = (betaz(i,j+1,k) - betaz(i,j-1,k)) / (2*dx(2)) + dbeta(1,3) = (betax(i,j,k+1) - betax(i,j,k-1)) / (2*dx(3)) + dbeta(2,3) = (betay(i,j,k+1) - betay(i,j,k-1)) / (2*dx(3)) + dbeta(3,3) = (betaz(i,j,k+1) - betaz(i,j,k-1)) / (2*dx(3)) + + ! d/dt gamma_ij = - 2 alpha K_ij + ! + gamma_kj d_i beta^k + ! + gamma_kj d_i beta^k + ! + beta^k d_k gamma_ij + + do a=1,3 + do b=1,3 + kk(a,b) = - gama_dot(a,b) + do c=1,3 + kk(a,b) = kk(a,b) + gama(k,j) * dbeta(k,i) & + & + gama(i,k) * dbeta(k,j) & + & + beta(k) * dgama(i,j,k) + end do + kk(a,b) = kk(a,b) / (2*alfa) + end do + end do + + kxx(i,j,k) = kk(1,1) + kxy(i,j,k) = kk(1,2) + kxz(i,j,k) = kk(1,3) + kyy(i,j,k) = kk(2,2) + kyz(i,j,k) = kk(2,3) + kzz(i,j,k) = kk(3,3) + + end do + end do + end do + + call CartSymGN (ierr, cctkGH, "ADMBase::curv") + if (ierr /= 0) call CCTK_WARN (0, "internal error") + ierr = Boundary_SelectGroupForBC (cctkGH, CCTK_ALL_FACES, bndwidth, & + int(options,ik), "ADMBase::curv", str_boundary) + if (ierr /= 0) call CCTK_WARN (0, "internal error") + +end subroutine CalcK |