From 7b342b357a587275f3cc4c047ebc717dcaeb9e7a Mon Sep 17 00:00:00 2001 From: miguel Date: Tue, 17 Apr 2001 09:51:14 +0000 Subject: Cleaning up and fixing the sign convention. It works! git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@40 a678b1cf-93e1-4b43-a69d-d43939e66649 --- src/brilldata.F | 72 +++++++++++++++++++++++----------------- src/brillq.F | 12 +++++-- src/finishbrilldata.F | 53 +++++++++++++++++------------- src/phif.F | 18 ++++++---- src/setupbrilldata2D.F | 89 ++++++++++++++++++++++++++++++++------------------ src/setupbrilldata3D.F | 47 ++++++++++++++++++++------ 6 files changed, 190 insertions(+), 101 deletions(-) (limited to 'src') diff --git a/src/brilldata.F b/src/brilldata.F index b811ee4..04c5eeb 100644 --- a/src/brilldata.F +++ b/src/brilldata.F @@ -1,3 +1,12 @@ +/*@@ + @file brilldata.F + @date + @author Carsten Gundlach (Cactus 4, Miguel Alcubierre) + @desc + Construct Brill wave initial data. + @enddesc + @version $Header$ +@@*/ #include "cctk.h" #include "cctk_Parameters.h" @@ -15,24 +24,35 @@ DECLARE_CCTK_FUNCTIONS integer ipsi,iMcoeff,iNcoeff + integer metpsi_index(7) integer ierr CCTK_REAL AbsTol(3),RelTol(3) +c Get indices for metric. + + call CCTK_VarIndex(metpsi_index(1), "einstein::gxx") + call CCTK_VarIndex(metpsi_index(2), "einstein::gxy") + call CCTK_VarIndex(metpsi_index(3), "einstein::gxz") + call CCTK_VarIndex(metpsi_index(4), "einstein::gyy") + call CCTK_VarIndex(metpsi_index(5), "einstein::gyz") + call CCTK_VarIndex(metpsi_index(6), "einstein::gzz") + call CCTK_VarIndex(metpsi_index(7), "einstein::psi") + c Get indices for grid functions. - call CCTK_VarIndex(ipsi , "idbrilldata::brillpsi") - if (ipsi .lt. 0) then + call CCTK_VarIndex(ipsi,"idbrilldata::brillpsi") + if (ipsi.lt.0) then call CCTK_WARN(0,"Grid variable index for iphi not found") end if - call CCTK_VarIndex(iMcoeff, "idbrilldata::brillMlinear") - if (iMcoeff .lt. 0) then + call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear") + if (iMcoeff.lt.0) then call CCTK_WARN(0,"Grid variable index for Mcoeff not found") end if - call CCTK_VarIndex(iNcoeff, "idbrilldata::brillNsource") - if (iNcoeff .lt. 0) then + call CCTK_VarIndex(iNcoeff,"idbrilldata::brillNsource") + if (iNcoeff.lt.0) then call CCTK_WARN(0,"Grid variable index for Ncoeff not found") end if @@ -54,39 +74,31 @@ c Tolerances for elliptic solve. RelTol(2)= -1 RelTol(3)= -1 -c Elliptic solver. - - if (axisym.eq.1) then - - call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf") - call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff") - - if (CCTK_EQUALS(brill_solver,"sor")) then - call Ell_LinFlatSolver(ierr,cctkGH,ipsi, - . iMcoeff,iNcoeff,AbsTol,RelTol,"sor") - end if - - if (CCTK_EQUALS(brill_solver,"petsc")) then - call Ell_LinFlatSolver(ierr,cctkGH,ipsi, - . iMcoeff,iNcoeff,AbsTol,RelTol,"petsc") - end if +c Boundaries. - if (CCTK_EQUALS(brill_solver,"bam")) then - call Ell_LinFlatSolver(ierr,cctkGH,ipsi, - . iMcoeff,iNcoeff,AbsTol,RelTol,"bam") - end if + call Ell_SetRealKey(ierr,1.0d0,"EllLinConfMetric::Bnd::Robin::inf") + call Ell_SetIntKey(ierr,1,"EllLinConfMetric::Bnd::Robin::falloff") - else +c Elliptic solver. + if (CCTK_EQUALS(brill_solver,"sor")) then + call Ell_LinConfMetricSolver(ierr,cctkGH, + . metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor") + else if (CCTK_EQUALS(brill_solver,"petsc")) then + call Ell_LinConfMetricSolver(ierr,cctkGH, + . metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"petsc") + else if (CCTK_EQUALS(brill_solver,"bam")) then + call Ell_LinConfMetricSolver(ierr,cctkGH, + . metpsi_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"bam") end if c Check for errors. - if (ierr .eq. ELL_SUCCESS) then + if (ierr.eq.ELL_SUCCESS) then call CCTK_INFO("Leaving elliptic solver: solve successful") - else if (ierr .eq. ELL_NOCONVERGENCE) then + else if (ierr.eq.ELL_NOCONVERGENCE) then call CCTK_INFO("Leaving elliptic solver: solver failed to converge") - else if (ierr .eq. ELL_NOSOLVER) then + else if (ierr.eq.ELL_NOSOLVER) then call CCTK_INFO("Elliptic solver not found") else call CCTK_INFO("Leaving elliptic solver: solve failed") diff --git a/src/brillq.F b/src/brillq.F index c7010f6..8d62064 100644 --- a/src/brillq.F +++ b/src/brillq.F @@ -1,11 +1,19 @@ +/*@@ + @file brilldata.F + @date + @author Carsten Gundlach, Miguel Alcubierre. + @desc + Construct Brill wave q function. + @enddesc + @version $Header$ +@@*/ + #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" function brillq(rho,z,phi) -C Authors: Carsten Gundlach, Miguel Alcubierre. -C C Calculates the function q that appear in the conformal C metric for Brill waves: C diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F index ee49b57..36060e9 100644 --- a/src/finishbrilldata.F +++ b/src/finishbrilldata.F @@ -1,3 +1,12 @@ +/*@@ + @file finishbrilldata.F + @date + @author Carsten Gundlach (Cactus 4, Miguel Alcubierre) + @desc + Reconstruct metric from conformal factor. + @enddesc + @version $Header$ +@@*/ #include "cctk.h" #include "cctk_Parameters.h" @@ -22,22 +31,22 @@ external brillq -C Numbers. +c Numbers. zero = 0.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. rhofudge = brill_rhofudge -C Replace flat metric left over from elliptic solve by -C Brill metric calculated from q and Psi. +c Replace flat metric left over from elliptic solve by +c Brill metric calculated from q and Psi. do k=1,nz do j=1,ny @@ -55,31 +64,31 @@ C Brill metric calculated from q and Psi. psi4 = brillpsi(i,j,k)**4 e2q = dexp(2.d0*brillq(rho1,z1,phi)) -C Fudge division by rho^2 on axis. (Physically, y^/rho^2, -C x^2/rho^2 and xy/rho^2 are of course regular.) -C Transform Brills form of the physical metric to Cartesian -C coordinates via -C -C e^2q (drho^2 + dz^2) + rho^2 dphi^2 = -C e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2 -C -C The individual coefficients can be read off as +c Fudge division by rho^2 on axis. (Physically, y^/rho^2, +c x^2/rho^2 and xy/rho^2 are of course regular.) +c Transform Brills form of the physical metric to Cartesian +c coordinates via +c +c e^2q (drho^2 + dz^2) + rho^2 dphi^2 = +c e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2 +c +c The individual coefficients can be read off as if (rho1 .gt. rhofudge) then - gzz(i,j,k) = psi4*e2q gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2) gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2) + gzz(i,j,k) = psi4*e2q gxy(i,j,k) = - psi4*(1.d0 - e2q)*x1*y1/rho2 else -C This fudge assumes that q = O(rho^2) near the axis. Which -C it should be, or the data will be singular. +c This fudge assumes that q = O(rho^2) near the axis. Which +c it should be, or the data will be singular. - gzz(i,j,k) = psi4 - gxx(i,j,k) = psi4 - gyy(i,j,k) = psi4 + gxx(i,j,k) = psi4 + gyy(i,j,k) = psi4 + gzz(i,j,k) = psi4 gxy(i,j,k) = zero end if @@ -88,12 +97,12 @@ C it should be, or the data will be singular. end do end do -C In any case, +c In any case, gxz = zero gyz = zero -C Vanishing extrinsic curvature completes the Cauchy data. +c Vanishing extrinsic curvature completes the Cauchy data. kxx = zero kyy = zero diff --git a/src/phif.F b/src/phif.F index c24f2ca..3d9bf79 100644 --- a/src/phif.F +++ b/src/phif.F @@ -1,17 +1,23 @@ +/*@@ + @file phif.F + @date October 1998 + @author Miguel Alcubierre + @desc + Function to find angle phi from coordinates {x,y}. + @enddesc + @version $Header$ +@@*/ + #include "cctk.h" CCTK_REAL function phif(x,y) -c Author: Miguel Alcubierre, October 1998. -c -c Function to find angle phi from coordinates {x,y}. - implicit none CCTK_REAL x,y - real*8 zero,half,one,two,pi + CCTK_REAL zero,half,one,two,pi -c Define numbers. +c Numbers. zero = 0.0D0 half = 0.5D0 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 diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F index 0f6abfa..7b289b1 100644 --- a/src/setupbrilldata3D.F +++ b/src/setupbrilldata3D.F @@ -1,3 +1,12 @@ +/*@@ + @file setupbrilldata3D.F + @date October 1998 + @author Miguel Alcubierre + @desc + Set up non-axisymmetric Brill data for elliptic solve. + @enddesc + @version $Header$ +@@*/ #include "cctk.h" #include "cctk_Parameters.h" @@ -7,13 +16,11 @@ subroutine setupbrilldata3D(CCTK_FARGUMENTS) -c Author: Miguel Alcubierre, October 1998. -c c Set up 3D Brill data for elliptic solve. The elliptic c equation we need to solve is: c c __ -c \/ psi = psi R / 8 +c \/ psi - psi R / 8 = 0 c c c c c where: @@ -26,9 +33,9 @@ c c c c The Ricci scalar for the conformal metric turns out to be: c -c / -2q 2 2 -2 2 2 \ -c R = 2 | e ( d q + d q ) + rho ( 3 (d q) + 2 d q ) | -c c \ z rho phi phi / +c / -2q 2 2 -2 2 2 \ +c R = - 2 | e ( d q + d q ) + rho ( 3 (d q) + 2 d q ) | +c c \ z rho phi phi / implicit none @@ -39,6 +46,9 @@ c c \ z rho phi phi / integer CCTK_Equals integer i,j,k integer nx,ny,nz + integer ierr + + integer, dimension(3) :: sym CCTK_REAL x1,y1,z1,rho1,rho2 CCTK_REAL phi,phif,e2q @@ -105,7 +115,19 @@ c Set up conformal metric. gxz = zero gyz = zero -c Set up coefficient Mlinear = - (1/8) Rc. +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 @@ -126,7 +148,7 @@ c Find M using centered differences over a small c interval. if (rho1.gt.rhofudge) then - brillMlinear(i,j,k) = - 0.25D0/e2q + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) . + brillq(rho1+eps,z1,phi) @@ -134,7 +156,7 @@ c interval. . - 4.0D0*brillq(rho1,z1,phi)) . / eps**2 else - brillMlinear(i,j,k) = - 0.25D0/e2q + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) . + two*brillq(eps,z1,phi) @@ -148,7 +170,7 @@ c must always be true otherwise the function c is not regular on the axis. if (rho1.gt.rhofudge) then - brillMlinear(i,j,k) = brillMlinear(i,j,k) - 0.25D0/rho2 + brillMlinear(i,j,k) = brillMlinear(i,j,k) + 0.25D0/rho2 . *(three*0.25D0*(brillq(rho1,z1,phi+eps) . - brillq(rho1,z1,phi-eps))**2 . + two*(brillq(rho1,z1,phi+eps) @@ -164,6 +186,11 @@ 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 -- cgit v1.2.3