diff options
Diffstat (limited to 'src/setupbrilldata2D.F')
| -rw-r--r-- | src/setupbrilldata2D.F | 89 |
1 files changed, 58 insertions, 31 deletions
diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F index 709370d..a7b35ff 100644 --- a/src/setupbrilldata2D.F +++ b/src/setupbrilldata2D.F @@ -1,3 +1,12 @@ +/*@@ + @file setupbrilldata2D.F + @date + @author Carsten Gundlach (Cactus 4, Miguel Alcubierre) + @desc + Set up axisymmetric Brill data for elliptic solve. + @enddesc + @version $Header$ +@@*/ #include "cctk.h" #include "cctk_Parameters.h" @@ -7,20 +16,18 @@ subroutine setupbrilldata2D(CCTK_FARGUMENTS) -C Author: Carsten Gundlach. -C -C Set up axisymmetric Brill data for elliptic solve. The elliptic -C equation we need to solve is: -C -C __ 2 2 -C \/ psi = psi ( d q + d q ) / 4 -C f z rho -C -C where: -C -C __ -C \/ = Flat space Laplacian. -C f +c Set up axisymmetric Brill data for elliptic solve. The elliptic +c equation we need to solve is: +c +c __ 2 2 +c \/ psi + psi ( d q + d q ) / 4 = 0 +c f z rho +c +c where: +c +c __ +c \/ = Flat space Laplacian. +c f implicit none @@ -29,6 +36,9 @@ C f integer i,j,k integer nx,ny,nz + integer ierr + + integer, dimension(3) :: sym CCTK_REAL x1,y1,z1,rho1 CCTK_REAL brillq,eps @@ -36,26 +46,26 @@ C f external brillq -C Numbers. +c Numbers. zero = 0.0D0 one = 1.0D0 -C Set up grid size. +c Set up grid size. nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) -C Parameters. +c Parameters. eps = brill_eps -C Initialize psi. +c Initialize psi. brillpsi = one -C Initialize metric. +c Initialize metric. psi = one @@ -67,7 +77,19 @@ C Initialize metric. gxz = zero gyz = zero -C Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho). +c Define the symmetries for the brill GFs. + + sym(1) = 1 + sym(2) = 1 + sym(3) = 1 + + call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillMlinear') + call SetCartSymVN(ierr,cctkGH,sym,'idbrilldata::brillNsource') + +c Set up coefficient arrays for elliptic solve. +c Notice that the Cactus conventions are: +c __ +c \/ psi + Mlinear*psi + Nsource = 0 do k=1,nz do j=1,ny @@ -79,25 +101,30 @@ C Set up coefficient Mlinear = - (1/4) (q,zz + q,rhorho). rho1 = dsqrt(x1**2 + y1**2) -C Centered derivatives. Note that here we may be calling brillq -C with a small negative rho, but that should be ok as long as -C brillq is even in rho - physically it must be, or the data -C will not be regular on the axis. +c Centered derivatives. Note that here we may be calling brillq +c with a small negative rho, but that should be ok as long as +c brillq is even in rho - physically it must be, or the data +c will not be regular on the axis. - brillMlinear(i,j,k) = - 0.25d0 - $ *(brillq(rho1,z1+eps,zero) - $ + brillq(rho1,z1-eps,zero) - $ + brillq(rho1+eps,z1,zero) - $ + brillq(rho1-eps,z1,zero) - $ - 4.d0*brillq(rho1,z1,zero))/eps**2 + brillMlinear(i,j,k) = 0.25d0 + . *(brillq(rho1,z1+eps,zero) + . + brillq(rho1,z1-eps,zero) + . + brillq(rho1+eps,z1,zero) + . + brillq(rho1-eps,z1,zero) + . - 4.d0*brillq(rho1,z1,zero))/eps**2 end do end do end do -C Set coefficient Nsource = 0. +c Set coefficient Nsource = 0. brillNsource = zero +c Synchronization and boundaries. + + call CCTK_SyncGroup(cctkGH,'idbrilldata::brillelliptic') + call ADM_Boundary(cctkGH,'idbrilldata::brillelliptic') + return end |