From a45a380d0c1c2bdf45d8fd92f44c00b586a032db Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 30 Jul 2024 16:46:29 -0300 Subject: [PATCH] SOURCES: Consider -B.div(B) terms for isothermal KEPES too. The -B.div(B) source term in the momentum equation was dropped for the isothermal MHD case. Apparently, its lack causes some numerical instabilities related to the accumulation of the divergence of B. Therefore, take it into account for the isothermal case too. Signed-off-by: Grzegorz Kowal --- sources/sources.F90 | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/sources/sources.F90 b/sources/sources.F90 index 201a456..be5da5e 100644 --- a/sources/sources.F90 +++ b/sources/sources.F90 @@ -571,21 +571,16 @@ module sources ! call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:)) -! update the momentum component increments (only in the adiabatic case for -! KEPES), i.e. +! update the momentum component increments, i.e. ! d/dt (ρv) + ∇.F = - (∇.B)B ! - if (glm_type < 3 .or. ien > 0) then - - pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) & + pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) & - db(:,:,:) * pdata%q(ibx,:,:,:) - pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) & + pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) & - db(:,:,:) * pdata%q(iby,:,:,:) - pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) & + pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) & - db(:,:,:) * pdata%q(ibz,:,:,:) - end if ! ien > 0 - ! add the HEGLM-MHD and KEPES source terms ! if (glm_type > 1) then