From e0dc2af4862d5ddb874328bd097f7f516231dd8c Mon Sep 17 00:00:00 2001 From: rhaas Date: Thu, 15 Sep 2011 16:49:53 +0000 Subject: add Multipatch support to GRHydro * not all features of GRHydro are supported yet, in particular only the HLLE solver supports Mulitpatch yet. Original commit by Christian Reisswig and Christian Ott git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@273 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Con2Prim.F90 | 114 +++++++++++++++++++++++------------------------ 1 file changed, 57 insertions(+), 57 deletions(-) (limited to 'src/GRHydro_Con2Prim.F90') diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index d3c9732..3c51302 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -101,11 +101,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) epsnegative = .false. - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) + gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),& + gbc(i,j,k),gcc(i,j,k)) !!$ if (det < 0.e0) then !!$ call CCTK_WARN(0, "nan produced (det negative)") @@ -135,7 +135,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 - vel(i,j,k,:) = 0.d0 + lvel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 if(evolve_temper.ne.0) then @@ -148,7 +148,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),keyerr,anyerr) else - ! w_lorentz=1, so the expression for tau reduces to: keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) @@ -157,6 +156,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) endif + ! w_lorentz=1, so the expression for tau reduces to: tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) cycle @@ -165,15 +165,15 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) if(evolve_temper.eq.0) then call Con2Prim_pt(GRHydro_eos_handle, 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), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), & + lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) else call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& - scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(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),temperature(i,j,k),y_e(i,j,k),& + scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),lvel(i,j,k,1),& + lvel(i,j,k,2),lvel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & @@ -186,7 +186,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) local_perc_ptol = GRHydro_perc_ptol*100.0d0 call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(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),& + lvel(i,j,k,1),lvel(i,j,k,2), lvel(i,j,k,3),eps(i,j,k),& temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & @@ -237,7 +237,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 - vel(i,j,k,:) = 0.d0 + lvel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -259,8 +259,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) !$OMP END CRITICAL call Con2Prim_ptPolytype(GRHydro_polytrope_handle, & 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), & + tau(i,j,k), rho(i,j,k), lvel(i,j,k,1), lvel(i,j,k,2), & + lvel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & uxx, uxy, uxz, uyy, uyz, uzz, det, x(i,j,k), y(i,j,k), & z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) #endif @@ -1048,18 +1048,18 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - 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)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) epsnegative = .false. @@ -1251,11 +1251,11 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) + gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),& + gbc(i,j,k),gcc(i,j,k)) if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -1277,8 +1277,8 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) call Con2Prim_ptPolytype(GRHydro_eos_handle, 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), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), & + lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) @@ -1600,18 +1600,18 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - 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)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(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) @@ -1686,18 +1686,18 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS) do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - 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)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(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) @@ -1878,7 +1878,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'velocities: ',& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3) + lvel(i,j,k,1),lvel(i,j,k,2),lvel(i,j,k,3) call CCTK_WARN(0,warnline) !$OMP END CRITICAL -- cgit v1.2.3