1 Texas Evolution code(T2) development

Randy Correll Mijan Huq Scott Klasky Richard Matzner January 11, 1996

1. Equations used in code.

2. Modules in DAGH and HPF.

3. Derivative routines.

4. Outer boundary.

5. Gauge conditions.

6. Future Plans.

2 Equations used in code (1)

The evolution equations that the TCODE solves are as follows:

∂_{t}g_{ij} = −2αK_{ij} + D_{i}β_{j} + D_{j}β_{i}  (1)
∂_{t}K_{j}^{i} = L_{β}K_{j}^{i} − D^{i}D_{j}α + α (R_{j}^{i} + KK_{j}^{i}) (2)

where g_{ij} is the 3-metric, K_{j}^{i} is the extrinsic curvature, D_{i} is the covariant derivative compatible with g_{ij}, α is the lapse function, β^{i} is the shift vector and L_{β} is the Lie derivative along the vector field β. L_{β}K_{j}^{i} is given by the following equation:

L_{β}K_{j}^{i} = β^{k}D_{k}K_{j}^{i} + β^{i}D_{k}K_{j}^{k} − β_{j}g^{kl}D_{l}K_{k}^{i} (3)

3 Variable Storage

Grid Function Number of independent components
g_{ij}  6 × 2 independent components(two = n+1 timestep)
g^{ij}  6 independent components(required for various stages)
K_{j}^{i}  9 × 2 components (multiply by two for n+1 timestep)
Γ_{jk}^{i}  18 components (stored instead of g_{ij,k})
R_{ij}  6 components (Raised when needed in eqn. (2))
α  scalar function
β^{i}  3 components (β_{i} can be obtained using g_{ij})

Hence the lower limit on required variables for the solution to the evolution is 64 grid variables. Plans are in the work for writing the equations using K_{ij} thus eliminating 6 variables.

4 Modules for the connection coefficients

The calculation of the connection coefficients can be broken up into several stages. This obviates the necessity for storing the first derivatives of the 3-metric which saves us 18 grid functions.

The connection coefficients are given by

Γ_{jk}^{i} = \frac{1}{2}g^{il} (g_{lj,k} + g_{lk,j} − g_{jk,l})  (4)

The calculation can be broken up into three separate sets using 9 temporary variables.

There are 3 steps:

5 3 Steps for Connection coefficients

For the first term in eqn. (4)

Γ_{jk}^{i} = \frac{1}{2}g^{il}g_{lj,k}
= \frac{1}{2} (g^{i1}g_{1j,k} + g^{i2}g_{2j,k} + g^{i3}g_{3j,k})  (5)

(6)

We call these parts C1-C3. For the second term in eqn. (4)

Γ_{jk}^{i} = Γ_{jk}^{i} + \frac{1}{2}g^{il}g_{lk,j}
= Γ_{jk}^{i} + \frac{1}{2} (g^{i1}g_{1k,j} + g^{i2}g_{2k,j} + g^{i3}g_{3k,j})  (7)

(8)

We call these parts C4-C6. For the third term in eqn. (4)

Γ_{jk}^{i} = Γ_{jk}^{i} − \frac{1}{2}g^{il}g_{jk,l}
= Γ_{jk}^{i} − \frac{1}{2} (g^{i1}g_{jk,1} + g^{i2}g_{jk,2} + g^{i3}g_{jk,3})  (9)

We call these parts C7-C9.

6 Modules for part C1

call diff(nx,ny,nz,h,g11,t1,'x',m,d,nl)
call diff(nx,ny,nz,h,g11,t2,'y',m,d,nl)
call diff(nx,ny,nz,h,g11,t3,'z',m,d,nl)
call diff(nx,ny,nz,h,g12,t4,'x',m,d,nl)
call diff(nx,ny,nz,h,g12,t5,'y',m,d,nl)
call diff(nx,ny,nz,h,g12,t6,'z',m,d,nl)
call diff(nx,ny,nz,h,g13,t7,'x',m,d,nl)
call diff(nx,ny,nz,h,g13,t8,'y',m,d,nl)
call diff(nx,ny,nz,h,g13,t9,'z',m,d,nl)
call connect1(g11, g12, g13, g22, g23, g33,
*                gup11, gup12, gup13, gup22, gup23, gup33,
*                c111, c112, c113, c122, c123, c133,
*                c211, c212, c213, c222, c223, c233,
*                c311, c312, c313, c322, c323, c333,
*                t1,   t2  ,   t3,   t4,   t5,   t6,
*                t7,   t8  ,   t9, h, nx , ny , nz)
                     .
                     .
                     .
c(i,j,k) = c(i,j,k) + 0.5d0*gup(i,j) * t(j,k)

Similar looking modules for other parts.

7 Ricci tensor

The calculation of the Ricci tensor can be similarly broken up into several segments with the use of 9 temporary grid functions once more. The Ricci tensor is given by

R_{ij} = Γ_{ij,a}^{a} − Γ_{ia,j}^{a} + Γ_{ab}^{a}Γ_{ij}^{b} − Γ_{bj}^{a}Γ_{ia}^{b}  (10)

The evaluation of the Ricci tensor is broken as follows:

R_{ij} = Γ_{ij,a}^{a}
= Γ_{ij,1}^{1} + Γ_{ij,2}^{2} + Γ_{ij,3}^{3}  (11)
R_{ij} = R_{ij} − Γ_{ia,j}^{a}
= R_{ij} − (Γ_{i1,j}^{1} + Γ_{i2,j}^{2} + Γ_{i3,j}^{3} )  (12)
R_{ij} = R_{ij} + Γ_{ab}^{a}Γ_{ij}^{b} − Γ_{bj}^{a}Γ_{ia}^{b}  (13)

where in eqns. (10,11) terms like Γ_{ij,1}^{1} and Γ_{i1,j}^{1} are calculated using temporary arrays.

8 Ricci tensor parts 1 and 7

   call ricci1(c111, c112, c113, c122, c123, c133,
*              c211, c212, c213, c222, c223, c233,
*              c311, c312, c313, c322, c323, c333,
*              t1, t2, t3, t4, t5, t6, t7, t8, t9,
*              r11, r12, r13, r22, r23, r33,
*              m,d,nl, h, nx, ny, nz)
   .
   .
   call ricci7(c111, c112, c113, c122, c123, c133,
*              c211, c212, c213, c222, c223, c233,
*              c311, c312, c313, c322, c323, c333,
*              t1, t2, t3, t4, t5, t6, t7, t8, t9,
*              r11, r12, r13, r22, r23, r33,
*              m,d,nl, h, nx, ny, nz)
 return
 end
}

9 Ricci tensor parts 1 and 7

C      ricci1 contains:
call diff(nx, ny, nz, h, c111, c111d1, 'x',m,d,nl)
call diff(nx, ny, nz, h, c112, c112d1, 'x',m,d,nl)
call diff(nx, ny, nz, h, c113, c113d1, 'x',m,d,nl)
call diff(nx, ny, nz, h, c122, c122d1, 'x',m,d,nl)
call diff(nx, ny, nz, h, c123, c123d1, 'x',m,d,nl)
call diff(nx, ny, nz, h, c133, c133d1, 'x',m,d,nl)
r11 = c111d1
r12 = c112d1
r13 = c113d1
r22 = c122d1
r23 = c123d1
r33 = c133d1
C     ricci7 contains:
FORALL (j=1:3,i=1:3,i<=j)
   r(i,j) = r(i,j) +
.              ( c(1,1,1)+c(2,2,1)+c(3,3,1) ) * c(1,i,j) +
.              ( c(1,1,2)+c(2,2,2)+c(3,3,2) ) * c(2,i,j) +
.              ( c(1,1,3)+c(2,2,3)+c(3,3,3) ) * c(3,i,j)
END FORALL
   r(i,j) = r(i,j)
.            - (c(1,i,1) * c(1,1,j) + c(2,i,1) * c(1,2,j)
.              +c(3,i,1) * c(1,3,j) )
.            - (c(1,i,2) * c(2,1,j) + c(2,i,2) * c(2,2,j)
.              +c(3,i,2) * c(2,3,j) )
.            - (c(1,i,3) * c(3,1,j) + c(2,i,3) * c(3,2,j)
.              +c(3,i,3) * c(3,3,j) )

10 Evolution equation for the 3-metric

∂_{t}g_{ij} = −2αK_{ij} + D_{i}β_{j} + D_{j}β_{i}.  (14)

The evaluation of the second and third terms is carried out by once again using 9 temporary grid functions. Recall that β^{i}, the contravariant form of the shift, is stored. Since g_{ij} is compatible with D_{i} we can rewrite D_{i}β_{j} = g_{jk}D_{i}β^{k}. Hence we are left with calculating

D_{i}β_{j} = g_{jk} (∂_{i}β^{k} − Γ_{il}^{k}β^{l}) .  (15)

Here once again we use a 9 grid functions to store the components of ∂_{i}β^{k} first and then add Γ_{il}^{k}β^{l} to form eqn. (15). Then evaluating eqn. (1) becomes a matter of using the temporary variables for D_{i}β_{j}.

Use second order finite differencing for time derivative. 5 Modules for these routines

C in evolve_g1
      call diff(nx,ny,nz,h,beta1,t11,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,t12,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta1,t13,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t21,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t22,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta2,t23,'z',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t31,'x',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t32,'y',ml,dl,nl)
      call diff(nx,ny,nz,h,beta3,t33,'z',ml,dl,nl)

11 5 Modules for evolution eqn

C in evolve_g2
FORALL (i=1:3,j=1:3)
   t(i,j) = t(i,j) +
.              c(1,i,j) * beta(1) +
.              c(2,i,j) * beta(2) +
.              c(3,i,j) * beta(3)
END FORALL
FORALL (i=1:3,j=1:3)
   t(i,j) = t(i,1) * g(1,j) + t(i,2) * g(2,j) + t(i,3) * g(3,j)
END FORALL
C in evolve_g3
FORALL (j=1:3,i=1:3,i<=j)
   t(i,j) = t(i,j) + t(j,i)
   t(j,i) = t(i,j)
END FORALL
C in evolve_g4
FORALL (i=1:3,j=1:3)
   t(i,j) = t(i,j) - 2.0d0*alpha(1) * ( g(1,i) * k(1,j) +
.            g(2,i) * k(2,j) + g(3,i) * k(3,j)  )
END FORALL
C in evolve_g5
FORALL (j=1:3,i=1:3,i<=j)
   g(i,j,1) = 2.0d0*dt*t(i,j) + g(i,j,-1)
END FORALL

12 Evolution equation for K_{j}^{i}

∂_{t}K_{j}^{i} = L_{β}K_{j}^{i} − D^{i}D_{j}α + α (R_{j}^{i} + KK_{j}^{i})

We break this equation up into 5 parts(5 modules)

C evolve_k1
call diff(nx,ny,nz,h,alpha,temp1,'x',ml,dl,nl)
call diff(nx,ny,nz,h,alpha,temp2,'y',ml,dl,nl)
call diff(nx,ny,nz,h,alpha,temp3,'z',ml,dl,nl)
C evolve_k2
call diff(nx,ny,nz,h,temp1,t11,'x',ml,dl,nl)
call diff(nx,ny,nz,h,temp1,t12,'y',ml,dl,nl)
call diff(nx,ny,nz,h,temp1,t13,'z',ml,dl,nl)
call diff(nx,ny,nz,h,temp2,t22,'y',ml,dl,nl)
call diff(nx,ny,nz,h,temp2,t23,'z',ml,dl,nl)
call diff(nx,ny,nz,h,temp3,t33,'z',ml,dl,nl)
t21 = t12
t31 = t13
t32 = t23
FORALL (j=1:3,i=1:3,i<=j)
   t(i,j) = t(i,j) - c(1,i,j) *dalpha(1) -
.              c(2,i,j)*dalpha(2) - c(3,i,j) * dalpha(3)
   t(j,i) = t(i,j)
END FORALL
c   raise the first index of t{ab}
call rais02i1()
FORALL (i=1:3,j=1:3)
   t(i,j) = -t(i,j)
END FORALL

13 Evolution equation for K_{j}^{i} continued

! Evolve K3
*     Lie_beta(K^1_1)
call diff(nx,ny,nz,h,k11,temp1,'x',ml,dl,nl)
call diff(nx,ny,nz,h,k11,temp2,'y',ml,dl,nl)
call diff(nx,ny,nz,h,k11,temp3,'z',ml,dl,nl)
t11 = t11 + beta1 * temp1  + beta2 * temp2 + beta3 * temp3
call diff(nx,ny,nz,h,beta1,temp1,'x',ml,dl,nl)
call diff(nx,ny,nz,h,beta1,temp2,'y',ml,dl,nl)
call diff(nx,ny,nz,h,beta1,temp3,'z',ml,dl,nl)
t11 = t11 - k11 * temp1 - k21 * temp2 - k31 * temp3
call diff(nx,ny,nz,h,beta1,temp1,'y',ml,dl,nl)
call diff(nx,ny,nz,h,beta2,temp2,'y',ml,dl,nl)
call diff(nx,ny,nz,h,beta3,temp3,'y',ml,dl,nl)
t11 = t11 + k11 * temp1 + k12 * temp2 + k13 * temp3
.
C Evolve K4
trace_k(1) = k11 + k22 + k33
FORALL (i=1:3,j=1:3)
   t(i,j) = t(i,j) + alpha(1) *(
.              gup(i,1) * r(1,j) + gup(i,2) * r(2,j) +
.              gup(i,3) * r(3,j) + trace_k(1) * k(i,j) )
END FORALL
C Evolve K5
FORALL (i=1:3,j=1:3)
   k(i,j,1) = 2.0d0*dt*t(i,j) + k(i,j,-1)
END FORALL

14 Derivative routines, Vertex Centered

Spatial derivative routines contain routines for Inner boundary, Outer boundary and interior. Time derivatives, O(h^{2}) use a 3 level leap frog scheme. All derivative routines are second order for the first derivative in x, y, z directions.

To handle the hole, 2 non-distributed lists are created, m(nl),d(nl) where m(l) contains the i,j,k location of the “hole” points and d(l) contains either (1,-1,0) indicating either forward, backwards or centered.

The derivatives around the hole are determined to second order. Second derivatives, for example f_{xx}, are only first order accurate. The interior second derivatives are second order accurate. Routine looks as follows:

      subroutine diff(h,f,df,dir,m,d,nl)
      use params
      implicit none
      real*8 h
      real*8, dimension(nx,ny,nz) :: f,df
      real*8, dimension(3*nl)     :: d
      integer, dimension(3*nl)    :: m
      character dir
!HPF$ PROCESSORS PROC(nproc)
!HPF$ DISTRIBUTE f *(BLOCK,*,*) ONTO *PROC
!HPF$ ALIGN(:,:,:) WITH *f(:,:,:)         :: df
      call diff_int(h,f,df,dir)
      call diff_ob(h,f,df,dir)
      if (nl.gt.1) call diff_ib(h,f,df,dir,m,d,nl)
      return
      end

subroutine diffx_ib(h,f,df,m,d,nl)
use params
implicit none
integer nl
real*8, dimension(nx,ny,nz) :: f,df
real*8, dimension(3*nl)     :: d
integer, dimension(3*nl)    :: m
real*8 h,h2
integer i,j,k,l
h2 = 0.5d0 / h
do l = 1 , nl
   i  = m(l)
   j  = m(nl+l)
   k  = m(2*nl+l)
   d1 = d(l)
   d2 = d(nl+l)
   d3 = d(2*nl+l)
   da = abs(d1)
   df(i,j,k) = h2*(d1*(-3.0d0*f(i,j,k)+4.0d0*f(i+int(d1),j,k)
$                       -      f(i+2*int(d1),j,k)) +
$               (1.0d0-da)*(f(i+1,j,k)-f(i-1,j,k) ) )
end do
return
end

subroutine diffx_ob(h,f,df)
use params
implicit none
real*8, dimension(nx,ny,nz) :: f,df
real*8  h,h2,hi
integer i,j,k
!HPF$ PROCESSORS PROC(nproc)
!HPF$ DISTRIBUTE f *(BLOCK,*,*) ONTO *PROC
!HPF$ ALIGN(:,:,:) WITH *f(:,:,:)         :: df
hi = 1.0d0/h
h2 = 0.5d0 / h
FORALL (k=1:nz,j=1:ny)
   df(1,j,k) = h2 * ( -3.0d0*f(1,j,k)
.             + 4.0d0*f(1+1,j,k) - f(1+2,j,k) )
   df(nx,j,k) = h2 * (  3.0d0*f(nx,j,k)
.              - 4.0d0*f(nx-1,j,k) + f(nx-2,j,k) )
END FORALL
FORALL (i=2:nx-1,k=1:nz)
      df(i,1,k) = h2 *(f(i+1,1,k)-f(i-1,1,k) )
      df(i,ny,k) = h2 *(f(i+1,ny,k)-f(i-1,ny,k) )
END FORALL
FORALL (i=2:nx-1,j=1:ny)
      df(i,j,1) = h2 *(f(i+1,j,1)-f(i-1,j,1) )
      df(i,j,nz) = h2 *(f(i+1,j,nz)-f(i-1,j,nz) )
END FORALL
return
end

15 Outer boundary

• Extrapolation conditions have been looked at: R. Correll

• One sided differencing conditions have been used: Correll, et. al.

• Outer boundary is not “killing” code, it’s the inner boundary.

• Use some form of characteristic code when it becomes available for outer boundary (future)

16 Gauge conditions

The following gauge conditions have been imposed:

• Geometric

• Algebraic (in future)

• Quasi-Isotropic , Geodesic slicing ( ran to about 20M)

17 Future Plans

• Moving the hole on the domain, use some sort of Alcubierre and Schutz method(implicit scheme using ADI). ( Huq & Matzner)

• Apparent Horizon solver written in a Data Parallel language. (Huq, Klasky & Matzner)

• Causal differencing around the inner boundary.

• Stability of methods around the inner boundary is being looked at by

• AMR in DAGH and HPF : Haupt,Huq,Klasky,Parashar Correll