User: rhaas Date: 2011/11/29 11:17 AM
Modified: /trunk/src/ GRHydro_CalcUpdate.F90, GRHydro_ENOReconstruct_drv.F90, GRHydro_HLLEM.F90, GRHydro_PPMMReconstruct_drv.F90, GRHydro_Prim2ConM.F90, GRHydro_TVDReconstruct.F90
Log: expand region where fluxes are computed for constraint transport
File Changes:
Directory: /trunk/src/ ======================
File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_CalcUpdate.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -48,7 +48,7 @@ if (use_weighted_fluxes == 0) then
!$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r) - do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! mean we need to compute fluxes for one more point the GRHydro_stencil would indicate + 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
File [modified]: GRHydro_ENOReconstruct_drv.F90 Delta lines: +12 -6 =================================================================== --- trunk/src/GRHydro_ENOReconstruct_drv.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_ENOReconstruct_drv.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -179,9 +179,11 @@ endif
if (flux_direction == 1) then + ! 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) - do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 - do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),& rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& @@ -259,9 +261,11 @@ end do !$OMP END PARALLEL DO else if (flux_direction == 2) then + ! 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) - do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 - do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),& rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& @@ -339,9 +343,11 @@ end do !$OMP END PARALLEL DO else if (flux_direction == 3) then + ! 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) - do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 - do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints if (CCTK_EQUALS(recon_vars,"primitive")) then call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),& rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),&
File [modified]: GRHydro_HLLEM.F90 Delta lines: +5 -3 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_HLLEM.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -74,9 +74,11 @@ call CCTK_WARN(0, "Flux direction not x,y,z") end if
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil - do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset)
f1 = 0.d0 lamminus = 0.d0
File [modified]: GRHydro_PPMMReconstruct_drv.F90 Delta lines: +12 -6 =================================================================== --- trunk/src/GRHydro_PPMMReconstruct_drv.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_PPMMReconstruct_drv.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -181,9 +181,11 @@
!!$ PPM starts: if (flux_direction == 1) then + ! 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) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints if(clean_divergence.ne.0) then call SimplePPM_1dM(GRHydro_eos_handle,0,nx,CCTK_DELTA_SPACE(1),& rho(:,j,k),velx(:,j,k),vely(:,j,k),velz(:,j,k),& @@ -252,9 +254,11 @@ end if
else if (flux_direction == 2) then + ! 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) - do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints if(clean_divergence.ne.0) then call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& @@ -322,9 +326,11 @@ end if
else if (flux_direction == 3) then + ! 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) - do k = GRHydro_stencil, ny - GRHydro_stencil + 1 - do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints if(clean_divergence.ne.0) then call SimplePPM_1dM(GRHydro_eos_handle,0,nz,CCTK_DELTA_SPACE(3),& rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),&
File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +10 -6 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_Prim2ConM.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -57,9 +57,11 @@ CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
- 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 + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + 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)) @@ -275,11 +277,13 @@ CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+ ! 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) - 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 + 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))
File [modified]: GRHydro_TVDReconstruct.F90 Delta lines: +5 -3 =================================================================== --- trunk/src/GRHydro_TVDReconstruct.F90 2011-11-29 17:17:04 UTC (rev 306) +++ trunk/src/GRHydro_TVDReconstruct.F90 2011-11-29 17:17:07 UTC (rev 307) @@ -52,10 +52,12 @@ !!$ Initially all Riemann problems are NON-trivial
trivial_rp = .false. + ! 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,dupw,dloc,delta,ratio,hdelta) - do k = GRHydro_stencil, nz-GRHydro_stencil+1 - do j = GRHydro_stencil, ny-GRHydro_stencil+1 - do i = GRHydro_stencil, nx-GRHydro_stencil+1 + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints*(1-xoffset) if (GRHydro_enable_internal_excision /= 0 .and. & (hydro_excision_mask(i,j,k) .ne. 0)) then trivial_rp(i-xoffset, j-yoffset, k-zoffset) = .true.