From 7a689eb93032bc54ddd3d7935188baef49b0f759 Mon Sep 17 00:00:00 2001
From: bmundim <bmundim@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45>
Date: Mon, 3 Oct 2011 16:05:30 +0000
Subject: RIT dev: Implement local Lax-Friedrichs flux calculation (MHD only).

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@281 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
---
 src/GRHydro_HLLEM.F90         | 34 +++++++++++++++++++++++++---------
 src/GRHydro_PreLoop.F90       |  7 +++++++
 src/GRHydro_RiemannSolveM.F90 |  2 +-
 src/GRHydro_Scalars.F90       |  1 +
 4 files changed, 34 insertions(+), 10 deletions(-)

(limited to 'src')

diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
index 200575f..7e9c541 100644
--- a/src/GRHydro_HLLEM.F90
+++ b/src/GRHydro_HLLEM.F90
@@ -33,6 +33,7 @@
 
 subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
   USE GRHydro_EigenproblemM
+  USE GRHydro_Scalars
 
   implicit none
   
@@ -43,7 +44,7 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
   CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff
   CCTK_REAL, dimension(10) :: prim_p, prim_m
   CCTK_REAL, dimension(5) :: lamminus,lamplus
-  CCTK_REAL ::  charmin, charmax, charpm,avg_alp,avg_det, sdet
+  CCTK_REAL ::  charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet
   CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
        uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, &
        cs2_p, cs2_m, dpdeps_p, dpdeps_m
@@ -404,6 +405,8 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
           charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
                lamplus(4),lamplus(5),  lamminus(1),lamminus(2),lamminus(3),&
                lamminus(4),lamminus(5))
+
+          chartop = max(-charmin,charmax)
           
           charpm = charmax - charmin
 
@@ -413,8 +416,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
              
              qdiff(m) = cons_m(m) - cons_p(m)
              
-             f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
-                  charmax * charmin * qdiff(m)) / charpm
+            if (HLLE) then
+              f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+                       charmax * charmin * qdiff(m)) / charpm
+            else if (LLF) then
+              f1(m) = 0.5d0 * (fplus(m) + fminus(m) - chartop * qdiff(m)) 
+            end if 
              
           end do
 
@@ -422,17 +429,26 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
            psidcdiff = psidcm - psidcp
            select case(whichpsidcspeed)
              case(0)
-             psidcf = (charmax * psidcfp - charmin * psidcfm + &
-                  charmax * charmin * psidcdiff) / charpm
+               if (HLLE) then
+                 psidcf = (charmax * psidcfp - charmin * psidcfm + &
+                           charmax * charmin * psidcdiff) / charpm
+               else if (LLF) then
+                 psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff)
+               end if
              case(1)
              !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic?
-             psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + &
-                  1.d0 * (-1.d0) * psidcdiff) / 2.d0
+             psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + &
+                  1.d0 * (-1.d0) * psidcdiff) 
              case(2)
              charmax = setcharmax
              charmin = setcharmin
-             psidcf = (charmax * psidcfp - charmin * psidcfm + &
-                  charmax * charmin * psidcdiff) / charpm
+               if (HLLE) then
+                 psidcf = (charmax * psidcfp - charmin * psidcfm + &
+                           charmax * charmin * psidcdiff) / charpm
+               else if (LLF) then
+                 chartop = max(-charmin,charmax)
+                 psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff)
+               end if
            end select
           endif
 
diff --git a/src/GRHydro_PreLoop.F90 b/src/GRHydro_PreLoop.F90
index e0b3bc2..bd15d1d 100644
--- a/src/GRHydro_PreLoop.F90
+++ b/src/GRHydro_PreLoop.F90
@@ -44,6 +44,8 @@ subroutine GRHydro_Scalar_Setup(CCTK_ARGUMENTS)
   MC1 = .false.
   MC2 = .false.
   SUPERBEE = .false.
+  HLLE = .false.
+  LLF = .false.
 
   if (CCTK_EQUALS(tvd_limiter,"minmod")) then
     MINMOD = .true.    
@@ -92,6 +94,11 @@ subroutine GRHydro_Scalar_Setup(CCTK_ARGUMENTS)
     call CCTK_WARN(0, "Numerical Viscosity Mode not recognized!")
   end if
 
+  if (CCTK_EQUALS(riemann_solver,"HLLE")) then
+    HLLE = .true.
+  else if (CCTK_EQUALS(riemann_solver,"LLF")) then
+    LLF = .true.
+  end if
 
 end subroutine GRHydro_Scalar_Setup
 
diff --git a/src/GRHydro_RiemannSolveM.F90 b/src/GRHydro_RiemannSolveM.F90
index bd16778..25f2006 100644
--- a/src/GRHydro_RiemannSolveM.F90
+++ b/src/GRHydro_RiemannSolveM.F90
@@ -38,7 +38,7 @@ subroutine RiemannSolveM(CCTK_ARGUMENTS)
 
   CCTK_INT :: i,j,k
   
-  if (CCTK_EQUALS(riemann_solver,"HLLE")) then
+  if (CCTK_EQUALS(riemann_solver,"HLLE").or.CCTK_EQUALS(riemann_solver,"LLF")) then
      
      call GRHydro_HLLEM(CCTK_PASS_FTOF)
      
diff --git a/src/GRHydro_Scalars.F90 b/src/GRHydro_Scalars.F90
index ae03e02..e62671a 100644
--- a/src/GRHydro_Scalars.F90
+++ b/src/GRHydro_Scalars.F90
@@ -19,6 +19,7 @@
    LOGICAL :: MINMOD, MINMOD2, MINMOD3, MC1, MC2, SUPERBEE, PPM3, PPM4
    LOGICAL :: ANALYTICAL
    LOGICAL :: FAST
+   LOGICAL :: HLLE, LLF
 
    
  end module GRHydro_Scalars
-- 
cgit v1.2.3