From 531b43c566ac6a5eebc091201dd9d29193976d2e Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 13 Aug 2013 14:56:21 +0000 Subject: Move Evec calculation in seperate loops this avoids accessing unitinialized values in the ghost zones when computing the hydro varaible updates. This also saves time in running those loops since they are smaller. From: Philipp Moesta git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@571 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_CalcUpdate.F90 | 64 +++++++++++++++++++++++++++++----------------- 1 file changed, 40 insertions(+), 24 deletions(-) (limited to 'src') diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index 3b08360..eee7d14 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -47,10 +47,46 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) if (use_weighted_fluxes == 0) then + if(evolve_mhd.ne.0 .and. transport_constraints.ne.0) then + !$OMP PARALLEL DO PRIVATE(i,j,k,alp_l,alp_r,alp_tmp) + 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 + + alp_l = 0.5d0 * (alp(i,j,k) + & + alp(i-xoffset,j-yoffset,k-zoffset)) + alp_r = 0.5d0 * (alp(i,j,k) + & + alp(i+xoffset,j+yoffset,k+zoffset)) + + ! we have to first compute all components of v\crossB = E and + ! combine them in the last substep into Bconshs + ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) + if(flux_direction.eq.1) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + elseif(flux_direction.eq.2) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + elseif(flux_direction.eq.3) then + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + end if + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif + !$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 + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil alp_l = 0.5d0 * (alp(i,j,k) + & alp(i-xoffset,j-yoffset,k-zoffset)) @@ -74,27 +110,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) alp_r * tauflux(i,j,k)) * idx if(evolve_mhd.ne.0) then - if(transport_constraints.ne.0) then - ! we have to first compute all components of v\crossB = E and - ! combine them in the last substep into Bconshs - ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) - if(flux_direction.eq.1) then - alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) - Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) - alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) - Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) - elseif(flux_direction.eq.2) then - alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) - Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) - alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) - Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) - elseif(flux_direction.eq.3) then - alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) - Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) - alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) - Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) - end if - else + if(transport_constraints.ne.1) then Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Bconsxflux(i,j,k)) * idx -- cgit v1.2.3