From 85cb2e6f9c6191fc9889a1055017236c9a071192 Mon Sep 17 00:00:00 2001 From: bmundim Date: Sat, 25 Dec 2010 06:44:22 +0000 Subject: RIT MHD dev: Update initial data routines to properly call out pressures and other hydro quantities using EOS_Omni rather than now-defunct interfaces. Add C2P_polytype interface. Un-hardwire fixed values of gamma for all MHD Con2Prim routines. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@117 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- interface.ccl.omni | 84 ------------------ param.ccl | 7 ++ schedule.ccl | 16 ++++ src/CheckParam.c | 11 ++- src/GRHydro_C2P2C.F90 | 20 +---- src/GRHydro_C2P2CM.F90 | 29 +++--- src/GRHydro_Con2Prim.F90 | 21 +---- src/GRHydro_P2C2P.F90 | 25 ++---- src/GRHydro_P2C2PM.F90 | 34 +++---- src/GRHydro_P2C2PM_polytype.F90 | 181 ++++++++++++++++++++++++++++++++++++++ src/GRHydro_ReadConformalData.F90 | 15 +--- src/GRHydro_ShockTube.F90 | 4 +- src/GRHydro_ShockTubeM.F90 | 101 +++++++++++++++++++-- src/GRHydro_SimpleWave.F90 | 37 +++----- src/make.code.defn | 3 +- 15 files changed, 363 insertions(+), 225 deletions(-) delete mode 100644 interface.ccl.omni create mode 100644 src/GRHydro_P2C2PM_polytype.F90 diff --git a/interface.ccl.omni b/interface.ccl.omni deleted file mode 100644 index e8fe1ae..0000000 --- a/interface.ccl.omni +++ /dev/null @@ -1,84 +0,0 @@ -# Interface definition for thorn GRHydro_Init_Data -# $Header$ - -implements: GRHydro_init_data -inherits: GRHydro grid - -#USES INCLUDE: EOS_Base.inc -USES INCLUDE: SpaceMask.h -#USES INCLUDE: EOS_Base.h - - -SUBROUTINE SpatialDet(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 OUT det) - - -SUBROUTINE Prim2ConPoly(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 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 Prim2ConGen(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 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, \ - CCTK_REAL OUT rho, CCTK_REAL OUT velx, \ - CCTK_REAL OUT vely, CCTK_REAL OUT velz, \ - CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \ - CCTK_REAL OUT w_lorentz, CCTK_REAL IN uxx, \ - CCTK_REAL IN uxy, CCTK_REAL IN uxz, CCTK_REAL IN uyy, \ - CCTK_REAL IN uyz, CCTK_REAL IN uzz, CCTK_REAL IN det, \ - CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ - CCTK_REAL IN r, CCTK_REAL IN rho_min, \ - CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) - - -USES FUNCTION SpatialDet -USES FUNCTION Prim2ConPoly -USES FUNCTION Prim2ConGen -USES FUNCTION Con2PrimPoly - -protected: - -CCTK_REAL simple_wave_grid_functions TYPE=GF TIMELEVELS=1 tags='checkpoint="no"' -{ - simple_tmp - c_s -} "1D arrays for the simple-wave routine" - -CCTK_REAL simple_wave_scalars TYPE=scalar -{ - simple_rho_0 - simple_eps_0 -} "values at v=0" - -CCTK_REAL simple_wave_output TYPE=GF TIMELEVELS=1 tags='checkpoint="no"' -{ - simple_rho - simple_eps -# simple_entropy -} "output variables for the simple-wave routine" - -private: - -CCTK_INT GRHydro_init_data_reflevel type = SCALAR tags='checkpoint="no"' "Refinement level GRHydro is working on right now" diff --git a/param.ccl b/param.ccl index bc000b0..01fbe7f 100644 --- a/param.ccl +++ b/param.ccl @@ -22,6 +22,7 @@ EXTENDS KEYWORD initial_data "" "con2primtest" :: "Testing the con -> prim conversion" "con2prim2con_test" :: "Testing the con -> prim -> con conversion" "prim2con2prim_test" :: "Testing the prim -> con -> prim conversion" + "prim2con2prim_polytype_test" :: "Testing the prim -> con -> prim conversion - polytype version" "reconstruction_test" :: "Testing reconstruction" } @@ -41,6 +42,11 @@ KEYWORD shock_case "Simple, Sod's problem or other?" "Simple" :: "GRAstro_Hydro test case" "Sod" :: "Sod's problem" "Blast" :: "Strong blast wave" + "Balsara1" :: "Balsara Test #1" + "Balsara2" :: "Balsara Test #2" + "Balsara3" :: "Balsara Test #3" + "Balsara4" :: "Balsara Test #4" + "Balsara5" :: "Balsara Test #5" } "Sod" REAL shock_xpos "Position of shock plane: x" @@ -106,6 +112,7 @@ REAL Bz_init "Initial B-field in the z-dir" shares:GRHydro +USES real GRHydro_eos_rf_prec USES real GRHydro_rho_central USES real GRHydro_eps_min "" USES real GRHydro_perc_ptol "" diff --git a/schedule.ccl b/schedule.ccl index 26bdee4..a35dfdb 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -77,6 +77,22 @@ if (CCTK_Equals(initial_data,"prim2con2prim_test")) { } } +if (CCTK_Equals(initial_data,"prim2con2prim_polytype_test")) { + STORAGE:GRHydro_init_data_reflevel + schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE p2c2p_call + { + LANG: Fortran + } "Calculate current refinement level" + + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + { + schedule p2c2pM_polytype in HydroBase_Initial AS p2c2p_call + { + LANG: Fortran + } "Testing primitive to conservative to primitive - MHD polytype version" + } +} + if (CCTK_Equals(initial_data,"reconstruction_test")) { schedule GRHydro_reconstruction_test in HydroBase_Initial { diff --git a/src/CheckParam.c b/src/CheckParam.c index 31222eb..124db7f 100644 --- a/src/CheckParam.c +++ b/src/CheckParam.c @@ -23,6 +23,15 @@ void GRHydro_InitData_CheckParameters(CCTK_ARGUMENTS) { CCTK_PARAMWARN("That test not yet implemented in MHD!"); } - + + if(!CCTK_Equals(Bvec_evolution_method,"GRHydro") && + ((CCTK_Equals(shock_case,"Balsara1")) || + (CCTK_Equals(shock_case,"Balsara2")) || + (CCTK_Equals(shock_case,"Balsara3")) || + (CCTK_Equals(shock_case,"Balsara4")) || + (CCTK_Equals(shock_case,"Balsara5")))) + { + CCTK_PARAMWARN("That test requires MHD! Set Bvec_evolution_method=GRHYDRO!"); + } } diff --git a/src/GRHydro_C2P2C.F90 b/src/GRHydro_C2P2C.F90 index 1ce1347..a1cc36d 100644 --- a/src/GRHydro_C2P2C.F90 +++ b/src/GRHydro_C2P2C.F90 @@ -34,13 +34,6 @@ subroutine c2p2c(CCTK_ARGUMENTS) 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 @@ -52,7 +45,6 @@ subroutine c2p2c(CCTK_ARGUMENTS) CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -63,7 +55,6 @@ subroutine c2p2c(CCTK_ARGUMENTS) 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") @@ -105,16 +96,13 @@ subroutine c2p2c(CCTK_ARGUMENTS) epsnegative = .false. GRHydro_rho_min = 1.e-10 -#if USE_EOS_OMNI - call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,& + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,& + call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,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(*,*) 'C2P2C test: initial values.' diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90 index 12c1ed4..3bd61b6 100644 --- a/src/GRHydro_C2P2CM.F90 +++ b/src/GRHydro_C2P2CM.F90 @@ -34,13 +34,6 @@ subroutine c2p2cM(CCTK_ARGUMENTS) 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 @@ -49,11 +42,10 @@ subroutine c2p2cM(CCTK_ARGUMENTS) 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_REAL pmin, epsmin, local_gam CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -64,7 +56,6 @@ subroutine c2p2cM(CCTK_ARGUMENTS) 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") @@ -110,16 +101,18 @@ subroutine c2p2cM(CCTK_ARGUMENTS) epsnegative = .false. GRHydro_rho_min = 1.e-10 -#if USE_EOS_OMNI - call EOS_Omni_press(GRHydro_eos_handle,keytemp,n,& + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,& + call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,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 + + local_gam=0.0 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + 1.d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr) + local_gam = local_gam+1.0 + C2P_failed = 0 write(*,*) 'C2P2CM test: initial values.' @@ -136,7 +129,7 @@ subroutine c2p2cM(CCTK_ARGUMENTS) 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, & + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,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,& diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index d6950ab..6b9ba34 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -52,14 +52,6 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS) CCTK_INT C2P_failed logical epsnegative -#ifndef USE_EOS_OMNI -#ifdef _EOS_BASE_INC_ -#undef _EOS_BASE_INC_ -#endif -#include "EOS_Base.inc" -#endif - -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -70,7 +62,6 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS) real*8 :: xtemp = 0.0d0 real*8 :: xye = 0.0d0 ! end EOS Omni vars -#endif call CCTK_WARN(1,"For this test, remember to use a polytropic EoS and to set eos_gamma = 2.0 and eos_k = 100.0") @@ -107,16 +98,12 @@ subroutine GRHydro_con2primtest(CCTK_ARGUMENTS) epsnegative = .false. -#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_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,& + call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,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(*,*) 'Con2Prim test: converting to primitive variables' diff --git a/src/GRHydro_P2C2P.F90 b/src/GRHydro_P2C2P.F90 index 10258c3..ba13d12 100644 --- a/src/GRHydro_P2C2P.F90 +++ b/src/GRHydro_P2C2P.F90 @@ -34,13 +34,6 @@ subroutine p2c2p(CCTK_ARGUMENTS) 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 @@ -52,7 +45,6 @@ subroutine p2c2p(CCTK_ARGUMENTS) CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -63,7 +55,6 @@ subroutine p2c2p(CCTK_ARGUMENTS) 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") @@ -103,19 +94,15 @@ subroutine p2c2p(CCTK_ARGUMENTS) 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,& + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_send,eps_send,xtemp,xye,press_send,keyerr,anyerr) + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,& + call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,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' diff --git a/src/GRHydro_P2C2PM.F90 b/src/GRHydro_P2C2PM.F90 index c68b735..1abf1ca 100644 --- a/src/GRHydro_P2C2PM.F90 +++ b/src/GRHydro_P2C2PM.F90 @@ -34,13 +34,6 @@ subroutine p2c2pm(CCTK_ARGUMENTS) 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 @@ -49,11 +42,10 @@ subroutine p2c2pm(CCTK_ARGUMENTS) 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_REAL pmin, epsmin, local_gam CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -64,7 +56,6 @@ subroutine p2c2pm(CCTK_ARGUMENTS) 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") @@ -108,19 +99,19 @@ subroutine p2c2pm(CCTK_ARGUMENTS) 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,& + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_send,eps_send,xtemp,xye,press_send,keyerr,anyerr) + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) - call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,n,& + call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,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 + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + 1.0d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr) + local_gam=local_gam+1.0 + C2P_failed = 0 write(*,*) 'P2C2PM test: the primitive variables are' @@ -155,7 +146,8 @@ subroutine p2c2pm(CCTK_ARGUMENTS) 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, & + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,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,& diff --git a/src/GRHydro_P2C2PM_polytype.F90 b/src/GRHydro_P2C2PM_polytype.F90 new file mode 100644 index 0000000..a8d0e77 --- /dev/null +++ b/src/GRHydro_P2C2PM_polytype.F90 @@ -0,0 +1,181 @@ + /*@@ + @file GRHydro_P2C2PM_polytype.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_polytype(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer didit,i,j,k,nx,ny,nz + CCTK_REAL det, local_K + 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, local_gam, sc_send + CCTK_INT C2P_failed + logical epsnegative + +! 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 + + 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 + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_send,1.0d0,xtemp,xye,press_send,keyerr,anyerr) + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_send,xeps,xtemp,xye,press_send,eps_send,keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + 1.d0,1.d0,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress + + local_gam=log(press_send/local_K)/log(rho_send) + + 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 prim2conpolytypeM(GRHydro_polytrope_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(6,*)'local_K:',local_K + sc_send=local_K*dens_send + + write(*,*) 'P2C2PM test: getting the associated primitive variables.' + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_polytrope_handle,local_gam,dens_send,sx_send,sy_send,sz_send, & + sc_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_polytype diff --git a/src/GRHydro_ReadConformalData.F90 b/src/GRHydro_ReadConformalData.F90 index c492d3c..665b91e 100644 --- a/src/GRHydro_ReadConformalData.F90 +++ b/src/GRHydro_ReadConformalData.F90 @@ -43,18 +43,9 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS -!#include "EOS_Base.h" -#ifdef _EOS_BASE_INC_ -#undef _EOS_BASE_INC_ -#endif -#include "EOS_Base.inc" - - CCTK_INT :: i,j,k,handle,ierr -! CCTK_REAL :: eos_k, eos_gamma CCTK_REAL :: rho_min, det -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: poly_eoskey = 0 @@ -68,7 +59,6 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS) real*8 :: xye(1) = 0.0d0 poly_eoskey = EOS_Omni_GetHandle("2D_Polytrope") ! end EOS Omni vars -#endif ! only gxx has been read in; copy it into gyy and gzz as well gyy = gxx @@ -155,14 +145,11 @@ subroutine GRHydro_ReadConformalData(CCTK_ARGUMENTS) do i=1,cctk_lsh(1) do j=1,cctk_lsh(2) do k=1,cctk_lsh(3) -#ifdef USE_EOS_OMNI + call EOS_Omni_EpsFromPress(poly_eoskey,keytemp,rf_precision,& n,rho(i,j,k),xeps,xtemp,xye,& press(i,j,k),eps(i,j,k),keyerr,anyerr) -#else - eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_eos_handle,rho(i,j,k),press(i,j,k)) -#endif end do end do end do diff --git a/src/GRHydro_ShockTube.F90 b/src/GRHydro_ShockTube.F90 index 80d8968..baf891d 100644 --- a/src/GRHydro_ShockTube.F90 +++ b/src/GRHydro_ShockTube.F90 @@ -11,6 +11,7 @@ #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) @@ -124,8 +125,7 @@ subroutine GRHydro_shocktube(CCTK_ARGUMENTS) eps(i,j,k) = epsr end if - call SpatialDet(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) + 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 Prim2ConPoly(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90 index 1169fef..96b0e7c 100644 --- a/src/GRHydro_ShockTubeM.F90 +++ b/src/GRHydro_ShockTubeM.F90 @@ -54,6 +54,13 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + bvcxl = Bx_init + bvcyl = By_init + bvczl = Bz_init + bvcxr = Bx_init + bvcyr = By_init + bvczr = Bz_init + nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -110,17 +117,95 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) velzr = 0.d0 epsl = 1500.d0 epsr = 1.5d-2 - else + else if (CCTK_EQUALS(shock_case,"Balsara1")) 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 + bvcxl=0.5d0 + bvcxr=0.5d0 + bvcyl=1.d0 + bvcyr=-1.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 1.0/rhol + epsr = 0.1/rhor + else if (CCTK_EQUALS(shock_case,"Balsara2")) 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 + bvcxl=5.0d0 + bvcxr=5.0d0 + bvcyl=6.d0 + bvcyr=0.7d0 + bvczl=6.d0 + bvczr=0.7d0 + epsl = 1.5d0*30.0d0/rhol + epsr = 1.5d0*1.0d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara3")) 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 + bvcxl=10.0d0 + bvcxr=10.0d0 + bvcyl=7.d0 + bvcyr=0.7d0 + bvczl=7.d0 + bvczr=0.7d0 + epsl = 1.5d0*1000.0d0/rhol + epsr = 1.5d0*0.1d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara4")) then + rhol = 1.d0 + rhor = 1.d0 + velxl = 0.999d0 + velxr = -0.999d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + bvcxl=10.0d0 + bvcxr=10.0d0 + bvcyl=7.d0 + bvcyr=-7.d0 + bvczl=7.d0 + bvczr=-7.d0 + epsl = 1.5d0*0.1d0/rhol + epsr = 1.5d0*0.1d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara5")) then + rhol = 1.08d0 + rhor = 1.d0 + velxl = 0.4d0 + velxr = -0.45d0 + velyl = 0.3d0 + velyr = -0.2d0 + velzl = 0.2d0 + velzr = 0.2d0 + bvcxl=2.0d0 + bvcxr=2.0d0 + bvcyl=0.3d0 + bvcyr=-0.7d0 + bvczl=0.3d0 + bvczr=0.5d0 + epsl = 1.5d0*0.95d0/rhol + epsr = 1.5d0*1.0d0/rhor + 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 diff --git a/src/GRHydro_SimpleWave.F90 b/src/GRHydro_SimpleWave.F90 index 2991ff4..c13916b 100644 --- a/src/GRHydro_SimpleWave.F90 +++ b/src/GRHydro_SimpleWave.F90 @@ -46,17 +46,9 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS -#ifndef USE_EOS_OMNI -#ifdef _EOS_BASE_INC_ -#undef _EOS_BASE_INC_ -#endif -#include "EOS_Base.inc" -#endif - CCTK_INT :: i, j, k, nx, ny, nz CCTK_REAL :: dr, k1, k2, k3, k4, in_data, old_data, source_data, new_data, c_0, det, pi -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -68,7 +60,6 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS) real*8 :: xye(1) = 0.0d0 real*8 :: rf_precision = 1.0d-10 ! end EOS Omni vars -#endif call CCTK_INFO("Setting initial data for a simple wave as Anile Miller Motta") @@ -143,21 +134,19 @@ subroutine GRHydro_SimpleWave(CCTK_ARGUMENTS) ! atmosphere - if ( (rho(i,1,1) < GRHydro_rho_min).OR.(velx(i,1,1) < 0) ) then - rho(i,1,1) = rho_abs_min -! rho(i,1,1) = 1.0 !the value of rho_min for the initial data - eps(i,1,1) = rho_abs_min**(1.d0/3.d0) - velx(i,1,1) = 0.d0 - w_lorentz(i,1,1) = 1.d0 -#if USE_EOS_OMNI - xeps = 1.0d0 - call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,rf_precision,n,& - rho(i,1,1),xeps,xtemp,xye,press(i,1,1),keyerr,anyerr) -#else - press(i,1,1) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,1,1), 1.0d0) -#endif - ! polytrope only (initial data) - end if + if ( (rho(i,1,1) < GRHydro_rho_min).OR.(velx(i,1,1) < 0) ) then + rho(i,1,1) = rho_abs_min + ! rho(i,1,1) = 1.0 !the value of rho_min for the initial data + eps(i,1,1) = rho_abs_min**(1.d0/3.d0) + velx(i,1,1) = 0.d0 + w_lorentz(i,1,1) = 1.d0 + + xeps = 1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,rf_precision,n,& + rho(i,1,1),xeps,xtemp,xye,press(i,1,1),keyerr,anyerr) + + ! polytrope only (initial data) + end if ! write(*,*) 'p',i, x(i,1,1), rho(i,1,1)**(4.d0/3.d0)/press(i,1,1) diff --git a/src/make.code.defn b/src/make.code.defn index dc78381..09a09f2 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -13,7 +13,8 @@ SRCS = GRHydro_C2P2C.F90 \ GRHydro_SimpleWave.F90 \ CheckParam.c\ GRHydro_P2C2P.F90 \ - GRHydro_P2C2PM.F90 + GRHydro_P2C2PM.F90 \ + GRHydro_P2C2PM_polytype.F90 # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3