From 74fb1e6ea34d6e03a35ff6c158f455c39904bf5a Mon Sep 17 00:00:00 2001 From: bmundim Date: Sun, 2 May 2010 20:59:32 +0000 Subject: file/parameter string replacement from whisky to GRHydro git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@112 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_HLLE.F90 | 576 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 576 insertions(+) create mode 100644 src/GRHydro_HLLE.F90 (limited to 'src/GRHydro_HLLE.F90') diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90 new file mode 100644 index 0000000..d32f4a7 --- /dev/null +++ b/src/GRHydro_HLLE.F90 @@ -0,0 +1,576 @@ + /*@@ + @file GRHydro_HLLE.F90 + @date Sat Jan 26 01:40:14 2002 + @author Ian Hawke, Pedro Montero, Toni Font + @desc + The HLLE solver. Called from the wrapper function, so works in + all directions. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include "SpaceMask.h" + + /*@@ + @routine GRHydro_HLLEGeneral + @date Sat Jan 26 01:41:02 2002 + @author Ian Hawke, Pedro Montero, Toni Font + @desc + The HLLE solver. Sufficiently simple that its just one big routine. + Rewritten for the new EOS interface. + @enddesc + @calls + @calledby + @history + Altered from Cactus 3 routines originally written by Toni Font. + @endhistory + +@@*/ + +subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) + + USE GRHydro_Eigenproblem + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: i, j, k, m + CCTK_REAL, dimension(5) :: cons_p,cons_m,fplus,fminus,lamplus + CCTK_REAL, dimension(5) :: f1,lamminus + CCTK_REAL, dimension(5) :: qdiff + CCTK_REAL, dimension(6) :: prim_p, prim_m + CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det + CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, psi4h, & + cs2_p, cs2_m, dpdeps_p, dpdeps_m + + CCTK_INT :: type_bits, trivial, not_trivial + + integer tadmor + + if(CCTK_EQUALS(HLLE_type,"Tadmor")) then + tadmor = 1 + else + tadmor = 0 + endif + + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & + &"trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", & + &"trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & + &"trivial") + else + 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 + + f1 = 0.d0 + lamminus = 0.d0 + lamplus = 0.d0 + cons_p = 0.d0 + cons_m = 0.d0 + fplus = 0.d0 + fminus = 0.d0 + qdiff = 0.d0 + +!!$ Set the left (p for plus) and right (m_i for minus, i+1) states + + cons_p(1) = densplus(i,j,k) + cons_p(2) = sxplus(i,j,k) + cons_p(3) = syplus(i,j,k) + cons_p(4) = szplus(i,j,k) + cons_p(5) = tauplus(i,j,k) + + cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) + + prim_p(1) = rhoplus(i,j,k) + prim_p(2) = velxplus(i,j,k) + prim_p(3) = velyplus(i,j,k) + prim_p(4) = velzplus(i,j,k) + prim_p(5) = epsplus(i,j,k) + + prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) + prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) + + prim_p(6) = pressplus(i,j,k) + cs2_p = eos_cs2_p(i,j,k) + dpdeps_p = eos_dpdeps_p(i,j,k) + + prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) + cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset) + dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset) + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. + + if (shift_state .ne. 0) then + if (flux_direction == 1) then + avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & + betax(i,j,k)) + else if (flux_direction == 2) then + avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & + betay(i,j,k)) + else if (flux_direction == 3) then + avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & + betaz(i,j,k)) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + else + avg_beta = 0.d0 + end if + + avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) + + gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) + gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) + gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) + gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) + gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) + gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) + + if (conformal_state .eq. 0) then + psi4h = 1.d0 + else + psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 + end if + + call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det) + +!!$ If the Riemann problem is trivial, just calculate the fluxes from the +!!$ left state and skip to the next cell + + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then + + if (flux_direction == 1) then + call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& + f1(1),f1(2),f1(3),f1(4),f1(5),& + prim_m(2),prim_m(6), & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& + f1(1),f1(3),f1(4),f1(2),f1(5),& + prim_m(3),prim_m(6), & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& + f1(1),f1(4),f1(2),f1(3),f1(5), & + prim_m(4),prim_m(6), & + avg_det,avg_alp,avg_beta) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + else !!! The end of this branch is right at the bottom of the routine + + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, & + psi4h*gyyh, psi4h*gyzh, psi4h*gzzh) + + if (flux_direction == 1) then + usendh = uxxh + else if (flux_direction == 2) then + usendh = uyyh + else if (flux_direction == 3) then + usendh = uzzh + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + +!!$ Calculate the jumps in the conserved variables + + qdiff(1) = cons_m(1) - cons_p(1) + qdiff(2) = cons_m(2) - cons_p(2) + qdiff(3) = cons_m(3) - cons_p(3) + qdiff(4) = cons_m(4) - cons_p(4) + qdiff(5) = cons_m(5) - cons_p(5) + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + call eigenvalues_general(& + prim_m(2),prim_m(3),prim_m(4),cs2_m, & + lamminus,& + psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(2),prim_p(3),prim_p(4),cs2_p, & + lamplus,& + psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + usendh,avg_alp,avg_beta) + call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& + fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),& + prim_p(2),prim_p(6), & + avg_det,avg_alp,avg_beta) + call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& + fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& + prim_m(2),prim_m(6), & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call eigenvalues_general(& + prim_m(3),prim_m(4),prim_m(2),cs2_m, & + lamminus,& + psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& + psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(3),prim_p(4),prim_p(2),cs2_p, & + lamplus,& + psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& + psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + usendh,avg_alp,avg_beta) + call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& + fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),& + prim_p(3),prim_p(6), & + avg_det,avg_alp,avg_beta) + call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& + fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& + prim_m(3),prim_m(6), & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call eigenvalues_general(& + prim_m(4),prim_m(2),prim_m(3),cs2_m, & + lamminus,& + psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& + psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(4),prim_p(2),prim_p(3),cs2_p, & + lamplus,& + psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& + psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + usendh,avg_alp,avg_beta) + call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& + fplus(1),fplus(4),fplus(2),fplus(3),fplus(5), & + prim_p(4),prim_p(6), & + avg_det,avg_alp,avg_beta) + call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& + fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & + prim_m(4),prim_m(6), & + avg_det,avg_alp,avg_beta) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + if(tadmor.eq.0) then + +!!$ Find minimum and maximum wavespeeds + + charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charpm = charmax - charmin + +!!$ Calculate flux by standard formula + + do m = 1,5 + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + + end do + + else + ! Tadmor's semi-discrete scheme: JcP 160, 241 (2000) + + charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), & + abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), & + abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5))) + + do m = 1,5 + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* & + qdiff(m) + + end do + + + end if + + + + end if !!! The end of the SpaceMask check for a trivial RP. + + densflux(i, j, k) = f1(1) + sxflux(i, j, k) = f1(2) + syflux(i, j, k) = f1(3) + szflux(i, j, k) = f1(4) + tauflux(i, j, k) = f1(5) + + end do + end do + end do + +end subroutine GRHydro_HLLEGeneral + + +subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) + + USE GRHydro_Eigenproblem + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: i, j, k, m + CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1 + CCTK_REAL, dimension(5) :: lamminus,lamplus + CCTK_REAL, dimension(number_of_tracers) :: qdiff + CCTK_REAL, dimension(6) :: prim_p, prim_m + CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det + CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, psi4h, & + cs2_p, cs2_m, dpdeps_p, dpdeps_m + + integer tadmor + + if(CCTK_EQUALS(HLLE_type,"Tadmor")) then + tadmor = 1 + else + tadmor = 0 + endif + + + 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 + + f1 = 0.d0 + lamminus = 0.d0 + lamplus = 0.d0 + cons_p = 0.d0 + cons_m = 0.d0 + fplus = 0.d0 + fminus = 0.d0 + qdiff = 0.d0 + +!!$ Set the left (p for plus) and right (m_i for minus, i+1) states + + cons_p(:) = cons_tracerplus(i,j,k,:) + cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + + prim_p(1) = rhoplus(i,j,k) + prim_p(2) = velxplus(i,j,k) + prim_p(3) = velyplus(i,j,k) + prim_p(4) = velzplus(i,j,k) + prim_p(5) = epsplus(i,j,k) + + prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) + prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) + + prim_p(6) = pressplus(i,j,k) + cs2_p = eos_cs2_p(i,j,k) + dpdeps_p = eos_dpdeps_p(i,j,k) + + prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) + cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset) + dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset) + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. + + if (shift_state .ne. 0) then + if (flux_direction == 1) then + avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & + betax(i,j,k)) + else if (flux_direction == 2) then + avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & + betay(i,j,k)) + else if (flux_direction == 3) then + avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & + betaz(i,j,k)) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + else + avg_beta = 0.d0 + end if + + avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) + + gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) + gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) + gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) + gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) + gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) + gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) + + if (conformal_state .eq. 0) then + psi4h = 1.d0 + else + psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 + end if + + call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det) + +!!$ If the Riemann problem is trivial, just calculate the fluxes from the +!!$ left state and skip to the next cell + + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, & + psi4h*gyyh, psi4h*gyzh, psi4h*gzzh) + + if (flux_direction == 1) then + usendh = uxxh + else if (flux_direction == 2) then + usendh = uyyh + else if (flux_direction == 3) then + usendh = uzzh + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + call eigenvalues_general(& + prim_m(2),prim_m(3),prim_m(4),cs2_m, & + lamminus,& + psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(2),prim_p(3),prim_p(4),cs2_p, & + lamplus,& + psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& + psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else if (flux_direction == 2) then + call eigenvalues_general(& + prim_m(3),prim_m(4),prim_m(2),cs2_m, & + lamminus,& + psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& + psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(3),prim_p(4),prim_p(2),cs2_p, & + lamplus,& + psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& + psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else if (flux_direction == 3) then + call eigenvalues_general(& + prim_m(4),prim_m(2),prim_m(3),cs2_m, & + lamminus,& + psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& + psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues_general(& + prim_p(4),prim_p(2),prim_p(3),cs2_p, & + lamplus,& + psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& + psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + if(tadmor.eq.0) then + +!!$ Find minimum and maximum wavespeeds + + charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + + charpm = charmax - charmin + +!!$ Calculate flux by standard formula + + do m = 1,number_of_tracers + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + end do + + else + ! Tadmor's semi-descrite scheme: JcP 160, 241 (2000) + + charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), & + abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), & + abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5))) + + do m = 1,number_of_tracers + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* & + qdiff(m) + + end do + + + end if + + + cons_tracerflux(i,j,k,:) = f1(:) + + end do + end do + end do + +end subroutine GRHydro_HLLE_TracerGeneral + -- cgit v1.2.3