User: rhaas Date: 2012/08/09 01:27 AM
Modified: /trunk/ param.ccl /trunk/src/ GRHydro_AdvectedLoopM.F90, GRHydro_AlfvenWaveM.F90
Log: Fix to oblique case for Advected loop, new implementation of Alfven wave test including used-defined pressure (to alter wavespeed) and oblique case
From: Josh Faber jfaberrit@jfaber3.local
File Changes:
Directory: /trunk/src/ ======================
File [modified]: GRHydro_AdvectedLoopM.F90 Delta lines: +3 -3 =================================================================== --- trunk/src/GRHydro_AdvectedLoopM.F90 2012-08-09 06:26:57 UTC (rev 149) +++ trunk/src/GRHydro_AdvectedLoopM.F90 2012-08-09 06:27:00 UTC (rev 150) @@ -133,9 +133,9 @@ dz = CCTK_DELTA_SPACE(3)
!!$ Note that the 3D test wasn't deviced to be used with AMR! - range_x = (cctk_gsh(1)-1)*dx - range_y = (cctk_gsh(2)-1)*dy - range_z = (cctk_gsh(3)-1)*dz + range_x = (cctk_gsh(1)-2*cctk_nghostzones(1))*dx + range_y = (cctk_gsh(2)-2*cctk_nghostzones(2))*dy + range_z = (cctk_gsh(3)-2*cctk_nghostzones(3))*dz
range_d = sqrt(range_z**2+range_x**2) diaglength = range_x*range_z/range_d
File [modified]: GRHydro_AlfvenWaveM.F90 Delta lines: +53 -30 =================================================================== --- trunk/src/GRHydro_AlfvenWaveM.F90 2012-08-09 06:26:57 UTC (rev 149) +++ trunk/src/GRHydro_AlfvenWaveM.F90 2012-08-09 06:27:00 UTC (rev 150) @@ -65,6 +65,7 @@ CCTK_REAL :: dx,dy,dz CCTK_REAL :: range_x,range_y,range_z,range_d CCTK_REAL :: cos_theta, sin_theta + CCTK_REAL :: diaglength,xnew,vparallel,vperp,Bparallel,Bperp CCTK_REAL :: Bvecx_d, Bvecz_d CCTK_REAL :: velx_d, velz_d CCTK_REAL :: det @@ -77,10 +78,13 @@
!!$pressure, density, B^x, specific internal energy and enthalpy: rhoval = 1.0d0 - pressval = 1.0d0 + pressval = alfvenwave_pressure Bxval = 1.0d0 epsval = pressval/(gam-1.0d0)/rhoval hval = 1.0d0 + epsval + pressval/rhoval + +!!$ DZ: rho=P=Bxval=AA = 1 +!!$ Using DZ parameters, epsval=1.5, hval=3.5
!!$ Alfven Wave Amplitude: AA=1.0d0 @@ -91,6 +95,14 @@ t3 = 0.5d0*(1.0d0+sqrt(1.0d0-t2**2)) valf = sqrt(Bxval**2/t1/t3)
+ write(*,*)'Alfven velocity:',valf + +!!$ Using DZ parameters with P=1.0, we have: +!!$ t1=5.5 +!!$ t2=4/11 +!!$ t3=0.96577 +!!$ valf=0.43389 + !!$ Vx value: vxval=0.0d0
@@ -103,45 +115,56 @@ dz = CCTK_DELTA_SPACE(3)
!!$ Note that the 3D test wasn't deviced to be used with AMR! - range_x = (cctk_gsh(1)-1)*dx - range_y = (cctk_gsh(2)-1)*dy - range_z = (cctk_gsh(3)-1)*dz + range_x = (cctk_gsh(1)-2*cctk_nghostzones(1))*dx + range_y = (cctk_gsh(2)-2*cctk_nghostzones(2))*dy + range_z = (cctk_gsh(3)-2*cctk_nghostzones(3))*dz
- range_d = sqrt(range_z**2+range_x**2) - - cos_theta = range_z/range_d - sin_theta = range_x/range_d - !!$ Alfven wave number - wnbr = 2.0d0*pi/range_x
do i=1,nx do j=1,ny do k=1,nz - - rho(i,j,k)=rhoval + + rho(i,j,k)=rhoval press(i,j,k)=pressval eps(i,j,k)=epsval - velx(i,j,k)=vxval - vely(i,j,k)=-valf*AA*cos(wnbr*x(i,j,k)) - velz(i,j,k)=-valf*AA*sin(wnbr*x(i,j,k)) - Bvecx(i,j,k)=Bxval - Bvecy(i,j,k)=AA*Bxval*cos(wnbr*x(i,j,k)) - Bvecz(i,j,k)=AA*Bxval*sin(wnbr*x(i,j,k)) + + if (CCTK_EQUALS(alfvenwave_type,"1D")) then
+ wnbr = 2.0d0*pi/range_x
-!!$ if (CCTK_EQUALS(advectedloop_type,"3D")) then -!!$ Bvecx_d = cos_theta*Bvecx(i,j,k)+sin_theta*Bvecz(i,j,k) -!!$ Bvecz_d = cos_theta*Bvecz(i,j,k)-sin_theta*Bvecx(i,j,k) -!!$ velx_d = cos_theta*velx(i,j,k)+sin_theta*velz(i,j,k) -!!$ velz_d = cos_theta*velz(i,j,k)-sin_theta*velx(i,j,k) -!!$ -!!$ Bvecx(i,j,k) = Bvecx_d -!!$ Bvecz(i,j,k) = Bvecz_d -!!$ velx(i,j,k) = velx_d -!!$ velz(i,j,k) = velz_d -!!$ end if - + velx(i,j,k)=vxval + vely(i,j,k)=-valf*AA*cos(wnbr*x(i,j,k)) + velz(i,j,k)=-valf*AA*sin(wnbr*x(i,j,k)) + Bvecx(i,j,k)=Bxval + Bvecy(i,j,k)=AA*Bxval*cos(wnbr*x(i,j,k)) + Bvecz(i,j,k)=AA*Bxval*sin(wnbr*x(i,j,k)) + + else if (CCTK_EQUALS(alfvenwave_type,"2D")) then + + diaglength=range_x*range_y/range_d + range_d = sqrt(range_x**2+range_y**2) + cos_theta = range_y/range_d + sin_theta = range_x/range_d + wnbr = 2.0d0*pi/diaglength + + xnew = cos_theta*x(i,j,k)+sin_theta*y(i,j,k) + + vparallel=vxval + vperp=-valf*AA*cos(wnbr*xnew) + velx(i,j,k)=vparallel*cos_theta-vperp*sin_theta + vely(i,j,k)=vparallel*sin_theta+vperp*cos_theta + velz(i,j,k)=-valf*AA*sin(wnbr*xnew) + Bparallel=Bxval + Bperp=AA*Bxval*cos(wnbr*xnew) + Bvecx(i,j,k)=Bparallel*cos_theta-Bperp*sin_theta + Bvecy(i,j,k)=Bparallel*sin_theta+Bperp*cos_theta + Bvecz(i,j,k)=AA*Bxval*sin(wnbr*xnew) + + else + call CCTK_WARN(0,"Alfven wave case not recognized!") + end if + 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: +14 -0 =================================================================== --- trunk/param.ccl 2012-08-09 06:26:57 UTC (rev 149) +++ trunk/param.ccl 2012-08-09 06:27:00 UTC (rev 150) @@ -191,6 +191,20 @@ "Numeric" :: "Finite difference approximation of the derivatives applied" } "Exact"
+#################################3 +# Alfven Wave test +################################3 +KEYWORD alfvenwave_type "1-dimensional or 2-dimensional?" +{ + "1D" :: "1-dimensional" + "2D" :: "2-dimensional (in x-y plane)" +} "1D" + +CCTK_REAL alfvenwave_pressure "P_gas for the Alfven wave" +{ + (0:* :: "positive" +} 1.0 + ################################################################################## # BONDI PARAMETERS: (black hole mass specified by parameters from the "Exact" thorn) ##################################################################################