diff options
Diffstat (limited to 'src/GRHydro_CalcUpdate.F90')
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 64 |
1 files changed, 40 insertions, 24 deletions
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 |