From e14b402798458759073298a0c31a7bd5f18e7774 Mon Sep 17 00:00:00 2001 From: bmundim Date: Fri, 31 Dec 2010 19:07:40 +0000 Subject: RIT MHD development: Add #include "cctk_Functions.h" and DECLARE_CCTK_FUNCTIONS to GRHydro_C2P2CM.F90. Solve several issues of rank mismatch in the EOS Omni interface. Initialize the velocities at GRHydro_C2P2CM.F90 since GRHydro_Con2PrimM_pt.c uses past values as initial guess for its Newton-Raphson procedure. These unitialized values were causing C2P2CM to fail from time to time, since it would access memory positions with very large numerical values, i.e. v>>1. C2P2CM test always passes now. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@118 ac85fae7-cede-4708-beff-ae01c7fa1c26 --- doc/documentation.tex | 2 +- interface.ccl | 17 ++++++++++++++ src/GRHydro_C2P2CM.F90 | 58 +++++++++++++++++++++++++++++----------------- src/GRHydro_ShockTubeM.F90 | 26 ++++++++++----------- 4 files changed, 68 insertions(+), 35 deletions(-) diff --git a/doc/documentation.tex b/doc/documentation.tex index 6383fd4..bee5e35 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -118,7 +118,7 @@ For testing purposes, this routine sets all the points to the values of the atmo \subsection{Simple wave} \label{sec:simple-wave} -This routine stes initial data for a simple wave with sinusoidal initial function for the velocity, +This routine testes initial data for a simple wave with sinusoidal initial function for the velocity, as described in Anile, Miller, Motta, {\it Formation and damping of relativistic strong shocks},Phys. Fluids {\bf 26}, 1450 (1983). diff --git a/interface.ccl b/interface.ccl index b8f7af7..4f3f7b5 100644 --- a/interface.ccl +++ b/interface.ccl @@ -79,12 +79,29 @@ SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) +void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ + CCTK_REAL INOUT tau, CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ + CCTK_REAL INOUT w_lorentz, \ + CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ + CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ + CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ + CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ + CCTK_REAL INOUT det, \ + CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + CCTK_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + USES FUNCTION SpatialDet USES FUNCTION Prim2ConPoly USES FUNCTION Prim2ConGen USES FUNCTION Prim2ConPolyM USES FUNCTION Prim2ConGenM USES FUNCTION Con2PrimPoly +USES FUNCTION Con2PrimGenM CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name) USES FUNCTION EOS_Omni_GetHandle diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90 index 3bd61b6..c382037 100644 --- a/src/GRHydro_C2P2CM.F90 +++ b/src/GRHydro_C2P2CM.F90 @@ -10,6 +10,7 @@ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" +#include "cctk_Functions.h" /*@@ @routine c2p2cM @@ -29,10 +30,12 @@ subroutine c2p2cM(CCTK_ARGUMENTS) + implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS integer didit,i,j,k,nx,ny,nz CCTK_REAL det @@ -42,19 +45,18 @@ 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, local_gam - CCTK_INT C2P_failed - logical epsnegative + CCTK_REAL :: pmin(1), epsmin(1), local_gam(1) + CCTK_REAL C2P_failed + CCTK_INT epsnegative + + CCTK_REAL :: rhoval(1), epsval(1) + CCTK_REAL :: gamval ! 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 + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1) + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0 ! end EOS Omni vars call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS") @@ -83,6 +85,12 @@ subroutine c2p2cM(CCTK_ARGUMENTS) uxy = 0.0d0 uxz = 0.0d0 uyz = 0.0d0 + +! 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 @@ -98,22 +106,27 @@ subroutine c2p2cM(CCTK_ARGUMENTS) press_send = 6.666666666666667d-7 w_lorentz_send = 1.0d0 - epsnegative = .false. +! epsnegative = .false. + epsnegative = 0 - GRHydro_rho_min = 1.e-10 + GRHydro_rho_min = 1.0d-10 + rhoval(1) = GRHydro_rho_min + epsval(1) = 1.0d0 call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) + rhoval,epsval,xtemp,xye,pmin,keyerr,anyerr) call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) + rhoval,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) - local_gam=0.0 + local_gam(1)=0.0d0 + rhoval(1) = 1.0d0 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 + rhoval,epsval,xtemp,xye,local_gam,keyerr,anyerr) + local_gam(1) = local_gam(1)+1.0 + gamval = local_gam(1) - C2P_failed = 0 + C2P_failed = 0.0d0 write(*,*) 'C2P2CM test: initial values.' write(*,*) ' conservative variables: ' @@ -129,7 +142,8 @@ subroutine c2p2cM(CCTK_ARGUMENTS) write(*,*) ' Bvecz : ',bvcz_send write(*,*) 'C2P2CM test: getting the associated primitive variables.' - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,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, & + call Con2PrimGenM(GRHydro_eos_handle,gamval,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,& @@ -149,9 +163,11 @@ subroutine c2p2cM(CCTK_ARGUMENTS) write(*,*) ' Bvecx : ',bvcx_send write(*,*) ' Bvecy : ',bvcy_send write(*,*) ' Bvecz : ',bvcz_send + write(*,*) ' C2P_failed : ',C2P_failed 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, & +! call prim2conM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, & + call Prim2ConGenM(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) diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90 index 96b0e7c..186cf69 100644 --- a/src/GRHydro_ShockTubeM.F90 +++ b/src/GRHydro_ShockTubeM.F90 @@ -118,22 +118,22 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) epsl = 1500.d0 epsr = 1.5d-2 else if (CCTK_EQUALS(shock_case,"Balsara1")) then - rhol = 1.d0 + rhol = 1.0d0 rhor = 0.125d0 - velxl = 0.d0 - velxr = 0.d0 - velyl = 0.d0 - velyr = 0.d0 - velzl = 0.d0 - velzr = 0.d0 + velxl = 0.0d0 + velxr = 0.0d0 + velyl = 0.0d0 + velyr = 0.0d0 + velzl = 0.0d0 + velzr = 0.0d0 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 + bvcyl=1.0d0 + bvcyr=-1.0d0 + bvczl=0.0d0 + bvczr=0.0d0 + epsl = 1.0d0/rhol + epsr = 0.1d0/rhor else if (CCTK_EQUALS(shock_case,"Balsara2")) then rhol = 1.d0 rhor = 1.d0 -- cgit v1.2.3