Commits April 2011

commits@lists.einsteintoolkit.org

12 participants
208 discussions

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 239)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 07:47 PM Modified: /branches/divclean/src/ GRHydro_SourceM.F90 Log: Fix Bcons sources. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_SourceM.F90 Delta lines: +26 -6 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 00:37:56 UTC (rev 238) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 00:47:43 UTC (rev 239) @@ -85,6 +85,7 @@ !!CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet + CCTK_REAL :: gdg_x, gdg_y, gdg_z !! g^{ik} d_k g_{ij} CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar @@ -463,23 +464,42 @@ dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + & 2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz) + !!$ g^{ik} d_k g_{li} + gdg_x = uxx*dx_gxx + uxy*dy_gxx + uxz*dz_gxx + & + uxy*dx_gxy + uyy*dy_gxy + uyz*dz_gxy + & + uxz*dx_gxz + uyz*dy_gxz + uzz*dz_gxz + + gdg_y = uxx*dx_gxy + uxy*dy_gxy + uxz*dz_gxy + & + uxy*dx_gyy + uyy*dy_gyy + uyz*dz_gyy + & + uxz*dx_gyz + uyz*dy_gyz + uzz*dz_gyz + + gdg_z = uxx*dx_gxz + uxy*dy_gxz + uxz*dz_gxz + & + uxy*dx_gyz + uyy*dy_gyz + uyz*dz_gyz + & + uxz*dx_gzz + uyz*dy_gzz + uzz*dz_gzz + bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_betax + & Bconsy(i,j,k)*dy_betax + Bconsz(i,j,k)*dz_betax ) + & psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + & - psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & - uxy*dy_det_bydet + uxz*dz_det_bydet ) + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & + uxy*dy_det_bydet + uxz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxx*gdg_x + uxy*gdg_y + & + uxz*gdg_z ) bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_betay + & Bconsy(i,j,k)*dy_betay + Bconsz(i,j,k)*dz_betay ) + & psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + & - psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & - uyy*dy_det_bydet + uyz*dz_det_bydet ) + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & + uyy*dy_det_bydet + uyz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxy*gdg_x + uyy*gdg_y + & + uyz*gdg_z ) bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_betaz + & Bconsy(i,j,k)*dy_betaz + Bconsz(i,j,k)*dz_betaz ) + & psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + & - psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & - uyz*dy_det_bydet + uzz*dz_det_bydet ) + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & + uyz*dy_det_bydet + uzz*dz_det_bydet ) - & + psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*gdg_x + uyz*gdg_y + & + uzz*gdg_z ) Bconsrhsx(i,j,k) = bvcx_source Bconsrhsy(i,j,k) = bvcy_source

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 238)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 07:37 PM Modified: /branches/divclean/src/ GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c, GRHydro_InterfacesM.h Log: Change interface for Con2Prim routines to reflect the new Bcons field. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +40 -47 =================================================================== --- branches/divclean/src/GRHydro_Con2PrimM.F90 2011-04-27 23:05:18 UTC (rev 237) +++ branches/divclean/src/GRHydro_Con2PrimM.F90 2011-04-28 00:37:56 UTC (rev 238) @@ -129,10 +129,6 @@ if(evolve_Y_e.ne.0) then Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) endif - - Bvec(i,j,k,1) = Bcons(i,j,k,1)/sqrt(det) - Bvec(i,j,k,2) = Bcons(i,j,k,2)/sqrt(det) - Bvec(i,j,k,3) = Bcons(i,j,k,3)/sqrt(det) if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then @@ -166,11 +162,11 @@ if(evolve_temper.eq.0) then call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& epsnegative,GRHydro_C2P_failed(i,j,k)) else stop "Please implement finite T con2prim routine in MHD part!" @@ -242,13 +238,13 @@ 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, & - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& epsnegative,GRHydro_C2P_failed(i,j,k)) - + #endif end if @@ -297,7 +293,7 @@ uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,sc + CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,scminus,scplus CCTK_INT :: epsnegative character(len=100) warnline @@ -390,11 +386,11 @@ call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) if (epsnegative .ne. 0) then @@ -410,16 +406,17 @@ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) - sc = local_K*densminus(i,j,k) + scminus = local_K*densminus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densminus(i,j,k), & - sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), sc, & - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus, & + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) + end if if (epsminus(i,j,k) .lt. 0.0d0) then @@ -446,13 +443,13 @@ epsnegative = 0 call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) - + if (epsnegative .ne. 0) then !$OMP CRITICAL call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!') @@ -466,15 +463,15 @@ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) - sc = local_K*densplus(i,j,k) + scplus = local_K*densplus(i,j,k) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densplus(i,j,k), & - sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), sc,& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) end if @@ -602,10 +599,6 @@ Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) endif - Bvec(i,j,k,1) = Bcons(i,j,k,1)/sqrt(det) - Bvec(i,j,k,2) = Bcons(i,j,k,2)/sqrt(det) - Bvec(i,j,k,3) = Bcons(i,j,k,3)/sqrt(det) - xrho=10.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& 1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) @@ -617,14 +610,14 @@ sc = local_K*dens(i,j,k) 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,& - rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & - Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,& - epsnegative,GRHydro_C2P_failed(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(i,j,k),& + vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end do end do end do @@ -730,20 +723,20 @@ call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), & sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,& - rhominus(i,j,k),& - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),& + Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & - Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,& epsnegative,GRHydro_C2P_failed(i,j,k)) call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), & sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& - rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),& + Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & - Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,& epsnegative,GRHydro_C2P_failed(i,j,k)) end do end do File [modified]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: +17 -8 =================================================================== --- branches/divclean/src/GRHydro_Con2PrimM_polytype_pt.c 2011-04-27 23:05:18 UTC (rev 237) +++ branches/divclean/src/GRHydro_Con2PrimM_polytype_pt.c 2011-04-28 00:37:56 UTC (rev 238) @@ -103,17 +103,18 @@ CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *sc_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -184,17 +185,18 @@ CCTK_REAL *dens_in, CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, CCTK_REAL *sc_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -228,12 +230,19 @@ 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," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_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 ); @@ -248,10 +257,6 @@ 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); @@ -267,6 +272,10 @@ 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; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: lg.Bsq = File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: +17 -10 =================================================================== --- branches/divclean/src/GRHydro_Con2PrimM_pt.c 2011-04-27 23:05:18 UTC (rev 237) +++ branches/divclean/src/GRHydro_Con2PrimM_pt.c 2011-04-28 00:37:56 UTC (rev 238) @@ -108,18 +108,17 @@ 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 *tau_in, + CCTK_REAL *tau_in, CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval); @@ -189,18 +188,19 @@ 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 *tau_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_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_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, - CCTK_REAL *bsq, CCTK_INT *epsnegative, CCTK_REAL *retval) @@ -231,12 +231,19 @@ 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," *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 ); @@ -251,10 +258,6 @@ 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); @@ -271,6 +274,10 @@ 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; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: lg.Bsq = File [modified]: GRHydro_InterfacesM.h Delta lines: +12 -8 =================================================================== --- branches/divclean/src/GRHydro_InterfacesM.h 2011-04-27 23:05:18 UTC (rev 237) +++ branches/divclean/src/GRHydro_InterfacesM.h 2011-04-28 00:37:56 UTC (rev 238) @@ -7,17 +7,18 @@ local_gam, dens, & sx, sy, sz, & tau, & + Bconsx, Bconsy, Bconsz, & rho, & velx, vely, velz,& epsilon, pressure,& + Bx, By, Bz, & + bsq,& w_lorentz, & gxx, gxy, gxz, & gyy, gyz, gzz, & uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - Bx, By, Bz, & - bsq,& epsnegative, & retval) @@ -27,17 +28,18 @@ CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL tau + CCTK_REAL Bconsx, Bconsy, Bconsz CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq CCTK_REAL w_lorentz CCTK_REAL gxx, gxy, gxz CCTK_REAL gyy, gyz, gzz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det - CCTK_REAL Bx, By, Bz - CCTK_REAL bsq CCTK_INT epsnegative CCTK_REAL retval end subroutine GRHydro_Con2PrimM_pt @@ -46,17 +48,18 @@ dens, & sx, sy, sz, & sc, & + Bconsx, Bconsy, Bconsz, & rho, & velx, vely, velz,& epsilon, pressure,& + Bx, By, Bz, & + bsq,& w_lorentz, & gxx, gxy, gxz, & gyy, gyz, gzz, & uxx, uxy, uxz,& uyy, uyz, uzz,& det,& - Bx, By, Bz, & - bsq,& epsnegative, & retval) @@ -66,17 +69,18 @@ CCTK_REAL dens CCTK_REAL sx, sy, sz CCTK_REAL sc + CCTK_REAL Bconsx, Bconsy, Bconsz CCTK_REAL rho CCTK_REAL velx, vely, velz CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq CCTK_REAL w_lorentz CCTK_REAL gxx, gxy, gxz CCTK_REAL gyy, gyz, gzz CCTK_REAL uxx, uxy, uxz CCTK_REAL uyy, uyz, uzz CCTK_REAL det - CCTK_REAL Bx, By, Bz - CCTK_REAL bsq CCTK_INT epsnegative CCTK_REAL retval end subroutine GRHydro_Con2PrimM_Polytype_pt

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 237)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 06:05 PM Modified: /branches/divclean/src/ GRHydro_SourceM.F90 Log: Fix variable names. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_SourceM.F90 Delta lines: +7 -7 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 22:58:37 UTC (rev 236) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 23:05:18 UTC (rev 237) @@ -444,7 +444,7 @@ taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & - dx_shiftx + dy_shifty + dz_shiftz ) * psidc(i,j,k) + & + dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & Bconsx(i,j,k) * ( dx_alp - half*alp(i,j,k) * & ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + 2.d0*uxy*dx_gxy + & 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/ sqrtdet + & @@ -463,20 +463,20 @@ dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + & 2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz) - bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftx + & - Bconsy(i,j,k)*dy_shiftx + Bconsz(i,j,k)*dz_shiftx ) + & + bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_betax + & + Bconsy(i,j,k)*dy_betax + Bconsz(i,j,k)*dz_betax ) + & psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + & psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & uxy*dy_det_bydet + uxz*dz_det_bydet ) - bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_shifty + & - Bconsy(i,j,k)*dy_shifty + Bconsz(i,j,k)*dz_shifty ) + & + bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_betay + & + Bconsy(i,j,k)*dy_betay + Bconsz(i,j,k)*dz_betay ) + & psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + & psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & uyy*dy_det_bydet + uyz*dz_det_bydet ) - bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftz + & - Bconsy(i,j,k)*dy_shiftz + Bconsz(i,j,k)*dz_shiftz ) + & + bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_betaz + & + Bconsy(i,j,k)*dy_betaz + Bconsz(i,j,k)*dz_betaz ) + & psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + & psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & uyz*dy_det_bydet + uzz*dz_det_bydet )

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 236)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:58 PM Modified: /branches/divclean/src/ GRHydro_RegisterVars.cc Log: Remove Bvecrhs. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_RegisterVars.cc Delta lines: +1 -1 =================================================================== --- branches/divclean/src/GRHydro_RegisterVars.cc 2011-04-27 22:56:11 UTC (rev 235) +++ branches/divclean/src/GRHydro_RegisterVars.cc 2011-04-27 22:58:37 UTC (rev 236) @@ -63,7 +63,7 @@ register_evolved("GRHydro::scon", "GRHydro::srhs"); if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { - register_constrained("HydroBase::Bvec", "GRHydro::Bvecrhs"); + register_constrained("HydroBase::Bvec"); register_evolved("HydroBase::Bcons", "GRHydro::Bconsrhs"); if(clean_divergence) { register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs");

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 235)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:56 PM Modified: /branches/divclean/src/ GRHydro_ReconstructPoly.F90 Log: Move erroneous end do statements. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: +1 -1 =================================================================== --- branches/divclean/src/GRHydro_ReconstructPoly.F90 2011-04-27 22:37:19 UTC (rev 234) +++ branches/divclean/src/GRHydro_ReconstructPoly.F90 2011-04-27 22:56:11 UTC (rev 235) @@ -693,6 +693,7 @@ end if end do end do + end do !$OMP END PARALLEL DO else if (flux_direction == 2) then !$OMP PARALLEL DO PRIVATE(i, j) @@ -764,7 +765,6 @@ end do end do end do - end do !$OMP END PARALLEL DO else if (flux_direction == 3) then !$OMP PARALLEL DO PRIVATE(i, j)

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 72)
by knarf＠cct.lsu.edu 27 Apr '11

27 Apr '11

User: knarf Date: 2011/04/27 05:55 PM Modified: / ET.tex Log: make clear that we should put all references in one file, using one convention File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +1 -0 =================================================================== --- ET.tex 2011-04-27 18:57:16 UTC (rev 71) +++ ET.tex 2011-04-27 22:55:25 UTC (rev 72) @@ -1889,6 +1889,7 @@ project cfz-411-aa}. \bibliographystyle{amsplain-url} +% TODO: all references should be in manifest/einsteintoolkit \bibliography{manifest/einsteintoolkit,local_bibtex/ott_references,local_bibtex/references} \end{document}

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 125)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:37 PM Modified: /trunk/ interface.ccl Log: Change Prim2Con functions interface. File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +17 -14 =================================================================== --- trunk/interface.ccl 2011-04-27 22:26:38 UTC (rev 124) +++ trunk/interface.ccl 2011-04-27 22:37:50 UTC (rev 125) @@ -40,18 +40,19 @@ CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) SUBROUTINE Prim2ConPolyM(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 Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, 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) + 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 OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, 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 IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, CCTK_REAL OUT w_lorentz) - SUBROUTINE Prim2ConGenM(CCTK_INT IN handle, \ CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ @@ -59,12 +60,14 @@ 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 Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, 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) + CCTK_REAL OUT press, CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, 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, \

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/ (Rev. 234)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:37 PM Modified: /branches/divclean/ interface.ccl Log: fix typo. File Changes: Directory: /branches/divclean/ ============================== File [modified]: interface.ccl Delta lines: +1 -1 =================================================================== --- branches/divclean/interface.ccl 2011-04-27 22:20:36 UTC (rev 233) +++ branches/divclean/interface.ccl 2011-04-27 22:37:19 UTC (rev 234) @@ -424,7 +424,7 @@ Bconsxplus,Bconsyplus,Bconszplus,Bconsxminus,Bconsyminus,Bconszminus } "Conservative variables extended to the cell boundaries" -real GRHydro_MHD_prin_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' +real GRHydro_MHD_prim_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { Bvecxplus,Bvecyplus,Bveczplus,Bvecxminus,Bvecyminus,Bveczminus } "Primitive mhd variables extended to the cell boundaries"

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 95)
by knarf＠cct.lsu.edu 27 Apr '11

27 Apr '11

User: knarf Date: 2011/04/27 05:33 PM Modified: /trunk/ einsteintoolkit.th Log: make the hg-version of Carpet standard in the dev-version of the ET toolkit File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +8 -7 =================================================================== --- trunk/einsteintoolkit.th 2011-04-21 19:59:51 UTC (rev 94) +++ trunk/einsteintoolkit.th 2011-04-27 22:33:23 UTC (rev 95) @@ -288,13 +288,14 @@ # Carpet, the AMR driver !TARGET = $ARR -!TYPE = git -!AUTH_URL = carpetgit@carpetcode.org:carpet -!URL = git://carpetcode.org/carpet -# For experimental version, comment previous three lines and uncomment next two -# Note that you will need to get access first (ask Erik Schnetter) -#!TYPE = hg -#!URL = ssh://carpetmercurial@carpetcode.org/carpet +!TYPE = hg +!AUTH_URL = ssh://carpetmercurial@carpetcode.org/carpet +!URL = http://www.carpetcode.org/hg/carpet +# For the git version of Carpet, comment the previous three lines and uncomment +# the next three +#!TYPE = git +#!AUTH_URL = carpetgit@carpetcode.org:carpet +#!URL = git://carpetcode.org/carpet !CHECKOUT = Carpet/doc Carpet/Carpet Carpet/CarpetEvolutionMask

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 124)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:26 PM Added: /trunk/src/ GRHydro_CylindricalExplosionM.F90, GRHydro_MonopoleM.F90 Log: RIT MHD dev: 1) Komissarov's cylindrical explosion test 2) Magnetic monopole test File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_CylindricalExplosionM.F90 Delta lines: +164 -0 =================================================================== --- trunk/src/GRHydro_CylindricalExplosionM.F90 (rev 0) +++ trunk/src/GRHydro_CylindricalExplosionM.F90 2011-04-27 22:26:38 UTC (rev 124) @@ -0,0 +1,164 @@ + /*@@ + @file GRHydro_CylindricalExplosionM.F90 + @date Apr 22, 2011 + @author Scott Noble, Joshua Faber, Bruno Mundim + @desc + Cylindrical magnetized shocks. + @enddesc + @@*/ + +#include "cctk.h" +#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) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) + + + /*@@ + @routine GRHydro_cylindricalexplosionM + @date Fri Apr 22 14:51:37 EDT 2011 + @author Scott Noble, Joshua Faber, Bruno Mundim + @desc + Initial data for cylidrically symmetric shocks with + magnetic fields. Meant to recreate tests demonstrated + in Komissarov (1999). + @enddesc + @calls + @calledby + @history + Using GRHydro_ShockTubeM.F90 as a template. + @endhistory + +@@*/ + +subroutine GRHydro_cylindricalexplosionM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: direction, det + CCTK_REAL :: rhol, rhor, pressl, pressr + CCTK_REAL :: bvcxl,bvcyl,bvczl + CCTK_REAL :: tmp,tmp2,gam,r_inner,r_outer + + CCTK_REAL :: cyl_fr + +!!$Uses Bx_init, By_init and Bz_init to set magnetic field strength +!!$Original tests had Bx_init = 0.1, 1.0 with By_init=Bz_init = 0 + + bvcxl = Bx_init + bvcyl = By_init + bvczl = Bz_init + +!!$Adiabatic index for test + gam = (4.d0/3.d0) + +!!$Inner radius and outer radius + r_inner = 8.d-1 + r_outer = 1.d0 + +!!$Inner values + rhol = 1.d-4 + pressl = 3.d-5 + +!!$Outer values + rhor = 1.d-2 + pressr = 1.d0 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + +!!$direction represents the cylindrical radius here. + if (CCTK_EQUALS(shocktube_type,"xshock")) then + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (y(i,j,k)-shock_ypos)**2) + else if (CCTK_EQUALS(shocktube_type,"yshock")) then + direction = sqrt((y(i,j,k)-shock_ypos)**2+& + (z(i,j,k)-shock_zpos)**2) + else if (CCTK_EQUALS(shocktube_type,"zshock")) then + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (z(i,j,k)-shock_zpos)**2) + end if + + tmp = cyl_fr(direction,r_inner,r_outer,rhol,rhor) + tmp2 = cyl_fr(direction,r_inner,r_outer,pressl,pressr)/( (gam - 1.d0) * tmp ) + rho(i,j,k) = tmp + eps(i,j,k) = tmp2 + + velx(i,j,k) = 0.d0 + vely(i,j,k) = 0.d0 + velz(i,j,k) = 0.d0 + Bvecx(i,j,k)=bvcxl + Bvecy(i,j,k)=bvcyl + Bvecz(i,j,k)=bvczl + + 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 Prim2ConPolyM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGenM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Bvecrhs = 0.d0 + + + return + +end subroutine GRHydro_CylindricalExplosionM + + + +function cyl_fr(R,r_inner,r_outer,min_f,max_f) + + implicit none + + CCTK_REAL :: cyl_fr + CCTK_REAL :: R,r_inner,r_outer,min_f,max_f + + if(R .gt. r_outer) then + cyl_fr = min_f + else if( R .lt. r_inner) then + cyl_fr = max_f + else + cyl_fr = exp( ( log(max_f)*(r_outer - R) + log(min_f)*(R - r_inner) ) / (r_outer - r_inner) ) + end if + + return + +end function cyl_fr File [added]: GRHydro_MonopoleM.F90 Delta lines: +109 -0 =================================================================== --- trunk/src/GRHydro_MonopoleM.F90 (rev 0) +++ trunk/src/GRHydro_MonopoleM.F90 2011-04-27 22:26:38 UTC (rev 124) @@ -0,0 +1,109 @@ + /*@@ + @file GRHydro_MonopoleM.F90 + @date Sep 23, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + Initial data of the shock tube type - MHD version. + @enddesc + @@*/ + +#include "cctk.h" +#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) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) + + /*@@ + @routine GRHydro_MonopoleM + @date Sat Jan 26 02:53:49 2002 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + A monopole in space + @enddesc + @calls + @calledby + @history + Expansion and alteration of the test code from GRAstro_Hydro, + written by Mark Miller. + @endhistory + +@@*/ + +subroutine GRHydro_MonopoleM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: direction, det + CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & + velzl, velzr, epsl, epsr + CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3 + + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + + rho(i,j,k) = 1.0 + velx(i,j,k) = 0.0 + vely(i,j,k) = 0.0 + velz(i,j,k) = 0.0 + eps(i,j,k) = 0.1 + if(i.eq.nx/2+1.and.j.eq.ny/2+1.and.k.eq.nz/2+1) then + Bvecx(i,j,k)=0.1 + else + Bvecx(i,j,k)=0.0 + endif + Bvecy(i,j,k)=0.0 + Bvecz(i,j,k)=0.0 + + 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 Prim2ConPolyM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGenM(GRHydro_eos_handle,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, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Bvecrhs = 0.d0 + + + return + +end subroutine GRHydro_MonopoleM

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits April 2011