[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 571)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Tue Aug 13 09:56:21 CDT 2013
User: rhaas
Date: 2013/08/13 09:56 AM
Modified:
/trunk/src/
GRHydro_CalcUpdate.F90
Log:
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 <pmoesta at tapir.caltech.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_CalcUpdate.F90
Delta lines: +40 -24
===================================================================
--- trunk/src/GRHydro_CalcUpdate.F90 2013-08-13 14:56:16 UTC (rev 570)
+++ trunk/src/GRHydro_CalcUpdate.F90 2013-08-13 14:56:21 UTC (rev 571)
@@ -47,10 +47,46 @@
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 @@
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
More information about the Commits
mailing list