Commits December 2010

commits@lists.einsteintoolkit.org

5 participants
41 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 204)
by cott＠tapir.caltech.edu 31 Dec '10

31 Dec '10

User: cott Date: 2010/12/31 02:03 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: * update prim2con_hot pointwise function to conform to interface definition of EOS_Omni functions File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +9 -8 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2010-12-31 19:55:57 UTC (rev 203) +++ trunk/src/GRHydro_Prim2Con.F90 2010-12-31 20:03:19 UTC (rev 204) @@ -220,11 +220,12 @@ implicit none DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det - CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& - deps, dpress, w, vlowx, vlowy, vlowz - CCTK_REAL :: temp,ye, x, y, z + CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho(1), dvelx, dvely, dvelz,& + deps(1), dpress(1), w, vlowx, vlowy, vlowz + CCTK_REAL :: temp(1),ye(1), x, y, z CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk character(len=512) warnline @@ -278,11 +279,11 @@ vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - ddens = sqrt(det) * drho * w - dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx - dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy - dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz - dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens + ddens = sqrt(det) * drho(1) * w + dsx = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowx + dsy = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowy + dsz = sqrt(det) * (drho(1)*(1+deps(1))+dpress(1))*w*w * vlowz + dtau = sqrt(det) * ((drho(1)*(1+deps(1))+dpress(1))*w*w - dpress(1)) - ddens end subroutine prim2con_hot

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 203)
by cott＠tapir.caltech.edu 31 Dec '10

31 Dec '10

User: cott Date: 2010/12/31 01:55 PM Modified: /trunk/ schedule.ccl Log: * some more cleaning up in schedule.ccl -- we don't need separate schedule entries for the case with or without evolved tracers -- the tracers won't get synced if they are not evolved, because they don't have storage in that case. * add entropy to the list of variables to be synced File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +21 -96 =================================================================== --- trunk/schedule.ccl 2010-12-31 18:54:31 UTC (rev 202) +++ trunk/schedule.ccl 2010-12-31 19:55:57 UTC (rev 203) @@ -367,8 +367,6 @@ } if (CCTK_Equals(EoS_Change_type,"GammaKS")) { - - if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { schedule GRHydro_EoSChangeGammaK_Shibata AT CCTK_Initial AFTER HydroBase_Initial BEFORE GRHydro_IVP { LANG: Fortran @@ -381,24 +379,10 @@ SYNC: hydrobase::eps SYNC: hydrobase::vel SYNC: hydrobase::temperature + SYNC: hydrobase::entropy SYNC: hydrobase::Y_e SYNC: Y_e_con } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution" - } else { - schedule GRHydro_EoSChangeGammaK_Shibata AT CCTK_Initial AFTER HydroBase_Initial BEFORE GRHydro_IVP - { - LANG: Fortran - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: hydrobase::rho - SYNC: hydrobase::press - SYNC: hydrobase::eps - SYNC: hydrobase::vel - SYNC: hydrobase::temperature - } "Reset the hydro variables if the EoS Gamma and K change between ID and evolution" - } } } @@ -1031,86 +1015,27 @@ { } "GRHydro Boundary conditions group" -if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) +schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound { - if (evolve_tracer) - { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions" - } else { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: hydrobase::temperature - SYNC: hydrobase::Y_e - SYNC: Y_e_con - } "Select GRHydro boundary conditions" - } -} else { - if (evolve_tracer) - { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: GRHydro_cons_tracers - SYNC: GRHydro_tracers - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions - MHD version" - } else { - schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound - { - LANG: Fortran - OPTIONS: LEVEL - SYNC: dens - SYNC: tau - SYNC: scon - SYNC: w_lorentz - SYNC: HydroBase::rho - SYNC: HydroBase::press - SYNC: HydroBase::eps - SYNC: HydroBase::vel - SYNC: HydroBase::Bvec - SYNC: hydrobase::temperature - } "Select GRHydro boundary conditions - MHD version" - } -} + LANG: Fortran + OPTIONS: LEVEL + SYNC: dens + SYNC: tau + SYNC: scon + SYNC: w_lorentz + SYNC: HydroBase::rho + SYNC: HydroBase::press + SYNC: HydroBase::eps + SYNC: HydroBase::vel + SYNC: HydroBase::Bvec + SYNC: GRHydro_cons_tracers + SYNC: GRHydro_tracers + SYNC: hydrobase::temperature + SYNC: hydrobase::entropy + SYNC: hydrobase::Y_e + SYNC: Y_e_con + SYNC: HydroBase::Bvec +} "Select GRHydro boundary conditions" ############################################################ ### Compute first differences of rho for mesh refinement ###

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 118)
by bcmsma＠astro.rit.edu 31 Dec '10

31 Dec '10

User: bmundim Date: 2010/12/31 01:07 PM Modified: /trunk/ interface.ccl /trunk/doc/ documentation.tex /trunk/src/ GRHydro_C2P2CM.F90, GRHydro_ShockTubeM.F90 Log: 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. File Changes: Directory: /trunk/doc/ ====================== File [modified]: documentation.tex Delta lines: +1 -1 =================================================================== --- trunk/doc/documentation.tex 2010-12-25 06:44:22 UTC (rev 117) +++ trunk/doc/documentation.tex 2010-12-31 19:07:40 UTC (rev 118) @@ -118,7 +118,7 @@ \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). Directory: /trunk/src/ ====================== File [modified]: GRHydro_C2P2CM.F90 Delta lines: +37 -21 =================================================================== --- trunk/src/GRHydro_C2P2CM.F90 2010-12-25 06:44:22 UTC (rev 117) +++ trunk/src/GRHydro_C2P2CM.F90 2010-12-31 19:07:40 UTC (rev 118) @@ -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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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) File [modified]: GRHydro_ShockTubeM.F90 Delta lines: +13 -13 =================================================================== --- trunk/src/GRHydro_ShockTubeM.F90 2010-12-25 06:44:22 UTC (rev 117) +++ trunk/src/GRHydro_ShockTubeM.F90 2010-12-31 19:07:40 UTC (rev 118) @@ -118,22 +118,22 @@ 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 Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +17 -0 =================================================================== --- trunk/interface.ccl 2010-12-25 06:44:22 UTC (rev 117) +++ trunk/interface.ccl 2010-12-31 19:07:40 UTC (rev 118) @@ -79,12 +79,29 @@ 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

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 202)
by bcmsma＠astro.rit.edu 31 Dec '10

31 Dec '10

User: bmundim Date: 2010/12/31 12:54 PM Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_pt.c, GRHydro_PPMM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_RegisterGZ.cc, GRHydro_RegisterVars.cc, GRHydro_TmunuM.F90, make.code.defn Log: RIT MHD development: Revert previous commit changes to schedule.ccl Merge MHD and GRHydro routines: Boundaries, CalcUpdate, PPM, Reconstruct(poly), RegisterGZ, RegisterVars Solve a few race conditions arising in the "magnetic" routines. This is still being tested, since it passes tests with intel compilers but fails for gnu ones. Alias function GRHydro_Con2PrimM_pt for C2P2CM test. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: <Unavailable> File [modified]: GRHydro_CalcUpdate.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_PPMM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Prim2ConM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.F90 Delta lines: <Unavailable> File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: <Unavailable> File [modified]: GRHydro_RegisterGZ.cc Delta lines: <Unavailable> File [modified]: GRHydro_RegisterVars.cc Delta lines: <Unavailable> File [modified]: GRHydro_TmunuM.F90 Delta lines: <Unavailable> File [modified]: make.code.defn Delta lines: <Unavailable> Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 201)
by cott＠tapir.caltech.edu 30 Dec '10

30 Dec '10

User: cott Date: 2010/12/30 02:49 AM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: * update error handling in _hot routine: warn only starting from a given reflevel File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +60 -53 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2010-12-25 06:53:39 UTC (rev 200) +++ trunk/src/GRHydro_Con2Prim.F90 2010-12-30 08:49:37 UTC (rev 201) @@ -529,16 +529,18 @@ pold = max(1.d-15,xpress) ! error handling if(anyerr.ne.0) then - call CCTK_WARN(failinfomode,"EOS error in c2p 0") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(failinfomode,"EOS error in c2p 0") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + endif endif @@ -572,18 +574,19 @@ rho,epsilon,temp,ye,xpress,keyerr,anyerr) ! error handling if(anyerr.ne.0) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(failinfomode,"EOS error in c2p 1") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + endif endif - f = pold - xpress @@ -600,16 +603,18 @@ GRHydro_C2P_failed = 1 !$OMP CRITICAL - call CCTK_WARN(failinfomode, 'count > GRHydro_countmax! ') - write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,'(a20,3i5)') 'ijk location: ',ii,jj,kk - call CCTK_WARN(failinfomode,warnline) - write(warnline,'(a20,g16.7)') 'radius: ',r - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failinfomode,"Setting the point to atmosphere") + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(failinfomode, 'count > GRHydro_countmax! ') + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,'(a20,3i5)') 'ijk location: ',ii,jj,kk + call CCTK_WARN(failinfomode,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failinfomode,"Setting the point to atmosphere") + endif !$OMP END CRITICAL ! for safety, let's set the point to atmosphere @@ -661,18 +666,19 @@ rho,epsilon,temp,ye,xpress,keyerr,anyerr) ! error handling if(anyerr.ne.0) then - call CCTK_WARN(failinfomode,"EOS error in c2p 2") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(failinfomode,"EOS error in c2p 2") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + endif endif - f = pnew - xpress enddo @@ -710,18 +716,19 @@ ! error handling if(anyerr.ne.0) then - call CCTK_WARN(failinfomode,"EOS error in c2p 3") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(failinfomode,"EOS error in c2p 3") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + endif endif - f = pold - xpress

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 60)
by eric9＠mac.com 28 Dec '10

28 Dec '10

User: eric9 Date: 2010/12/28 11:13 AM Modified: /trunk/ einsteintoolkit.th Log: checkout GetComponents from GitHub instead File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +6 -0 =================================================================== --- trunk/einsteintoolkit.th 2010-12-13 17:01:05 UTC (rev 59) +++ trunk/einsteintoolkit.th 2010-12-28 17:13:56 UTC (rev 60) @@ -300,7 +300,13 @@ Carpet/LoopControl +!TARGET = $ROOT/bin +!TYPE = git +!URL = git://github.com/gridaphobe/CRL.git +!AUTH_URL = git@github.com:gridaphobe/CRL.git +!CHECKOUT = GetComponents + # Private thorns !TARGET = $ARR !TYPE = ignore

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 200)
by bcmsma＠astro.rit.edu 24 Dec '10

24 Dec '10

User: bmundim Date: 2010/12/25 12:53 AM Added: /trunk/src/ GRHydro_Con2PrimM_polytype_pt.c Log: Add missing file from previous commit. File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: +916 -0 =================================================================== --- trunk/src/GRHydro_Con2PrimM_polytype_pt.c (rev 0) +++ trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2010-12-25 06:53:39 UTC (rev 200) @@ -0,0 +1,916 @@ +/*********************************************************************************** + 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" + +/* 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 */ +static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq,rho_gm1; +static CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW; +static CCTK_REAL t2,t4,t7,t24,two_Bsq,t300,t400,s200,s100; + + +// Declarations: +static CCTK_REAL vsq_calc(CCTK_REAL W); + +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][1], + CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL) ); + +static void func_W( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], + CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); +static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ; + +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL) ); +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); + +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos); +/* pressure as a function of rho0 and u */ +static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) +{ + return((gammaeos - 1.)*u) ; +} + +/* Pressure as a function of rho0 and w = rho0 + u + p */ +static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w, CCTK_REAL gammaeos) +{ + return((gammaeos-1.)*(w - rho0)/gammaeos) ; +} + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( + CCTK_INT *handle, CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *sc_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + 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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_INT *epsnegative, + CCTK_REAL *retval); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_Polytype_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 = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + sc_in = K D, where K is the polytropic constant + 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_Polytype_pt) ( + CCTK_INT *handle, CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *sc_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + 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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_1d[1]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL dens, gammaeos; + CCTK_REAL QdotB,utsq,gamma_sq; + CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL rho_g, x_rho[1]; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_INT i,j, i_increase,ntries ; + + /* Assume ok initially: */ + *retval = 0.; + *epsnegative = 0; + + gammaeos = *gamma_eos; + +#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," *Sc = %26.16e \n", *sc_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," *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," *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," *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; + 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; + + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + 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)) ; + QdotBsq = QdotB*QdotB ; + + Qtsq = (usx * sx + usy * sy + usz * sz) ; + + D = dens; + + half_Bsq = 0.5*Bsq; + + t2 = D*D; + t4 = QdotBsq*t2; + t7 = Bsq*Bsq; + t24 = 1/t2; + two_Bsq = Bsq + Bsq; + t300 = QdotBsq*Bsq*t2; + t400 = Qtsq*t2; + + s200 = D*gammaeos*(*sc_in); + CCTK_REAL g_o_gm1 = (gammaeos/(gammaeos-1.)); + s100 = g_o_gm1*(*sc_in); + + /* 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 = D / gamma ; + CCTK_REAL gm1 = gammaeos-1.; + rho_gm1 = pow(rho0,gm1); + p = (*sc_in) * rho_gm1 / gamma; + u = p / gm1; + w = rho0 + u + p ; + + W_last = w*gammasq ; + + + // Make sure that W is large enough so that v^2 < 1 : + i_increase = 0; + while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) + - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq)) + && (i_increase < 10) ) { + W_last *= 10.; + i_increase++; + } + + W_for_gnr2 = W_for_gnr2_old = W_last; + rho_for_gnr2 = rho_for_gnr2_old = rho0; + + // Calculate W: + x_1d[0] = W_last; + + *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, func_W ) ; + + W = x_1d[0]; + + /* Problem with solver, so return denoting error before doing anything further */ + if( ((*retval) != 0.) || (W == FAIL_VAL) ) { + *retval = *retval*100.+1.; + return; + } + else{ + if(W <= 0. || W > W_TOO_BIG) { + *retval = 3.; + return; + } + } + + rho_g = x_rho[0] = rho_for_gnr2; + + ntries = 0; + while ( (*retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + rho_g *= 10.; + x_rho[0] = rho_g; + } + + rho_for_gnr2 = x_rho[0]; + + if( (*retval != 0) ) { + *retval = 10; + return; + } + + // Calculate v^2 : + rho0 = rho_for_gnr2; + rho_gm1 = pow(rho0,gm1); + + utsq = (D-rho0)*(D+rho0)/(rho0*rho0); + + gamma_sq = 1.+utsq; + gamma = sqrt(gamma_sq); + + // Calculate v^2: + if( vsq >= 1. ) { + *retval = 4.; + return; + } + + // Recover the primitive variables from the scalars and conserved variables: + + w = W / gamma_sq; + + // printf("doublecheck - S, rho, gamma: %e %e %e\n",*sc_in, rho_gm1,gamma); + + p = (*sc_in) * rho_gm1 / gamma; + + u = p / gm1; + + // 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+Bsq); + QdB_o_W = QdotB / W; + *bsq = Bsq * (1.-vsq) + 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(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,"gam = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + + +/**********************************************************************/ +/**************************************************************************** + vsq_calc(): + + -- evaluate v^2 (spatial, normalized velocity) from + W = \gamma^2 w + +****************************************************************************/ +static CCTK_REAL vsq_calc(CCTK_REAL W) +{ + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (Bsq + W); + Xsq = Bsq_W * Bsq_W; + + return( ( Wsq * Qtsq + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); +} + + +/******************************************************************** + + x1_of_x0(): + + -- calculates v^2 from W with some physical bounds checking; + -- asumes x0 is already physical + -- makes v^2 physical if not; + +*********************************************************************/ + +static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) +{ + CCTK_REAL x1,vsq; + CCTK_REAL dv = 1.e-15; + + vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive + + return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); + +} + +/******************************************************************** + + validate_x(): + + -- makes sure that x[0,1] have physical values, based upon + their definitions: + +*********************************************************************/ + +static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) +{ + + const CCTK_REAL dv = 1.e-15; + + /* Always take the absolute value of x[0] and check to see if it's too big: */ + x[0] = fabs(x[0]); + x[0] = (x[0] > W_TOO_BIG) ? x0[0] : x[0]; + + x[1] = (x[1] < 0.) ? 0. : x[1]; /* if it's too small */ + x[1] = (x[1] > 1.) ? (1. - dv) : x[1]; /* if it's too big */ + + return; + +} + +/************************************************************ + + oned_newton_raphson(): + + -- performs Newton-Rapshon method on an 2d system. + + -- inspired in part by Num. Rec.'s routine newt(); + +*****************************************************************/ +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_INT, CCTK_REAL) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, id, jd, i_extra, doing_extra; + CCTK_REAL dW,W,W_old; + + CCTK_INT keep_iterating, i_increase; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + //-fast for( id = 0; id < n ; id++) x_old[id] = x_orig[id] = x[id] ; + 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, n, gammaeos); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + + //-fast for( id = 0; id < n ; id++) { x_old[id] = x[id] ; } + x_old[0] = x[0] ; + + /* don't use line search : */ + //-fast for( id = 0; id < n ; id++) { x[id] += dx[id] ; } + x[0] += dx[0] ; + +// //METHOD specific: +// i_increase = 0; +// while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*Bsq ) - +// QdotBsq*(2.*x[0] + Bsq) ) <= x[0]*x[0]*(Qtsq-Bsq*Bsq)) +// && (i_increase < 10) ) { +// x[0] -= (1.*i_increase) * dx[0] / 10. ; +// i_increase++; +// } + + /****************************************/ + /* 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])) ) { +#if(DEBUG_CON2PRIMM) + fprintf(stderr,"\ngnr not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", + f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); +#endif + return(2); + } + + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + else { + return(1); + } + + return(0); + +} + + +/**********************************************************************/ +/********************************************************************************* + func_W() + + -- calculates the residuals, and Newton step for general_newton_raphson(); + -- for this method, x=W here; + + 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_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) +{ + CCTK_INT retval, ntries; + CCTK_REAL W, x_rho[1], rho, rho_g ; + CCTK_REAL t15, t200,t1000 ; + + W = x[0]; + W_for_gnr2_old = W_for_gnr2; + W_for_gnr2 = W; + + // get rho from NR: + rho_g = x_rho[0] = rho_for_gnr2; + + ntries = 0; + while ( (retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + rho_g *= 10.; + x_rho[0] = rho_g; + } + +#if(DEBUG_CON2PRIMM) + if( x_rho[0] <= 0. ) { + fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fflush(stderr); + } + + if( retval ) { + fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fflush(stderr); + } +#endif + + rho_for_gnr2_old = rho_for_gnr2; + rho = rho_for_gnr2 = x_rho[0]; + + CCTK_REAL gm1 = gammaeos-1.; + + rho_gm1 = pow(rho,gm1); + drho_dW = -rho*rho/( -rho_gm1*s200 + W*rho); + + // t6 = rho*rho; + // t100 = (D-rho)*(D+rho); // t2 - t6; + t15 = -(D-rho)*(D+rho); // t6-t2 + t200 = W + two_Bsq; + t1000 = rho*drho_dW; + resid[0] = (t300+(t4+t4+(t400+t15*(t7+(t200)*W))*W)*W)*t24; + jac[0][0] = 2*(t4+(t400+t15*t7+(3.0*t15*Bsq+t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*t24; + + dx[0] = -resid[0]/jac[0][0]; + + *df = - resid[0]*resid[0]; + *f = -0.5*(*df); + + +// fprintf(stdout,"QdotBsq = %28.18e ; \n",QdotBsq ); +// fprintf(stdout,"Sc = %28.18e ; \n",Sc ); +// fprintf(stdout,"Bsq = %28.18e ; \n",Bsq ); +// fprintf(stdout,"Qtsq = %28.18e ; \n",Qtsq ); +// fprintf(stdout,"Dc = %28.18e ; \n",D ); +// fprintf(stdout,"drhodW = %28.18e ; \n",drho_dW ); +// fprintf(stdout,"W = %28.18e ; \n",W ); +// fprintf(stdout,"rho = %28.18e ; \n",rho ); +// fprintf(stdout,"resid_W = %28.18e ; \n",resid[0] ); +// fprintf(stdout,"jac_W = %28.18e ; \n",jac[0][0]); +// fprintf(stdout,"deriv1 %g %g %g %g \n",W,resid[0],jac[0][0],dx[0]); + + return; + +} + + +/***********************************************************/ +/********************************************************************** + + gnr2() + + -- used to calculate rho from W + +*****************************************************************/ +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, id,jd, i_extra, doing_extra; + CCTK_REAL dW,W,W_old; + + CCTK_INT keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + //-fast for( id = 0; id < n ; id++) x_old[id] = x_orig[id] = x[id] ; + x_old[0] = x_orig[0] = x[0] ; + + n_iter = 0; + + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + //-fast errx = 0.; + //-fast for( id = 0; id < n ; id++) { x_old[id] = x[id] ; } + x_old[0] = x[0] ; + + /* Make the newton step: */ + //-fast for( id = 0; id < n ; id++) { x[id] += dx[id] ; } + x[0] += dx[0] ; + + /* Calculate the convergence criterion */ + //-fast for( id = 0; id < n ; id++) { errx += (x[id]==0.) ? fabs(dx[id]) : fabs(dx[id]/x[id]); } + //-fast errx /= 1.*n; + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + /* Make sure that the new x[] is physical : */ + // METHOD specific: + 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_TOL2) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 0; + } + + if( doing_extra == 1 ) i_extra++ ; + + // See if we've done the extra iterations, or have done too many iterations: + if( ((fabs(errx) <= NEWT_TOL2)&&(doing_extra == 0)) || + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { +#if(DEBUG_CON2PRIMM) + fprintf(stderr,"\ngnr2 not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", + f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); +#endif + return(2); + } + + // Return in different ways depending on whether a solution was found: + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + else { + return(1); + } + + return(0); + +} + +/********************************************************************************* + func_rho(): + + -- residual/jacobian routine to calculate rho from W via the polytrope: + + W = ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho + + 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[]; + *********************************************************************************/ +// for the isentropic version: eq. (27) +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) +{ + + CCTK_REAL A, B, C, rho, W, B0; + CCTK_REAL t40,t14; + + CCTK_REAL gm1 = gammaeos-1.; + + rho = x[0]; + W = W_for_gnr2; + + rho_gm1 = t40 = pow(rho,gm1); + + resid[0] = (rho*W+(-t40*s100-D)*D); + t14 = t40/rho; // rho^(g-2) + jac[0][0] = -t14*s200 + W; + // drho_dW = -rho/jac[0][0]; + + dx[0] = -resid[0]/jac[0][0]; + *df = - resid[0]*resid[0]; + *f = -0.5*(*df); + + // fprintf(stdout,"deriv3 %g %g %g %g %g \n",rho,W,resid[0],jac[0][0],dx[0]); +// fprintf(stdout,"Dc := %28.18e ; \n",D); +// fprintf(stdout,"Sc := %28.18e ; \n",Sc); +// fprintf(stdout,"rho := %28.18e ; \n",rho); +// fprintf(stdout,"W := %28.18e ; \n",W); +// fprintf(stdout,"resid_rho := %28.18e ; \n",resid[0] ); +// fprintf(stdout,"jac_rho := %28.18e ; \n",jac[0][0] ); + + return; + +} + + +/********************************************************************** + ********************************************************************** + + The following routines specify the equation of state. All routines + above here should be indpendent of EOS. If the user wishes + to use another equation of state, the below functions must be replaced + by equivalent routines based upon the new EOS. + + ********************************************************************** +**********************************************************************/ + +/**********************************************************************/ +/********************************************************************** + eos_info(): + + -- returns with all the EOS-related values needed; + **********************************************************************/ +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos) +{ + register double ftmp,gtmp; + + ftmp = 1. - vsq; + gtmp = sqrt(ftmp); + + CCTK_REAL gam_m1_o_gam = (gammaeos-1.)/gammaeos; + + *dpdw = gam_m1_o_gam * ftmp ; + *dpdvsq = gam_m1_o_gam * ( 0.5 * D/gtmp - W ) ; + + return( gam_m1_o_gam * ( W * ftmp - D * gtmp ) ); // p + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 199)
by bcmsma＠astro.rit.edu 24 Dec '10

24 Dec '10

User: bmundim Date: 2010/12/25 12:53 AM Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_pt.c, GRHydro_InterfacesM.h, GRHydro_RegisterVarsM.cc, make.code.defn Log: RIT MHD dev: Fix to GRHydro_RegisterVarsM.cc and allows MHD to work without Y_E bing evolved. Add C2P polytype interface. Change schedule.ccl to represent a more detailed version of choosing what gridfucntions get sync'd in Boundaries depending on what is being evolved. Un-hardwire fixed values of gamma for all MHD Con2Prim routines. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_InterfacesM.h Delta lines: <Unavailable> File [modified]: GRHydro_RegisterVarsM.cc Delta lines: <Unavailable> File [modified]: make.code.defn Delta lines: <Unavailable> Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 117)
by bcmsma＠astro.rit.edu 24 Dec '10

24 Dec '10

User: bmundim Date: 2010/12/25 12:44 AM Added: /trunk/src/ GRHydro_P2C2PM_polytype.F90 Removed: /trunk/ interface.ccl.omni Modified: /trunk/ param.ccl, schedule.ccl /trunk/src/ CheckParam.c, GRHydro_C2P2C.F90, GRHydro_C2P2CM.F90, GRHydro_Con2Prim.F90, GRHydro_P2C2P.F90, GRHydro_P2C2PM.F90, GRHydro_ReadConformalData.F90, GRHydro_ShockTube.F90, GRHydro_ShockTubeM.F90, GRHydro_SimpleWave.F90, make.code.defn Log: 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. File Changes: Directory: /trunk/src/ ====================== File [modified]: CheckParam.c Delta lines: +10 -1 =================================================================== --- trunk/src/CheckParam.c 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/CheckParam.c 2010-12-25 06:44:22 UTC (rev 117) @@ -23,6 +23,15 @@ { 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!"); + } } File [modified]: GRHydro_C2P2C.F90 Delta lines: +4 -16 =================================================================== --- trunk/src/GRHydro_C2P2C.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_C2P2C.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -34,13 +34,6 @@ 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 @@ CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -63,7 +55,6 @@ 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 @@ 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.' File [modified]: GRHydro_C2P2CM.F90 Delta lines: +11 -18 =================================================================== --- trunk/src/GRHydro_C2P2CM.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_C2P2CM.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -34,13 +34,6 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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,& File [modified]: GRHydro_Con2Prim.F90 Delta lines: +4 -17 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_Con2Prim.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -52,14 +52,6 @@ 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 @@ 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 @@ 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' File [modified]: GRHydro_P2C2P.F90 Delta lines: +6 -19 =================================================================== --- trunk/src/GRHydro_P2C2P.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_P2C2P.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -34,13 +34,6 @@ 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 @@ CCTK_INT C2P_failed logical epsnegative -#if USE_EOS_OMNI ! begin EOS Omni vars integer :: n = 1 integer :: keytemp = 0 @@ -63,7 +55,6 @@ 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 @@ 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' File [modified]: GRHydro_P2C2PM.F90 Delta lines: +13 -21 =================================================================== --- trunk/src/GRHydro_P2C2PM.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_P2C2PM.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -34,13 +34,6 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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,& File [added]: GRHydro_P2C2PM_polytype.F90 Delta lines: +181 -0 =================================================================== --- trunk/src/GRHydro_P2C2PM_polytype.F90 (rev 0) +++ trunk/src/GRHydro_P2C2PM_polytype.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -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 File [modified]: GRHydro_ReadConformalData.F90 Delta lines: +1 -14 =================================================================== --- trunk/src/GRHydro_ReadConformalData.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_ReadConformalData.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -43,18 +43,9 @@ 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 @@ 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 @@ 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 File [modified]: GRHydro_ShockTube.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_ShockTube.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_ShockTube.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -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 @@ 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),& File [modified]: GRHydro_ShockTubeM.F90 Delta lines: +93 -8 =================================================================== --- trunk/src/GRHydro_ShockTubeM.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_ShockTubeM.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -54,6 +54,13 @@ 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 @@ 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 File [modified]: GRHydro_SimpleWave.F90 Delta lines: +13 -24 =================================================================== --- trunk/src/GRHydro_SimpleWave.F90 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/GRHydro_SimpleWave.F90 2010-12-25 06:44:22 UTC (rev 117) @@ -46,17 +46,9 @@ 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 @@ 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 @@ ! 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) File [modified]: make.code.defn Delta lines: +2 -1 =================================================================== --- trunk/src/make.code.defn 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/src/make.code.defn 2010-12-25 06:44:22 UTC (rev 117) @@ -13,7 +13,8 @@ GRHydro_SimpleWave.F90 \ CheckParam.c\ GRHydro_P2C2P.F90 \ - GRHydro_P2C2PM.F90 + GRHydro_P2C2PM.F90 \ + GRHydro_P2C2PM_polytype.F90 # Subdirectories containing source files SUBDIRS = Directory: /trunk/ ================== File [removed]: interface.ccl.omni Delta lines: +0 -84 =================================================================== --- trunk/interface.ccl.omni 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/interface.ccl.omni 2010-12-25 06:44:22 UTC (rev 117) @@ -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" File [modified]: param.ccl Delta lines: +7 -0 =================================================================== --- trunk/param.ccl 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/param.ccl 2010-12-25 06:44:22 UTC (rev 117) @@ -22,6 +22,7 @@ "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 @@ "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 @@ shares:GRHydro +USES real GRHydro_eos_rf_prec USES real GRHydro_rho_central USES real GRHydro_eps_min "" USES real GRHydro_perc_ptol "" File [modified]: schedule.ccl Delta lines: +16 -0 =================================================================== --- trunk/schedule.ccl 2010-11-21 03:29:35 UTC (rev 116) +++ trunk/schedule.ccl 2010-12-25 06:44:22 UTC (rev 117) @@ -77,6 +77,22 @@ } } +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 {

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 198)
by schnetter＠cct.lsu.edu 22 Dec '10

22 Dec '10

User: eschnett Date: 2010/12/22 08:47 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: Loop in most prim2con routines over all grid points, not only the evolved ones. In this way, the conserved variables are defined on all grid points. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +18 -18 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2010-12-22 06:52:38 UTC (rev 197) +++ trunk/src/GRHydro_Prim2Con.F90 2010-12-23 02:47:41 UTC (rev 198) @@ -318,9 +318,9 @@ if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, det) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) 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)) @@ -337,9 +337,9 @@ !$OMP END PARALLEL DO else !$OMP PARALLEL DO PRIVATE(i, j, det) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) 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)) @@ -360,9 +360,9 @@ if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) enddo enddo @@ -551,9 +551,9 @@ CCTK_REAL :: det !$OMP PARALLEL DO PRIVATE(i, j, det) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) 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)) @@ -571,9 +571,9 @@ if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) enddo enddo @@ -764,9 +764,9 @@ ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall) - do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 - do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 - do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + do k = 1,cctk_lsh(3) + do j = 1,cctk_lsh(2) + do i = 1,cctk_lsh(1) gxxpt = gxx(i,j,k) gxypt = gxy(i,j,k)

1 0

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits December 2010