From 86a32a01b174be105ec341c0a055d9e059f3db3f Mon Sep 17 00:00:00 2001 From: rhaas Date: Fri, 11 Jan 2013 15:04:08 +0000 Subject: GRHydro_InitData: initialize c2p2c values as parameters and update routine to test inversions in the gr context. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@192 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- param.ccl | 91 +++++++++++++++++++++++++++++++++++++++++++++++++ src/GRHydro_C2P2CM.F90 | 92 +++++++++++++++++++++++++++++--------------------- 2 files changed, 145 insertions(+), 38 deletions(-) diff --git a/param.ccl b/param.ccl index 67e51e1..bb9dad6 100644 --- a/param.ccl +++ b/param.ccl @@ -126,6 +126,8 @@ BOOLEAN attenuate_atmosphere "Attenuate the velocity in the atmosphere" { } "no" +# Initial magnetic field used by different tests + REAL Bx_init "Initial B-field in the x-dir" { *:* :: "Anything" @@ -141,6 +143,95 @@ REAL Bz_init "Initial B-field in the z-dir" *:* :: "Anything" } 0.0 +# Initial primitive values: +REAL rho_init "Initial rest mass density" +{ + (0:* :: "Anything positive." +} 1.0d-6 + +REAL velx_init "Initial x velocity" +{ + *:* :: "Anything." +} 1.0d-1 + +REAL vely_init "Initial y velocity" +{ + *:* :: "Anything." +} 1.0d-1 + +REAL velz_init "Initial z velocity" +{ + *:* :: "Anything." +} 1.0d-1 + +REAL eps_init "Initial specific internal energy" +{ + (0:* :: "Anything positive." +} 1.0d-6 + +REAL press_init "Initial pressure" +{ + (0:* :: "Anything positive." +} 6.666666666666667d-7 + +# Initial conservative values used by c2p tests: +REAL dens_init "Initial conserved mass density" +{ + (0:* :: "Anything positive." +} 1.29047362 + +REAL sx_init "Initial x component of conserved momentum density" +{ + *:* :: "Anything." +} 0.166666658 + +REAL sy_init "Initial y component of conserved momentum density" +{ + *:* :: "Anything." +} 0.166666658 + +REAL sz_init "Initial z component of conserved momentum density" +{ + *:* :: "Anything." +} 0.166666658 + +REAL tau_init "Initial conserved total energy density" +{ + (0:* :: "Anything positive." +} 0.484123939 + +# Initial values for 3-metric components. Default to euclidian 3-metric +REAL gxx_init "Initial xx metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 1.0 + +REAL gxy_init "Initial xy metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 0.0 + +REAL gxz_init "Initial xz metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 0.0 + +REAL gyy_init "Initial yy metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 1.0 + +REAL gyz_init "Initial yz metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 0.0 + +REAL gzz_init "Initial zz metric componenent" +{ + *:* :: "Anything, but be carefull to set a positive definite 3-metric!" +} 1.0 + + KEYWORD monopole_type "Which kind of monopole?" { "Point" :: "Single point with Bx /= 0" diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90 index c167504..9394005 100644 --- a/src/GRHydro_C2P2CM.F90 +++ b/src/GRHydro_C2P2CM.F90 @@ -11,6 +11,8 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" + /*@@ @routine c2p2cM @@ -37,13 +39,13 @@ subroutine c2p2cM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: det + CCTK_REAL :: det, sdet, invdet 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 :: bconsx_send, bconsy_send, bconsz_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 :: press_send,w_lorentz_send CCTK_REAL :: bvcx_send, bvcy_send, bvcz_send, b2_send CCTK_REAL :: C2P_failed CCTK_INT :: epsnegative @@ -58,50 +60,64 @@ subroutine c2p2cM(CCTK_ARGUMENTS) call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS") - 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 + gxx_send = gxx_init + gxy_send = gxy_init + gxz_send = gxz_init + gyy_send = gyy_init + gyz_send = gyz_init + gzz_send = gzz_init - det = 1.0d0 - - uxx = 1.0d0 - uyy = 1.0d0 - uzz = 1.0d0 - uxy = 0.0d0 - uxz = 0.0d0 - uyz = 0.0d0 + det = SPATIAL_DETERMINANT(gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send) + sdet = sqrt(det) + invdet = 1.d0 / det + uxx = (-gyz_send**2 + gyy_send*gzz_send)*invdet + uxy = (gxz_send*gyz_send - gxy_send*gzz_send)*invdet + uyy = (-gxz_send**2 + gxx_send*gzz_send)*invdet + uxz = (-gxz_send*gyy_send + gxy_send*gyz_send)*invdet + uyz = (gxy_send*gxz_send - gxx_send*gyz_send)*invdet + uzz = (-gxy_send**2 + gxx_send*gyy_send)*invdet ! Initialize the velocity as GRHydro_Con2PrimM_pt may use these ! values as initial guess for its Newton-Raphson procedure. - velx_send = 0.1d0 - vely_send = 0.1d0 - velz_send = 0.1d0 - - dens_send = 1.29047362d0 - sx_send = 0.166666658d0 - sy_send = 0.166666658d0 - sz_send = 0.166666658d0 - tau_send = 0.484123939d0 +! velx_send = 0.1d0 +! vely_send = 0.1d0 +! velz_send = 0.1d0 +! +! dens_send = 1.29047362d0 +! sx_send = 0.166666658d0 +! sy_send = 0.166666658d0 +! sz_send = 0.166666658d0 +! tau_send = 0.484123939d0 + + velx_send = velx_init + vely_send = vely_init + velz_send = velz_init + + dens_send = dens_init + sx_send = sx_init + sy_send = sy_init + sz_send = sz_init + tau_send = tau_init bvcx_send = Bx_init bvcy_send = By_init bvcz_send = Bz_init - bconsx_send = Bx_init - bconsy_send = By_init - bconsz_send = Bz_init - - rho_send=1.0d-6 - eps_send = 1.0d-6 - press_send = 6.666666666666667d-7 - w_lorentz_send = 1.0d0 + bconsx_send = sdet*Bx_init + bconsy_send = sdet*By_init + bconsz_send = sdet*Bz_init + + rho_send = rho_init + eps_send = eps_init + press_send = press_init + w_lorentz_send = sqrt(1.0d0-(gxx_send*velx_send**2+gyy_send*vely_send**2+ & + gzz_send*velz_send**2 & + +2.0d0*(gxy_send*velx_send*vely_send+ & + gxz_send*velx_send*velz_send+ & + gyz_send*vely_send*velz_send & + ) & + ) & + ) + w_lorentz_send = 1.0d0/w_lorentz_send epsnegative = 0 -- cgit v1.2.3