#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" subroutine solve_iter_IVP(CCTK_FARGUMENTS,pass) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS INTEGER CCTK_Equals INTEGER, DIMENSION(3) :: sym integer pass integer i,j,k,bbox_save(6) INTEGER ierr CCTK_REAL one,two,divw,pi,bam_robin_limit_orig CCTK_REAL bam_robin_power_orig integer counter CCTK_REAL vlowx,vlowy,vlowz,ww_lorentz CCTK_REAL dx, dy, dz CCTK_REAL my_elliptic_tol integer octant, bitant, quadrant octant = (CCTK_Equals(domain,"octant")) bitant = (CCTK_Equals(domain,"bitant")) quadrant = (CCTK_Equals(domain,"quadrant")) c if first pass, use smaller tolerence for elliptic solvers if((pass.eq.1).and.(ivp_2pass.ne.0)) then my_elliptic_tol = ivp_elliptic_tol_pass1 else my_elliptic_tol = ivp_elliptic_tol endif one = 1.0d0 two = 2.0d0 pi = 3.141592653589793D0 c NOW, COMPUTE A^TWIDLE^2 c use kxx, etc. as temporary variables for A^Twidle kxx = 0.0d0 kxy = 0.0d0 kxz = 0.0d0 kyy = 0.0d0 kyz = 0.0d0 kzz = 0.0d0 c dx settings dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) do k=2,cctk_lsh(3)-1 do j=2,cctk_lsh(2)-1 do i=2,cctk_lsh(1)-1 divw = ( (sqrt(ivp_det(i+1,j,k))*ivp_wx(i+1,j,k) - . sqrt(ivp_det(i-1,j,k))*ivp_wx(i-1,j,k))/(2.0d0*dx) + . (sqrt(ivp_det(i,j+1,k))*ivp_wy(i,j+1,k) - . sqrt(ivp_det(i,j-1,k))*ivp_wy(i,j-1,k))/(2.0d0*dy) + . (sqrt(ivp_det(i,j,k+1))*ivp_wz(i,j,k+1) - . sqrt(ivp_det(i,j,k-1))*ivp_wz(i,j,k-1))/(2.0d0*dz))/ . sqrt(ivp_det(i,j,k)) kxx(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uxx(i+1,j,k)-ivp_uxx(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uxx(i,j+1,k)-ivp_uxx(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uxx(i,j,k+1)-ivp_uxx(i,j,k-1))/(2.0d0*dz) + . ivp_uxx(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) + . ivp_uxy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) + . ivp_uxz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxx(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) + . ivp_uxy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) + . ivp_uxz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uxx(i,j,k) kyy(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uyy(i+1,j,k)-ivp_uyy(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uyy(i,j+1,k)-ivp_uyy(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uyy(i,j,k+1)-ivp_uyy(i,j,k-1))/(2.0d0*dz) + . ivp_uxy(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) + . ivp_uyy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) + . ivp_uyz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxy(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) + . ivp_uyy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) + . ivp_uyz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uyy(i,j,k) kzz(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uzz(i+1,j,k)-ivp_uzz(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uzz(i,j+1,k)-ivp_uzz(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uzz(i,j,k+1)-ivp_uzz(i,j,k-1))/(2.0d0*dz) + . ivp_uxz(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) + . ivp_uyz(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) + . ivp_uzz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxz(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) + . ivp_uyz(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) + . ivp_uzz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uzz(i,j,k) kxy(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uxy(i+1,j,k)-ivp_uxy(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uxy(i,j+1,k)-ivp_uxy(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uxy(i,j,k+1)-ivp_uxy(i,j,k-1))/(2.0d0*dz) + . ivp_uxx(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) + . ivp_uxy(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) + . ivp_uxz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxy(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) + . ivp_uyy(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) + . ivp_uyz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uxy(i,j,k) kxz(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uxz(i+1,j,k)-ivp_uxz(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uxz(i,j+1,k)-ivp_uxz(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uxz(i,j,k+1)-ivp_uxz(i,j,k-1))/(2.0d0*dz) + . ivp_uxx(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) + . ivp_uxy(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) + . ivp_uxz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxz(i,j,k)* (ivp_wx(i+1,j,k) -ivp_wx(i-1,j,k ))/(2.0d0*dx) + . ivp_uyz(i,j,k)* (ivp_wx(i,j+1,k) -ivp_wx(i,j-1,k ))/(2.0d0*dy) + . ivp_uzz(i,j,k)* (ivp_wx(i,j,k+1) -ivp_wx(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uxz(i,j,k) kyz(i,j,k) = - . ivp_wx(i,j,k) * (ivp_uyz(i+1,j,k)-ivp_uyz(i-1,j,k))/(2.0d0*dx) - . ivp_wy(i,j,k) * (ivp_uyz(i,j+1,k)-ivp_uyz(i,j-1,k))/(2.0d0*dy) - . ivp_wz(i,j,k) * (ivp_uyz(i,j,k+1)-ivp_uyz(i,j,k-1))/(2.0d0*dz) + . ivp_uxy(i,j,k)* (ivp_wz(i+1,j,k) -ivp_wz(i-1,j,k ))/(2.0d0*dx) + . ivp_uyy(i,j,k)* (ivp_wz(i,j+1,k) -ivp_wz(i,j-1,k ))/(2.0d0*dy) + . ivp_uyz(i,j,k)* (ivp_wz(i,j,k+1) -ivp_wz(i,j,k-1 ))/(2.0d0*dz) + . ivp_uxz(i,j,k)* (ivp_wy(i+1,j,k) -ivp_wy(i-1,j,k ))/(2.0d0*dx) + . ivp_uyz(i,j,k)* (ivp_wy(i,j+1,k) -ivp_wy(i,j-1,k ))/(2.0d0*dy) + . ivp_uzz(i,j,k)* (ivp_wy(i,j,k+1) -ivp_wy(i,j,k-1 ))/(2.0d0*dz) - . (2.0d0/3.0d0) * divw * ivp_uyz(i,j,k) c now, add in astar to get A^twidle kxx(i,j,k) = kxx(i,j,k) + ivp_astarxx(i,j,k) kxy(i,j,k) = kxy(i,j,k) + ivp_astarxy(i,j,k) kxz(i,j,k) = kxz(i,j,k) + ivp_astarxz(i,j,k) kyy(i,j,k) = kyy(i,j,k) + ivp_astaryy(i,j,k) kyz(i,j,k) = kyz(i,j,k) + ivp_astaryz(i,j,k) kzz(i,j,k) = kzz(i,j,k) + ivp_astarzz(i,j,k) enddo enddo enddo call CCTK_SyncGroup(cctkGH, "einstein::curv") c should apply inner boundary conditions here bbox_save = cctk_bbox if(octant.ne.0) then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz') sym(1) = -1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy') sym(1) = -1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz') sym(1) = 1 sym(2) = -1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz') endif if(quadrant.ne.0)then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(5) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz') sym(1) = -1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy') sym(1) = -1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz') endif if(bitant.ne.0)then cctk_bbox(1) = 0 cctk_bbox(2) = 0 cctk_bbox(3) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxx') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyy') call SetCartSymVN(ierr,cctkGH,sym,'einstein::kzz') sym(1) = -1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxy') sym(1) = -1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kxz') sym(1) = 1 sym(2) = -1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'einstein::kyz') endif cctk_bbox = bbox_save c add A^twidle^2 term to (phi**-7) source ivp_bam_m7 = & gxx**2*kxx**2+4.d0*gxx*gxy*kxx*kxy+2.d0*gxy**2*kxy**2+2.d0*gxx*g & yy*kxy**2+4.d0*gxx*gxz*kxx*kxz+4.d0*gxy*gxz*kxy*kxz+4.d0*gxx*gyz & *kxy*kxz+2.d0*gxz**2*kxz**2+2.d0*gxx*gzz*kxz**2+2.d0*gxy**2*kxx* & kyy+4.d0*gxy*gyy*kxy*kyy+4.d0*gxy*gyz*kxz*kyy+gyy**2*kyy**2+4.d0 & *gxy*gxz*kxx*kyz+4.d0*gxz*gyy*kxy*kyz+4.d0*gxy*gyz*kxy*kyz+4.d0* & gxz*gyz*kxz*kyz+4.d0*gxy*gzz*kxz*kyz+4.d0*gyy*gyz*kyy*kyz+2.d0*g & yz**2*kyz**2+2.d0*gyy*gzz*kyz**2+2.d0*gxz**2*kxx*kzz+4.d0*gxz*gy & z*kxy*kzz+4.d0*gxz*gzz*kxz*kzz+2.d0*gyz**2*kyy*kzz+4.d0*gyz*gzz* & kyz*kzz+gzz**2*kzz**2 ivp_bam_m7 = ivp_bam_m7 / 8.0d0 c add scalar curvature contributions Mlinear = - ivp_scalarc / 8.0d0 c add trace of extrinsic curvature contributions ivp_bam_p5 = - (trK**2) / 12.0d0 Nsource = 0.0d0 c now set up vector equations ivp_wx_source = 0.0d0 ivp_wy_source = 0.0d0 ivp_wz_source = 0.0d0 do k=2,cctk_lsh(3)-1 do j=2,cctk_lsh(2)-1 do i=2,cctk_lsh(1)-1 ivp_wx_source(i,j,k) = (-2.0d0/3.0d0) * (ivp_psi(i,j,k)**6) * . ( ivp_uxx(i,j,k)*(trK(i+1,j,k)-trK(i-1,j,k))/(2.0d0*dx) + . ivp_uxy(i,j,k)*(trK(i,j+1,k)-trK(i,j-1,k))/(2.0d0*dy) + . ivp_uxz(i,j,k)*(trK(i,j,k+1)-trK(i,j,k-1))/(2.0d0*dz)) ivp_wy_source(i,j,k) = (-2.0d0/3.0d0) * (ivp_psi(i,j,k)**6) * . ( ivp_uxy(i,j,k)*(trK(i+1,j,k)-trK(i-1,j,k))/(2.0d0*dx) + . ivp_uyy(i,j,k)*(trK(i,j+1,k)-trK(i,j-1,k))/(2.0d0*dy) + . ivp_uyz(i,j,k)*(trK(i,j,k+1)-trK(i,j,k-1))/(2.0d0*dz)) ivp_wz_source(i,j,k) = (-2.0d0/3.0d0) * (ivp_psi(i,j,k)**6) * . ( ivp_uxz(i,j,k)*(trK(i+1,j,k)-trK(i-1,j,k))/(2.0d0*dx) + . ivp_uyz(i,j,k)*(trK(i,j+1,k)-trK(i,j-1,k))/(2.0d0*dy) + . ivp_uzz(i,j,k)*(trK(i,j,k+1)-trK(i,j,k-1))/(2.0d0*dz)) enddo enddo enddo if (assume_astar_divfree.eq.0) then call add_div_astar_vec_source(CCTK_FARGUMENTS) endif ivp_bam_m3 = 0.0d0 if (mahc_evolution.ne.0) then ivp_bam_m3 = 2.0d0 * pi * . ( (rho + rho*eps + press) * w_lorentz**2 - press) ivp_wx_source = ivp_wx_source - 8.0d0 * pi * . ( (rho + rho*eps + press) * w_lorentz**2 * velx) * . ( ivp_psi ** (ivp_matterdens_pow + 8.0d0)) ivp_wy_source = ivp_wy_source - 8.0d0 * pi * . ( (rho + rho*eps + press) * w_lorentz**2 * vely) * . ( ivp_psi ** (ivp_matterdens_pow + 8.0d0)) ivp_wz_source = ivp_wz_source - 8.0d0 * pi * . ( (rho + rho*eps + press) * w_lorentz**2 * velz) * . ( ivp_psi ** (ivp_matterdens_pow + 8.0d0)) endif 1000 continue call CCTK_SyncGroup(cctkGH, "ivp::ivp_wi_group") c should apply inner boundary conditions here bbox_save = cctk_bbox if(octant.ne.0) then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::Mlinear') call SetCartSymVN(ierr,cctkGH,sym,'ivp::Nsource') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx_source') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy_source') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz_source') endif if(quadrant.ne.0)then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(5) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::Mlinear') call SetCartSymVN(ierr,cctkGH,sym,'ivp::Nsource') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx_source') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy_source') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz_source') endif if(bitant.ne.0)then cctk_bbox(1) = 0 cctk_bbox(2) = 0 cctk_bbox(3) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::Mlinear') call SetCartSymVN(ierr,cctkGH,sym,'ivp::Nsource') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx_source') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy_source') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz_source') endif cctk_bbox = bbox_save c*********************************************************************** c solve the equations! c*********************************************************************** c do not change when do_ivp = yes c and do_ivp_for_HAM = 'no' if(do_ivp_for_HAM.eq.0) then ivp_psi = 1.0d0 else c set up outer boundary for Hamiltonian constraint ******************** c case 1: when you use BAM Robin BC (BAM_Elliptic) if(CCTK_Equals(ivp_robin_method, "BAM").eq.1) then c we are messing with parameter database here, so be careful c to leave it in its original form! write(*,*)'Something needs to be done with these' write(*,*)'BAM parameters' STOP c bam_robin_limit_orig = bam_robin_limit c bam_robin_power_orig = bam_robin_power c bam_robin_limit = 1.0d0 c bam_robin_power = 1.0d0 c write(*,*) 'IVP: setting robin parameters to ', c . bam_robin_limit, bam_robin_power c case 2: fix elseif(CCTK_Equals(ivp_robin_method,"fix").ne.0) then ivp_psi=1.0 endif c end of if(CCTKEquals ...ivp_robin_method...BAM) c----------------------------------------------------------------------- c solve hamiltonian constraint c----------------------------------------------------------------------- if(CCTK_Equals(ivp_ham_solver,"bam_nonlin").ne.0) then call ivp_bam_nonlin(cctkGH,my_elliptic_tol) elseif(CCTK_Equals(ivp_ham_solver,"linear").ne.0) then call linear_ham_solver(CCTK_FARGUMENTS,my_elliptic_tol) else write(*,*) 'Error in IVP: must choose ham_solver to be ' write(*,*) ' either bam_nonlin or linear!' STOP endif c endif of if (do_ivp_for_HAM) .. else .. endif c set up outer boundary for momentum constraints ********************** c do not solve when time symmetric c do not solve when do_ivp_for_MOM = no if ((ivp_ts.ne.0).or.(do_ivp_for_MOM.eq.0)) then ivp_wx = 0.0d0 ivp_wy = 0.0d0 ivp_wz = 0.0d0 else c set up outer boundary for momentum constraints ********************** if(CCTK_Equals(ivp_robin_method,"BAM").eq.1) then write(*,*)'Something needs to be done with these' write(*,*)'BAM parameters' STOP c bam_robin_limit = 0.0d0 c bam_robin_power = ivp_wvec_falloff_power c write(*,*) 'IVP: call ivp_bam_veclaplace w/robin parameters ', c . bam_robin_limit,bam_robin_power endif if(CCTK_Equals(ivp_robin_method,"fix").ne.0) then ivp_wx = 0.0d0 ivp_wy = 0.0d0 ivp_wz = 0.0d0 endif c----------------------------------------------------------------------- c solve momentum constraints c----------------------------------------------------------------------- call ivp_bam_veclaplace(cctkGH,my_elliptic_tol) endif c end of if (do_ivp_for_MOM) .. else .. c now, set the parameter database to its original state! if(CCTK_Equals(ivp_robin_method,"BAM").eq.1) then write(*,*)'Something needs to be done with these' write(*,*)'BAM parameters' STOP c bam_robin_limit = bam_robin_limit_orig c bam_robin_power = bam_robin_power_orig c write(*,*) 'IVP: setting robin parameters back to', c . bam_robin_limit,bam_robin_power endif c should apply inner boundary conditions here bbox_save = cctk_bbox if(octant.ne.0) then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz') endif if(quadrant.ne.0)then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(5) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz') endif if(bitant.ne.0)then cctk_bbox(1) = 0 cctk_bbox(2) = 0 cctk_bbox(3) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wx') sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wy') sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_wz') endif cctk_bbox = bbox_save return end