diff options
| -rw-r--r-- | param.ccl | 38 | ||||
| -rw-r--r-- | src/AHFinder.F | 46 | ||||
| -rw-r--r-- | src/AHFinder_gau.F | 31 | ||||
| -rw-r--r-- | src/AHFinder_int.F | 36 |
4 files changed, 121 insertions, 30 deletions
@@ -81,6 +81,36 @@ REAL dhole2_zmax "Upper z bound of excised region 2" : :: "Any value" } 0.0 +REAL dhole3_xmin "Lower x bound of excised region 2" +{ +: :: "Any value" +} 0.0 + +REAL dhole3_ymin "Lower y bound of excised region 2" +{ +: :: "Any value" +} 0.0 + +REAL dhole3_zmin "Lower z bound of excised region 2" +{ +: :: "Any value" +} 0.0 + +REAL dhole3_xmax "Upper x bound of excised region 2" +{ +: :: "Any value" +} 0.0 + +REAL dhole3_ymax "Upper y bound of excised region 2" +{ +: :: "Any value" +} 0.0 + +REAL dhole3_zmax "Upper z bound of excised region 2" +{ +: :: "Any value" +} 0.0 + ############################## ### PRIVATE PARAMETERS ### @@ -479,9 +509,9 @@ REAL ahf_shiftcoeff "Coefficient for shift" shares: IO -EXTENDS STRING outdir "Name of IO output directory" +STRING outdir "" { -} "." +} ### FROM GRID ### @@ -492,6 +522,10 @@ KEYWORD domain "" { } +KEYWORD bitant_plane "" +{ +} + ### FROM EINSTEIN ### diff --git a/src/AHFinder.F b/src/AHFinder.F index 178ab05..e6ceb06 100644 --- a/src/AHFinder.F +++ b/src/AHFinder.F @@ -499,6 +499,14 @@ else if (CCTK_Equals(domain,"quadrant").eq.1) then if (xc.eq.zero) refx = .true. if (yc.eq.zero) refy = .true. + else if (CCTK_Equals(domain,"bitant").eq.1) then + if (CCTK_Equals(bitant_plane,"xy").eq.1) then + if (zc.eq.zero) refz = .true. + else if (CCTK_Equals(bitant_plane,"xz").eq.1) then + if (yc.eq.zero) refy = .true. + else + if (xc.eq.zero) refx = .true. + end if end if ! Find {stepx,stepy,stepz}. @@ -678,6 +686,36 @@ rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc-zmn) end if +! Bitant. + + else if (CCTK_Equals(domain,"bitant").eq.1) then + + if (CCTK_Equals(bitant_plane,"xy").eq.1) then + + if (zc.eq.zero) then + rmx = min(xmx-xc,ymx-yc,xc-xmn,yc-ymn,zmx) + else + rmx = min(xmx-xc,ymx-yc,xc-xmn,yc-ymn,zmx-zc,zc) + end if + + else if (CCTK_Equals(bitant_plane,"xz").eq.1) then + + if (yc.eq.zero) then + rmx = min(xmx-xc,zmx-zc,xc-xmn,zc-zmn,ymx) + else + rmx = min(xmx-xc,zmx-zc,xc-xmn,zc-zmn,ymx-yc,yc) + end if + + else + + if (xc.eq.zero) then + rmx = min(ymx-yc,zmx-zc,yc-ymn,zc-zmn,xmx) + else + rmx = min(ymx-yc,zmx-zc,yc-ymn,zc-zmn,xmx-xc,xc) + end if + + end if + ! Full. else @@ -1536,7 +1574,13 @@ ! Third hole. -! This is not here since the parameters do not exist yet! + if (dhole3_xmin.gt.(xc-0.5D0*rhor)) dhole3_xmin = xc-0.5D0*rhor + if (dhole3_ymin.gt.(yc-0.5D0*rhor)) dhole3_ymin = yc-0.5D0*rhor + if (dhole3_zmin.gt.(zc-0.5D0*rhor)) dhole3_zmin = zc-0.5D0*rhor + + if (dhole3_xmax.lt.(xc+0.5D0*rhor)) dhole3_xmax = xc+0.5D0*rhor + if (dhole3_ymax.lt.(yc+0.5D0*rhor)) dhole3_ymax = yc+0.5D0*rhor + if (dhole3_zmax.lt.(zc+0.5D0*rhor)) dhole3_zmax = zc+0.5D0*rhor end if diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F index efe359c..521c54c 100644 --- a/src/AHFinder_gau.F +++ b/src/AHFinder_gau.F @@ -28,21 +28,16 @@ logical firstgau logical firstcal(4) - - integer gxx_index,gyy_index,gzz_index,gxy_index,gxz_index,gyz_index - integer ahfgauss_index - integer x_index,y_index,z_index integer handle,sum_handle - integer j,k,l,m,n,p - CCTK_INT i + integer i,j,k,l,m,n,p integer npoints integer error1,error2,ierror CCTK_REAL LEGEN CCTK_REAL theta,phi,xp,yp,zp,rp CCTK_REAL cost,sint,cosp,sinp - CCTK_REAL dtheta,dphi,dtp,idtheta,idphi + CCTK_REAL dtheta,dphi,dtp,idtheta,idphi,phistart CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx CCTK_REAL xmn_l,ymn_l,zmn_l,xmx_l,ymx_l,zmx_l CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp @@ -80,6 +75,7 @@ data firstcal(2) / .true. / data firstcal(3) / .true. / data firstcal(4) / .true. / + save firstcal save firstgau @@ -105,6 +101,8 @@ ! dphi Grid spacing in phi. ! dtp dtheta*dphi. ! +! phistart Origin for phi (normally 0, but for some symmetries -pi/2). +! ! idtheta 1/(2 dtheta) ! idphi 1/(2 dphi) ! @@ -284,7 +282,9 @@ call CCTK_CoordRange(ierror,cctkGH,ymn,ymx,-1,"y","cart3d"); call CCTK_CoordRange(ierror,cctkGH,zmn,zmx,-1,"z","cart3d"); -! Find {dtheta,dphi,dtp}. +! Find {phistart,dtheta,dphi,dtp}. + + phistart = zero dtheta = pi/dble(ntheta) @@ -305,7 +305,7 @@ dphi = 0.25D0*dphi else if (refx) then dphi = half*dphi - call CCTK_INFO('This combination of symmetries has not been properly implemented yet!') + phistart = - half*pi else if (refy) then dphi = half*dphi end if @@ -509,7 +509,7 @@ ! Find {theta,phi}. theta = dtheta*dble(i) - phi = dphi*dble(j) + phi = dphi*dble(j) + phistart ! Find sines and cosines. @@ -637,7 +637,7 @@ ! Interpolator handle. - call CCTK_InterpHandle (handle,"simple") + call CCTK_InterpHandle(handle,"simple") ! Local origin of spatial coordinates. @@ -648,9 +648,10 @@ ! Interpolation. call CCTK_Interp(ierror,cctkGH,handle,npoints,3,7,7, - . nx,ny,nz,xa,ya,za,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, - . CCTK_VARIABLE_REAL,xmn_l,ymn_l,zmn_l, - . dx,dy,dz,gxx,gyy,gzz,gxy,gxz,gyz,ahfgauss, + . nx,ny,nz,xa,ya,za, + . CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + . xmn_l,ymn_l,zmn_l,dx,dy,dz, + . gxx,gyy,gzz,gxy,gxz,gyz,ahfgauss, . CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, . CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, . CCTK_VARIABLE_REAL, @@ -672,7 +673,7 @@ ! Find {theta,phi}. theta = dtheta*dble(i) - phi = dphi*dble(j) + phi = dphi*dble(j) + phistart ! Find {rp}. diff --git a/src/AHFinder_int.F b/src/AHFinder_int.F index 81a952e..a28c198 100644 --- a/src/AHFinder_int.F +++ b/src/AHFinder_int.F @@ -42,15 +42,17 @@ CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx CCTK_REAL xmn_l,ymn_l,zmn_l,xmx_l,ymx_l,zmx_l - CCTK_REAL dtheta,dphi,dtp,idtheta,idphi + CCTK_REAL dtheta,dphi,dtp,idtheta,idphi,phistart CCTK_REAL det,idet CCTK_REAL zero,quarter,half,one,two,three,four,pi CCTK_REAL sina,cosa,intw,grad,sigma,h0 CCTK_REAL aux,aux1,aux2 - CCTK_REAL gupij(3,3) CCTK_REAL tempv(0:lmax),tempm(lmax,lmax) - CCTK_REAL, dimension(3) :: origin,maxval,minval + + CCTK_REAL, dimension(3) :: origin,maxval,minval + CCTK_REAL, dimension(3,3) :: gupij + CCTK_REAL, allocatable, dimension(:,:) :: theta,phi CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za,da,exp,gradn CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz @@ -64,7 +66,7 @@ ! Variables to be saved for next call. save xmn,ymn,zmn,xmx,ymx,zmx - save dtheta,dphi,dtp,idtheta,idphi + save dtheta,dphi,dtp,idtheta,idphi,phistart save npt,npp,l_ntheta,l_nphi,theta0,phi0 ! Description of variables: @@ -140,6 +142,8 @@ ! idtheta 1/(2 dtheta) ! idphi 1/(2 dphi) ! +! phistart Origin for phi (normally 0, but for some symmetries -pi/2). +! ! da Area element. @@ -157,6 +161,11 @@ pi = acos(-one) + +! ***************************************************************** +! *** TOTAL NUMBER OF PROCESSORS AND LOCAL PROCESSOR NUMBER *** +! ***************************************************************** + myproc = CCTK_MyProc(cctkGH) nprocs = CCTK_nProcs(cctkGH) @@ -190,9 +199,11 @@ call CCTK_CoordRange(ierror,cctkGH,zmn,zmx,-1,"z","cart3d") -! ********************************* -! *** FIND dtheta,dphi,dtp *** -! ********************************* +! ****************************************** +! *** FIND phistart,dtheta,dphi,dtp *** +! ****************************************** + + phistart = zero dtheta = pi/dble(ntheta) @@ -212,7 +223,8 @@ if (refx.and.refy) then dphi = quarter*dphi else if (refx) then - call CCTK_INFO('This combination of symmetries has not been properly implemented yet!') + dphi = half*dphi + phistart = - half*pi else if (refy) then dphi = half*dphi end if @@ -344,8 +356,8 @@ i = mod(myproc,npt) j = myproc/npt - theta0 = i*l_ntheta - phi0 = j*l_nphi + theta0 = i*l_ntheta + phi0 = j*l_nphi ! Now correct all this if the total number of grid points ! in a given direction was not exactly divisible by the @@ -477,7 +489,7 @@ ! Find {theta,phi}. theta(i,j) = dtheta*dble(i+theta0) - phi(i,j) = dphi*dble(j+phi0) + phi(i,j) = dphi*dble(j+phi0) + phistart ! Find sines and cosines. @@ -644,7 +656,7 @@ . CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, . CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL) -! For H of C flows, we need the interpolated expansion, the +! For H or C flows, we need the interpolated expansion, the ! norm of the gradient of the horizon function, and the mask. else |