Commits May 2012

commits@lists.einsteintoolkit.org

11 participants
183 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 332)
by roland.haas＠physics.gatech.edu 14 May '12

14 May '12

User: rhaas Date: 2012/05/14 10:37 AM Modified: /trunk/src/ GRHydro_Reconstruct.F90 Log: GRHydro: correct atmosphere handling after primitve reconstruction I have fixed a bug in GRHydro's atmosphere handling after reconstruction. After reconstruction, the reconstructed plus and minus face values are tested for whether they drop below atmosphere level and are then reset. In particular, this means that plus and minus face values for cell i don't generally coincide with the minus/plus face values of the corresponding cells i-1 and i+1. If the reconstruction previously found that the Riemann problem was trivial between, say rhominus(i+1) and rhoplus(i), when changing rhoplus(i), it is not anymore! However, the mask "trivial_rp" testing for a trivial Riemann problem is not changed! This leads to inconsistent behavior at the boundary of the atmosphere. The symptom was a drift in the center of mass of a static and perfectly symmetric TOV star. Patch by Christian Reisswig. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Reconstruct.F90 Delta lines: +33 -8 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2012-05-14 15:31:15 UTC (rev 331) +++ trunk/src/GRHydro_Reconstruct.F90 2012-05-14 15:37:35 UTC (rev 332) @@ -44,6 +44,7 @@ CCTK_INT :: i,j,k CCTK_REAL :: local_min_tracer + CCTK_INT :: type_bits, not_trivial ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -83,12 +84,30 @@ end if end if + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemX", & + &"not_trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemY", & + &"not_trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(not_trivial, "Hydro_RiemannProblemZ", & + &"not_trivial") + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) - if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & - rhominus(i,j,k).lt.GRHydro_rho_min) then + 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(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & + rhominus(i,j,k).lt.GRHydro_rho_min) then + rhoplus(i,j,k) = rho(i,j,k) rhominus(i,j,k) = rho(i,j,k) epsplus(i,j,k) = eps(i,j,k) @@ -107,11 +126,16 @@ tracerminus(i,j,k,:) = tracer(i,j,k,:) end where end if - enddo - enddo + ! Riemann problem might not be trivial anymore!!!! + SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial) + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial) + enddo enddo - !$OMP END PARALLEL DO + enddo + !$OMP END PARALLEL DO + + if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then @@ -133,3 +157,4 @@ end subroutine Reconstruction +

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 331)
by roland.haas＠physics.gatech.edu 14 May '12

14 May '12

User: rhaas Date: 2012/05/14 10:31 AM Added: /trunk/src/ GRHydro_Analysis.F90 Log: GRHydro: add file missing from r330 File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_Analysis.F90 Delta lines: +91 -0 =================================================================== --- trunk/src/GRHydro_Analysis.F90 (rev 0) +++ trunk/src/GRHydro_Analysis.F90 2012-05-14 15:31:15 UTC (rev 331) @@ -0,0 +1,91 @@ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" +#include "SpaceMask.h" + +subroutine GRHydro_Analysis_Init(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i,j,k + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + divB(i,j,k) = 0.0d0 + end do + end do + end do + !$OMP END PARALLEL DO + +end subroutine GRHydro_Analysis_Init + +subroutine GRHydro_CalcDivB(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i,j,k,itracer + CCTK_REAL :: idx, idy, idz, Bcons_l1, Bcons_r1, Bcons_l2, Bcons_r2, Bcons_l3, Bcons_r3 + + idx = 1.d0 / CCTK_DELTA_SPACE(1) + idy = 1.d0 / CCTK_DELTA_SPACE(2) + idz = 1.d0 / CCTK_DELTA_SPACE(3) + + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,Bcons_l1,Bcons_r1,Bcons_l2,Bcons_r2,Bcons_l3,Bcons_r3) + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + + if(evolve_mhd.ne.0) then + if(track_divB.ne.0) then + if(transport_constraints.ne.0) then + ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) + divB(i,j,k) = divB(i,j,k) + & + 0.25d0*(Bcons(i+1,j,k,1)-Bcons(i,j,k,1)+ & + Bcons(i+1,j+1,k,1)-Bcons(i,j+1,k,1)+ & + Bcons(i+1,j,k+1,1)-Bcons(i,j,k+1,1)+ & + Bcons(i+1,j+1,k+1,1)-Bcons(i,j+1,k+1,1))*idx + & + 0.25d0*(Bcons(i,j+1,k,2)-Bcons(i,j,k,2)+ & + Bcons(i+1,j+1,k,2)-Bcons(i+1,j,k,2)+ & + Bcons(i,j+1,k+1,2)-Bcons(i,j,k+1,2)+ & + Bcons(i+1,j+1,k+1,2)-Bcons(i+1,j,k+1,2))*idy + & + 0.25d0*(Bcons(i,j,k+1,3)-Bcons(i,j,k,3)+ & + Bcons(i+1,j,k+1,3)-Bcons(i+1,j,k,3)+ & + Bcons(i,j+1,k+1,3)-Bcons(i,j+1,k,3)+ & + Bcons(i+1,j+1,k+1,3)-Bcons(i+1,j+1,k,3))*idz + else + Bcons_l1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i-1,j,k,1)) + Bcons_l2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j-1,k,2)) + Bcons_l3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k-1,3)) + Bcons_r1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i+1,j,k,1)) + Bcons_r2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j+1,k,2)) + Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k+1,3)) + !write(*,*) "Bcons_r1 = ", Bcons_r1 + divB(i,j,k) = divB(i,j,k) + (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + !divB(i,j,k) = 0.0d0 + endif + endif + endif + + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end subroutine GRHydro_CalcDivB

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 330)
by roland.haas＠physics.gatech.edu 14 May '12

14 May '12

User: rhaas Date: 2012/05/14 10:04 AM Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_CalcUpdate.F90, GRHydro_HLLEM.F90, GRHydro_SourceM.F90, make.code.defn Log: GRHydro: fixes for GRMHD All by Philipp Moesta. 1) Fix parity of psidc and divb 2) Fix a wrong index in the source terms of scon 3) Fix wrong indices of derivatives of space-time metric in the source of the divergence cleaning scalar. 4) Calculate divergence of B in MoL PseudoEvolution and set its Prolongation="Restrict". 5) Correct the source terms and fluxes for the Bfield and the divergence cleaning field when having a non-flat space-time. 6) Make sure alpha factors match between UpdateCalculation and fluxes definition. 7) Include 1/sqrt(detg) factor in calculation of \epsilon^{\muijk} in the cross product to obtain the Bfield form the vector potential. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +12 -12 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/src/GRHydro_CalcUpdate.F90 2012-05-14 15:04:12 UTC (rev 330) @@ -39,7 +39,7 @@ DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k,itracer - CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r, alp_tmp + CCTK_REAL :: idx, alp_l, alp_r, Bcons_l, Bcons_r, alp_tmp idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction) @@ -47,7 +47,7 @@ if (use_weighted_fluxes == 0) then - !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,alp_tmp,Bvec_l,Bvec_r) + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,alp_tmp,Bcons_l,Bcons_r) do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil @@ -119,11 +119,11 @@ Bcons(i+xoffset,j+1-xoffset,k+1 ,flux_direction)-Bcons(i ,j+zoffset ,k+1-zoffset,flux_direction)+ & Bcons(i+1 ,j+1 ,k+1 ,flux_direction)-Bcons(i+1-xoffset,j+1-yoffset,k+1-zoffset,flux_direction))*idx else - Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) - Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) - divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx + Bcons_l = 0.5d0 * (Bcons(i,j,k,flux_direction) + & + Bcons(i-xoffset,j-yoffset,k-zoffset,flux_direction)) + Bcons_r = 0.5d0 * (Bcons(i,j,k,flux_direction) + & + Bcons(i+xoffset,j+yoffset,k+zoffset,flux_direction)) + divB(i,j,k) = divB(i,j,k) + (Bcons_l - Bcons_r ) * idx endif endif @@ -276,11 +276,11 @@ psidcflux(i,j,k)) * idx endif if(track_divB.ne.0) then - Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) - Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) - divB(i,j,k) = divB(i,j,k) + ( Bvec_l - Bvec_r ) * idx + Bcons_l = 0.5d0 * (Bcons(i,j,k,flux_direction) + & + Bcons(i-xoffset,j-yoffset,k-zoffset,flux_direction)) + Bcons_r = 0.5d0 * (Bcons(i,j,k,flux_direction) + & + Bcons(i+xoffset,j+yoffset,k+zoffset,flux_direction)) + divB(i,j,k) = divB(i,j,k) + ( Bcons_l - Bcons_r ) * idx endif endif File [modified]: GRHydro_HLLEM.F90 Delta lines: +37 -36 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/src/GRHydro_HLLEM.F90 2012-05-14 15:04:12 UTC (rev 330) @@ -226,10 +226,10 @@ vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz - psidcf = avg_alp*cons_p(6)/sdet-psidcp*avg_betax + f1(6)=f1(6) + 1.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + psidcf = cons_p(6)/sdet-psidcp*avg_betax/avg_alp endif else if (flux_direction == 2) then @@ -239,10 +239,10 @@ vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz - psidcf = avg_alp*cons_p(7)/sdet-psidcp*avg_betay + f1(6)=f1(6) + 1.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + psidcf = cons_p(7)/sdet-psidcp*avg_betay/avg_alp endif else if (flux_direction == 3) then @@ -252,10 +252,10 @@ vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz - psidcf = avg_alp*cons_p(8)/sdet-psidcp*avg_betaz + f1(6)=f1(6) + 1.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + psidcf = cons_p(8)/sdet-psidcp*avg_betaz/avg_alp endif else @@ -315,14 +315,14 @@ avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - cons_m(6)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - cons_m(6)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - cons_m(6)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz - psidcfp = avg_alp*cons_p(6)/sdet-avg_betax*psidcp - psidcfm = avg_alp*cons_m(6)/sdet-avg_betax*psidcm + fminus(6)=fminus(6) + 1.0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + psidcfp = cons_p(6)/sdet-avg_betax*psidcp/avg_alp + psidcfm = cons_m(6)/sdet-avg_betax*psidcm/avg_alp endif else if (flux_direction == 2) then @@ -349,14 +349,14 @@ avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - cons_m(7)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - cons_m(7)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - cons_m(7)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz - psidcfp = avg_alp*cons_p(7)/sdet-avg_betay*psidcp - psidcfm = avg_alp*cons_m(7)/sdet-avg_betay*psidcm + fminus(6)=fminus(6) + 1.0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + psidcfp = cons_p(7)/sdet-avg_betay*psidcp/avg_alp + psidcfm = cons_m(7)/sdet-avg_betay*psidcm/avg_alp endif else if (flux_direction == 3) then @@ -383,14 +383,14 @@ avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - cons_m(8)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - cons_m(8)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - cons_m(8)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz - psidcfp = avg_alp*cons_p(8)/sdet-avg_betaz*psidcp - psidcfm = avg_alp*cons_m(8)/sdet-avg_betaz*psidcm + fminus(6)=fminus(6) + 1.0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + psidcfp = cons_p(8)/sdet-avg_betaz*psidcp/avg_alp + psidcfm = cons_m(8)/sdet-avg_betaz*psidcm/avg_alp endif else @@ -428,6 +428,7 @@ end do if(clean_divergence.ne.0) then + psidcdiff = psidcm - psidcp select case(whichpsidcspeed) case(0) File [modified]: GRHydro_SourceM.F90 Delta lines: +20 -10 =================================================================== --- trunk/src/GRHydro_SourceM.F90 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/src/GRHydro_SourceM.F90 2012-05-14 15:04:12 UTC (rev 330) @@ -395,7 +395,10 @@ halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +& txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz - sx_source = t00*& + + + + sx_source = t00*& (halfshiftdgx - alp(i,j,k)*dx_alp) + halfTdgx + & t0x*(shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) +& t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) +& @@ -403,7 +406,7 @@ rhohstarW2*invalp*(vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz) -& bt*(bxlow*dx_betax + bylow*dx_betay + bzlow*dx_betaz) - sy_source = t00*& + sy_source = t00*& (halfshiftdgy - alp(i,j,k)*dy_alp) + halfTdgy + & t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) +& t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) +& @@ -411,8 +414,8 @@ rhohstarW2*invalp*(vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz) -& bt*(bxlow*dy_betax + bylow*dy_betay + bzlow*dy_betaz) - sz_source = t00*& - (halfshiftdgz - alp(i,j,k)*dz_alp) + halfTdgy + & + sz_source = t00*& + (halfshiftdgz - alp(i,j,k)*dz_alp) + halfTdgz + & t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) +& t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) +& t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) +& @@ -423,20 +426,21 @@ srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source - taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then + psidcrhs(i,j,k) = -1.d0 * (kap_dc*alp(i,j,k) + & dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & - Bconsx(i,j,k) * ( dx_alp - half*alp(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 + & + 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/sqrtdet + & Bconsy(i,j,k) * (dy_alp - half*alp(i,j,k) * & ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & - 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + & + 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/sqrtdet + & Bconsz(i,j,k) * (dz_alp - half*alp(i,j,k) * & - ( uxx*dz_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & - 2.d0*uxz*dz_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + ( uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + 2.d0*uxy*dz_gxy + & + 2.d0*uxz*dz_gxz + 2.d0*uyz*dz_gyz ) )/sqrtdet !!$ g^{jk} d_i g_{kj} = d_i (g) / det dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + & @@ -511,6 +515,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -525,6 +530,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -539,6 +545,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -553,6 +560,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -567,6 +575,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do @@ -581,6 +590,7 @@ Bcons(i,j,k,1) = -1d100 Bcons(i,j,k,2) = -1d100 Bcons(i,j,k,3) = -1d100 + psidc(i,j,k) = -1d100 end do end do end do File [modified]: make.code.defn Delta lines: +1 -0 =================================================================== --- trunk/src/make.code.defn 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/src/make.code.defn 2012-05-14 15:04:12 UTC (rev 330) @@ -4,6 +4,7 @@ # Source files in this directory SRCS = Utils.F90 \ + GRHydro_Analysis.F90 \ GRHydro_Boundaries.F90 \ GRHydro_CalcBcom.F90 \ GRHydro_CalcUpdate.F90 \ Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +2 -2 =================================================================== --- trunk/interface.ccl 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/interface.ccl 2012-05-14 15:04:12 UTC (rev 330) @@ -358,7 +358,7 @@ #real w_lorentz type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Lorentz factor" -real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Psi parameter for divergence cleaning" +real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 tensorparity=-1 jacobian="inverse_jacobian" interpolator="matter"' "Psi parameter for divergence cleaning" real densrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for dens" real taurhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for tau" @@ -367,7 +367,7 @@ real psidcrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for psidc" -real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint" +real divB type = GF Timelevels = 1 tags='Prolongation="Restrict" checkpoint="no" tensorparity=-1' "Magnetic field constraint" real bcom[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" tensorparity=-1 interpolator="matter"' "b^i: comoving contravariant magnetic field 4-vector spatial components" real bcom0 type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "b^0 component of the comoving contravariant magnetic field 4-vector" File [modified]: schedule.ccl Delta lines: +19 -0 =================================================================== --- trunk/schedule.ccl 2012-05-13 19:22:38 UTC (rev 329) +++ trunk/schedule.ccl 2012-05-14 15:04:12 UTC (rev 330) @@ -1162,3 +1162,22 @@ { } "Calculate the stress-energy tensor" +# Calculate divergence in PseudoEvolution + +if(track_divB) +{ + # bit of a misnomer. Currently only contains divB + schedule group GRHydroAnalysis IN MoL_PseudoEvolution + { + } "Calculate analysis quantities" + + schedule GRHydro_Analysis_Init IN GRHydroAnalysis + { + LANG: Fortran + } "Initialize variables" + + schedule GRHydro_CalcDivB IN GRHydroAnalysis AFTER GRHydro_Analysis_Init + { + LANG: Fortran + } "Calculate divB" +}

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 180)
by knarf＠cct.lsu.edu 14 May '12

14 May '12

User: knarf Date: 2012/05/14 09:05 AM Modified: /trunk/ einsteintoolkit.th Log: add ADMMass File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +1 -0 =================================================================== --- trunk/einsteintoolkit.th 2012-05-10 15:52:47 UTC (rev 179) +++ trunk/einsteintoolkit.th 2012-05-14 14:05:23 UTC (rev 180) @@ -268,6 +268,7 @@ !AUTH_URL = https://svn.aei.mpg.de/numrel/$1/$2/trunk !URL = http://svn.aei.mpg.de/numrel/$1/$2/trunk !CHECKOUT = +AEIThorns/ADMMass AEIThorns/AEILocalInterp AEIThorns/PunctureTracker AEIThorns/SystemStatistics

1 0

New changesets in Carpet Mercurial
by Carpet mercurial repositories 13 May '12

13 May '12

changeset: 3551:0749f85a5257 user: Roland Haas <rhaas(a)caltech.edu> date: Sun Apr 22 19:13:11 2012 -0700 summary: CarpetReduce: output calling routine when warning about reductions changeset: 3552:47e0c5bd5d07 user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Thu May 10 16:49:18 2012 -0700 summary: CarpetIOScalar: correct column labels changeset: 3553:82b61aa7e667 tag: tip user: Roland Haas <roland.haas(a)physics.gatech.edu> date: Thu May 10 17:24:01 2012 -0700 summary: CarpetIOScalar: add test cases for all combinations of

1 0

New changesets in Carpet Mercurial
by Carpet mercurial repositories 13 May '12

13 May '12

changeset: 3550:636bbb1a96ff tag: tip user: Christian Reisswig <reisswig(a)tapir.caltech.edu> date: Sun May 13 21:11:59 2012 -0700 summary: CarpetLib: Fill overlap zones via prolongation during regrid.

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 329)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:22 PM Modified: /trunk/ interface.ccl Log: GRHydro: correct vector type of Bcons to match Bvec File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +2 -2 =================================================================== --- trunk/interface.ccl 2012-05-13 19:16:05 UTC (rev 328) +++ trunk/interface.ccl 2012-05-13 19:22:38 UTC (rev 329) @@ -345,9 +345,9 @@ real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "generalized momenta" -real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorparity=-1 tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable" +real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" tensorparity=-1 tensorweight=+1.0 jacobian="jacobian" interpolator="matter"' "B-field conservative variable" -real Evec[3] type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Electric field at edges" +real Evec[3] type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" tensorweight=+1.0 jacobian="jacobian" interpolator="matter"' "Electric field at edges" real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 328)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:16 PM Modified: /trunk/src/ GRHydro_RegisterVars.cc Log: GRHydro: Correct syntax error Detected by gcc 4.7. Patch by Erik Schnetter File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_RegisterVars.cc Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_RegisterVars.cc 2012-05-13 19:13:54 UTC (rev 327) +++ trunk/src/GRHydro_RegisterVars.cc 2012-05-13 19:16:05 UTC (rev 328) @@ -34,7 +34,7 @@ // Main function called by Cactus to register variables with MoL -extern "C"void GRHydro_Register(CCTK_ARGUMENTS) +extern "C" void GRHydro_Register(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS;

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 327)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:13 PM Modified: /trunk/src/ GRHydro_CalcUpdate.F90 Log: GRHydro: OpenMP parallelize constraint transport loop File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +2 -0 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2012-05-13 19:13:31 UTC (rev 326) +++ trunk/src/GRHydro_CalcUpdate.F90 2012-05-13 19:13:54 UTC (rev 327) @@ -293,6 +293,7 @@ if (transport_constraints.ne.0 .and. flux_direction.eq.1) then ! HACK: x direction is last ! FIXME: I think one could wrap all of this into a single do loop and remove the ! Evec storage + !$OMP PARALLEL DO PRIVATE(i,j,k) do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil @@ -311,6 +312,7 @@ enddo enddo enddo + !$OMP END PARALLEL DO end if return

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 326)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:13 PM Modified: /trunk/src/ GRHydro_CalcUpdate.F90 Log: GRHydro: fix OpenMP data race condition in UpdateCalculation for constraint transport File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2012-05-13 19:12:48 UTC (rev 325) +++ trunk/src/GRHydro_CalcUpdate.F90 2012-05-13 19:13:31 UTC (rev 326) @@ -47,7 +47,7 @@ if (use_weighted_fluxes == 0) then - !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r) + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,alp_tmp,Bvec_l,Bvec_r) do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil

1 0

← Newer
1
...
8
9
10
11
12
13
14
...
19
Older →

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits May 2012