gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

SOURCES: Consider -B.div(B) terms for isothermal KEPES too.

Browse Source 
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 <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2024-07-30 16:46:29 -03:00
parent e339cf15e8
commit a45a380d0c
1 changed files with 4 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
sources/sources.F90
View File

@ -571,21 +571,16 @@ module sources
! !
call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:)) call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))
! update the momentum component increments (only in the adiabatic case for ! update the momentum component increments, i.e.
! KEPES), i.e.
! d/dt (ρv) + ∇.F = - (∇.B)B ! 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,:,:,:) - db(:,:,:) * pdata%q(ibx,:,:,:)
pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) & pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) &
- db(:,:,:) * pdata%q(iby,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) & pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) &
- db(:,:,:) * pdata%q(ibz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
end if ! ien > 0
! add the HEGLM-MHD and KEPES source terms ! add the HEGLM-MHD and KEPES source terms
! !
if (glm_type > 1) then if (glm_type > 1) then