From fb25489487a676605510f5ba338f30802d4a7dbb Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Wed, 6 Mar 2024 21:21:59 -0300
Subject: [PATCH] SCHEMES: Fix ROE solvers for MHD.

The tearing test problem was numerically unstable in the adiabatic case.
Setting bty to 1 in the case of br beeing zero solved the issue. Make it
consistent for the isothermal ROE solver as well. Additionaly, one typo
was fixed in the calculation of the left eigenvectors.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 sources/schemes.F90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/sources/schemes.F90 b/sources/schemes.F90
index ea13213..958ac5d 100644
--- a/sources/schemes.F90
+++ b/sources/schemes.F90
@@ -1969,7 +1969,7 @@ module schemes
         bty  = qi(iby) / br
         btz  = qi(ibz) / br
       else
-        bty  = 0.0d+00
+        bty  = 1.0d+00
         btz  = 0.0d+00
       end if
 
@@ -3170,7 +3170,7 @@ module schemes
         btz = qi(ibz) / br
       else
         br  = 0.0d+00
-        bty = 0.0d+00
+        bty = 1.0d+00
         btz = 0.0d+00
       end if
 
@@ -3306,18 +3306,18 @@ module schemes
       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,5) =   f1
       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,5) =   f1
       lvec(7,6) =   lvec(1,6)
       lvec(7,7) =   lvec(1,7)
 
@@ -3327,18 +3327,18 @@ module schemes
       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,5) =   f1
       lvec(3,6) = - f3 * bty - f1 * qi(iby)
-      lvec(3,6) = - f3 * btz - f1 * qi(ibz)
+      lvec(3,7) = - 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,5) =   f1
       lvec(5,6) =   lvec(3,6)
       lvec(5,7) =   lvec(3,7)
 
@@ -3359,10 +3359,10 @@ module schemes
 
       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,5) = - f1
       lvec(4,6) =   f1 * qi(iby)
       lvec(4,7) =   f1 * qi(ibz)