[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 124)
bcmsma at astro.rit.edu
bcmsma at astro.rit.edu
Wed Apr 27 17:26:39 CDT 2011
User: bmundim
Date: 2011/04/27 05:26 PM
Added:
/trunk/src/
GRHydro_CylindricalExplosionM.F90, GRHydro_MonopoleM.F90
Log:
RIT MHD dev:
1) Komissarov's cylindrical explosion test
2) Magnetic monopole test
File Changes:
Directory: /trunk/src/
======================
File [added]: GRHydro_CylindricalExplosionM.F90
Delta lines: +164 -0
===================================================================
--- trunk/src/GRHydro_CylindricalExplosionM.F90 (rev 0)
+++ trunk/src/GRHydro_CylindricalExplosionM.F90 2011-04-27 22:26:38 UTC (rev 124)
@@ -0,0 +1,164 @@
+ /*@@
+ @file GRHydro_CylindricalExplosionM.F90
+ @date Apr 22, 2011
+ @author Scott Noble, Joshua Faber, Bruno Mundim
+ @desc
+ Cylindrical magnetized shocks.
+ @enddesc
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "GRHydro_Macros.h"
+
+#define velx(i,j,k) vel(i,j,k,1)
+#define vely(i,j,k) vel(i,j,k,2)
+#define velz(i,j,k) vel(i,j,k,3)
+#define sx(i,j,k) scon(i,j,k,1)
+#define sy(i,j,k) scon(i,j,k,2)
+#define sz(i,j,k) scon(i,j,k,3)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+
+
+ /*@@
+ @routine GRHydro_cylindricalexplosionM
+ @date Fri Apr 22 14:51:37 EDT 2011
+ @author Scott Noble, Joshua Faber, Bruno Mundim
+ @desc
+ Initial data for cylidrically symmetric shocks with
+ magnetic fields. Meant to recreate tests demonstrated
+ in Komissarov (1999).
+ @enddesc
+ @calls
+ @calledby
+ @history
+ Using GRHydro_ShockTubeM.F90 as a template.
+ @endhistory
+
+@@*/
+
+subroutine GRHydro_cylindricalexplosionM(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
+
+ CCTK_INT :: i, j, k, nx, ny, nz
+ CCTK_REAL :: direction, det
+ CCTK_REAL :: rhol, rhor, pressl, pressr
+ CCTK_REAL :: bvcxl,bvcyl,bvczl
+ CCTK_REAL :: tmp,tmp2,gam,r_inner,r_outer
+
+ CCTK_REAL :: cyl_fr
+
+!!$Uses Bx_init, By_init and Bz_init to set magnetic field strength
+!!$Original tests had Bx_init = 0.1, 1.0 with By_init=Bz_init = 0
+
+ bvcxl = Bx_init
+ bvcyl = By_init
+ bvczl = Bz_init
+
+!!$Adiabatic index for test
+ gam = (4.d0/3.d0)
+
+!!$Inner radius and outer radius
+ r_inner = 8.d-1
+ r_outer = 1.d0
+
+!!$Inner values
+ rhol = 1.d-4
+ pressl = 3.d-5
+
+!!$Outer values
+ rhor = 1.d-2
+ pressr = 1.d0
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ do i=1,nx
+ do j=1,ny
+ do k=1,nz
+
+!!$direction represents the cylindrical radius here.
+ if (CCTK_EQUALS(shocktube_type,"xshock")) then
+ direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+ (y(i,j,k)-shock_ypos)**2)
+ else if (CCTK_EQUALS(shocktube_type,"yshock")) then
+ direction = sqrt((y(i,j,k)-shock_ypos)**2+&
+ (z(i,j,k)-shock_zpos)**2)
+ else if (CCTK_EQUALS(shocktube_type,"zshock")) then
+ direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+ (z(i,j,k)-shock_zpos)**2)
+ end if
+
+ tmp = cyl_fr(direction,r_inner,r_outer,rhol,rhor)
+ tmp2 = cyl_fr(direction,r_inner,r_outer,pressl,pressr)/( (gam - 1.d0) * tmp )
+ rho(i,j,k) = tmp
+ eps(i,j,k) = tmp2
+
+ velx(i,j,k) = 0.d0
+ vely(i,j,k) = 0.d0
+ velz(i,j,k) = 0.d0
+ Bvecx(i,j,k)=bvcxl
+ Bvecy(i,j,k)=bvcyl
+ Bvecz(i,j,k)=bvczl
+
+ 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
+ call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
+ gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+ tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+ velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+ eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+ else
+ call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
+ gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+ tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+ velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+ eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+ end if
+ enddo
+ enddo
+ enddo
+
+ densrhs = 0.d0
+ srhs = 0.d0
+ taurhs = 0.d0
+ Bvecrhs = 0.d0
+
+
+ return
+
+end subroutine GRHydro_CylindricalExplosionM
+
+
+
+function cyl_fr(R,r_inner,r_outer,min_f,max_f)
+
+ implicit none
+
+ CCTK_REAL :: cyl_fr
+ CCTK_REAL :: R,r_inner,r_outer,min_f,max_f
+
+ if(R .gt. r_outer) then
+ cyl_fr = min_f
+ else if( R .lt. r_inner) then
+ cyl_fr = max_f
+ else
+ cyl_fr = exp( ( log(max_f)*(r_outer - R) + log(min_f)*(R - r_inner) ) / (r_outer - r_inner) )
+ end if
+
+ return
+
+end function cyl_fr
File [added]: GRHydro_MonopoleM.F90
Delta lines: +109 -0
===================================================================
--- trunk/src/GRHydro_MonopoleM.F90 (rev 0)
+++ trunk/src/GRHydro_MonopoleM.F90 2011-04-27 22:26:38 UTC (rev 124)
@@ -0,0 +1,109 @@
+ /*@@
+ @file GRHydro_MonopoleM.F90
+ @date Sep 23, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+ @desc
+ Initial data of the shock tube type - MHD version.
+ @enddesc
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "GRHydro_Macros.h"
+
+#define velx(i,j,k) vel(i,j,k,1)
+#define vely(i,j,k) vel(i,j,k,2)
+#define velz(i,j,k) vel(i,j,k,3)
+#define sx(i,j,k) scon(i,j,k,1)
+#define sy(i,j,k) scon(i,j,k,2)
+#define sz(i,j,k) scon(i,j,k,3)
+#define Bvecx(i,j,k) Bvec(i,j,k,1)
+#define Bvecy(i,j,k) Bvec(i,j,k,2)
+#define Bvecz(i,j,k) Bvec(i,j,k,3)
+
+ /*@@
+ @routine GRHydro_MonopoleM
+ @date Sat Jan 26 02:53:49 2002
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+ @desc
+ A monopole in space
+ @enddesc
+ @calls
+ @calledby
+ @history
+ Expansion and alteration of the test code from GRAstro_Hydro,
+ written by Mark Miller.
+ @endhistory
+
+@@*/
+
+subroutine GRHydro_MonopoleM(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
+
+ CCTK_INT :: i, j, k, nx, ny, nz
+ CCTK_REAL :: direction, det
+ CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, &
+ velzl, velzr, epsl, epsr
+ CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr
+ CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3
+
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ do i=1,nx
+ do j=1,ny
+ do k=1,nz
+
+ rho(i,j,k) = 1.0
+ velx(i,j,k) = 0.0
+ 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
+ else
+ Bvecx(i,j,k)=0.0
+ 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
+ call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
+ gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+ tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+ velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+ eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+ else
+ call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
+ gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+ tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),&
+ velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+ eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+ end if
+
+ enddo
+ enddo
+ enddo
+
+ densrhs = 0.d0
+ srhs = 0.d0
+ taurhs = 0.d0
+ Bvecrhs = 0.d0
+
+
+ return
+
+end subroutine GRHydro_MonopoleM
More information about the Commits
mailing list