User: bmundim
Date: 2010/10/08 03:29 PM
Added:
/trunk/src/
GRHydro_C2P2CM.F90, GRHydro_Macros.h, GRHydro_P2C2P.F90, GRHydro_P2C2PM.F90, GRHydro_ShockTubeM.F90
Modified:
/trunk/
interface.ccl, param.ccl, schedule.ccl
/trunk/src/
CheckParam.c, make.code.defn
Log:
Add the magnetized routine versions for prim2con2prim
and vice versa.
Add the magnetized version for the ShockTube routine.
File Changes:
Directory: /trunk/src/
======================
File [modified]: CheckParam.c
Delta lines: +13 -0
===================================================================
--- trunk/src/CheckParam.c 2010-08-18 18:42:06 UTC (rev 113)
+++ trunk/src/CheckParam.c 2010-10-08 20:29:50 UTC (rev 114)
@@ -11,5 +11,18 @@
{
CCTK_PARAMWARN("You have to set 'HydroBase::timelevels to at least 2");
}
+
+ if(CCTK_Equals(Bvec_evolution_method,"GRHydro") &&
+ ((CCTK_Equals(initial_hydro,"ony_atmo")) ||
+ (CCTK_Equals(initial_hydro,"read_conformal")) ||
+ (CCTK_Equals(initial_hydro,"simple_wave")) ||
+ (CCTK_Equals(initial_data,"con2primtest")) ||
+ (CCTK_Equals(initial_data,"reconstruction_test")) ||
+ (CCTK_Equals(shocktube_type,"diagshock")) ||
+ (CCTK_Equals(shocktube_type,"sphere"))))
+ {
+ CCTK_PARAMWARN("That test not yet implemented in MHD!");
+ }
+
}
File [added]: GRHydro_C2P2CM.F90
Delta lines: +182 -0
===================================================================
--- trunk/src/GRHydro_C2P2CM.F90 (rev 0)
+++ trunk/src/GRHydro_C2P2CM.F90 2010-10-08 20:29:50 UTC (rev 114)
@@ -0,0 +1,182 @@
+ /*@@
+ @file GRHydro_C2P2CM.F90
+ @date Sep 23, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ A test of the conservative <--> primitive variable exchange
+ @enddesc
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+ @routine c2p2cM
+ @date Sep 23, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ Testing the conservative <--> primitive variable transformations.
+ The values before and after should match.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+@@*/
+
+subroutine c2p2cM(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+ integer didit,i,j,k,nx,ny,nz
+ CCTK_REAL det
+ CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
+ CCTK_REAL gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send
+ CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send
+ CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
+ CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
+ CCTK_REAL bvcx_send, bvcy_send, bvcz_send, b2_send
+ CCTK_REAL pmin, epsmin
+ CCTK_INT C2P_failed
+ logical epsnegative
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+! end EOS Omni vars
+#endif
+
+ call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ x_send = 0.0d0
+ y_send = 0.0d0
+ z_send = 0.0d0
+ r_send = 0.0d0
+
+ gxx_send = 1.0d0
+ gyy_send = 1.0d0
+ gzz_send = 1.0d0
+ gxy_send = 0.0d0
+ gxz_send = 0.0d0
+ gyz_send = 0.0d0
+
+ det = 1.0d0
+
+ uxx = 1.0d0
+ uyy = 1.0d0
+ uzz = 1.0d0
+ uxy = 0.0d0
+ uxz = 0.0d0
+ uyz = 0.0d0
+
+ dens_send = 1.29047362d0
+ sx_send = 0.166666658d0
+ sy_send = 0.166666658d0
+ sz_send = 0.166666658d0
+ tau_send = 0.484123939d0
+
+ bvcx_send = Bx_init
+ bvcy_send = By_init
+ bvcz_send = Bz_init
+
+ eps_send = 1.0d-6
+ press_send = 6.666666666666667d-7
+ w_lorentz_send = 1.0d0
+
+ epsnegative = .false.
+
+ GRHydro_rho_min = 1.e-10
+#if USE_EOS_OMNI
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+#else
+ pmin = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
+ epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
+#endif
+ C2P_failed = 0
+
+ write(*,*) 'C2P2CM test: initial values.'
+ write(*,*) ' conservative variables: '
+ write(*,*) ' dens: ',dens_send
+ write(*,*) ' sx : ',sx_send
+ write(*,*) ' sy : ',sy_send
+ write(*,*) ' sz : ',sz_send
+ write(*,*) ' tau : ',tau_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ write(*,*) 'C2P2CM test: getting the associated primitive variables.'
+ call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,dens_send,sx_send,sy_send,sz_send, &
+ tau_send,rho_send,velx_send,vely_send,velz_send, &
+ eps_send,press_send,w_lorentz_send, &
+ gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
+ uxx,uxy,uxz,uyy,uyz,uzz,det,&
+ bvcx_send,bvcy_send,bvcz_send,b2_send,&
+ epsnegative,C2P_failed)
+
+ write(*,*) 'C2P2CM test: the primitive variables are'
+ write(*,*) ' primitive variables: '
+ write(*,*) ' rho : ',rho_send
+ write(*,*) ' velx : ',velx_send
+ write(*,*) ' vely : ',vely_send
+ write(*,*) ' velz : ',velz_send
+ write(*,*) ' press : ',press_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ write(*,*) 'C2P2CM test: converting back to conserved variables.'
+ call prim2conM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+ dens_send, sx_send, sy_send, sz_send, tau_send, bvcx_send, bvcy_send, bvcz_send, rho_send, &
+ velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send)
+
+ write(*,*) 'C2P2CM test: the conserved variables are'
+ write(*,*) ' conservative variables: '
+ write(*,*) ' dens: ',dens_send
+ write(*,*) ' sx : ',sx_send
+ write(*,*) ' sy : ',sy_send
+ write(*,*) ' sz : ',sz_send
+ write(*,*) ' tau : ',tau_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ STOP
+
+ return
+
+end subroutine c2p2cM
File [added]: GRHydro_Macros.h
Delta lines: +12 -0
===================================================================
--- trunk/src/GRHydro_Macros.h (rev 0)
+++ trunk/src/GRHydro_Macros.h 2010-10-08 20:29:50 UTC (rev 114)
@@ -0,0 +1,12 @@
+#define SPATIAL_DETERMINANT(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_) \
+ (-(gxz_)**2*(gyy_) + 2*(gxy_)*(gxz_)*(gyz_) - (gxx_)*(gyz_)**2 - (gxy_)**2*(gzz_) \
+ + (gxx_)*(gyy_)*(gzz_))
+
+#define DOTP(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x1_,y1_,z1_,x2_,y2_,z2_) \
+ ( (gxx_)*(x1_)*(x2_)+(gyy_)*(y1_)*(y2_)+(gzz_)*(z1_)*(z2_)+ \
+ (gxy_)*( (x1_)*(y2_)+(y1_)*(x2_) )+(gxz_)*( (x1_)*(z2_)+(z1_)*(x2_) )+\
+ (gyz_)*( (y1_)*(z2_)+(z1_)*(y2_) ) )
+
+#define DOTP2(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x_,y_,z_) \
+ ( (gxx_)*(x_)**2+(gyy_)*(y_)**2+(gzz_)*(z_)**2+ \
+ 2.0*( (gxy_)*(x_)*(y_)+(gxz_)*(x_)*(z_)+(gyz_)*(y_)*(z_) ) )
File [added]: GRHydro_P2C2P.F90
Delta lines: +167 -0
===================================================================
--- trunk/src/GRHydro_P2C2P.F90 (rev 0)
+++ trunk/src/GRHydro_P2C2P.F90 2010-10-08 20:29:50 UTC (rev 114)
@@ -0,0 +1,167 @@
+ /*@@
+ @file GRHydro_P2C2P.F90
+ @date Sep 25, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ A test of the conservative <--> primitive variable exchange
+ @enddesc
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+ @routine p2c2p
+ @date Sep 25, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ Testing the conservative <--> primitive variable transformations.
+ The values before and after should match.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+@@*/
+
+subroutine p2c2p(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+ integer didit,i,j,k,nx,ny,nz
+ CCTK_REAL det
+ CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
+ CCTK_REAL gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send
+ CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send
+ CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
+ CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
+ CCTK_REAL pmin, epsmin
+ CCTK_INT C2P_failed
+ logical epsnegative
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+! end EOS Omni vars
+#endif
+
+ call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ x_send = 0.0d0
+ y_send = 0.0d0
+ z_send = 0.0d0
+ r_send = 0.0d0
+
+ gxx_send = 1.0d0
+ gyy_send = 1.0d0
+ gzz_send = 1.0d0
+ gxy_send = 0.0d0
+ gxz_send = 0.0d0
+ gyz_send = 0.0d0
+
+ det = 1.0d0
+
+ uxx = 1.0d0
+ uyy = 1.0d0
+ uzz = 1.0d0
+ uxy = 0.0d0
+ uxz = 0.0d0
+ uyz = 0.0d0
+
+ rho_send = 1.29047362d0
+ velx_send = 0.166666658d0
+ vely_send = 0.166666658d0
+ velz_send = 0.166666658d0
+ eps_send = 0.484123939d0
+
+ w_lorentz_send = 1.d0/sqrt(1.0d0-velx_send*velx_send-vely_send*vely_send-velz_send*velz_send)
+
+ epsnegative = .false.
+
+ GRHydro_rho_min = 1.e-10
+#if USE_EOS_OMNI
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+ rho_send,1.0d0,xtemp,eps_send,press_send,keyerr,anyerr)
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+#else
+ press_send = EOS_Pressure(GRHydro_eos_handle, rho_send, eps_send)
+ pmin = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
+ epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
+#endif
+ C2P_failed = 0
+
+ write(*,*) 'P2C2P test: the primitive variables are'
+ write(*,*) ' primitive variables: '
+ write(*,*) ' rho : ',rho_send
+ write(*,*) ' velx : ',velx_send
+ write(*,*) ' vely : ',vely_send
+ write(*,*) ' velz : ',velz_send
+ write(*,*) ' press : ',press_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+
+ write(*,*) 'P2C2P test: converting back to conserved variables.'
+ call prim2con(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+ dens_send, sx_send, sy_send, sz_send, tau_send, rho_send, &
+ velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send)
+
+ write(*,*) 'P2C2P test: initial values.'
+ write(*,*) ' conservative variables: '
+ write(*,*) ' dens: ',dens_send
+ write(*,*) ' sx : ',sx_send
+ write(*,*) ' sy : ',sy_send
+ write(*,*) ' sz : ',sz_send
+ write(*,*) ' tau : ',tau_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+
+ write(*,*) 'P2C2P test: getting the associated primitive variables.'
+ call Con2Prim_pt(GRHydro_eos_handle,dens_send,sx_send,sy_send,sz_send, &
+ tau_send,rho_send,velx_send,vely_send,velz_send, &
+ eps_send,press_send,w_lorentz_send, &
+ uxx,uxy,uxz,uyy,uyz,uzz,det,x_send,y_send,z_send,r_send,&
+ epsnegative,GRHydro_rho_min, pmin, epsmin, GRHydro_init_data_reflevel,C2P_failed)
+
+ write(*,*) 'P2C2P test: the primitive variables are'
+ write(*,*) ' primitive variables: '
+ write(*,*) ' rho : ',rho_send
+ write(*,*) ' velx : ',velx_send
+ write(*,*) ' vely : ',vely_send
+ write(*,*) ' velz : ',velz_send
+ write(*,*) ' press : ',press_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+
+ STOP
+
+ return
+
+end subroutine p2c2p
File [added]: GRHydro_P2C2PM.F90
Delta lines: +183 -0
===================================================================
--- trunk/src/GRHydro_P2C2PM.F90 (rev 0)
+++ trunk/src/GRHydro_P2C2PM.F90 2010-10-08 20:29:50 UTC (rev 114)
@@ -0,0 +1,183 @@
+ /*@@
+ @file GRHydro_P2C2PM.F90
+ @date Sep 25, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ A test of the conservative <--> primitive variable exchange
+ @enddesc
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+ @routine p2c2pm
+ @date Sep 25, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Luca Baiotti
+ @desc
+ Testing the conservative <--> primitive variable transformations.
+ The values before and after should match.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+@@*/
+
+subroutine p2c2pm(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+ integer didit,i,j,k,nx,ny,nz
+ CCTK_REAL det
+ CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz
+ CCTK_REAL gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send
+ CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send
+ CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send
+ CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send
+ CCTK_REAL bvcx_send,bvcy_send,bvcz_send,b2_send
+ CCTK_REAL pmin, epsmin
+ CCTK_INT C2P_failed
+ logical epsnegative
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+! end EOS Omni vars
+#endif
+
+ call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ x_send = 0.0d0
+ y_send = 0.0d0
+ z_send = 0.0d0
+ r_send = 0.0d0
+
+ gxx_send = 1.0d0
+ gyy_send = 1.0d0
+ gzz_send = 1.0d0
+ gxy_send = 0.0d0
+ gxz_send = 0.0d0
+ gyz_send = 0.0d0
+
+ det = 1.0d0
+
+ uxx = 1.0d0
+ uyy = 1.0d0
+ uzz = 1.0d0
+ uxy = 0.0d0
+ uxz = 0.0d0
+ uyz = 0.0d0
+
+ rho_send = 1.29047362d0
+ velx_send = 0.166666658d0
+ vely_send = 0.166666658d0
+ velz_send = 0.166666658d0
+ eps_send = 0.484123939d0
+
+ bvcx_send = Bx_init
+ bvcy_send = By_init
+ bvcz_send = Bz_init
+
+ w_lorentz_send = 1.d0/sqrt(1.0d0-velx_send*velx_send-vely_send*vely_send-velz_send*velz_send)
+
+ epsnegative = .false.
+
+ GRHydro_rho_min = 1.e-10
+#if USE_EOS_OMNI
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+ rho_send,1.0d0,xtemp,eps_send,press_send,keyerr,anyerr)
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+
+ call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,&
+ GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+#else
+ press_send = EOS_Pressure(GRHydro_eos_handle, rho_send, eps_send)
+ pmin = EOS_Pressure(GRHydro_eos_handle, GRHydro_rho_min, 1.0d0)
+ epsmin = EOS_SpecificIntEnergy(GRHydro_eos_handle, GRHydro_rho_min, pmin)
+#endif
+ C2P_failed = 0
+
+ write(*,*) 'P2C2PM test: the primitive variables are'
+ write(*,*) ' primitive variables: '
+ write(*,*) ' rho : ',rho_send
+ write(*,*) ' velx : ',velx_send
+ write(*,*) ' vely : ',vely_send
+ write(*,*) ' velz : ',velz_send
+ write(*,*) ' press : ',press_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ write(*,*) 'P2C2PM test: converting back to conserved variables.'
+ call prim2conM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+ dens_send, sx_send, sy_send, sz_send, tau_send, bvcx_send, bvcy_send, bvcz_send, rho_send, &
+ velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send)
+
+ write(*,*) 'P2C2PM test: initial values.'
+ write(*,*) ' conservative variables: '
+ write(*,*) ' dens: ',dens_send
+ write(*,*) ' sx : ',sx_send
+ write(*,*) ' sy : ',sy_send
+ write(*,*) ' sz : ',sz_send
+ write(*,*) ' tau : ',tau_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ write(*,*) 'P2C2PM test: getting the associated primitive variables.'
+ call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,dens_send,sx_send,sy_send,sz_send, &
+ tau_send,rho_send,velx_send,vely_send,velz_send, &
+ eps_send,press_send,w_lorentz_send, &
+ gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
+ uxx,uxy,uxz,uyy,uyz,uzz,det,&
+ bvcx_send,bvcy_send,bvcz_send,b2_send,&
+ epsnegative,C2P_failed)
+
+ write(*,*) 'P2C2PM test: the primitive variables are'
+ write(*,*) ' primitive variables: '
+ write(*,*) ' rho : ',rho_send
+ write(*,*) ' velx : ',velx_send
+ write(*,*) ' vely : ',vely_send
+ write(*,*) ' velz : ',velz_send
+ write(*,*) ' press : ',press_send
+ write(*,*) ' eps : ',eps_send
+ write(*,*) ' W : ',w_lorentz_send
+ write(*,*) ' Bvecx : ',bvcx_send
+ write(*,*) ' Bvecy : ',bvcy_send
+ write(*,*) ' Bvecz : ',bvcz_send
+
+ STOP
+
+ return
+
+end subroutine p2c2pm
File [added]: GRHydro_ShockTubeM.F90
Delta lines: +174 -0
===================================================================
--- trunk/src/GRHydro_ShockTubeM.F90 (rev 0)
+++ trunk/src/GRHydro_ShockTubeM.F90 2010-10-08 20:29:50 UTC (rev 114)
@@ -0,0 +1,174 @@
+ /*@@
+ @file GRHydro_ShockTubeM.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_shocktubeM
+ @date Sat Jan 26 02:53:49 2002
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+ @desc
+ Initial data for shock tubes, parallel to
+ a coordinate axis. Either Sods problem or the standard shock tube.
+ @enddesc
+ @calls
+ @calledby
+ @history
+ Expansion and alteration of the test code from GRAstro_Hydro,
+ written by Mark Miller.
+ @endhistory
+
+@@*/
+
+subroutine GRHydro_shocktubeM(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
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ do i=1,nx
+ do j=1,ny
+ do k=1,nz
+ if (CCTK_EQUALS(shocktube_type,"diagshock")) then
+ direction = x(i,j,k) - shock_xpos + &
+ y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos
+ else if (CCTK_EQUALS(shocktube_type,"xshock")) then
+ direction = x(i,j,k) - shock_xpos
+ else if (CCTK_EQUALS(shocktube_type,"yshock")) then
+ direction = y(i,j,k) - shock_ypos
+ else if (CCTK_EQUALS(shocktube_type,"zshock")) then
+ direction = z(i,j,k) - shock_zpos
+ else if (CCTK_EQUALS(shocktube_type,"sphere")) then
+ direction = sqrt((x(i,j,k)-shock_xpos)**2+&
+ (y(i,j,k)-shock_ypos)**2+&
+ (z(i,j,k)-shock_zpos)**2)-shock_radius
+ end if
+ if (CCTK_EQUALS(shock_case,"Simple")) then
+ rhol = 10.d0
+ rhor = 1.d0
+ velxl = 0.d0
+ velxr = 0.d0
+ velyl = 0.d0
+ velyr = 0.d0
+ velzl = 0.d0
+ velzr = 0.d0
+ epsl = 2.d0
+ epsr = 1.d-6
+ else if (CCTK_EQUALS(shock_case,"Sod")) then
+ rhol = 1.d0
+ rhor = 0.125d0
+ velxl = 0.d0
+ velxr = 0.d0
+ velyl = 0.d0
+ velyr = 0.d0
+ velzl = 0.d0
+ velzr = 0.d0
+ epsl = 1.5d0
+ epsr = 1.2d0
+!!$This line only for polytrope, k=1
+!!$ epsr = 0.375d0
+ else if (CCTK_EQUALS(shock_case,"Blast")) then
+ rhol = 1.d0
+ rhor = 1.d0
+ velxl = 0.d0
+ velxr = 0.d0
+ velyl = 0.d0
+ velyr = 0.d0
+ velzl = 0.d0
+ velzr = 0.d0
+ epsl = 1500.d0
+ epsr = 1.5d-2
+ else
+ call CCTK_WARN(0,"Shock case not recognized")
+ end if
+
+ bvcxl = Bx_init
+ bvcyl = By_init
+ bvczl = Bz_init
+ bvcxr = Bx_init
+ bvcyr = By_init
+ bvczr = Bz_init
+
+ if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.&
+ ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then
+ rho(i,j,k) = rhol
+ velx(i,j,k) = velxl
+ vely(i,j,k) = velyl
+ velz(i,j,k) = velzl
+ eps(i,j,k) = epsl
+ Bvecx(i,j,k)=bvcxl
+ Bvecy(i,j,k)=bvcyl
+ Bvecz(i,j,k)=bvczl
+ else
+ rho(i,j,k) = rhor
+ velx(i,j,k) = velxr
+ vely(i,j,k) = velyr
+ velz(i,j,k) = velzr
+ eps(i,j,k) = epsr
+ Bvecx(i,j,k)=bvcxr
+ Bvecy(i,j,k)=bvcyr
+ Bvecz(i,j,k)=bvczr
+ 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
+ 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_shocktubeM
File [modified]: make.code.defn
Delta lines: +5 -1
===================================================================
--- trunk/src/make.code.defn 2010-08-18 18:42:06 UTC (rev 113)
+++ trunk/src/make.code.defn 2010-10-08 20:29:50 UTC (rev 114)
@@ -3,13 +3,17 @@
# Source files in this directory
SRCS = GRHydro_C2P2C.F90 \
+ GRHydro_C2P2CM.F90 \
GRHydro_Con2Prim.F90 \
GRHydro_ReconstructTest.F90 \
GRHydro_ShockTube.F90 \
+ GRHydro_ShockTubeM.F90 \
GRHydro_Only_Atmo.F90 \
GRHydro_ReadConformalData.F90 \
GRHydro_SimpleWave.F90 \
- CheckParam.c
+ CheckParam.c\
+ GRHydro_P2C2P.F90 \
+ GRHydro_P2C2PM.F90
# Subdirectories containing source files
SUBDIRS =
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +28 -0
===================================================================
--- trunk/interface.ccl 2010-08-18 18:42:06 UTC (rev 113)
+++ trunk/interface.ccl 2010-10-08 20:29:50 UTC (rev 114)
@@ -39,6 +39,32 @@
CCTK_REAL IN velz, CCTK_REAL IN epsilon, \
CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
+SUBROUTINE Prim2ConPolyM(CCTK_INT IN handle, \
+ CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \
+ CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
+ CCTK_REAL IN det, CCTK_REAL OUT dens, \
+ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
+ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
+ CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \
+ CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \
+ CCTK_REAL IN vely, \
+ CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \
+ CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
+
+
+SUBROUTINE Prim2ConGenM(CCTK_INT IN handle, \
+ CCTK_REAL IN gxx, CCTK_REAL IN gxy, \
+ CCTK_REAL IN gxz, CCTK_REAL IN gyy, \
+ CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
+ CCTK_REAL IN det, CCTK_REAL OUT dens, \
+ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
+ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
+ CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \
+ CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \
+ CCTK_REAL IN vely, \
+ CCTK_REAL IN velz, CCTK_REAL IN epsilon, \
+ CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz)
+
SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \
CCTK_REAL OUT sx, CCTK_REAL OUT sy, \
CCTK_REAL OUT sz, CCTK_REAL OUT tau, \
@@ -56,6 +82,8 @@
USES FUNCTION SpatialDet
USES FUNCTION Prim2ConPoly
USES FUNCTION Prim2ConGen
+USES FUNCTION Prim2ConPolyM
+USES FUNCTION Prim2ConGenM
USES FUNCTION Con2PrimPoly
CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name)
File [modified]: param.ccl
Delta lines: +18 -0
===================================================================
--- trunk/param.ccl 2010-08-18 18:42:06 UTC (rev 113)
+++ trunk/param.ccl 2010-10-08 20:29:50 UTC (rev 114)
@@ -4,6 +4,7 @@
shares:HydroBase
USES CCTK_INT timelevels
+USES KEYWORD Bvec_evolution_method
EXTENDS KEYWORD initial_hydro ""
{
@@ -20,6 +21,7 @@
# "shocktube" :: "Shock tube initial data for GRHydro"
"con2primtest" :: "Testing the con -> prim conversion"
"con2prim2con_test" :: "Testing the con -> prim -> con conversion"
+ "prim2con2prim_test" :: "Testing the prim -> con -> prim conversion"
"reconstruction_test" :: "Testing reconstruction"
}
@@ -86,6 +88,22 @@
{
} "no"
+REAL Bx_init "Initial B-field in the x-dir"
+{
+ *:* :: "Anything"
+} 0.0
+
+REAL By_init "Initial B-field in the y-dir"
+{
+ *:* :: "Anything"
+} 0.0
+
+REAL Bz_init "Initial B-field in the z-dir"
+{
+ *:* :: "Anything"
+} 0.0
+
+
shares:GRHydro
USES real GRHydro_rho_central
File [modified]: schedule.ccl
Delta lines: +45 -6
===================================================================
--- trunk/schedule.ccl 2010-08-18 18:42:06 UTC (rev 113)
+++ trunk/schedule.ccl 2010-10-08 20:29:50 UTC (rev 114)
@@ -6,10 +6,20 @@
} "Check parameters"
if (CCTK_Equals(initial_hydro,"shocktube")) {
- schedule GRHydro_shocktube in HydroBase_Initial
+ if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
- LANG: Fortran
- } "Shocktube initial data"
+ schedule GRHydro_shocktubeM in HydroBase_Initial
+ {
+ LANG: Fortran
+ } "Shocktube initial data - MHD version"
+
+ } else {
+ schedule GRHydro_shocktube in HydroBase_Initial
+ {
+ LANG: Fortran
+ } "Shocktube initial data"
+
+ }
}
if (CCTK_Equals(initial_data,"con2primtest")) {
@@ -27,15 +37,44 @@
if (CCTK_Equals(initial_data,"con2prim2con_test")) {
STORAGE:GRHydro_init_data_reflevel
- schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE c2p2c
+ schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE c2p2c_call
{
LANG: Fortran
} "Calculate current refinement level"
- schedule c2p2c in HydroBase_Initial
+ if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
+ schedule c2p2cM in HydroBase_Initial AS c2p2c_call
+ {
+ LANG: Fortran
+ } "Testing conservative to primitive to conservative - MHD version"
+ } else {
+ schedule c2p2c in HydroBase_Initial AS c2p2c_call
+ {
+ LANG: Fortran
+ } "Testing conservative to primitive to conservative"
+ }
+}
+
+if (CCTK_Equals(initial_data,"prim2con2prim_test")) {
+ STORAGE:GRHydro_init_data_reflevel
+ schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE p2c2p_call
+ {
LANG: Fortran
- } "Testing conservative to primitive to conservative"
+ } "Calculate current refinement level"
+
+ if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
+ {
+ schedule p2c2pM in HydroBase_Initial AS p2c2p_call
+ {
+ LANG: Fortran
+ } "Testing primitive to conservative to primitive - MHD version"
+ } else {
+ schedule p2c2p in HydroBase_Initial AS p2c2p_call
+ {
+ LANG: Fortran
+ } "Testing primitive to conservative to primitive"
+ }
}
if (CCTK_Equals(initial_data,"reconstruction_test")) {
