diff options
Diffstat (limited to 'src/GRHydro_Prim2Con.F90')
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 295 |
1 files changed, 175 insertions, 120 deletions
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 805e486..6d03e05 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -50,6 +50,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k + integer :: keytemp CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& g11r,g12r,g13r,g22r,g23r,g33r,avg_detr CCTK_REAL :: xtemp(1) @@ -122,71 +123,118 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) end do !$OMP END PARALLEL DO else - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& - !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & - !$OMP g11r,g12r,g13r,g22r,g23r,g33r) - 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 + if(reconstruct_temper.ne.0) then + keytemp = 1 + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) + 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 + + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) - g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) - g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) - g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) - g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) - g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) - g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) - g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) - g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) - g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) - g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) - g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + + call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel,& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11l,g12l,g13l,g22l,& + g23l,g33l, & + avg_detl,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),& + tempminus(i,j,k),y_e_minus(i,j,k)) + + call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel, & + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11r,g12r,g13r,& + g22r,g23r,g33r,& + avg_detr, 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),tempplus(i,j,k), & + y_e_plus(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + else + keytemp = 0 + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) + 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 + + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) - - ! 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 prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& - r(i,j,k),& - g11l,g12l,g13l,g22l,& - g23l,g33l, & - avg_detl,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),& - 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 prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & - cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& - r(i,j,k),& - g11r,g12r,g13r,& - g22r,g23r,g33r,& - avg_detr, 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),xtemp, & - y_e_plus(i,j,k)) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) + call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel,& + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11l,g12l,g13l,g22l,& + g23l,g33l, & + avg_detl,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),& + xtemp,y_e_minus(i,j,k)) + + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) + call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel, & + cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + r(i,j,k),& + g11r,g12r,g13r,& + g22r,g23r,g33r,& + avg_detr, 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),xtemp, & + y_e_plus(i,j,k)) + tempminus(i,j,k) = xtemp(1) + end do end do end do - end do - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO + endif endif end subroutine primitive2conservative @@ -245,7 +293,7 @@ subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & end subroutine prim2con -subroutine prim2con_hot(handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & +subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & temp,ye) @@ -264,9 +312,9 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & character(len=512) warnline ! begin EOS Omni vars - integer :: n,keytemp,anyerr,keyerr(1) + integer :: n,keytemp,lkeytemp,anyerr,keyerr(1) real*8 :: temp0(1) - n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0 + n = 1;lkeytemp=keytemp;anyerr = 0;keyerr(1) = 0 ! end EOS Omni vars temp0 = temp @@ -276,66 +324,70 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & ye = max(min(ye,GRHydro_Y_e_max),GRHydro_Y_e_min) - call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + call EOS_Omni_press(handle,lkeytemp,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,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r - 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) - !$OMP END CRITICAL - endif + if(reconstruct_temper.ne.0) then + 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. - if(GRHydro_eos_hot_prim2con_warn.ne.0) then - !$OMP CRITICAL - call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") - write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r - 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 - endif - 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,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r - 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 + 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 + lkeytemp=1 + call EOS_Omni_press(handle,lkeytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + lkeytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r + 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) + !$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. + if(GRHydro_eos_hot_prim2con_warn.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r + 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 + endif + lkeytemp=1 + temp = temp0 + call EOS_Omni_press(handle,lkeytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + lkeytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(i8,3i5,1P10E15.6)") cctk_iteration,ii,jj,kk,x,y,z,r + 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 endif @@ -387,6 +439,7 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k + CCTK_INT :: keytemp CCTK_REAL :: det character(len=512) :: warnline @@ -433,6 +486,7 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) end do !$OMP END PARALLEL DO else + keytemp = 0 !$OMP PARALLEL DO PRIVATE(i, j, k, det) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) @@ -441,7 +495,8 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,cctk_iteration,& + call prim2con_hot(GRHydro_eos_handle,keytemp,& + GRHydro_reflevel,cctk_iteration,& i,j,k,& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),g11(i,j,k),& g12(i,j,k),g13(i,j,k),& |