diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/ParamCheck.c | 31 | ||||
-rw-r--r-- | src/Startup.c | 4 | ||||
-rw-r--r-- | src/brilldata.F | 30 | ||||
-rw-r--r-- | src/brillq.F | 69 | ||||
-rw-r--r-- | src/finishbrilldata.F | 6 | ||||
-rw-r--r-- | src/setupbrilldata2D.F | 58 | ||||
-rw-r--r-- | src/setupbrilldata3D.F | 118 |
7 files changed, 144 insertions, 172 deletions
diff --git a/src/ParamCheck.c b/src/ParamCheck.c index fe55173..7a04cec 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -52,7 +52,9 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS) if( ! CCTK_Equals(metric_type, "physical") && ! CCTK_Equals(metric_type, "static conformal")) { - CCTK_PARAMWARN("Unknown ADMBase::metric_type - known types are \"physical\" and \"static conformal\""); + CCTK_VParamWarn("Unknown ADMBase::metric_type (%s): " + "known types are \"physical\" and \"static conformal\"", + metric_type); } if (CCTK_Equals(metric_type, "static conformal")) @@ -72,6 +74,33 @@ void IDBrillData_ParamChecker(CCTK_ARGUMENTS) { CCTK_INFO(" ... using physical metric"); } + + if (CCTK_Equals(initial_data,"brilldata")) + { + CCTK_INFO(" ... using full 3D solver"); + } + else if (CCTK_Equals(initial_data,"brilldata2D")) + { + CCTK_INFO(" ... using reduced 2D solver"); + if (brill3d_d != 0.0) + { + CCTK_WARN(0,"But 3D parameter brill3d_d set for the 2D solver !?!"); + } + } + + if (CCTK_Equals(q_function,"exp")) + { + CCTK_INFO(" ... Exponential brill function"); + } + else if (CCTK_Equals(q_function,"eppley")) + { + CCTK_INFO(" ... generalised Eppley Brill function"); + } + else if (CCTK_Equals(q_function,"exp")) + { + CCTK_INFO(" ... Gundlach-Holz Brill function"); + } + } /******************************************************************** diff --git a/src/Startup.c b/src/Startup.c index 8660cbe..a6d4c13 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -18,7 +18,9 @@ void BrillData_InitSymBound(CCTK_ARGUMENTS) sym[1] = 1; sym[2] = 1; - SetCartSymVN(cctkGH, sym,"IDBrillData::brillpsi"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillpsi"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillMlinear"); + SetCartSymVN(cctkGH,sym,"IDBrillData::brillNsource"); return; } diff --git a/src/brilldata.F b/src/brilldata.F index 7dc01e0..3c91b97 100644 --- a/src/brilldata.F +++ b/src/brilldata.F @@ -46,22 +46,33 @@ c Get indices for grid functions. call CCTK_VarIndex(iMcoeff,"idbrilldata::brillMlinear") if (iMcoeff.lt.0) then - call CCTK_WARN(0,"Grid variable index for Mcoeff not found") + call CCTK_WARN(0,"Grid variable index for brillMlinear not found") end if call CCTK_VarIndex(iNcoeff,"idbrilldata::brillNsource") if (iNcoeff.lt.0) then - call CCTK_WARN(0,"Grid variable index for Ncoeff not found") + call CCTK_WARN(0,"Grid variable index for brillNsource not found") end if c Set up background metric and coefficients for linear solve. - if (axisym.eq.1) then + if (CCTK_EQUALS(initial_data,"brilldata2D")) then call setupbrilldata2D(CCTK_ARGUMENTS) else call setupbrilldata3D(CCTK_ARGUMENTS) end if + if (output_coeffs .eq. 1) then + call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillMlinear","Debug_brillMlinear") + call CCTK_OutputVarAs(ierr,cctkGH,"IDBrillData::brillNsource","Debug_brillNsource") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxx","Debug_gxx") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxy","Debug_gxy") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gxz","Debug_gxz") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyy","Debug_gyy") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gyz","Debug_gyz") + call CCTK_OutputVarAs(ierr,cctkGH,"ADMBase::gzz","Debug_gzz") + end if + c Tolerances for elliptic solve. AbsTol(1)= brill_thresh @@ -74,14 +85,17 @@ c Tolerances for elliptic solve. c Boundaries. - call Ell_SetStrKey(ierr,"yes","EllLinConfMetric::Bnd::Robin") - call Ell_SetRealKey(ierr,1.0D0,"EllLinConfMetric::Bnd::Robin::inf") - call Ell_SetIntKey(ierr,1,"EllLinConfMetric::Bnd::Robin::falloff") + call Ell_SetStrKey(ierr,"yes","EllLinMetric::Bnd::Robin") + call Ell_SetRealKey(ierr,1.0D0,"EllLinMetric::Bnd::Robin::inf") + call Ell_SetIntKey(ierr,1,"EllLinMetric::Bnd::Robin::falloff") + +c Elliptic solver -c Elliptic solver. +c Just in case some solvers do not use the standard interface + conformal_state = 0 if (CCTK_EQUALS(brill_solver,"sor")) then - call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit"); + call Ell_SetIntKey(ierr,sor_maxit,"Ell::SORmaxit") call Ell_LinMetricSolver(ierr,cctkGH, . metric_index,ipsi,iMcoeff,iNcoeff,AbsTol,RelTol,"sor") else if (CCTK_EQUALS(brill_solver,"petsc")) then diff --git a/src/brillq.F b/src/brillq.F index bbd9220..df6a65d 100644 --- a/src/brillq.F +++ b/src/brillq.F @@ -50,93 +50,86 @@ c q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ] implicit none DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS logical firstcall integer qtype - integer b,c CCTK_REAL brillq,rho,z,phi - CCTK_REAL a,r0,sr,rho0,srho,sz - CCTK_REAL d,e,m,n,phi0 data firstcall /.true./ save firstcall,qtype - save a,b,c,r0,sr,rho0,srho,sz - save d,e,m,n,phi0 c Get parameters at first call. if (firstcall) then - qtype = brill_q - - a = brill_a - b = brill_b - c = brill_c - - r0 = brill_r0 - sr = brill_sr - rho0 = brill_rho0 - srho = brill_srho - sz = brill_sz - - if (axisym.eq.0) then - d = brill_d - e = brill_e - m = brill_m - n = brill_n - phi0 = brill_phi0 + if (CCTK_EQUALS(q_function,"exp")) then + qtype = 0 + else if (CCTK_EQUALS(q_function,"eppley")) then + qtype = 1 + else if (CCTK_EQUALS(q_function,"gundlach")) then + qtype = 2 + else + call CCTK_WARN(0,"Brill wave data type not recognised") end if firstcall = .false. end if + if (rho < 0) then + rho = -rho + end if + c Calculate q(rho,z) from a choice of forms. -c brill_q = 0. +c "exp" brill data if (qtype.eq.0) then - brillq = a*dabs(rho)**(2 + b) - . / dexp((rho - rho0)**2)/(rho**2 + z**2) + brillq = exp_a*dabs(rho)**(2 + exp_b) + . / dexp((rho - exp_rho0)**2)/(rho**2 + z**2) - if (sz.ne.0.0D0) then - brillq = brillq/dexp((z/sz)**2) + if (exp_sigmaz.ne.0.0D0) then + brillq = brillq/dexp((z/exp_sigmaz)**2) end if else if (qtype.eq.1) then -c brill_q = 1. This includes Eppleys choice of q. +c "Eppley" brill data. This includes Eppleys choice of q. c But note that q(Eppley) = 2q(Cactus). - brillq = a*(rho/srho)**b - . / ( 1.0D0 + ((rho**2 + z**2 - r0**2)/sr**2)**(c/2) ) + brillq = eppley_a*(rho/eppley_sigmarho)**eppley_b + & / ( 1.0D0 + ((rho**2 + z**2 - eppley_r0**2)/ + & eppley_sigmar**2)**(eppley_c/2) ) else if (qtype.eq.2) then -c brill_q = 2. This includes my (Carstens) notion of what a +c "Gundlach" brill data This includes my (Carstens) notion of what a c smooth "pure quadrupole" q should be, plus the choice of c Holz et al. - brillq = a*(rho/srho)**b - . / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2)) + brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b + & / dexp(((rho**2 + z**2 - gundlach_r0**2)/ + & gundlach_sigmar**2)**(gundlach_c/2)) else c Unknown type for q function. - call CCTK_WARN(0,"Unknown type of Brill function q") + call CCTK_WARN(0,"brillq: Unknown type of Brill function q") end if c If desired, multiply with a phi dependent factor. - if (axisym.eq.0) then - brillq = brillq*(1.0D0 + d*rho**m*cos(n*(phi-phi0))**2 - . / (1.0D0 + e*rho**m)) + if (CCTK_EQUALS(initial_data,"brilldata")) then + brillq = brillq*(1.0D0 + brill3d_d*rho**brill3d_m* + & cos(brill3D_n*(phi-brill3d_phi0))**2 + & / (1.0D0 + brill3d_e*rho**brill3d_m)) end if return diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F index 9af9fbf..467c28d 100644 --- a/src/finishbrilldata.F +++ b/src/finishbrilldata.F @@ -55,7 +55,7 @@ c Brill metric calculated from q and Psi. y1 = y(i,j,k) z1 = z(i,j,k) - rho2 = x1**2 + y1**2 + rho2 = x1*x1 + y1*y1 rho1 = dsqrt(rho2) phi = phif(x1,y1) @@ -74,8 +74,8 @@ c The individual coefficients can be read off as if (rho1.gt.rhofudge) then - gxx(i,j,k) = (e2q + (one - e2q)*y1**2/rho2) - gyy(i,j,k) = (e2q + (one - e2q)*x1**2/rho2) + gxx(i,j,k) = (e2q + (one - e2q)*y1*y1/rho2) + gyy(i,j,k) = (e2q + (one - e2q)*x1*x1/rho2) gzz(i,j,k) = e2q gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F index f84e825..f2e5f13 100644 --- a/src/setupbrilldata2D.F +++ b/src/setupbrilldata2D.F @@ -33,11 +33,8 @@ c f DECLARE_CCTK_PARAMETERS integer i,j,k - integer nx,ny,nz integer ierr - integer, dimension(3) :: sym - CCTK_REAL x1,y1,z1,rho1 CCTK_REAL brillq,eps CCTK_REAL zp,zm,rhop,rhom @@ -50,53 +47,34 @@ c Numbers. zero = 0.0D0 one = 1.0D0 -c Set up grid size. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - c Epsilon for finite differencing. eps = cctk_delta_space(1) -c Initialize psi. - - brillpsi = one - -c Initialize metric. - - gxx = one - gyy = one - gzz = one - - gxy = zero - gxz = zero - gyz = zero - -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 - do i=1,nx + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) x1 = x(i,j,k) y1 = y(i,j,k) z1 = z(i,j,k) - rho1 = dsqrt(x1**2 + y1**2) + rho1 = dsqrt(x1*x1 + y1*y1) + + brillpsi(i,j,k) = one + + gxx(i,j,k) = one + gyy(i,j,k) = one + gzz(i,j,k) = one + + gxy(i,j,k) = zero + gxz(i,j,k) = zero + gyz(i,j,k) = zero 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 @@ -116,17 +94,15 @@ c will not be regular on the axis. . + brillq(rhom,z1,zero) . - 4.0D0*brillq(rho1,z1,zero))/eps**2 + brillNsource(i,j,k) = zero + end do end do end do -c Set coefficient Nsource = 0. - - brillNsource = zero - c Synchronization and boundaries. - call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') return diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F index ecd0fd3..97cf28b 100644 --- a/src/setupbrilldata3D.F +++ b/src/setupbrilldata3D.F @@ -43,22 +43,13 @@ c c \ z rho phi phi / DECLARE_CCTK_FUNCTIONS 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 CCTK_REAL brillq,rhofudge,eps CCTK_REAL zero,one,two,three -c Set up grid size. - - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - c Define numbers zero = 0.0D0 @@ -74,78 +65,44 @@ c Epsilon for finite differencing. eps = cctk_delta_space(1) -c Initialize psi. - - brillpsi = one - - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) x1 = x(i,j,k) y1 = y(i,j,k) z1 = z(i,j,k) - rho2 = x1**2 + y1**2 + rho2 = x1*x1 + y1*y1 rho1 = dsqrt(rho2) - phi = phif(x1,y1) e2q = dexp(two*brillq(rho1,z1,phi)) - if (rho1.gt.rhofudge) then - gxx(i,j,k) = e2q + (one - e2q)*y1**2/rho2 - gyy(i,j,k) = e2q + (one - e2q)*x1**2/rho2 - gzz(i,j,k) = e2q - gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 - else - gxx(i,j,k) = one - gyy(i,j,k) = one - gzz(i,j,k) = one - gxy(i,j,k) = zero - end if - - end do - end do - end do - - gxz = zero - gyz = zero - -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 - do i=1,nx - - x1 = x(i,j,k) - y1 = y(i,j,k) - z1 = z(i,j,k) - - rho2 = x1**2 + y1**2 - rho1 = dsqrt(rho2) +c Initialise psi - phi = phif(x1,y1) + brillpsi(i,j,k) = one - e2q = dexp(two*brillq(rho1,z1,phi)) +c Set up metric and coefficient arrays for elliptic solve. +c Notice that the Cactus conventions are: +c __ +c \/ psi + Mlinear*psi + Nsource = 0 c Find M using centered differences over a small c interval. +c Here we assume that for very small rho, the +c phi derivatives are essentially zero. This +c must always be true otherwise the function +c is not regular on the axis. + if (rho1.gt.rhofudge) then + gxx(i,j,k) = e2q + (one - e2q)*y1*y1/rho2 + gyy(i,j,k) = e2q + (one - e2q)*x1*x1/rho2 + gzz(i,j,k) = e2q + gxy(i,j,k) = - (one - e2q)*x1*y1/rho2 + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) @@ -154,8 +111,20 @@ c interval. . - 4.0D0*brillq(rho1,z1,phi)) . / eps**2 + 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) + . - two*brillq(rho1,z1,phi) + . + brillq(rho1,z1,phi-eps)))/eps**2 + else + gxx(i,j,k) = one + gyy(i,j,k) = one + gzz(i,j,k) = one + gxy(i,j,k) = zero + brillMlinear(i,j,k) = 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) @@ -165,32 +134,21 @@ c interval. end if -c Here I assume that for very small rho, the -c phi derivatives are essentially zero. This -c must always be true otherwise the function -c is not regular on the axis. + gxz(i,j,k) = zero + gyz(i,j,k) = zero - if (rho1.gt.rhofudge) then - 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) - . - two*brillq(rho1,z1,phi) - . + brillq(rho1,z1,phi-eps)))/eps**2 - end if +c Set coefficient Nsource = 0 + brillNsource(i,j,k) = zero end do end do end do -c Set coefficient Nsource = 0. - - brillNsource = zero c Synchronization and boundaries. - call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') - call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CCTK_SyncGroup(ierr,cctkGH,'idbrilldata::brillelliptic') +c call CartSymGN(ierr,cctkGH,'idbrilldata::brillelliptic') return end |