diff options
Diffstat (limited to 'src/GRHydro_MonopoleM.F90')
| -rw-r--r-- | src/GRHydro_MonopoleM.F90 | 44 |
1 files changed, 40 insertions, 4 deletions
diff --git a/src/GRHydro_MonopoleM.F90 b/src/GRHydro_MonopoleM.F90 index 6a6ea32..b8bee8c 100644 --- a/src/GRHydro_MonopoleM.F90 +++ b/src/GRHydro_MonopoleM.F90 @@ -56,7 +56,7 @@ subroutine GRHydro_MonopoleM(CCTK_ARGUMENTS) velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3 - + CCTK_REAL :: rr2,rg2 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -71,11 +71,47 @@ subroutine GRHydro_MonopoleM(CCTK_ARGUMENTS) vely(i,j,k) = 0.0 velz(i,j,k) = 0.0 eps(i,j,k) = 0.1 - if(i.eq.nx/2+1.and.j.eq.ny/2+1.and.k.eq.nz/2+1) then - Bvecx(i,j,k)=0.1 + if(CCTK_EQUALS(monopole_type,"Point")) then + if(i.eq.nx/2+1.and.j.eq.ny/2+1.and.k.eq.nz/2+1) then + Bvecx(i,j,k)=0.1 + else + Bvecx(i,j,k)=0.0 + endif + else if(CCTK_EQUALS(monopole_type,"Gauss")) then + rr2=x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2 + rg2=R_Gauss*R_Gauss + if(rr2.lt.rg2) then + Bvecx(i,j,k) = exp(-1.0*rr2/rg2)-1.0/exp(1.0) + else + Bvecx(i,j,k) = 0.0 + endif + else if(CCTK_EQUALS(monopole_type,"1dalt")) then + rr2=x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2 + rg2=R_Gauss*R_Gauss + if(rr2.lt.rg2) then + Bvecx(i,j,k) = exp(-1.0*rr2/rg2)-1.0/exp(1.0) + else + Bvecx(i,j,k) = 0.0 + endif + if(mod(i+j,2).eq.0)Bvecx(i,j,k)=-1.0*Bvecx(i,j,k) + else if(CCTK_EQUALS(monopole_type,"2dalt")) then + rr2=x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2 + rg2=R_Gauss*R_Gauss + if(rr2.lt.rg2) then + Bvecx(i,j,k) = exp(-1.0*rr2/rg2)-1.0/exp(1.0) + Bvecy(i,j,k) = exp(-1.0*rr2/rg2)-1.0/exp(1.0) + else + Bvecx(i,j,k) = 0.0 + Bvecy(i,j,k) = 0.0 + endif + if(mod(i+j,2).eq.0)then + Bvecx(i,j,k)=-1.0*Bvecx(i,j,k) + Bvecy(i,j,k)=-1.0*Bvecy(i,j,k) + endif else - Bvecx(i,j,k)=0.0 + call CCTK_WARN(0,"Unrecognized monopole type!!!") endif + Bvecy(i,j,k)=0.0 Bvecz(i,j,k)=0.0 |