[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 149)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2012/08/09 01:26 AM

Modified:
 /trunk/
  param.ccl
 /trunk/src/
  GRHydro_MonopoleM.F90

Log:
 Slight tweaks to the monopole routine options (patch by Josh)
 
 From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_MonopoleM.F90
Delta lines: +26 -5
===================================================================
--- trunk/src/GRHydro_MonopoleM.F90	2012-08-09 06:26:55 UTC (rev 148)
+++ trunk/src/GRHydro_MonopoleM.F90	2012-08-09 06:26:57 UTC (rev 149)
@@ -71,9 +71,14 @@
            vely(i,j,k) = 0.0
            velz(i,j,k) = 0.0
            eps(i,j,k) = 0.1
+
+           Bvecx(i,j,k)=0.0
+           Bvecy(i,j,k)=0.0
+           Bvecz(i,j,k)=0.0
+
            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
+                 Bvecx(i,j,k)=Monopole_point_Bx
               else
                  Bvecx(i,j,k)=0.0
               endif
@@ -93,7 +98,7 @@
               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)
+              if(mod(i+j+k,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
@@ -104,17 +109,33 @@
                  Bvecx(i,j,k) = 0.0
                  Bvecy(i,j,k) = 0.0
               endif
+              if(mod(i+j+k,2).eq.0)then
+                 Bvecx(i,j,k)=-1.0*Bvecx(i,j,k)
+!!$  Only vary one component, for different character
+!!$                 Bvecy(i,j,k)=-1.0*Bvecy(i,j,k)
+              endif
+           else if(CCTK_EQUALS(monopole_type,"3dalt")) 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)
+                 Bvecz(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
+                 Bvecz(i,j,k) = 0.0
+              endif
+!!$  Different spatial pattern!
               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)
+                 Bvecz(i,j,k)=-1.0*Bvecz(i,j,k)
               endif
            else
               call CCTK_WARN(0,"Unrecognized monopole type!!!")
            endif
 
-           Bvecy(i,j,k)=0.0
-           Bvecz(i,j,k)=0.0
-
            det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
 
         if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then

Directory: /trunk/
==================

File [modified]: param.ccl
Delta lines: +6 -0
===================================================================
--- trunk/param.ccl	2012-08-09 06:26:55 UTC (rev 148)
+++ trunk/param.ccl	2012-08-09 06:26:57 UTC (rev 149)
@@ -147,6 +147,7 @@
   "Gauss" :: "Gaussian w/radius R_Gauss"
   "1dalt" :: "1-d alternating"
   "2dalt" :: "2-d alternating"
+  "3dalt" :: "3-d alternating"
 } "Point"
 
 CCTK_REAL R_Gauss "Radius for a Gaussian monopole"
@@ -154,6 +155,11 @@
   0:* :: "Any positive number"
 } 1.0
 
+CCTK_REAL Monopole_point_Bx "Pointlike Monopole Bx value"
+{
+  *:* :: "Any number"
+} 1.0
+
 # For cylindrical explosion test:
 CCTK_REAL cyl_r_inner "Inner Radius"
 {