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

[rhaas at tapir.caltech.edu]

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

Modified:
 /trunk/src/
  GRHydro_AdvectedLoopM.F90

Log:
 Revised parameters for advected loop test
 
 From: Josh Faber <jfaberrit at jafsmamacbook.main.ad.rit.edu>

File Changes:

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

File [modified]: GRHydro_AdvectedLoopM.F90
Delta lines: +60 -12
===================================================================
--- trunk/src/GRHydro_AdvectedLoopM.F90	2012-08-09 06:26:42 UTC (rev 141)
+++ trunk/src/GRHydro_AdvectedLoopM.F90	2012-08-09 06:26:44 UTC (rev 142)
@@ -67,7 +67,7 @@
   CCTK_REAL :: range_x,range_y,range_z,range_d
   CCTK_REAL :: cos_theta, sin_theta, tan_theta
   CCTK_REAL :: Bvecx_d, Bvecy_d, Bvecz_d
-  CCTK_REAL :: x_d, y_d, z_d
+  CCTK_REAL :: x_d, y_d, z_d,diaglength
   CCTK_REAL :: dx_d, dy_d, dz_d, dx_x, dz_x
 
 !!$Adiabatic index for test:
@@ -83,19 +83,47 @@
   rhoval   = 1.0d0
   pressval = 3.0d0
 
+  if (CCTK_EQUALS(advectedloop_type,"2D")) then	
+
+
 !!$ Vx, Vy and Vz values:
-  if (CCTK_EQUALS(advectedloop_case,"V^z/=0")) then
-   vxval=0.2d0/sqrt(6.0d0)
-   vyval=0.5d0*vxval
-   vzval=vyval
-  else if (CCTK_EQUALS(advectedloop_case,"V^z=0")) then
-   vxval=0.2d0/sqrt(6.0d0)
-   vyval=0.5d0*vxval
-   vzval=0.0d0
-  else
-   call CCTK_WARN(0,"V^z component case not recognized!")
-  end if 
+    if (CCTK_EQUALS(advectedloop_case,"V^z/=0")) then
+
+!!$vxval=0.2d0/sqrt(6.0d0)
+!!$ This new choice yields a crossing time of exactly t=24
+!!$ assuming -1<x<1; -0.5<y<0.5
+     vxval=1.d0/12.d0
+     vyval=0.5d0*vxval
+     vzval=vyval
+
+    else if (CCTK_EQUALS(advectedloop_case,"V^z=0")) then
+!!$   vxval=0.2d0/sqrt(6.0d0)
+     vxval=1.d0/12.d0
+     vyval=0.5d0*vxval
+     vzval=0.0d0
+    else
+     call CCTK_WARN(0,"V^z component case not recognized!")
+    end if 
     
+  else if (CCTK_EQUALS(advectedloop_type,"3D")) then	
+
+     vxval=0.2d0*sqrt(2.0d0)
+     vyval=0.2d0
+
+    if (CCTK_EQUALS(advectedloop_case,"V^z/=0")) then
+
+!!$vxval=0.2d0/sqrt(6.0d0)
+!!$ This new choice yields a crossing time of exactly t=5
+!!$ assuming -1<x<1; -0.5<y<0.5
+     vzval=0.1
+
+   else if (CCTK_EQUALS(advectedloop_case,"V^z=0")) then
+     vzval=0.0d0
+    else
+     call CCTK_WARN(0,"V^z component case not recognized!")
+    end if 
+
+  endif
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
   nz = cctk_lsh(3)
@@ -110,7 +138,10 @@
   range_z = (cctk_gsh(3)-1)*dz
 
   range_d = sqrt(range_z**2+range_x**2)
+  diaglength = range_x*range_z/range_d
 
+!!$ For 3-d case, the grid is going to be assumed to be a cube
+
   cos_theta = range_z/range_d
   sin_theta = range_x/range_d
   tan_theta = sin_theta/cos_theta
@@ -162,16 +193,33 @@
 
            else if (CCTK_EQUALS(advectedloop_type,"3D")) then
 
+
+!!$  tangential velocity should be parallel to (1,1,1) plus a normal component (-1,0,1)
+
              velx(i,j,k)=cos_theta*vxval-sin_theta*vzval
              vely(i,j,k)=vyval
              velz(i,j,k)=cos_theta*vzval+sin_theta*vxval
 
              Bvecz_d=0.0d0
 
+!!$ x_d = (x+z)/sqrt(2) => x_d=0 is equivalent to x+z=0
+
              x_d = cos_theta*x(i,j,k)+sin_theta*z(i,j,k)
              y_d = y(i,j,k)
              z_d = cos_theta*z(i,j,k)-sin_theta*x(i,j,k)
           
+!!$ need to make x_d periodic!
+
+       	     if(x_d.gt.1.5*diaglength) then
+	        x_d=x_d-2.0*diaglength
+             else if (x_d.gt.0.5*diaglength .and. x_d.lt.1.5*diaglength) then
+	        x_d=x_d-diaglength
+	     else if(x_d.lt.-1.5*diaglength) then
+	        x_d=x_d+2.0*diaglength
+             else if (x_d.lt.(-0.5*diaglength) .and. x_d.gt.(-1.5*diaglength)) then
+	        x_d=x_d+diaglength
+             endif
+	
              radius = sqrt(x_d**2+y_d**2)
   
              if (CCTK_EQUALS(advectedloop_delA,"Exact")) then