************************************************************************ * Computes the residual of the hamiltonian and momentum constraint * equations. * Randall R. Correll * Mijan F. Huq * Scott A. Klasky * Richard A. Matzner * * COPYRIGHT 1995-1996 ************************************************************************ subroutine constraint(nx,ny,nz,dx,dy,dz,mask, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ hamnorm,mom1norm,mom2norm,mom3norm, $ c111,c112,c113,c122,c123,c133, $ c211,c212,c213,c222,c223,c233, $ c311,c312,c313,c322,c323,c333, $ r11,r12,r13,r22,r23,r33, $ gup11,gup12,gup13,gup22,gup23,gup33, $ t1,t2,t3,t4,t5,t6,t7,t8,t9,temp1, $ h,m,d,nl) implicit none integer nx, ny, nz, nl real*8 dx,dy,dz,h real*8 mask(nx,ny,nz) c extrinsic curvature: K^i_j, Kold^i_j and tr(k) real*8 k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz), & k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz), & k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz) real*8 r11(nx,ny,nz), r12(nx,ny,nz), r13(nx,ny,nz), $ r22(nx,ny,nz), r23(nx,ny,nz), r33(nx,ny,nz) real*8 gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz), $ gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz) real*8 hamnorm,mom1norm,mom2norm,mom3norm, n2i c connection coefficients: C^i_jk real*8 c111(nx,ny,nz), c112(nx,ny,nz), c113(nx,ny,nz), & c122(nx,ny,nz), c123(nx,ny,nz), c133(nx,ny,nz), & c211(nx,ny,nz), c212(nx,ny,nz), c213(nx,ny,nz), & c222(nx,ny,nz), c223(nx,ny,nz), c233(nx,ny,nz), & c311(nx,ny,nz), c312(nx,ny,nz), c313(nx,ny,nz), & c322(nx,ny,nz), c323(nx,ny,nz), c333(nx,ny,nz) real*8 t1(nx,ny,nz), t2(nx,ny,nz), t3(nx,ny,nz), $ t4(nx,ny,nz), t5(nx,ny,nz), t6(nx,ny,nz), $ t7(nx,ny,nz), t8(nx,ny,nz), t9(nx,ny,nz), $ temp1(nx,ny,nz) real*8 d(3*nl) integer m(3*nl) c local variables integer i, j, k real*8 sum hamnorm = 0.0d0 mom1norm = 0.0d0 mom2norm = 0.0d0 mom3norm = 0.0d0 sum = 0.0d0 do k = 1, nz do j = 1, ny do i = 1, nx c compute hamiltonian error according to the template below: c R + (trK)^2 - K^a_b K^b_a = 0 t7(i,j,k) = k11(i,j,k)+k22(i,j,k)+k33(i,j,k) temp1(i,j,k) = mask(i,j,k)*( $ (gup11(i,j,k)*r11(i,j,k) $ +gup12(i,j,k)*r12(i,j,k) * 2.0e0 $ +gup13(i,j,k)*r13(i,j,k) * 2.0e0 $ +gup22(i,j,k)*r22(i,j,k) $ +gup23(i,j,k)*r23(i,j,k) * 2.0e0 $ +gup33(i,j,k)*r33(i,j,k)) & + t7(i,j,k)*t7(i,j,k) & + ( k11(i,j,k)*k11(i,j,k) & + k12(i,j,k)*k21(i,j,k) & + k13(i,j,k)*k31(i,j,k) & + k21(i,j,k)*k12(i,j,k) & + k22(i,j,k)*k22(i,j,k) & + k23(i,j,k)*k32(i,j,k) & + k31(i,j,k)*k13(i,j,k) & + k32(i,j,k)*k23(i,j,k) & + k33(i,j,k)*k33(i,j,k) ) & ) end do end do end do c compute momentum error according to the template below: c D_b(K^b_a) - D_a(k^b_b) = 0 c moma(i,j,k) = mask(i,j,k)*( c & + k1ad1(i,j,k) + k2ad2(i,j,k) + k3ad3(i,j,k) c & + ( c111(i,j,k)*k1a(i,j,k) c & + c112(i,j,k)*k2a(i,j,k) c & + c113(i,j,k)*k3a(i,j,k) c & + c212(i,j,k)*k1a(i,j,k) c & + c222(i,j,k)*k2a(i,j,k) c & + c223(i,j,k)*k3a(i,j,k) c & + c313(i,j,k)*k1a(i,j,k) c & + c323(i,j,k)*k2a(i,j,k) c & + c333(i,j,k)*k3a(i,j,k) ) c & - ( c1a1(i,j,k)*k11(i,j,k) c & + c2a1(i,j,k)*k12(i,j,k) c & + c3a1(i,j,k)*k13(i,j,k) c & + c1a2(i,j,k)*k21(i,j,k) c & + c2a2(i,j,k)*k22(i,j,k) c & + c3a2(i,j,k)*k23(i,j,k) c & + c1a3(i,j,k)*k31(i,j,k) c & + c2a3(i,j,k)*k32(i,j,k) c & + c3a3(i,j,k)*k33(i,j,k) ) c & - ( k11da(i,j,k) + k22da(i,j,k) + k33da(i,j,k) ) c & ) call diff(nx,ny,nz,dx,k11,t1,'x',m,d,nl) call diff(nx,ny,nz,dx,k21,t2,'y',m,d,nl) call diff(nx,ny,nz,dx,k31,t3,'z',m,d,nl) call diff(nx,ny,nz,dx,k11,t4,'x',m,d,nl) call diff(nx,ny,nz,dx,k22,t5,'x',m,d,nl) call diff(nx,ny,nz,dx,k33,t6,'x',m,d,nl) do k = 1, nz do j = 1, ny do i = 1, nx t7(i,j,k) = mask(i,j,k)*( & + t1(i,j,k) + t2(i,j,k) + t3(i,j,k) & + ( c111(i,j,k)*k11(i,j,k) & + c112(i,j,k)*k21(i,j,k) & + c113(i,j,k)*k31(i,j,k) & + c212(i,j,k)*k11(i,j,k) & + c222(i,j,k)*k21(i,j,k) & + c223(i,j,k)*k31(i,j,k) & + c313(i,j,k)*k11(i,j,k) & + c323(i,j,k)*k21(i,j,k) & + c333(i,j,k)*k31(i,j,k) ) & - ( c111(i,j,k)*k11(i,j,k) & + c211(i,j,k)*k12(i,j,k) & + c311(i,j,k)*k13(i,j,k) & + c112(i,j,k)*k21(i,j,k) & + c212(i,j,k)*k22(i,j,k) & + c312(i,j,k)*k23(i,j,k) & + c113(i,j,k)*k31(i,j,k) & + c213(i,j,k)*k32(i,j,k) & + c313(i,j,k)*k33(i,j,k) ) & - ( t4(i,j,k) + t5(i,j,k) + t6(i,j,k) ) & ) end do end do end do call diff(nx,ny,nz,dx,k12,t1,'x',m,d,nl) call diff(nx,ny,nz,dx,k22,t2,'y',m,d,nl) call diff(nx,ny,nz,dx,k32,t3,'z',m,d,nl) call diff(nx,ny,nz,dx,k11,t4,'y',m,d,nl) call diff(nx,ny,nz,dx,k22,t5,'y',m,d,nl) call diff(nx,ny,nz,dx,k33,t6,'y',m,d,nl) do k = 1, nz do j = 1, ny do i = 1, nx t8(i,j,k) = mask(i,j,k)*( & + t1(i,j,k) + t2(i,j,k) + t3(i,j,k) & + ( c111(i,j,k)*k12(i,j,k) & + c112(i,j,k)*k22(i,j,k) & + c113(i,j,k)*k32(i,j,k) & + c212(i,j,k)*k12(i,j,k) & + c222(i,j,k)*k22(i,j,k) & + c223(i,j,k)*k32(i,j,k) & + c313(i,j,k)*k12(i,j,k) & + c323(i,j,k)*k22(i,j,k) & + c333(i,j,k)*k32(i,j,k) ) & - ( c112(i,j,k)*k11(i,j,k) & + c212(i,j,k)*k12(i,j,k) & + c312(i,j,k)*k13(i,j,k) & + c122(i,j,k)*k21(i,j,k) & + c222(i,j,k)*k22(i,j,k) & + c322(i,j,k)*k23(i,j,k) & + c123(i,j,k)*k31(i,j,k) & + c223(i,j,k)*k32(i,j,k) & + c323(i,j,k)*k33(i,j,k) ) & - ( t4(i,j,k) + t5(i,j,k) + t6(i,j,k) ) & ) end do end do end do call diff(nx,ny,nz,dx,k13,t1,'x',m,d,nl) call diff(nx,ny,nz,dx,k23,t2,'y',m,d,nl) call diff(nx,ny,nz,dx,k33,t3,'z',m,d,nl) call diff(nx,ny,nz,dx,k11,t4,'z',m,d,nl) call diff(nx,ny,nz,dx,k22,t5,'z',m,d,nl) call diff(nx,ny,nz,dx,k33,t6,'z',m,d,nl) do k = 1, nz do j = 1, ny do i = 1, nx t9(i,j,k) = mask(i,j,k)*( & + t1(i,j,k) + t2(i,j,k) + t3(i,j,k) & + ( c111(i,j,k)*k13(i,j,k) & + c112(i,j,k)*k23(i,j,k) & + c113(i,j,k)*k33(i,j,k) & + c212(i,j,k)*k13(i,j,k) & + c222(i,j,k)*k23(i,j,k) & + c223(i,j,k)*k33(i,j,k) & + c313(i,j,k)*k13(i,j,k) & + c323(i,j,k)*k23(i,j,k) & + c333(i,j,k)*k33(i,j,k) ) & - ( c113(i,j,k)*k11(i,j,k) & + c213(i,j,k)*k12(i,j,k) & + c313(i,j,k)*k13(i,j,k) & + c123(i,j,k)*k21(i,j,k) & + c223(i,j,k)*k22(i,j,k) & + c323(i,j,k)*k23(i,j,k) & + c133(i,j,k)*k31(i,j,k) & + c233(i,j,k)*k32(i,j,k) & + c333(i,j,k)*k33(i,j,k) ) & - ( t4(i,j,k) + t5(i,j,k) + t6(i,j,k) ) & ) end do end do end do c compute L2-norm of residuals c do k = 1, nz c do j = 1, ny c do i = 1, nx c c hamnorm = hamnorm + temp1(i,j,k)*temp1(i,j,k) c mom1norm = mom1norm + t7(i,j,k)*t7(i,j,k) c mom2norm = mom2norm + t8(i,j,k)*t8(i,j,k) c mom3norm = mom3norm + t9(i,j,k)*t9(i,j,k) c sum = sum + mask(i,j,k) c c end do c end do c end do c c hamnorm = sqrt(hamnorm/sum) c mom1norm = sqrt(mom1norm/sum) c mom2norm = sqrt(mom2norm/sum) c mom3norm = sqrt(mom3norm/sum) hamnorm = n2i(nx, temp1) mom1norm = n2i(nx, t7) mom2norm = n2i(nx, t8) mom3norm = n2i(nx, t9) return end