Commits March 2013

commits@lists.einsteintoolkit.org

9 participants
84 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 204)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ interface.ccl, param.ccl /trunk/src/ GRHydro_C2P2CM.F90 Log: Finish the implementation of con2prim using the entropy equation. * Tested successfully pointwisely. It still needs to be tested on evolution. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_C2P2CM.F90 Delta lines: +87 -2 =================================================================== --- trunk/src/GRHydro_C2P2CM.F90 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/src/GRHydro_C2P2CM.F90 2013-03-28 01:49:03 UTC (rev 204) @@ -44,10 +44,12 @@ 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 :: entropy_send, entropycons_send; CCTK_REAL :: rho_send,velx_send,vely_send,velz_send,eps_send CCTK_REAL :: press_send,w_lorentz_send CCTK_REAL :: bvcx_send, bvcy_send, bvcz_send, b2_send CCTK_REAL :: C2P_failed + CCTK_REAL :: BdotB, Bdotv, b2, beta_mag CCTK_INT :: epsnegative ! begin EOS Omni vars @@ -119,6 +121,7 @@ ) w_lorentz_send = 1.0d0/w_lorentz_send + epsnegative = 0 xrho = 1.0d-10 @@ -136,8 +139,46 @@ xrho,epsval,xtemp,xye,local_gam,keyerr,anyerr) local_gam = local_gam + 1.0 + if(use_c2p_with_entropy_eqn.eq.1)then + entropy_send = ( local_gam(1) - 1.0d0 ) * eps_send * & + rho_send**(2.0d0 - local_gam(1)) + entropycons_send = sdet * w_lorentz_send * entropy_send + endif + C2P_failed = 0.0d0 + write(*,*) 'C2P2CM test: metric values.' + write(*,*) ' gxx: ', gxx_send + write(*,*) ' gxy: ', gxy_send + write(*,*) ' gxz: ', gxz_send + write(*,*) ' gyy: ', gyy_send + write(*,*) ' gyz: ', gyz_send + write(*,*) ' gzz: ', gzz_send + write(*,*) ' uxx: ', uxx + write(*,*) ' uxy: ', uxy + write(*,*) ' uxz: ', uxz + write(*,*) ' uyy: ', uyy + write(*,*) ' uyz: ', uyz + write(*,*) ' uzz: ', uzz + write(*,*) ' det: ', det + + write(*,*) 'C2P2CM test: initial primitive guess values.' + 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 + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropy : ',entropy_send + endif + write(*,*) ' W : ',w_lorentz_send + write(*,*) ' Bvecx : ',bvcx_send + write(*,*) ' Bvecy : ',bvcy_send + write(*,*) ' Bvecz : ',bvcz_send + write(*,*) ' C2P_failed : ',C2P_failed + write(*,*) 'C2P2CM test: initial values.' write(*,*) ' conservative variables: ' write(*,*) ' dens: ',dens_send @@ -148,19 +189,57 @@ write(*,*) ' Bconsx : ',bconsx_send write(*,*) ' Bconsy : ',bconsy_send write(*,*) ' Bconsz : ',bconsz_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropycons : ',entropycons_send + endif write(*,*) ' eps : ',eps_send write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send write(*,*) ' Bvecy : ',bvcy_send write(*,*) ' Bvecz : ',bvcz_send + +! Is this point magnetically dominated? + BdotB = gxx_send*bvcx_send**2+gyy_send*bvcy_send**2+gzz_send*bvcz_send**2& + +2.0*(gxy_send*bvcx_send*bvcy_send+gxz_send*bvcx_send*bvcy_send+ & + gyz_send*bvcy_send*bvcz_send) + Bdotv = gxx_send*bvcx_send*velx_send+gyy_send*bvcy_send*vely_send & + +gzz_send*bvcz_send*velz_send & + +gxy_send*(bvcx_send*vely_send+bvcy_send*velx_send) & + +gxz_send*(bvcx_send*velz_send+bvcz_send*velx_send) & + +gyz_send*(bvcy_send*velz_send+bvcz_send*vely_send) + + b2 = BdotB/w_lorentz_send**2 + Bdotv**2 + + beta_mag = 2.0*press_send/b2; + + write(*,*) ' BdotB : ',BdotB + write(*,*) ' Bdotv : ',Bdotv + write(*,*) ' b2 : ',b2 + write(*,*) 'beta_mag : ',beta_mag - write(*,*) 'C2P2CM test: getting the associated primitive variables.' - call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + if(use_c2p_with_entropy_eqn.eq.0)then + write(*,*) 'C2P2CM test: getting the associated primitive variables.' + call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & tau_send,bconsx_send,bconsy_send,bconsz_send,xtemp(1),xye(1),rho_send,velx_send,vely_send,velz_send, & eps_send,press_send,bvcx_send,bvcy_send,bvcz_send,b2_send,w_lorentz_send, & gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,& uxx,uxy,uxz,uyy,uyz,uzz,det,& epsnegative,C2P_failed) + else + write(*,*) 'C2P2CMee test: getting the associated primitive variables.' + call Con2PrimGenMee(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + tau_send,bconsx_send,bconsy_send,bconsz_send,entropycons_send,xtemp(1),xye(1),rho_send,velx_send,vely_send,velz_send, & + eps_send,press_send,bvcx_send,bvcy_send,bvcz_send,b2_send,w_lorentz_send, & + gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,& + uxx,uxy,uxz,uyy,uyz,uzz,det,& + epsnegative,C2P_failed) + endif + + if(use_c2p_with_entropy_eqn.eq.1)then + entropy_send = ( local_gam(1) - 1.0d0 ) * eps_send * & + rho_send**(2.0d0 - local_gam(1)) + entropycons_send = sdet * w_lorentz_send * entropy_send + endif write(*,*) 'C2P2CM test: the primitive variables are' write(*,*) ' primitive variables: ' @@ -170,6 +249,9 @@ write(*,*) ' velz : ',velz_send write(*,*) ' press : ',press_send write(*,*) ' eps : ',eps_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropy : ',entropy_send + endif write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send write(*,*) ' Bvecy : ',bvcy_send @@ -191,6 +273,9 @@ write(*,*) ' Bconsx : ',bconsx_send write(*,*) ' Bconsy : ',bconsy_send write(*,*) ' Bconsz : ',bconsz_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropycons : ',entropycons_send + endif write(*,*) ' eps : ',eps_send write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +22 -0 =================================================================== --- trunk/interface.ccl 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/interface.ccl 2013-03-28 01:49:03 UTC (rev 204) @@ -131,6 +131,27 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ + CCTK_REAL INOUT prec, 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 Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT entropycons, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, 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 Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + 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_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + void FUNCTION Con2PrimPolyM(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 sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ @@ -157,6 +178,7 @@ USES FUNCTION Con2PrimPoly USES FUNCTION Con2PrimGen USES FUNCTION Con2PrimGenM +USES FUNCTION Con2PrimGenMee USES FUNCTION Con2PrimPolyM CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name) File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/param.ccl 2013-03-28 01:49:03 UTC (rev 204) @@ -188,6 +188,11 @@ } 6.666666666666667d-7 # Initial conservative values used by c2p tests: + +CCTK_BOOLEAN use_c2p_with_entropy_eqn "Use the con2prim routine that uses the entropy equation instead of the energy equation" +{ +} no + REAL dens_init "Initial conserved mass density" { (0:* :: "Anything positive."

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 203)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ param.ccl Log: GRHydro_InitData: Allow n_p=0 again for poloidal magnetic field From: Tanja Bode <tanja.bode(a)physics.gatech.edu> File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -1 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:00 UTC (rev 202) +++ trunk/param.ccl 2013-03-28 01:49:02 UTC (rev 203) @@ -424,7 +424,7 @@ CCTK_INT poloidal_n_p "Vector potential strength" { - (0:* :: "Any positive integer." + 0:* :: "Any positive integer." } 3 CCTK_REAL poloidal_P_cut "Pressure used to confine the B field inside a star"

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 202)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: damp hydro quantities to atmosphere inside of rKS < 1 M From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +13 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:59 UTC (rev 201) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:00 UTC (rev 202) @@ -1238,6 +1238,19 @@ Bvecz(i) = bondi_bmag*SQR(M)*z[i]/sqrt(det)/CUBE(r[i]); } + // damp everything down to atmosphere inside of 1M + if(rspher < M) { + const CCTK_REAL smooth = 0.5*(1+tanh(tan(M_PI*(rspher/M-0.5)))); + rho[i] *= smooth; + velx(i) *= smooth; + vely(i) *= smooth; + velz(i) *= smooth; + eps[i] *= smooth; + Bvecx(i) *= smooth; + Bvecy(i) *= smooth; + Bvecz(i) *= smooth; + } + Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], gxz[i],gyy[i],gyz[i],gzz[i], det, &dens[i],&sx(i),&sy(i),&sz(i),

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 201)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: compute w_lorentz from uT and alp instead of vel later on prim2con will of course compute the grid function using vel From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +4 -12 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:57 UTC (rev 200) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:59 UTC (rev 201) @@ -1225,21 +1225,13 @@ vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * inv_alp; velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * inv_alp; - w_lorentz[i] = 1./sqrt(1-(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ - 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i))) ); - if(isnan(w_lorentz[i])) { - CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, - "nan in w_lorentz at (%g,%g,%g): v2 = %g", - x[i],y[i],z[i],(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ - 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i)))); - } - SpatialDet(gxx[i],gxy[i],gxz[i],gyy[i],gyz[i],gzz[i],&det); if(bvec_method == BVEC_TRANSFORM) { - Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; - Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; - Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; + const CCTK_REAL w = ucon[TT] * alp[i]; + Bvecx(i) = w*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; + Bvecy(i) = w*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; + Bvecz(i) = w*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; } else { Bvecx(i) = bondi_bmag*SQR(M)*x[i]/sqrt(det)/CUBE(r[i]); Bvecy(i) = bondi_bmag*SQR(M)*y[i]/sqrt(det)/CUBE(r[i]);

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 200)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: remove rho_guess from find_bondi_solution was used to call the newton solver twice with the very same arguments. Not sure what the idea was. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +1 -12 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:56 UTC (rev 199) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:57 UTC (rev 200) @@ -647,7 +647,7 @@ const int ntries = 10000; int itry; - CCTK_REAL rhotmp, rho_guess; + CCTK_REAL rhotmp; CCTK_REAL dr,ur; @@ -662,17 +662,12 @@ } else { // rhotmp = (sqrt(Qdot) - 1.) * (gamma_eos - 1.) / ( gamma_eos * K ); - // rho_guess = pow( rhotmp , (1./(gamma_eos - 1.)) ); if(r < rs) { ur = pow(r,-0.5) ; } else { ur = 0.5*pow(r,-1.5) ; } *rho = Mdot / (4.*M_PI * r * r * ur); } } - // safe guess value for multiple tries - rho_guess = *rho; - - // set global variables needed by residual function: r_sol = r ; @@ -681,12 +676,6 @@ retval = gnr_bondi( rho, NEWT_DIM_B, bondi_resid); - // first try guess if failure - if( retval ) { - *rho = rho_guess; - retval = gnr_bondi( rho, NEWT_DIM_B, bondi_resid); - } - // If we were unsure about the guess and solver fails, then creep from known solution to desired point: if( retval ) {

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 199)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_InitDat: clear psidc and c2p_failed in ID routine this is mostly important when we freeze the evolution afterwards From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +4 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:55 UTC (rev 198) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:56 UTC (rev 199) @@ -1267,6 +1267,10 @@ eps[i],&press[i], Bvecx(i),Bvecy(i),Bvecz(i), &w_lorentz[i]); + + GRHydro_C2P_failed[i] = 0.; + if(psidc) + psidc[i] = 0.; }

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 198)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_BondiM.c Log: GRHydro_InitData: set Bvec according to UIUC formula From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +46 -26 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:55 UTC (rev 198) @@ -901,6 +901,8 @@ -- driver routine for the Bondi solution; ***********************************************************************************/ +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max); + void GRHydro_BondiM(CCTK_ARGUMENTS) { GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); @@ -914,17 +916,18 @@ } -void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - CCTK_REAL det; + CCTK_REAL det, inv_alp; CCTK_REAL rhotmp, utmp, vtmp, rspher, xpos[NDIM], x_spher[NDIM], ucon[NDIM], bcon[NDIM]; int retval; CCTK_REAL *r_bondi, *logr_bondi, *rho_bondi, *u_bondi, *v_bondi; CCTK_REAL dlogr,logrmin; CCTK_REAL rmin_bondi,rmax_bondi; -// CCTK_INT GRHydro_reflevel=0; + int bvec_method; + enum {BVEC_DIRECT, BVEC_TRANSFORM}; int nbondi; @@ -946,14 +949,26 @@ if( CCTK_EQUALS(bondi_coordinates,"Boyer-Lindquist") ) { coord_type = COORD_BOYERLINDQUIST ; - } - if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { coord_type = COORD_KERRSCHILD ; - } - if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { coord_type = COORD_ISOTROPIC ; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown coordinate type '%s'", bondi_coordinates); } + if( CCTK_EQUALS(bondi_Bvec_method, "direct")) { + bvec_method = BVEC_DIRECT; + } else if( CCTK_EQUALS(bondi_Bvec_method, "transform")) { + bvec_method = BVEC_TRANSFORM; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown bvec setup method '%s'", bondi_Bvec_method); + } + + + /* xyz location of the black hole : */ i = 0; x_cen = bh_bondi_pos_x[i] ; @@ -1215,11 +1230,11 @@ } } - det = 1./alp[i]; /* temp var */ + inv_alp = 1./alp[i]; - velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * det; - vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * det; - velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * det; + velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * inv_alp; + vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * inv_alp; + velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * inv_alp; w_lorentz[i] = 1./sqrt(1-(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i))) ); @@ -1230,24 +1245,29 @@ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i)))); } - - Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; - Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; - Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; - SpatialDet(gxx[i],gxy[i],gxz[i],gyy[i],gyz[i],gzz[i],&det); - Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], - gxz[i],gyy[i],gyz[i],gzz[i], - det, &dens[i],&sx(i),&sy(i),&sz(i), - &tau[i], - &Bconsx(i),&Bconsy(i),&Bconsz(i), - rho[i], - velx(i),vely(i),velz(i), - eps[i],&press[i], - Bvecx(i),Bvecy(i),Bvecz(i), - &w_lorentz[i]); + if(bvec_method == BVEC_TRANSFORM) { + Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; + Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; + Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; + } else { + Bvecx(i) = bondi_bmag*SQR(M)*x[i]/sqrt(det)/CUBE(r[i]); + Bvecy(i) = bondi_bmag*SQR(M)*y[i]/sqrt(det)/CUBE(r[i]); + Bvecz(i) = bondi_bmag*SQR(M)*z[i]/sqrt(det)/CUBE(r[i]); + } + Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], + gxz[i],gyy[i],gyz[i],gzz[i], + det, &dens[i],&sx(i),&sy(i),&sz(i), + &tau[i], + &Bconsx(i),&Bconsy(i),&Bconsz(i), + rho[i], + velx(i),vely(i),velz(i), + eps[i],&press[i], + Bvecx(i),Bvecy(i),Bvecz(i), + &w_lorentz[i]); + } free(r_bondi); free(logr_bondi); free(rho_bondi); free(u_bondi); free(v_bondi); Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +6 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/param.ccl 2013-03-28 01:48:55 UTC (rev 198) @@ -395,6 +395,12 @@ (0:* :: "positive" } 1.0 +CCTK_KEYWORD bondi_Bvec_method "how to compute the magnetic field vector" +{ + "direct" :: "directly from Cartesian metric" + "transform" :: "transform Schwarzschild solution to Kerr Schild" +} "direct" + CCTK_BOOLEAN bondi_evolve_only_annulus "reset to initial data outside of bondi_freeze_inner_radius and bondi_freeze_outer_radius" { } "no"

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 197)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/ param.ccl, schedule.ccl /trunk/src/ GRHydro_BondiM.c Log: GRHydro: add parameters to set Bondi solution only in certain range. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +17 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:53 UTC (rev 197) @@ -84,6 +84,7 @@ //some math helpers #define SQR(x) ((x)*(x)) +#define CUBE(x) ((x)*(x)*(x)) static CCTK_REAL Mdot, rs, vs_sq, vs, cs_sq, cs, rhos, hs, K, Qdot, gamma_eos, r_sol; @@ -902,6 +903,19 @@ ***********************************************************************************/ void GRHydro_BondiM(CCTK_ARGUMENTS) { + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); +} + +void GRHydro_BondiM_Range(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS; + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, bondi_freeze_inner_radius, + bondi_freeze_outer_radius); +} + + +void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +{ DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; CCTK_REAL det; @@ -1134,6 +1148,9 @@ rspher = x_spher[RR]; + if(rspher > range_min && rspher < range_max) + continue; + /* Find nearest point in the Bondi solution : */ j = (int) ( 0.5 + (log10(rspher) - logrmin) / dlogr ) ; Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +14 -0 =================================================================== --- trunk/param.ccl 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/param.ccl 2013-03-28 01:48:53 UTC (rev 197) @@ -395,6 +395,20 @@ (0:* :: "positive" } 1.0 +CCTK_BOOLEAN bondi_evolve_only_annulus "reset to initial data outside of bondi_freeze_inner_radius and bondi_freeze_outer_radius" +{ +} "no" + +CCTK_REAL bondi_freeze_inner_radius "reset to initial at radii below this" +{ + *:* :: "any value" +} -1. + +CCTK_REAL bondi_freeze_outer_radius "reset to initial at radii above this" +{ + *:* :: "any value" +} 1e300 + # For Poloidal Magnetic field test: CCTK_REAL poloidal_A_b "Vector potential strength" File [modified]: schedule.ccl Delta lines: +7 -0 =================================================================== --- trunk/schedule.ccl 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/schedule.ccl 2013-03-28 01:48:53 UTC (rev 197) @@ -250,6 +250,13 @@ } "setup GRHydro vars for the magnetized Bondi solution" } +if(bondi_evolve_only_annulus) { +SCHEDULE GRHydro_BondiM_Range IN HydroBase_Con2Prim BEFORE Con2Prim +{ + LANG: C +} "force analytic solution outside anulus" +} + if (CCTK_EQUALS(initial_Bvec, "poloidalmagfield") || CCTK_EQUALS(initial_Avec, "poloidalmagfield")) { # SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_INITIAL AFTER IN HydroBase_Initial AFTER rnsid_init AFTER TOV_Initial_Data after CCCC_StarMapper_InitialData

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 501)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Modified: /trunk/src/ GRHydro_TVDReconstruct_drv.F90 Log: GRHydro: support reconstruct_Wv in TVD reconstructor From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_TVDReconstruct_drv.F90 Delta lines: +107 -8 =================================================================== --- trunk/src/GRHydro_TVDReconstruct_drv.F90 2013-03-28 01:47:14 UTC (rev 500) +++ trunk/src/GRHydro_TVDReconstruct_drv.F90 2013-03-28 01:47:16 UTC (rev 501) @@ -50,6 +50,8 @@ ! save memory when MP is not used ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz TARGET vel, lvel TARGET Bvec, lBvec @@ -67,14 +69,32 @@ CCTK_INT :: ierr + ! variables used when reconstruct_Wv == true + CCTK_REAL :: inv_w ! 1/Lorentz factor + CCTK_REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: wvel ! vel*w_lorentz + CCTK_REAL :: agxx,agxy,agxz,agyy,agyz,agzz ! metric components + ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc vup => lvel Bprim => lBvec else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz vup => vel Bprim => Bvec end if @@ -83,6 +103,29 @@ if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with trivial_rp") end if + + !!$ reconstruct w^i = w_lorentz*v^i to ensure slower than light + !!$ speeds? + if (reconstruct_Wv .ne. 0) then + ! all velocity like quantities are now w_lorentz*velocity. We will convert + ! back to ordinary velocity at the end, using w_lorentz = sqrt(1 + g_{ij} + ! w^i w^j). + allocate(wvel(cctk_ash(1),cctk_ash(2),cctk_ash(3),3),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with wvel") + end if + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + wvel(i,j,k,1) = vup(i,j,k,1) * w_lorentz(i,j,k) + wvel(i,j,k,2) = vup(i,j,k,2) * w_lorentz(i,j,k) + wvel(i,j,k,3) = vup(i,j,k,3) * w_lorentz(i,j,k) + end do + end do + end do + !$OMP END PARALLEL DO + end if call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & @@ -184,14 +227,69 @@ call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) + if (reconstruct_Wv.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + ! divide out the Loretnz factor obtained from w_lorentz = + ! sqrt(1+g_{ij} w^i w^j) for both the + ! plus and minus quantities this should by construction ensure + ! that any Lorentz factor calculated from them later on is + ! physical (ie. > 1.d0) + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i,j,k,inv_w,agxx,agxy,agxz,agyy,agyz,agzz) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints*(1-xoffset) + agxx = 0.5d0*( g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset) ) + agxy = 0.5d0*( g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset) ) + agxz = 0.5d0*( g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset) ) + agyy = 0.5d0*( g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset) ) + agyz = 0.5d0*( g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset) ) + agzz = 0.5d0*( g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset) ) + inv_w = 1d0/sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) & + + agyy*velyminus(i,j,k)*velyminus(i,j,k) & + + agzz*velzminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) & + + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) ) + velxminus(i,j,k) = velxminus(i,j,k)*inv_w + velyminus(i,j,k) = velyminus(i,j,k)*inv_w + velzminus(i,j,k) = velzminus(i,j,k)*inv_w + + agxx = 0.5d0*( g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset) ) + agxy = 0.5d0*( g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset) ) + agxz = 0.5d0*( g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset) ) + agyy = 0.5d0*( g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset) ) + agyz = 0.5d0*( g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset) ) + agzz = 0.5d0*( g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset) ) + inv_w = 1d0/sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) & + + agyy*velyplus(i,j,k)*velyplus(i,j,k) & + + agzz*velzplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) & + + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) ) + velxplus(i,j,k) = velxplus(i,j,k)*inv_w + velyplus(i,j,k) = velyplus(i,j,k)*inv_w + velzplus(i,j,k) = velzplus(i,j,k)*inv_w + end do + end do + end do + !$OMP END PARALLEL DO + else + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + end if + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) if(evolve_mhd.ne.0) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & Bprim(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) @@ -270,6 +368,7 @@ !!$ TVD ends. deallocate(trivial_rp) + ! Fortran 90 will automatically deallocate wvel for us end subroutine GRHydro_TVDReconstruct_drv

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 500)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Added: /trunk/src/ GRHydro_Con2PrimM_ptee.c Modified: /trunk/ interface.ccl /trunk/src/ GRHydro_Con2PrimM.F90, make.code.defn Log: Finish the implementation of con2prim using the entropy equation. * Tested successfully pointwisely. It still needs to be tested on evolution. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +34 -10 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:14 UTC (rev 500) @@ -334,7 +334,11 @@ Bvecz_tmp = Bprim(i,j,k,3) keytemp = 0 - + !Watch out for the values returned to b2. Here b2 is the Bprim^2 + !while inside the point-wise con2prim routines it is the square + !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten + !in this routine, but we may need to find a better notation + !avoid future confusions. call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & @@ -347,7 +351,13 @@ epsnegative,GRHydro_C2P_failed(i,j,k)) if(evolve_entropy.ne.0) then - entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + if(GRHydro_C2P_failed(i,j,k).ne.0) then + !Use previous time step for rho: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + else + !Use the current correct value of rho returned by con2prim: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho_tmp + endif endif if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -430,14 +440,28 @@ Bvecy_tmp = Bprim(i,j,k,2) Bvecz_tmp = Bprim(i,j,k,3) - call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & - dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & - Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& - velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& - Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & - epsnegative,GRHydro_C2P_failed(i,j,k)) + if(evolve_entropy.ne.0) then + call GRHydro_Con2PrimM_ptee(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), & + entropycons(i,j,k), xye(1), & + xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + else + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end if rho(i,j,k) = rho_tmp press(i,j,k) = press_tmp File [added]: GRHydro_Con2PrimM_ptee.c Delta lines: +592 -0 =================================================================== --- trunk/src/GRHydro_Con2PrimM_ptee.c (rev 0) +++ trunk/src/GRHydro_Con2PrimM_ptee.c 2013-03-28 01:47:14 UTC (rev 500) @@ -0,0 +1,592 @@ +/*********************************************************************************** + Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, + and Luca Del Zanna. + + PVS_GRMHD + + This file was derived from PVS_GRMHD. The authors of PVS_GRMHD include + Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna. + PVS_GRMHD is available under the GPL from: + http://rainman.astro.uiuc.edu/codelib/ + + You are morally obligated to cite the following paper in his/her + scientific literature that results from use of this file: + + [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, + Astrophysical Journal, 641, 626. + + PVS_GRMHD is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + PVS_GRMHD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with PVS_GRMHD; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + If the user has any questions, please direct them to Scott C. Noble at + scn(a)astro.rit.edu . + +***********************************************************************************/ + + + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <string.h> +#include <math.h> +#include <float.h> +#include <complex.h> + +#include "cctk.h" +#include "cctk_Parameters.h" + +/* Set this to be 1 if you want debug output */ +#define DEBUG_CON2PRIMM (0) + + +/* Adiabatic index used for the state equation */ + +#define MAX_NEWT_ITER (30) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ +#define NEWT_TOL (1.0e-10) /* Min. of tolerance allowed for Newton-Raphson iterations */ +#define MIN_NEWT_TOL (1.0e-10) /* Max. of tolerance allowed for Newton-Raphson iterations */ +#define EXTRA_NEWT_ITER (2) + +#define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ +#define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ + +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed + to always be smaller than this. This + is used to detect solver failures */ + +#define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ + +/************************************************** + The following functions assume a Gamma-law EOS: +***************************************************/ + +/* Local Globals */ +struct LocGlob { + CCTK_REAL Bsq, QdotBsq, Qtsq, Qdotn, D, half_Bsq, Sc, g_o_gm1, + W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW ; +} ; + +// Declarations: + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) ; + +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_pt(): + ----------------------------- + + -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess. + + -- Uses the 2D method of Noble et al. (2006): + -- Solves for two independent variables (W,v^2) via a 2D + Newton-Raphson method + -- Can be used (in principle) with a general equation of state. + + -- Minimizes two residual functions using a homemade Newton-Raphson routine. + -- It is homemade so that it can catch exceptions and handle them correctly, plus it is + optimized for this particular problem. + + -- Note that the notation used herein is that of Noble et al. (2006) except for the argument + list. + + +INPUT: (using GRHydro variable defintions) + + s[x,y,z] = scons[0,1,2] = \alpha \sqrt(\gamma) T^0_i + dens, tau = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + tau = \alpha^2 \sqrt(\gamma) T^{00} - D + g[x,y,z][x,y,x] = spatial metric corresponding to \gamma + u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] + det = sqrt(\gamma) + B[x,y,z] = Bvec[0,1,2] + bsq = b^\mu b_\mu + + epsnegative = (integer) + = 0 if rho and epsilon are positive + != 0 otherwise + + + -- (users should set B[x,y,z] = 0 for hydrodynamic runs) + + +OUTPUT: (using GRHydro variable defintions) + rho, eps = as defined in GRHydro, primitive variables + vel[x,y,z] = as defined in GRHydro, primitive variables + + +RETURN VALUE: of retval = (i*100 + j) where + i = 0 -> Newton-Raphson solver either was not called (yet or not used) + or returned successfully; + 1 -> Newton-Raphson solver did not converge to a solution with the + given tolerances; + 2 -> Newton-Raphson procedure encountered a numerical divergence + (occurrence of "nan" or "+/-inf" ; + + j = 0 -> success + 1 -> failure: some sort of failure in Newton-Raphson; + 2 -> failure: unphysical vsq = v^2 value at initial guess; + 3 -> failure: W<0 or W>W_TOO_BIG + 4 -> failure: v^2 > 1 + ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) + +**********************************************************************************/ + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_1d[1]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL tau, dens, gammaeos; + CCTK_REAL QdotB; + CCTK_REAL rho0,u,p,w,gammasq,gamma,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL rho_gm1; + CCTK_REAL utsq; + + DECLARE_CCTK_PARAMETERS; + + struct LocGlob lg; + + gammaeos = *gamma_eos; + + /* Assume ok initially: */ + // BCM: But let the driver function take care of its initialization + //*retval = 0.; + *epsnegative = 0; + +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *entropycons = %26.16e \n", *entropycons_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *pressure = %26.16e \n", *pressure ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif + + /* First undensitize all conserved variables : */ + sx = ( *sx_in) * inv_sqrt_detg; + sy = ( *sy_in) * inv_sqrt_detg; + sz = ( *sz_in) * inv_sqrt_detg; + tau = ( *tau_in) * inv_sqrt_detg; + dens = (*dens_in) * inv_sqrt_detg; + + usx = (*uxx)*sx + (*uxy)*sy + (*uxz)*sz; + usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; + usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + + lg.Sc = (*entropycons_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + lg.g_o_gm1 = gammaeos/(gammaeos-1.0); + + lg.Bsq = + (*gxx) * (*Bx) * (*Bx) + + (*gyy) * (*By) * (*By) + + (*gzz) * (*Bz) * (*Bz) + + 2*( + (*gxy) * (*Bx) * (*By) + + (*gxz) * (*Bx) * (*Bz) + + (*gyz) * (*By) * (*Bz) ); + + QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + lg.QdotBsq = QdotB*QdotB ; + + lg.Qdotn = -(tau + dens) ; + + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; + + lg.D = dens; + + lg.half_Bsq = 0.5*lg.Bsq; + + /* calculate W from last timestep and use for guess */ + vsq = + (*gxx) * (*velx) * (*velx) + + (*gyy) * (*vely) * (*vely) + + (*gzz) * (*velz) * (*velz) + + 2*( + (*gxy) * (*velx) * (*vely) + + (*gxz) * (*velx) * (*velz) + + (*gyz) * (*vely) * (*velz) ); + + if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { + vsq = fabs(vsq); + } + if(vsq < 0. || vsq > 1. ) { + *retval = 2.; + return; + } + + gammasq = 1. / (1. - vsq); + gamma = sqrt(gammasq); + + // Always calculate rho from D and gamma so that using D in EOS remains consistent + // i.e. you don't get positive values for dP/d(vsq) . + rho0 = lg.D / gamma ; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + +// W_last = w*gammasq ; + + // Calculate W and vsq: + x_1d[0] = rho0; +// *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; + *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ; + rho0 = x_1d[0]; + + /* Problem with solver, so return denoting error before doing anything further */ + if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) { + *retval = *retval*100.+1.; + return; + } + else{ + if( rho0 > W_TOO_BIG) { + *retval = 3.; + return; + } + } + + // Calculate v^2: + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); + gammasq = 1.+utsq; + gamma = sqrt(gammasq); + + if( utsq < 0. ) { + *retval = 4.; + return; + } + + + // Recover the primitive variables from the scalars and conserved variables: + rho0 = lg.D / gamma; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + W = w * gammasq; + + // User may want to handle this case differently, e.g. do NOT return upon + // a negative rho/u, calculate v^i so that rho/u can be floored by other routine: + if( (rho0 <= 0.) || (u <= 0.) ) { + *epsnegative = 1; + return; + } + + *rho = rho0; + *epsilon = u / rho0; + *w_lorentz = gamma; + *pressure = p ; + + g_o_WBsq = 1./(W+lg.Bsq); + QdB_o_W = QdotB / W; + *bsq = lg.Bsq / gammasq + QdB_o_W*QdB_o_W; + + *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; + *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; + *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } + + +#if(DEBUG_CON2PRIMM) + fprintf(stdout,"rho = %26.16e \n",*rho ); + fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); + fprintf(stdout,"pressure = %26.16e \n",*pressure ); + fprintf(stdout,"w_lorentz = %26.16e \n",*w_lorentz); + fprintf(stdout,"bsq = %26.16e \n",*bsq ); + fprintf(stdout,"velx = %26.16e \n",*velx ); + fprintf(stdout,"vely = %26.16e \n",*vely ); + fprintf(stdout,"velz = %26.16e \n",*velz ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + +/**********************************************************************/ +/************************************************************ + + general_newton_raphson(): + + -- performs Newton-Rapshon method on an arbitrary system. + + -- inspired in part by Num. Rec.'s routine newt(); + + Arguements: + + -- x[] = set of independent variables to solve for; + -- n = number of independent variables and residuals; + -- funcd = name of function that calculates residuals, etc.; + +*****************************************************************/ +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, i_extra, doing_extra; + CCTK_REAL W,W_old; + + int keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + + W = W_old = 0.; + + n_iter = 0; + + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + /* don't use line search : */ + x[0] += dx[0] ; + + /****************************************/ + /* Calculate the convergence criterion */ + /****************************************/ + + /* For the new criterion, always look at error in "W" : */ + // METHOD specific: + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + x[0] = fabs(x[0]); + + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) || + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { + return(2); + } + + + if( fabs(errx) > MIN_NEWT_TOL){ + return(1); + } + if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + + return(0); + +} + +/**********************************************************************/ +/********************************************************************************* + func_rho(): + + -- residual/jacobian routine to calculate rho Qtsq equation with + the definition of W + W = ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho + substituted in. + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + + CCTK_REAL t1,t2; + CCTK_REAL t12; + CCTK_REAL t17; + CCTK_REAL t3, t100,rhosq,t200,rho,W,dWdrho,dvsqdrho,vsq; + + rho = x[0]; + rhosq = rho*rho; + t200 = 1./(lgp->D*lgp->D); + t100 = lgp->g_o_gm1*lgp->Sc*pow(rho,(gammaeos-1.)); + W = lgp->D*( lgp->D + t100 ) / rho ; + dWdrho = lgp->D * ( -lgp->D + t100*(gammaeos-2.) ) / rhosq; + t1 = W*W; + t2 = lgp->Bsq+W; + // t3 = pow(Bsq+W,2.0); + t3 = t2*t2; + vsq = (lgp->D-rho)*(lgp->D+rho)*t200; + dvsqdrho = -2*rho*t200; + resid[0] = t1*(lgp->Qtsq-vsq*t3)+lgp->QdotBsq*(t2+W); + t12 = lgp->Bsq*lgp->Bsq; + t17 = dWdrho*vsq; + jac[0][0] = 2*lgp->QdotBsq*dWdrho + +((lgp->Qtsq-vsq*t12)*2*dWdrho+(-6*t17*lgp->Bsq-dvsqdrho*t12 + +(-2*dvsqdrho*lgp->Bsq-4*t17-dvsqdrho*W)*W)*W)*W; + + dx[0] = -resid[0]/jac[0][0]; + *f = 0.5*resid[0]*resid[0]; + *df = -2. * (*f); + + return; + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM File [modified]: make.code.defn Delta lines: +1 -0 =================================================================== --- trunk/src/make.code.defn 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/src/make.code.defn 2013-03-28 01:47:14 UTC (rev 500) @@ -46,6 +46,7 @@ GRHydro_Tmunu.F90 \ GRHydro_Con2PrimM.F90 \ GRHydro_Con2PrimM_pt.c \ + GRHydro_Con2PrimM_ptee.c \ GRHydro_Con2PrimM_polytype_pt.c \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +22 -0 =================================================================== --- trunk/interface.ccl 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/interface.ccl 2013-03-28 01:47:14 UTC (rev 500) @@ -84,6 +84,27 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ + CCTK_REAL INOUT prec, 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 Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT entropycons, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, 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 Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + 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_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + #void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,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 Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ @@ -207,6 +228,7 @@ #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenMee WITH GRHydro_Con2PrimM_ptee LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran

1 0

← Newer
1
2
3
4
5
6
7
8
9
Older →

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits March 2013