/*@@ @file GRHydro_CalcUpdateM.F90 @date Oct 10, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc Calculates the update terms given the fluxes. Moved to here so that @enddesc @@*/ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Functions.h" #include "SpaceMask.h" /*@@ @routine UpdateCalculationM @date Oct 10, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc Calculates the update terms from the fluxes. @enddesc @calls @calledby @history Moved out of the Riemann solver routines to make the FishEye / weighted flux calculation easier. @endhistory @@*/ subroutine UpdateCalculationM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k,itracer CCTK_REAL :: idx, alp_l, alp_r CCTK_INT :: type_bits, atmosphere, not_atmosphere call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& "in_atmosphere") call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& "not_in_atmosphere") idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction) if (CCTK_EQUALS(method_type, "RSA FV")) then if (use_weighted_fluxes == 0) then !$OMP PARALLEL DO PRIVATE(i,j,itracer,alp_l, alp_r) 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 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)) densrhs(i,j,k) = densrhs(i,j,k) + & (alp_l * densflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * densflux(i,j,k)) * idx srhs(i,j,k,1) = srhs(i,j,k,1) + & (alp_l * sxflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * sxflux(i,j,k)) * idx srhs(i,j,k,2) = srhs(i,j,k,2) + & (alp_l * syflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * syflux(i,j,k)) * idx srhs(i,j,k,3) = srhs(i,j,k,3) + & (alp_l * szflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * szflux(i,j,k)) * idx taurhs(i,j,k) = taurhs(i,j,k) + & (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * tauflux(i,j,k)) * idx Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + & (alp_l * Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Bvecxflux(i,j,k)) * idx Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + & (alp_l * Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Bvecyflux(i,j,k)) * idx Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + & (alp_l * Bveczflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Bveczflux(i,j,k)) * idx if(clean_divergence.ne.0) then psidcrhs(i,j,k) = psidcrhs(i,j,k) + & (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * psidcflux(i,j,k)) * idx endif if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers cons_tracerrhs(i,j,k,itracer) = cons_tracerrhs(i,j,k,itracer) + & (alp_l * cons_tracerflux(i-xoffset,j-yoffset,k-zoffset,itracer) - & alp_r * cons_tracerflux(i,j,k,itracer)) * idx enddo end if if (evolve_Y_e .ne. 0) then Y_e_con_rhs(i,j,k) = Y_e_con_rhs(i,j,k) + & (alp_l * Y_e_con_flux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Y_e_con_flux(i,j,k)) * idx end if if (wk_atmosphere .eq. 1) then if ( (atmosphere_mask(i,j,k) .eq. 1) .or. & (SpaceMask_CheckStateBitsF90(space_mask,i,j,k,type_bits,atmosphere)) ) then !!$ We are in the atmosphere so the momentum flux must vanish srhs(i,j,k,:) = 0.d0 if ( ( (atmosphere_mask(i-1,j ,k ) .eq. 1) .and. & (atmosphere_mask(i+1,j ,k ) .eq. 1) .and. & (atmosphere_mask(i ,j-1,k ) .eq. 1) .and. & (atmosphere_mask(i ,j+1,k ) .eq. 1) .and. & (atmosphere_mask(i ,j ,k-1) .eq. 1) .and. & (atmosphere_mask(i ,j ,k+1) .eq. 1) & ) .or. & ( (SpaceMask_CheckStateBitsF90(space_mask,i-1,j ,k ,type_bits,atmosphere)) .and. & (SpaceMask_CheckStateBitsF90(space_mask,i+1,j ,k ,type_bits,atmosphere)) .and. & (SpaceMask_CheckStateBitsF90(space_mask,i ,j-1,k ,type_bits,atmosphere)) .and. & (SpaceMask_CheckStateBitsF90(space_mask,i ,j+1,k ,type_bits,atmosphere)) .and. & (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k-1,type_bits,atmosphere)) .and. & (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k+1,type_bits,atmosphere)) & ) & ) then !!$ All neighbours are also atmosphere so all rhs vanish densrhs(i,j,k) = 0.d0 taurhs(i,j,k) = 0.d0 !!$ !!$ We should still evolve the B-field in the atmosphere !!$ end if end if end if enddo enddo enddo !$OMP END PARALLEL DO else call CCTK_WARN(0, "Not supported") !!$ 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 !!$ !!$ 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)) !!$ !!$ densrhs(i,j,k) = densrhs(i,j,k) + & !!$ (alp_l * & !!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * & !!$ &densflux(i-xoffset,j-yoffset,k-zoffset) - & !!$ alp_r * & !!$ &cell_surface(i,j,k,flux_direction) * & !!$ &densflux(i,j,k)) * idx / cell_volume(i,j,k) !!$ sxrhs(i,j,k) = sxrhs(i,j,k) + & !!$ (alp_l * & !!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * & !!$ &sxflux(i-xoffset,j-yoffset,k-zoffset) - & !!$ alp_r * & !!$ &cell_surface(i,j,k,flux_direction) * & !!$ &sxflux(i,j,k)) * idx / cell_volume(i,j,k) !!$ syrhs(i,j,k) = syrhs(i,j,k) + & !!$ (alp_l * & !!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * & !!$ &syflux(i-xoffset,j-yoffset,k-zoffset) - & !!$ alp_r * & !!$ &cell_surface(i,j,k,flux_direction) * & !!$ &syflux(i,j,k)) * idx / cell_volume(i,j,k) !!$ szrhs(i,j,k) = szrhs(i,j,k) + & !!$ (alp_l * & !!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * & !!$ &szflux(i-xoffset,j-yoffset,k-zoffset) - & !!$ alp_r * & !!$ &cell_surface(i,j,k,flux_direction) * & !!$ &szflux(i,j,k)) * idx / cell_volume(i,j,k) !!$ taurhs(i,j,k) = taurhs(i,j,k) + & !!$ (alp_l * & !!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * & !!$ &tauflux(i-xoffset,j-yoffset,k-zoffset) - & !!$ alp_r * & !!$ &cell_surface(i,j,k,flux_direction) * & !!$ &tauflux(i,j,k)) * idx / cell_volume(i,j,k) !!$ !!$ enddo !!$ enddo !!$ enddo end if else if (CCTK_EQUALS(method_type, "Flux split FD")) then 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 densrhs(i,j,k) = densrhs(i,j,k) + & (densflux(i-xoffset,j-yoffset,k-zoffset) - & densflux(i,j,k)) * idx srhs(i,j,k,1) = srhs(i,j,k,1) + & (sxflux(i-xoffset,j-yoffset,k-zoffset) - & sxflux(i,j,k)) * idx srhs(i,j,k,2) = srhs(i,j,k,2) + & (syflux(i-xoffset,j-yoffset,k-zoffset) - & syflux(i,j,k)) * idx srhs(i,j,k,3) = srhs(i,j,k,3) + & (szflux(i-xoffset,j-yoffset,k-zoffset) - & szflux(i,j,k)) * idx taurhs(i,j,k) = taurhs(i,j,k) + & (tauflux(i-xoffset,j-yoffset,k-zoffset) - & tauflux(i,j,k)) * idx Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + & (Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - & Bvecxflux(i,j,k)) * idx Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + & (Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - & Bvecyflux(i,j,k)) * idx Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + & (Bveczflux(i-xoffset,j-yoffset,k-zoffset) - & Bveczflux(i,j,k)) * idx if(clean_divergence.ne.0) then psidcrhs(i,j,k) = psidcrhs(i,j,k) + & (psidcflux(i-xoffset,j-yoffset,k-zoffset) - & psidcflux(i,j,k)) * idx endif enddo enddo enddo end if return end subroutine UpdateCalculationM