From 9df45a235ce39bfd46d83369f2c4ae859088abf9 Mon Sep 17 00:00:00 2001 From: rhaas Date: Mon, 27 Aug 2012 19:19:17 +0000 Subject: GRHydro: Changes to make GRHydro MHD work with nuclear/hot equation of state: 1) Switch to EOSOmni pointwise C2P routine and modify where necessary. 2) Modify Con2PrimM.F90 to allow for the evolution of temperature and adjust the wrapper routine. 3) Create EigenProblemM_hot pointwise routine and call that from HLLEM.F90 when temperature is evolved. Additionally adjust HLLEM where necessary. 4) Adjust InterfacesM.h to incorporate the newly created functions. 5) Fix a loop problem (not taking into account constraint transport) in PPM reconstruction of Y_e 6) Introduce Prim2ConM_hot and call this pointwise routine from Prim2ConM.F90 when temperature is evolved. Additionally also make this routine available to initial data routine in GRHydro_InitData 7) Adjust loops in GRHydro_PoloidalMagFieldM.F90 to not set boundary points it cannot set but instead call boundary group afterwards! Pay attention as this will not work with boundary conditions set to "none" in MHD case anymore but is the correct thing to do. 8) Allow StarMapper to extend HydroBase::initial_hydro = "starmapper". 9) Smaller fixes. From: Philipp Moesta git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@410 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Prim2ConM.F90 | 271 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 270 insertions(+), 1 deletion(-) (limited to 'src/GRHydro_Prim2ConM.F90') diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index 59acf0d..a6aef7d 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -56,9 +56,29 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + !CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,& + ! g11r,g12r,g13r,g22r,g23r,g33r + CCTK_REAL :: xtemp(1) + character(len=256) NaN_WarnLine + !TARGET gxx, gxy, gxz, gyy, gyz, gzz + + !CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + !g11 => gxx + !g12 => gxy + !g13 => gxz + !g22 => gyy + !g23 => gyz + !g33 => gzz + ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces + + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) @@ -99,6 +119,79 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO + + else + + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) + + !$OMP CRITICAL + if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then + write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k) + call CCTK_WARN(1, NaN_WarnLine) + endif + !$OMP END CRITICAL + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(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), w_lorentzminus(i, j, k), xtemp, y_e_minus(i,j,k)) + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i+xoffset,j+yoffset,k+zoffset)) + + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j,k),tauplus(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), & + w_lorentzplus(i,j,k), xtemp, y_e_plus(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif + end subroutine primitive2conservativeM @@ -117,6 +210,152 @@ end subroutine primitive2conservativeM @@*/ +subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & + dBconsx, dBconsy, dBconsz, & + drho, dvelx, dvely, dvelz,& + deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz + CCTK_REAL :: temp(1),ye(1), yein(1), x, y, z + CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk + CCTK_REAL :: w + CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz + character(len=256) NaN_WarnLine + character(len=512) warnline + +! begin EOS Omni vars + integer :: n, keytemp, anyerr, keyerr(1) + ! real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) + real*8 :: temp0(1) + n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 + !xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 +! end EOS Omni vars + + temp0 = temp + w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) + +!!$ BEGIN: Check for NaN value + if (w .ne. w) then + !$OMP CRITICAL + write(NaN_WarnLine,'(a100,6g15.6)') 'NaN produced in sqrt(): (gxx,gxy,gxz,gyy,gyz,gzz)', gxx, gxy, gxz, gyy, gyz, gzz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (x,y,z)', x, y, z + call CCTK_WARN(1, NaN_WarnLine) + !write(NaN_WarnLine,'(a100,g15.6)') 'NaN produced in sqrt(): w', w + !call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL + endif +!!$ END: Check for NaN value + + ye = max(min(ye,GRHydro_Y_e_max),GRHydro_Y_e_min) + + yein = ye + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then + ! in this case (coarse grid error that is hopefully restricted + ! away), we use the average temperature between cells and call + ! the EOS with keytemp=1 + keytemp=1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye,yein + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + endif + else + ! This is a way of recovering even on finer refinement levels: + ! Use the average temperature at the interface instead of the + ! reconstructed specific internal energy. + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + keytemp=1 + temp = temp0 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + endif + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) + ab0=w*Bdotv + b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 + blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & + w*Bdotv*vlowx + blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & + w*Bdotv*vlowy + blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & + w*Bdotv*vlowz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ab0*blowx) + dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + ab0*blowy) + dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + ab0*blowz) + dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + + dBconsx = sqrt(det)*dBvcx + dBconsy = sqrt(det)*dBvcy + dBconsz = sqrt(det)*dBvcz + + !!$OMP CRITICAL + !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w + !call CCTK_WARN(1, NaN_WarnLine) + !!$OMP END CRITICAL +end subroutine prim2conM_hot + + subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) @@ -202,7 +441,6 @@ end subroutine prim2conM @@*/ - subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) implicit none @@ -210,9 +448,13 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + CCTK_REAL :: xtemp(1) CCTK_INT :: i, j, k CCTK_REAL :: det + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,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 @@ -232,6 +474,31 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO + else + !$OMP PARALLEL DO PRIVATE(k,j,i,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 + + 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)) + + call prim2conM_hot(GRHydro_eos_handle,GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(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),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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),Bvecx(i,j,k), & + Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k), temperature(i,j,k), Y_e(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j, k) @@ -451,6 +718,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det + !$OMP PARALLEL DO PRIVATE(k,j,i,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 @@ -469,6 +737,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j, k) -- cgit v1.2.3