SOURCES: Add source terms related to the MHD KEPES solver.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-01-06 16:46:36 -03:00
parent 25ea3c3daf
commit 226998e80e

View File

@ -121,6 +121,9 @@ module sources
case("heglm", "HEGLM")
glm_type = 2
glm_name = "HEGLM"
case("kepes", "KEPES")
glm_type = 3
glm_name = "KEPES"
case default
glm_type = 0
glm_name = "none"
@ -556,7 +559,7 @@ module sources
dh(2) = ady(pdata%meta%level)
dh(3) = adz(pdata%meta%level)
! add the EGLM-MHD source terms
! add the EGLM-MHD and KEPES source terms
!
if (glm_type > 0) then
@ -571,23 +574,7 @@ module sources
du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
! update the energy equation
!
if (ien > 0 .and. ibp > 0) then
! calculate the gradient of divergence potential
!
call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
! add the divergence potential source term to the energy equation, i.e.
! d/dt E + .F = - B.(ψ)
!
du(ien,:,:,:) = du(ien,:,:,:) &
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
end if ! ien > 0
! add the HEGLM-MHD source terms
! add the HEGLM-MHD and KEPES source terms
!
if (glm_type > 1) then
@ -616,7 +603,41 @@ module sources
end if ! glm_type > 1
end if ! glmtype > 0
if (ibp > 0 .and. (ien > 0 .or. glm_type == 3)) then
! calculate the gradient of divergence correcting field
!
call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
db(:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
! update the divergence correcting field
! d/dt ψ+ .F = - (v.)ψ
!
du(ibp,:,:,:) = du(ibp,:,:,:) - db(:,:,:)
! update the energy equation
!
if (ien > 0) then
if (glm_type == 3) then
! add the divergence potential source term to the energy equation, i.e.
! d/dt E + .F = - ψ(v.)ψ
!
du(ien,:,:,:) = du(ien,:,:,:) - pdata%q(ibp,:,:,:) * db(:,:,:)
else
! add the divergence potential source term to the energy equation, i.e.
! d/dt E + .F = - B.(ψ)
!
du(ien,:,:,:) = du(ien,:,:,:) &
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
end if
end if ! ien > 0
end if ! glm_type == 3
end if ! glm_type > 0
! if anomalous resistivity is enabled
!