SCHEMES: Rewrite riemann_mhd_adi_roe().

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>

sources/schemes.F90

@@ -3250,129 +3250,294 @@ module schemes
   subroutine riemann_mhd_adi_roe(ql, qr, fi)
 
     use coordinates, only : nn => bcells
-    use equations  , only : nf, idn, ivx, ivz, imx, imy, imz, ipr, ien, &
+    use equations  , only : nf, adiabatic_index
+    use equations  , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien, &
                             ibx, iby, ibz, ibp
-    use equations  , only : prim2cons, fluxspeed, eigensystem_roe
+    use equations  , only : prim2cons, fluxspeed
 
     implicit none
 
     real(kind=8), dimension(:,:), intent(in)  :: ql, qr
     real(kind=8), dimension(:,:), intent(out) :: fi
 
-    integer      :: i, p
+    integer      :: i
     real(kind=8) :: sdl, sdr, sds
-    real(kind=8) :: xx, yy
+    real(kind=8) :: xx, yy, pml, pmr
-    real(kind=8) :: pml, pmr
+
+    real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs
+    real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty
+    real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm
+    real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb
+    real(kind=8) :: f1, f2, f3, f4, f5
 
     real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
-    real(kind=8), dimension( 2,nn) :: cl, cr
+    real(kind=8), dimension(nf   ) :: qi
-    real(kind=8), dimension(nf   ) :: qi, ci, al
+    real(kind=8), dimension(7    ) :: lm, al, df
-    real(kind=8), dimension(nf,nf) :: li, ri
+
+    logical                     , save :: first = .true.
+    integer     , dimension(7)  , save :: ivs
+    real(kind=8)                , save :: adi_m1, adi_m2, adi_m2d1
+    real(kind=8), dimension(7,7), save :: rvec, lvec
+!$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec)
 
 !-------------------------------------------------------------------------------
 !
+    if (first) then
+      adi_m1   = adiabatic_index - 1.0d+00
+      adi_m2   = adiabatic_index - 2.0d+00
+      adi_m2d1 = adi_m2 / adi_m1
+
+      rvec(:,:) = 0.0d+00
+      lvec(:,:) = 0.0d+00
+      rvec(4,1) = 1.0d+00
+
+      ivs = [ idn, imx, imy, imz, ipr, iby, ibz ]
+
+      first = .false.
+    end if
+
     call prim2cons(ql, ul)
     call prim2cons(qr, ur)
 
-    call fluxspeed(ql, ul, fl, cl)
+    call fluxspeed(ql, ul, fl)
-    call fluxspeed(qr, ur, fr, cr)
+    call fluxspeed(qr, ur, fr)
 
     do i = 1, nn
 
-! calculate variables for the Roe intermediate state
+      pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
-!
+      pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
       sdl = sqrt(ql(idn,i))
       sdr = sqrt(qr(idn,i))
       sds = sdl + sdr
 
-! prepare magnetic pressures
-!
-      pml = 0.5d+00 * sum(ql(ibx:ibz,i)**2)
-      pmr = 0.5d+00 * sum(qr(ibx:ibz,i)**2)
-
-! prepare the Roe intermediate state vector (eq. 11.60 in [2])
-!
       qi(idn)     = sdl * sdr
-      qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
+      qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
-      qi(ipr)     = ((ul(ien,i) + ql(ipr,i) + pml) / sdl                       &
+      qi(ipr)     = ((ul(ien,i) + ql(ipr,i) + pml) / sdl &
                   +  (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds
+      qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
       qi(ibx)     = ql(ibx,i)
-      qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
       qi(ibp)     = ql(ibp,i)
 
-! prepare coefficients
+      f1 = qr(iby,i) - ql(iby,i)
-!
+      f2 = qr(ibz,i) - ql(ibz,i)
-      xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2                               &
+      xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
-                                        + (ql(ibz,i) - qr(ibz,i))**2) / sds**2
+      yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
-      yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)
-
-! obtain eigenvalues and eigenvectors
-!
-      call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))
-
-! return upwind fluxes
-!
-      if (ci(1) >= 0.0d+00) then
-
-        fi(:,i) = fl(:,i)
-
-      else if (ci(nf) <= 0.0d+00) then
-
-        fi(:,i) = fr(:,i)
 
+      br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
+      if (br2 > 0.0d+00) then
+        br  = sqrt(br2)
+        bty = qi(iby) / br
+        btz = qi(ibz) / br
       else
+        br  = 0.0d+00
+        bty = 0.0d+00
+        btz = 0.0d+00
+      end if
 
-! prepare fluxes which do not change across the states
+      v2   = sum(qi(ivx:ivz) * qi(ivx:ivz))
-!
+      v2h  = 5.0d-01 * v2
-        fi(ibx,i) = fl(ibx,i)
+      ayy  = adi_m1 - adi_m2 * yy
-        fi(ibp,i) = fl(ibp,i)
+      sqty = sqrt(ayy)
+      ca2  = qi(ibx) * qi(ibx)
+      hp   = qi(ien) - (ca2 + br2) / qi(idn)
+      ca2  = ca2 / qi(idn)
+      f1   = adi_m1 * (hp - v2h)
+      f2   = adi_m2 * xx
+      if (f1 > f2) then
+        cc2 = f1 - f2
+        cc  = sqrt(cc2)
+        f1  = ayy * br2 / qi(idn)
+        f2  = ca2 + f1
+        f3  = f2 + cc2
+        f4  = f2 - cc2
+        f5  = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1)
+        cf2 = 5.0d-01 * (f3 + f5)
+        cs2 = cc2 * ca2 / cf2
+      else
+        cf2 = ca2 + ayy * br2 / qi(idn)
+        cc2 = 0.0d+00
+        cc  = 0.0d+00
+        cs2 = 0.0d+00
+      end if
 
-! calculate wave amplitudes α = L.ΔU
+      if (cf2 > cs2) then
-!
+        f2 = cf2 - cs2
-        al(:) = 0.0d+00
+        if (cc2 > cs2) then
-        do p = 1, nf
+          alf = sqrt((cc2 - cs2) / f2)
-          al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
-        end do
-
-! calculate the flux average
-!
-        fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
-        fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
-        fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
-        fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
-        fi(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i))
-        fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
-        fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))
-
-! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
-!
-        if (qi(ivx) >= 0.0d+00) then
-          do p = 1, nf
-            xx       = - 0.5d+00 * al(p) * abs(ci(p))
-
-            fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
-            fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
-            fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
-            fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
-            fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
-            fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
-            fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
-          end do
         else
-          do p = nf, 1, -1
+          alf = 0.0d+00
-            xx       = - 0.5d+00 * al(p) * abs(ci(p))
-
-            fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
-            fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
-            fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
-            fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
-            fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
-            fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
-            fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
-          end do
         end if
+        if (cf2 > cc2) then
+          als = sqrt((cf2 - cc2) / f2)
+        else
+          als = 0.0d+00
+        end if
+      else
+        alf = 1.0d+00
+        als = 1.0d+00
+      end if
 
-      end if ! sl < 0 < sr
+      cf = sign(sqrt(cf2), qi(ivx))
+      ca = sign(sqrt(ca2), qi(ivx))
+      cs = sign(sqrt(cs2), qi(ivx))
+
+      lm(1) = qi(ivx) + cf
+      lm(2) = qi(ivx) + ca
+      lm(3) = qi(ivx) + cs
+      lm(4) = qi(ivx)
+      lm(5) = qi(ivx) - cs
+      lm(6) = qi(ivx) - ca
+      lm(7) = qi(ivx) - cf
+
+      vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty
+
+      sqrtd    = sqrt(qi(idn))
+      qf       = cf * sign(alf, qi(ibx))
+      qs       = cs * sign(als, qi(ibx))
+      f1       = cc / sqrtd
+      af_prime = f1 * alf
+      as_prime = f1 * als
+      f1       = br / sqty
+      afpbb    = af_prime * f1
+      aspbb    = as_prime * f1
+
+      f1 = qs / sqty
+      f2 = qs * vbet
+      f3 = as_prime / sqty
+      rvec(1,1) = alf
+      rvec(1,2) = alf * lm(1)
+      rvec(1,3) = alf * qi(ivy) - f1 * bty
+      rvec(1,4) = alf * qi(ivz) - f1 * btz
+      rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb
+      rvec(1,6) = f3 * bty
+      rvec(1,7) = f3 * btz
+
+      rvec(7,1) = alf
+      rvec(7,2) = alf * lm(7)
+      rvec(7,3) = alf * qi(ivy) + f1 * bty
+      rvec(7,4) = alf * qi(ivz) + f1 * btz
+      rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb
+      rvec(7,6) = rvec(1,6)
+      rvec(7,7) = rvec(1,7)
+
+      f1 = qf / sqty
+      f2 = qf * vbet
+      f3 = - af_prime / sqty
+      rvec(3,1) = als
+      rvec(3,2) = als * lm(3)
+      rvec(3,3) = als * qi(ivy) + f1 * bty
+      rvec(3,4) = als * qi(ivz) + f1 * btz
+      rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb
+      rvec(3,6) = f3 * bty
+      rvec(3,7) = f3 * btz
+
+      rvec(5,1) = als
+      rvec(5,2) = als * lm(5)
+      rvec(5,3) = als * qi(ivy) - f1 * bty
+      rvec(5,4) = als * qi(ivz) - f1 * btz
+      rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb
+      rvec(5,6) = rvec(3,6)
+      rvec(5,7) = rvec(3,7)
+
+      f1 = sign(1.0d+00, qi(ivx))
+      f2 = sign(1.0d+00 / sqrtd, qi(ibx))
+      rvec(2,3) =   f1 * btz
+      rvec(2,4) = - f1 * bty
+      rvec(2,5) =   f1 * (qi(ivy) * btz - qi(ivz) * bty)
+      rvec(2,6) = - f2 * btz
+      rvec(2,7) =   f2 * bty
+
+      rvec(6,3) = - rvec(2,3)
+      rvec(6,4) = - rvec(2,4)
+      rvec(6,5) = - rvec(2,5)
+      rvec(6,6) =   rvec(2,6)
+      rvec(6,7) =   rvec(2,7)
+
+      rvec(4,2) = qi(ivx)
+      rvec(4,3) = qi(ivy)
+      rvec(4,4) = qi(ivz)
+      rvec(4,5) = v2h + adi_m2d1 * xx
+
+! === update the varying elements of the left eigenvectors matrix
+!
+      norm  = 2.0d+00 * cc2
+      f1    = sqty * br / norm
+      afpb  = af_prime * f1
+      aspb  = as_prime * f1
+      sqty  = sqty / norm
+      vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty
+
+      f1 = adi_m1 * alf / norm
+      f2 = qs * sqty
+      f3 = as_prime * qi(idn) * sqty
+      f4 = alf * cf / norm
+      f5 = qs * vqstr
+      lvec(1,1) =   f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb
+      lvec(1,5) =   f1
+      lvec(1,2) = - f1 * qi(ivx) + f4
+      lvec(1,3) = - f1 * qi(ivy) - f2 * bty
+      lvec(1,4) = - f1 * qi(ivz) - f2 * btz
+      lvec(1,6) =   f3 * bty - f1 * qi(iby)
+      lvec(1,7) =   f3 * btz - f1 * qi(ibz)
+
+      lvec(7,1) =   f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb
+      lvec(7,5) =   f1
+      lvec(7,2) = - f1 * qi(ivx) - f4
+      lvec(7,3) = - f1 * qi(ivy) + f2 * bty
+      lvec(7,4) = - f1 * qi(ivz) + f2 * btz
+      lvec(7,6) =   lvec(1,6)
+      lvec(7,7) =   lvec(1,7)
+
+      f1 = adi_m1 * als / norm
+      f2 = qf * sqty
+      f3 = af_prime * qi(idn) * sqty
+      f4 = als * cs / norm
+      f5 = qf * vqstr
+      lvec(3,1) =   f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb
+      lvec(3,5) =   f1
+      lvec(3,2) = - f1 * qi(ivx) + f4
+      lvec(3,3) = - f1 * qi(ivy) + f2 * bty
+      lvec(3,4) = - f1 * qi(ivz) + f2 * btz
+      lvec(3,6) = - f3 * bty - f1 * qi(iby)
+      lvec(3,6) = - f3 * btz - f1 * qi(ibz)
+
+      lvec(5,1) =   f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb
+      lvec(5,5) =   f1
+      lvec(5,2) = - f1 * qi(ivx) - f4
+      lvec(5,3) = - f1 * qi(ivy) - f2 * bty
+      lvec(5,4) = - f1 * qi(ivz) - f2 * btz
+      lvec(5,6) =   lvec(3,6)
+      lvec(5,7) =   lvec(3,7)
+
+      f1 = sign(5.0d-01, qi(ivx)) * bty
+      f2 = sign(5.0d-01, qi(ivx)) * btz
+      f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx))
+      lvec(2,1) =   qi(ivz) * f1 - qi(ivy) * f2
+      lvec(2,3) =   f2
+      lvec(2,4) = - f1
+      lvec(2,6) = - f2 * f3
+      lvec(2,7) =   f1 * f3
+
+      lvec(6,1) = - lvec(2,1)
+      lvec(6,3) = - lvec(2,3)
+      lvec(6,4) = - lvec(2,4)
+      lvec(6,6) =   lvec(2,6)
+      lvec(6,7) =   lvec(2,7)
+
+      f1 = 1.0d+00 / cc2
+      lvec(4,1) =   1.0d+00 - f1 * (v2h - adi_m2d1 * xx)
+      lvec(4,5) = - f1
+      lvec(4,2) =   f1 * qi(ivx)
+      lvec(4,3) =   f1 * qi(ivy)
+      lvec(4,4) =   f1 * qi(ivz)
+      lvec(4,6) =   f1 * qi(iby)
+      lvec(4,7) =   f1 * qi(ibz)
+
+      al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
+      df = matmul(al, rvec)
+      fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
+      fi(ibx,i) = fl(ibx,i)
+      fi(ibp,i) = fl(ibp,i)
 
     end do