Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 552)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_Reconstruct.F90, GRHydro_Source.F90 Log: GRHydro: fix compiler warnings * remove unused variable * pacify compiler warning about possibly uninitialized value not sure if one could somehow get the compiler to realize that the condition on which access is based is false in both case, ie that evolve_MHD is an effective constant. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Reconstruct.F90 Delta lines: +3 -0 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2013-07-06 18:10:53 UTC (rev 551) +++ trunk/src/GRHydro_Reconstruct.F90 2013-07-06 18:10:57 UTC (rev 552) @@ -108,6 +108,9 @@ else local_min_tracer = 0.0d0 end if + else + ! pacify compiler warning about unintialized local_min_tracer + local_min_tracer = 1d42 end if File [modified]: GRHydro_Source.F90 Delta lines: +2 -3 =================================================================== --- trunk/src/GRHydro_Source.F90 2013-07-06 18:10:53 UTC (rev 551) +++ trunk/src/GRHydro_Source.F90 2013-07-06 18:10:57 UTC (rev 552) @@ -77,7 +77,6 @@ CCTK_REAL :: db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc CCTK_REAL :: dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc CCTK_REAL :: da, db, dc, ida, idb, idc - CCTK_REAL :: psi4pt, da_psi4, db_psi4, dc_psi4 CCTK_REAL :: shiftshiftk, shiftka, shiftkb, shiftkc CCTK_REAL :: sumTK CCTK_REAL :: halfshiftdga, halfshiftdgb, halfshiftdgc @@ -186,7 +185,7 @@ !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& !$OMP localgaa,localgab,localgac,localgbb,localgbc,localgcc,& - !$OMP psi4pt,det,sqrtdet,rhoenthalpyW2,shifta,shiftb,shiftc,& + !$OMP det,sqrtdet,rhoenthalpyW2,shifta,shiftb,shiftc,& !$OMP da_betaa,da_betab,da_betac,db_betaa,db_betab,db_betac,& !$OMP dc_betaa,dc_betab,dc_betac,velashift,velbshift,velcshift,& !$OMP vlowa,vlowb,vlowc,t00,t0a,t0b,t0c,taa,tab,tac,tbb,tbc,tcc,& @@ -195,7 +194,7 @@ !$OMP da_gaa, da_gab, da_gac, da_gbb, da_gbc, da_gcc,& !$OMP db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc,& !$OMP dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc,& - !$OMP da_psi4,db_psi4,dc_psi4,shiftshiftk,shiftka,shiftkb,shiftkc,& + !$OMP shiftshiftk,shiftka,shiftkb,shiftkc,& !$OMP sumTK,halfshiftdga,halfshiftdgb,halfshiftdgc,& !$OMP halfTdga,halfTdgb,halfTdgc,invalp,invalp2) do k=1 + GRHydro_stencil,nc - GRHydro_stencil

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 551)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90 Log: GRHydro: Con2PrimM: Set points to atmosphere if in excised region. From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +44 -4 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2013-07-06 18:10:49 UTC (rev 550) +++ trunk/src/GRHydro_Con2PrimM.F90 2013-07-06 18:10:53 UTC (rev 551) @@ -63,7 +63,7 @@ CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1) CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold - CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum + CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum, dummy1, dummy2 CCTK_INT :: niter CCTK_INT :: epsnegative character(len=100) warnline @@ -167,14 +167,17 @@ !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, & !$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, & !$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, & - !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress) + !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress, dummy1, dummy2) do k = 1, nz do j = 1, ny do i = 1, nx + !do not compute if in atmosphere region! + if (atmosphere_mask(i,j,k) .gt. 0) cycle + !do not compute if in atmosphere or in excised region - if ((atmosphere_mask(i,j,k) .ne. 0) .or. & - (hydro_excision_mask(i,j,k) .ne. 0)) cycle + !if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + ! (hydro_excision_mask(i,j,k) .ne. 0)) cycle epsnegative = 0 @@ -201,6 +204,43 @@ endif + ! In excised region, set to atmosphere! + if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .gt. 0)) then + SET_ATMO_MIN(dens(i,j,k), sqrt(det)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) !GRHydro_rho_min + scon(i,j,k,:) = 0.d0 + vup(i,j,k,:) = 0.d0 + w_lorentz(i,j,k) = 1.d0 + + if(evolve_temper.ne.0) then + ! set hot atmosphere values + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k) + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + else + 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) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + endif + + !!$ tau does need to take into account the existing B-field + !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + + if(tau(i,j,k).le.sdet*b2*0.5d0)then + tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0 + endif + + cycle + + endif + if(evolve_Y_e.ne.0) then Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),& GRHydro_Y_e_min)

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 550)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ Con2Prim_fortran_interfaces.F90 Log: GRHydro: whitespace change, tab->8 spaces From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: Con2Prim_fortran_interfaces.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/Con2Prim_fortran_interfaces.F90 2013-07-06 18:10:45 UTC (rev 549) +++ trunk/src/Con2Prim_fortran_interfaces.F90 2013-07-06 18:10:49 UTC (rev 550) @@ -6,8 +6,8 @@ interface subroutine Con2Prim_pt(& - cctk_iteration,ii,jj,kk,& - handle, & + cctk_iteration,ii,jj,kk,& + handle, & dens, & sx, sy, sz, & tau, &

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 549)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_Con2PrimHot.F90 Log: GRHydro: small cosmetic change to Con2PrimHot From: Christian David Ott <cott(a)zwicky.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimHot.F90 Delta lines: +5 -1 =================================================================== --- trunk/src/GRHydro_Con2PrimHot.F90 2013-07-06 18:10:41 UTC (rev 548) +++ trunk/src/GRHydro_Con2PrimHot.F90 2013-07-06 18:10:45 UTC (rev 549) @@ -229,11 +229,15 @@ ! Now reset this point to minimum temperature. GRHydro_C2P_failed(i,j,k) = 0 !$OMP CRITICAL - write(warnline,"(A10,i7,4i3,1P10E15.6)") "reset T:",& + write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& cctk_iteration,GRHydro_Reflevel,& i,j,k,rho(i,j,k),& eps(i,j,k),temperature(i,j,k),y_e(i,j,k) call CCTK_WARN(1,warnline) + write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& + cctk_iteration,GRHydro_Reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) + call CCTK_WARN(1,warnline) !$OMP END CRITICAL temperature(i,j,k) = GRHydro_hot_atmo_temp keytemp = 1

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 548)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_PPM.F90 Log: GRHydro: remove unreachable code From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +0 -11 =================================================================== --- trunk/src/GRHydro_PPM.F90 2013-07-06 18:10:37 UTC (rev 547) +++ trunk/src/GRHydro_PPM.F90 2013-07-06 18:10:41 UTC (rev 548) @@ -315,17 +315,6 @@ if (.not.apply_enhanced_ppm) then ! Only for 1984 PPM scheme! if (ppm_detect .ne. 0) then - - if (apply_enhanced_ppm) then - ! make sure drho, d2rho and dmrho are computed everywhere! - do i = 2, nx-1 - drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) - d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1)) - end do - do i = 2, nx-1 - STEEP(rho, drho, dmrho) - end do - endif do i = 3, nx - 2 if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - &

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 547)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_CalcBcom.F90 Log: GRHydro: Don't compute Bcom twice. From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcBcom.F90 Delta lines: +5 -0 =================================================================== --- trunk/src/GRHydro_CalcBcom.F90 2013-07-06 18:10:32 UTC (rev 546) +++ trunk/src/GRHydro_CalcBcom.F90 2013-07-06 18:10:37 UTC (rev 547) @@ -40,6 +40,11 @@ call CCTK_WARN(0,"No storage for GRHydro::bcom") end if +! Only compute if Tmunu is not computed! (Bcom is also computed with Tmunu) + if (stress_energy_state .ne. 0) then + return + endif + !$OMP PARALLEL DO PRIVATE(i,j,k,velxlow, velylow, velzlow,& !$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum1,dum2,b2,bxlow,bylow,bzlow)

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 546)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_CalcBcom.F90 Log: GRHydro: Make "CalcBcom" compatible with multipatch by enforcing scheduling after local->global trafo of primitives. From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcBcom.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_CalcBcom.F90 2013-07-06 18:10:28 UTC (rev 545) +++ trunk/src/GRHydro_CalcBcom.F90 2013-07-06 18:10:32 UTC (rev 546) @@ -40,10 +40,6 @@ call CCTK_WARN(0,"No storage for GRHydro::bcom") end if - if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then - call CCTK_WARN(0,"MP not yet supported in GRHydro_CalcBcom") - end if - !$OMP PARALLEL DO PRIVATE(i,j,k,velxlow, velylow, velzlow,& !$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum1,dum2,b2,bxlow,bylow,bzlow) Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +31 -1 =================================================================== --- trunk/schedule.ccl 2013-07-06 18:10:28 UTC (rev 545) +++ trunk/schedule.ccl 2013-07-06 18:10:32 UTC (rev 546) @@ -1363,6 +1363,36 @@ } # sync_conserved_only + +#SCHEDULE GRHydro_TestSync in HydroBase_PostStep AFTER HydroBase_Boundaries BEFORE HydroBase_Con2Prim +#{ +# LANG: C +# OPTIONS: LEVEL +#} "Sync test variables" + +#SCHEDULE GROUP Bla in HydroBase_PostStep AFTER HydroBase_Boundaries AFTER HydroBase_Con2Prim +#{ +#} "bla" + +#SCHEDULE GRHydro_TestSync in Bla +#{ +# LANG: C +# OPTIONS: LEVEL +#} "Sync test variables" + +#SCHEDULE GRHydro_TestSync in CCTK_POSTRESTRICT AFTER Refluxing_CorrectState BEFORE MoL_PostStep +#{ +# LANG: C +# OPTIONS: LEVEL +#} "Sync test variables" + +#SCHEDULE GRHydro_TestSync in CCTK_POSTREGRID AFTER Refluxing_Reset +#{ +# LANG: C +# OPTIONS: LEVEL +#} "Sync test variables" + + ############################################################ ### Compute first differences of rho for mesh refinement ### ############################################################ @@ -1509,7 +1539,7 @@ # When Tmunu is not calculated, then the following # routine is used to evaluate bcom, if requested. if (calculate_bcom){ - schedule GRHydro_CalcBcom IN HydroBase_PostStep + schedule GRHydro_CalcBcom IN HydroBase_PostStep AFTER GRHydroTransformPrimToGlobalBasis { LANG: Fortran } "Compute comoving magnetic field, pressure, etc..."

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 545)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_PPM.F90, GRHydro_PPMReconstruct_drv.F90, GRHydro_ReconstructPoly.F90 Log: GRHydro: add option to switch to oPPM on a given refinement level When reconstruction method is enhanced PPM, this change allows using old PPM on one refinement level, whereas the rest of the reflevels will use enhanced PPM. It is only for hydro, but could be extended to MHD also. From: abdik <abdik(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +28 -18 =================================================================== --- trunk/src/GRHydro_PPM.F90 2013-07-06 18:10:24 UTC (rev 544) +++ trunk/src/GRHydro_PPM.F90 2013-07-06 18:10:28 UTC (rev 545) @@ -52,7 +52,7 @@ (iand(mask((i)),(type_bits)).eq.(state_bits)) -subroutine SimplePPM_1d(handle,poly,& +subroutine SimplePPM_1d(apply_enhanced_ppm,handle,poly,& nx,dx,rho,velx,vely,velz,eps,press,rhominus,& velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,& velzplus,epsplus,trivial_rp, hydro_excision_mask,& @@ -68,6 +68,8 @@ CCTK_REAL, parameter :: one = 1 + logical :: apply_enhanced_ppm + CCTK_INT :: handle,poly,nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps @@ -131,7 +133,7 @@ end if - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! This is the original PPM algorithm by Colella & Woodward 1984. !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 @@ -310,11 +312,11 @@ !!$So this is just dropped from eq. (3.2) of C&W. !!$We can get around this by just rescaling the constant k0 (ppm_k0 here). - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then ! Only for 1984 PPM scheme! if (ppm_detect .ne. 0) then - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then ! make sure drho, d2rho and dmrho are computed everywhere! do i = 2, nx-1 drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) @@ -470,7 +472,7 @@ end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then ! In 1984 PPM, flattening is applied before constraining parabolic profiles. if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. do i = 3, nx - 2 @@ -777,7 +779,7 @@ endif - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then do i = 3, nx - 2 @@ -852,7 +854,7 @@ endif - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 ! original Colella&Woodward monotonicity preservation @@ -1019,7 +1021,7 @@ end subroutine SimplePPM_1d !!!! routine for doing PPM to temperature -subroutine SimplePPM_temperature_1d(& +subroutine SimplePPM_temperature_1d(apply_enhanced_ppm,& nx,dx,velx,temperature,press,& tempminus,& tempplus, hydro_excision_mask) @@ -1031,6 +1033,8 @@ DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1.0d0 CCTK_INT :: nx @@ -1063,7 +1067,7 @@ end if - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! This is the original PPM algorithm by Colella & Woodward 1984. !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 @@ -1165,7 +1169,7 @@ do i = 2, nx-1 dpress(i) = press(i+1) - press(i-1) end do - + endif !!$ Zone flattening. See appendix of C&W, p. 197-8. @@ -1187,7 +1191,7 @@ end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then ! In 1984 PPM, flattening is applied before constraining parabolic profiles. if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3. do i = 3, nx - 2 @@ -1439,7 +1443,7 @@ endif - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then do i = 3, nx - 2 @@ -1469,7 +1473,7 @@ endif - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 ! original Colella&Woodward monotonicity preservation @@ -1535,7 +1539,8 @@ -subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & +subroutine SimplePPM_tracer_1d(apply_enhanced_ppm,& + nx,dx,rho,velx,vely,velz, & tracer,tracerminus,tracerplus,press) USE GRHydro_Scalars @@ -1544,6 +1549,8 @@ DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1 CCTK_INT :: nx @@ -1749,7 +1756,8 @@ end subroutine SimplePPM_tracer_1d -subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & +subroutine SimplePPM_ye_1d(apply_enhanced_ppm,& + nx,dx,rho,velx,vely,velz, & Y_e,Y_e_minus,Y_e_plus,press) USE GRHydro_Scalars @@ -1758,6 +1766,8 @@ DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1 CCTK_INT :: nx @@ -1779,7 +1789,7 @@ end do - if (use_enhanced_PPM .eq. 0) then + if (.not.apply_enhanced_PPM) then !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 !!$ This is the expression for an even grid. @@ -1900,7 +1910,7 @@ end if end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then if (PPM3) then do i = 3, nx - 2 flatten = tilde_flatten(i) @@ -2018,7 +2028,7 @@ - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then ! Constrain parabolic profiles if (PPM3) then do i = GRHydro_stencil, nx - GRHydro_stencil + 1 File [modified]: GRHydro_PPMReconstruct_drv.F90 Delta lines: +34 -12 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv.F90 2013-07-06 18:10:24 UTC (rev 544) +++ trunk/src/GRHydro_PPMReconstruct_drv.F90 2013-07-06 18:10:28 UTC (rev 545) @@ -69,6 +69,8 @@ CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + logical :: apply_enhanced_ppm + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then g11 => gaa g12 => gab @@ -118,6 +120,14 @@ ny = cctk_lsh(2) nz = cctk_lsh(3) + ! if use_enhanced_ppm, allow old PPM on one level + if (GRHydro_oppm_reflevel .eq. (-1) .or. & + GRHydro_reflevel .ne. GRHydro_oppm_reflevel) then + apply_enhanced_ppm = use_enhanced_ppm .ne. 0 + else + apply_enhanced_ppm = .false. + end if + !!$ Initialize variables that store reconstructed quantities: !$OMP PARALLEL DO PRIVATE(i,j,k) @@ -166,7 +176,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& @@ -192,7 +203,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_temperature_1d(nx,CCTK_DELTA_SPACE(1),velx(:,j,k), & + call SimplePPM_temperature_1d(apply_enhanced_ppm,& + nx,CCTK_DELTA_SPACE(1),velx(:,j,k), & temperature(:,j,k),press(:,j,k),& tempminus(:,j,k),tempplus(:,j,k), hydro_excision_mask) end do @@ -203,7 +215,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & tracer(:,j,k,:),tracerminus(:,j,k,:),tracerplus(:,j,k,:), & press(:,j,k)) @@ -215,7 +228,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + call SimplePPM_ye_1d(apply_enhanced_ppm,& + nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & press(:,j,k)) @@ -228,7 +242,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& @@ -254,7 +269,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_temperature_1d(ny,CCTK_DELTA_SPACE(2),vely(j,:,k), & + call SimplePPM_temperature_1d(apply_enhanced_ppm,& + ny,CCTK_DELTA_SPACE(2),vely(j,:,k), & temperature(j,:,k),press(j,:,k),& tempminus(j,:,k),tempplus(j,:,k), hydro_excision_mask) end do @@ -265,7 +281,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & vely(j,:,k),velz(j,:,k),velx(j,:,k), & tracer(j,:,k,:),tracerminus(j,:,k,:),tracerplus(j,:,k,:), & press(j,:,k)) @@ -277,7 +294,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + call SimplePPM_ye_1d(apply_enhanced_ppm,& + ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & vely(j,:,k),velz(j,:,k),velx(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & press(j,:,k)) @@ -290,7 +308,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& @@ -316,7 +335,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_temperature_1d(nz,CCTK_DELTA_SPACE(3),velz(j,k,:), & + call SimplePPM_temperature_1d(apply_enhanced_ppm,& + nz,CCTK_DELTA_SPACE(3),velz(j,k,:), & temperature(j,k,:),press(j,k,:),& tempminus(j,k,:),tempplus(j,k,:), hydro_excision_mask) end do @@ -328,7 +348,8 @@ do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velz(j,k,:),velx(j,k,:),vely(j,k,:), & tracer(j,k,:,:),tracerminus(j,k,:,:),tracerplus(j,k,:,:), & press(j,k,:)) @@ -340,7 +361,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + call SimplePPM_ye_1d(apply_enhanced_ppm,& + nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velz(j,k,:),velx(j,k,:),vely(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & press(j,k,:)) File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: +28 -9 =================================================================== --- trunk/src/GRHydro_ReconstructPoly.F90 2013-07-06 18:10:24 UTC (rev 544) +++ trunk/src/GRHydro_ReconstructPoly.F90 2013-07-06 18:10:28 UTC (rev 545) @@ -52,6 +52,8 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + logical :: apply_enhanced_ppm + integer :: nx, ny, nz, i, j, k, itracer logical, dimension(:,:,:), allocatable :: trivial_rp @@ -109,6 +111,14 @@ call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", & &"not_trivial") + ! if use_enhanced_ppm, allow old PPM on one level + if (GRHydro_oppm_reflevel .eq. (-1) .or. & + GRHydro_reflevel .ne. GRHydro_oppm_reflevel) then + apply_enhanced_ppm = use_enhanced_ppm .ne. 0 + else + apply_enhanced_ppm = .false. + end if + nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -329,7 +339,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,1,nx,CCTK_DELTA_SPACE(1),& rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),eps(:,j,k),& press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& @@ -356,7 +367,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & tracer(:,j,k,:),tracerminus(:,j,k,:),tracerplus(:,j,k,:), & press(:,j,k)) @@ -368,7 +380,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & press(:,j,k)) @@ -440,7 +453,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,1,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& @@ -467,7 +481,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & vely(j,:,k),velz(j,:,k),velx(j,:,k), & tracer(j,:,k,:),tracerminus(j,:,k,:),tracerplus(j,:,k,:), & press(j,:,k)) @@ -479,7 +494,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & velx(j,:,k),vely(j,:,k),velz(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & press(j,:,k)) @@ -551,7 +567,8 @@ !$OMP PARALLEL DO PRIVATE(i, j, k) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_1d(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& + call SimplePPM_1d(apply_enhanced_ppm,& + GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),eps(j,k,:),& press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& @@ -579,7 +596,8 @@ do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velz(j,k,:),velx(j,k,:),vely(j,k,:), & tracer(j,k,:,:),tracerminus(j,k,:,:),tracerplus(j,k,:,:), & press(j,k,:)) @@ -591,7 +609,8 @@ !$OMP PARALLEL DO PRIVATE(j, k) do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 - call SimplePPM_tracer_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & + call SimplePPM_tracer_1d(apply_enhanced_ppm,& + nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velx(j,k,:),vely(j,k,:),velz(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & press(j,k,:)) Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2013-07-06 18:10:24 UTC (rev 544) +++ trunk/param.ccl 2013-07-06 18:10:28 UTC (rev 545) @@ -252,6 +252,11 @@ { } no +int GRHydro_oppm_reflevel "Ref level where oPPM is used instead of ePPM (used with use_enhaced_ppm=yes)." +{ + -1:10 :: "0-10 (the reflevel number) or -1 (off)" +} -1 + real enhanced_ppm_C2 "Parameter for enhancecd ppm limiter. Default from McCorquodale & Colella 2011" STEERABLE=ALWAYS { *:* :: "must be greater than 1. According to Colella&Sekora 2008, enhanced ppm is insensitive to C in [1.25,5]"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 544)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_Prim2ConM.F90 Log: GRHydro: Make MHD version of Prim2Con compatible with multipatch. From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +208 -90 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2013-07-06 18:10:18 UTC (rev 543) +++ trunk/src/GRHydro_Prim2ConM.F90 2013-07-06 18:10:24 UTC (rev 544) @@ -13,15 +13,15 @@ #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 velx(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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) +#define Bvecx(i,j,k) Bprim(i,j,k,1) +#define Bvecy(i,j,k) Bprim(i,j,k,2) +#define Bvecz(i,j,k) Bprim(i,j,k,3) #define Bconsx(i,j,k) Bcons(i,j,k,1) #define Bconsy(i,j,k) Bcons(i,j,k,2) #define Bconsz(i,j,k) Bcons(i,j,k,3) @@ -50,56 +50,76 @@ implicit none + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS 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 :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_detr CCTK_REAL :: xtemp(1) character(len=256) NaN_WarnLine - !TARGET gxx, gxy, gxz, gyy, gyz, gzz + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: Bprim - !CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if - !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) + !$OMP g11l,g12l,g13l,g22l,g23l,g33l,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) 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)) + 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(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) - call prim2conM(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + call prim2conM(GRHydro_eos_handle, 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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & @@ -107,7 +127,7 @@ 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)) - call prim2conM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + call prim2conM(GRHydro_eos_handle, 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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & @@ -129,27 +149,27 @@ 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) + !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) 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)) + 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(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + 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 @@ -166,7 +186,7 @@ !$OMP END CRITICAL endif 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, & + i,j,k,x(i,j,k),y(i,j,k),z(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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & @@ -182,7 +202,7 @@ 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, & + i,j,k,x(i,j,k),y(i,j,k),z(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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& @@ -452,13 +472,46 @@ implicit none + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k CCTK_REAL :: det CCTK_REAL :: maxtau0 + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: Bprim + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if + + maxtau0 = -1.0d60 if(evolve_temper.ne.1) then @@ -468,12 +521,12 @@ 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)) + 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 prim2conM(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),& + call prim2conM(GRHydro_eos_handle,g11(i,j,k),& + g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(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),& @@ -509,13 +562,13 @@ 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)) + 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 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),& + g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(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),& @@ -565,39 +618,72 @@ implicit none + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS 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,avg_detl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: Bprim + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if + ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces - !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_detr) 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)) + 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(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) - call prim2conpolytypeM(GRHydro_eos_handle, gxxl,gxyl,gxzl,& - gyyl,gyzl,gzzl, & + call prim2conpolytypeM(GRHydro_eos_handle, 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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k),& @@ -606,8 +692,8 @@ 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)) - call prim2conpolytypeM(GRHydro_eos_handle, gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + call prim2conpolytypeM(GRHydro_eos_handle, 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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& @@ -740,22 +826,54 @@ implicit none + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k CCTK_REAL :: det + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: Bprim + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if !$OMP PARALLEL DO PRIVATE(k,j,i,det) 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)) + 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 prim2conpolytypeM(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),& + call prim2conpolytypeM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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),&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 543)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:10 PM Modified: /trunk/src/ GRHydro_Con2PrimHot.F90 Log: GRHydro: Fixed a problem with excision mask and atmopshere treatment in hot Con2Prim. From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimHot.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_Con2PrimHot.F90 2013-07-06 18:10:16 UTC (rev 542) +++ trunk/src/GRHydro_Con2PrimHot.F90 2013-07-06 18:10:18 UTC (rev 543) @@ -83,7 +83,7 @@ !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& - !$OMP warnline, dummy1, dummy2) + !$OMP warnline, dummy1, dummy2,reset_to_atmo) do k = 1, nz do j = 1, ny do i = 1, nx @@ -125,7 +125,7 @@ reset_to_atmo = 1 endif - if (reset_to_atmo .gt. 0 .or. hydro_excision_mask(i,j,k) .gt. 0) then + if (reset_to_atmo .gt. 0 .or. (GRHydro_enable_internal_excision /= 0 .and. hydro_excision_mask(i,j,k) .gt. 0)) then SET_ATMO_MIN(dens(i,j,k), sqrt(det)*GRHydro_rho_min, r(i,j,k)) SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) scon(i,j,k,:) = 0.d0

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits