[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 281) - Commits

3 Oct 2011

User: bmundim
Date: 2011/10/03 11:05 AM
Modified:
 /trunk/
  param.ccl
 /trunk/src/
  GRHydro_HLLEM.F90, GRHydro_PreLoop.F90, GRHydro_RiemannSolveM.F90, GRHydro_Scalars.F90
Log:
 RIT dev: Implement local Lax-Friedrichs flux calculation (MHD only).
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_HLLEM.F90
Delta lines: +25 -9
===================================================================

--- trunk/src/GRHydro_HLLEM.F90	2011-09-26 15:50:04 UTC (rev 280)
+++ trunk/src/GRHydro_HLLEM.F90	2011-10-03 16:05:30 UTC (rev 281)
@@ -33,6 +33,7 @@
subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
   USE GRHydro_EigenproblemM
+  USE GRHydro_Scalars
implicit none
@@ -43,7 +44,7 @@
   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 @@
           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 @@
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 @@
            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
File [modified]: GRHydro_PreLoop.F90
Delta lines: +7 -0
===================================================================
--- trunk/src/GRHydro_PreLoop.F90	2011-09-26 15:50:04 UTC (rev 280)
+++ trunk/src/GRHydro_PreLoop.F90	2011-10-03 16:05:30 UTC (rev 281)
@@ -44,6 +44,8 @@
   MC1 = .false.
   MC2 = .false.
   SUPERBEE = .false.
+  HLLE = .false.
+  LLF = .false.
if (CCTK_EQUALS(tvd_limiter,"minmod")) then
     MINMOD = .true.    
@@ -92,6 +94,11 @@
     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
File [modified]: GRHydro_RiemannSolveM.F90
Delta lines: +1 -1
===================================================================
--- trunk/src/GRHydro_RiemannSolveM.F90	2011-09-26 15:50:04 UTC (rev 280)
+++ trunk/src/GRHydro_RiemannSolveM.F90	2011-10-03 16:05:30 UTC (rev 281)
@@ -38,7 +38,7 @@
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)
File [modified]: GRHydro_Scalars.F90
Delta lines: +1 -0
===================================================================
--- trunk/src/GRHydro_Scalars.F90	2011-09-26 15:50:04 UTC (rev 280)
+++ trunk/src/GRHydro_Scalars.F90	2011-10-03 16:05:30 UTC (rev 281)
@@ -19,6 +19,7 @@
    LOGICAL :: MINMOD, MINMOD2, MINMOD3, MC1, MC2, SUPERBEE, PPM3, PPM4
    LOGICAL :: ANALYTICAL
    LOGICAL :: FAST
+   LOGICAL :: HLLE, LLF
end module GRHydro_Scalars
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +1 -0
===================================================================
--- trunk/param.ccl	2011-09-26 15:50:04 UTC (rev 280)
+++ trunk/param.ccl	2011-10-03 16:05:30 UTC (rev 281)
@@ -203,6 +203,7 @@
   "Roe"		:: "Standard Roe solver"
   "Marquina"	:: "Marquina flux"
   "HLLE"	:: "HLLE"
+  "LLF"	:: "Local Lax-Friedrichs (MHD only at the moment)"
 } "HLLE"

    