fix a bug in artificial viscosity (quokka-astro#1605)

### Description
In the current implementation of artificial viscosity, the viscosity
term is added in the wrong direction. This PR fixes the bug by using the
variable `F_canonical` to store intermediate fluxes and then applying
the correct permutation.

### Related issues
Are there any GitHub issues that are fixed by this pull request? Add a
link to them here.

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [ ] I have added a link to any related issues (if applicable; see
above).
- [ ] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [ ] I have added tests for any new physics that this PR adds to the
code.
- [x] *(For quokka-astro org members)* I have manually triggered the GPU
tests with the magic comment `/azp run`.

Author: mxchen39

src/hydro/hydro_system.hpp

@@ -1410,13 +1410,6 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 			x1FSpds(i, j, k, 1) = fspd_p;
 		}
 
-		quokka::valarray<double, nHydroScalars_> F = F_canonical;
-
-		// permute momentum components according to flux direction DIR
-		F[velN_index] = F_canonical[x1Momentum_index];
-		F[velV_index] = F_canonical[x2Momentum_index];
-		F[velW_index] = F_canonical[x3Momentum_index];
-
 		// add artificial viscosity
 		// following Colella & Woodward (1984), eq. (4.2)
 		const double div_v = AMREX_D_TERM(du, +0.5 * (dvl + dvr), +0.5 * (dwl + dwr));
@@ -1441,7 +1434,14 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 			}
 		}
 
-		F = F + viscosity * (U_L - U_R);
+		F_canonical = F_canonical + viscosity * (U_L - U_R);
+
+		quokka::valarray<double, nHydroScalars_> F = F_canonical;
+
+		// permute momentum components according to flux direction DIR
+		F[velN_index] = F_canonical[x1Momentum_index];
+		F[velV_index] = F_canonical[x2Momentum_index];
+		F[velW_index] = F_canonical[x3Momentum_index];
 
 		// set energy fluxes to zero if EOS is isothermal
 		if constexpr (HydroSystem<problem_t>::is_eos_isothermal()) {