Grzegorz Kowal
a45a380d0c
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>
777 lines
23 KiB
Fortran
777 lines
23 KiB
Fortran
!!******************************************************************************
|
||
!!
|
||
!! This file is part of the AMUN source code, a program to perform
|
||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||
!! adaptive mesh.
|
||
!!
|
||
!! Copyright (C) 2008-2024 Grzegorz Kowal <grzegorz@amuncode.org>
|
||
!!
|
||
!! This program is free software: you can redistribute it and/or modify
|
||
!! it under the terms of the GNU General Public License as published by
|
||
!! the Free Software Foundation, either version 3 of the License, or
|
||
!! (at your option) any later version.
|
||
!!
|
||
!! This program is distributed in the hope that it will be useful,
|
||
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
!! GNU General Public License for more details.
|
||
!!
|
||
!! You should have received a copy of the GNU General Public License
|
||
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
!!
|
||
!!******************************************************************************
|
||
!!
|
||
!! module: SOURCES
|
||
!!
|
||
!! This modules adds source terms.
|
||
!!
|
||
!!******************************************************************************
|
||
!
|
||
module sources
|
||
|
||
implicit none
|
||
|
||
! interfaces for the extra source terms subroutine
|
||
!
|
||
abstract interface
|
||
subroutine update_extra_sources_iface(t, dt, pdata)
|
||
use blocks, only : block_data
|
||
implicit none
|
||
real(kind=8) , intent(in) :: t, dt
|
||
type(block_data), pointer , intent(inout) :: pdata
|
||
end subroutine
|
||
end interface
|
||
|
||
! pointer to the extra source terms subroutine
|
||
!
|
||
procedure(update_extra_sources_iface), pointer, save :: &
|
||
update_extra_sources => null()
|
||
|
||
! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
|
||
!
|
||
integer , save :: glm_type = 0
|
||
character(len=32), save :: glm_name = "none"
|
||
|
||
! viscosity coefficient
|
||
!
|
||
real(kind=8) , save :: viscosity = 0.0d+00
|
||
|
||
! resistivity coefficient
|
||
!
|
||
real(kind=8) , save :: resistivity = 0.0d+00
|
||
real(kind=8) , save :: anomalous = 0.0d+00
|
||
real(kind=8) , save :: jcrit = 1.0d+00
|
||
|
||
private
|
||
|
||
public :: initialize_sources, finalize_sources, print_sources
|
||
public :: update_sources, update_extra_sources
|
||
public :: viscosity, resistivity
|
||
|
||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||
!
|
||
contains
|
||
!
|
||
!===============================================================================
|
||
!!
|
||
!!*** PUBLIC SUBROUTINES *****************************************************
|
||
!!
|
||
!===============================================================================
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INITIALIZE_SOURCES:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine initializes module SOURCES.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine initialize_sources(verbose, status)
|
||
|
||
use parameters, only : get_parameter
|
||
|
||
implicit none
|
||
|
||
logical, intent(in) :: verbose
|
||
integer, intent(out) :: status
|
||
|
||
character(len=8) :: tglm = "none"
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
! get the type of the GLM source terms
|
||
!
|
||
call get_parameter("glm_source_terms", tglm)
|
||
|
||
! set the glm_type variable to correct value
|
||
!
|
||
select case(trim(tglm))
|
||
case("eglm", "EGLM")
|
||
glm_type = 1
|
||
glm_name = "EGLM"
|
||
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"
|
||
end select
|
||
|
||
! get viscosity coefficient
|
||
!
|
||
call get_parameter("viscosity", viscosity)
|
||
|
||
if (viscosity < 0.0d+00) then
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "Negative viscosity coefficient!"
|
||
end if
|
||
status = 1
|
||
end if
|
||
|
||
! get resistivity coefficients
|
||
!
|
||
call get_parameter("resistivity", resistivity)
|
||
|
||
if (resistivity < 0.0d+00) then
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "Negative resistivity coefficient!"
|
||
end if
|
||
status = 1
|
||
end if
|
||
|
||
call get_parameter("anomalous", anomalous)
|
||
|
||
if (anomalous < 0.0d+00) then
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "Negative anomalous resistivity coefficient!"
|
||
end if
|
||
status = 1
|
||
end if
|
||
|
||
call get_parameter("jcrit", jcrit)
|
||
|
||
if (jcrit <= 0.0d+00) then
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "Non-positive critical current density coefficient!"
|
||
end if
|
||
status = 1
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine initialize_sources
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FINALIZE_SOURCES:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine releases memory used by the module.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine finalize_sources(status)
|
||
|
||
implicit none
|
||
|
||
integer, intent(out) :: status
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine finalize_sources
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRINT_SOURCES:
|
||
! ------------------------
|
||
!
|
||
! Subroutine prints module parameters.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine print_sources(verbose)
|
||
|
||
use equations, only : magnetized
|
||
use helpers , only : print_section, print_parameter
|
||
|
||
implicit none
|
||
|
||
logical, intent(in) :: verbose
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (verbose) then
|
||
|
||
call print_section(verbose, "Source terms")
|
||
call print_parameter(verbose, "viscosity" , viscosity )
|
||
if (magnetized) then
|
||
call print_parameter(verbose, "resistivity" , resistivity)
|
||
if (anomalous > 0.0d+00) then
|
||
call print_parameter(verbose, "anomalous" , anomalous )
|
||
call print_parameter(verbose, "jcrit" , jcrit )
|
||
end if
|
||
call print_parameter(verbose, "glm source terms", glm_name )
|
||
end if
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine print_sources
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine UPDATE_SOURCES:
|
||
! -------------------------
|
||
!
|
||
! Subroutine adds source terms to the array of increments dU.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! t, dt - the time and time increment;
|
||
! pdata - the pointer to a data block;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine update_sources(t, dt, pdata)
|
||
|
||
use blocks , only : block_data
|
||
use coordinates, only : nn => bcells
|
||
use coordinates, only : ax, adx, ay, ady, adz
|
||
#if NDIMS == 3
|
||
use coordinates, only : az
|
||
#endif /* NDIMS == 3 */
|
||
use equations , only : magnetized
|
||
use equations , only : inx, iny, inz
|
||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
|
||
use equations , only : ibx, iby, ibz, ibp
|
||
use gravity , only : gravity_enabled, gravitational_acceleration
|
||
use helpers , only : print_message
|
||
use operators , only : divergence, gradient, laplace, curl
|
||
use workspace , only : resize_workspace, work, work_in_use
|
||
|
||
implicit none
|
||
|
||
real(kind=8) , intent(in) :: t, dt
|
||
type(block_data), pointer , intent(inout) :: pdata
|
||
|
||
logical, save :: first = .true.
|
||
|
||
integer :: status
|
||
integer :: i, j, k
|
||
real(kind=8) :: fc, gc
|
||
real(kind=8) :: gx, gy, gz
|
||
real(kind=8) :: dvxdx, dvxdy, dvxdz, divv
|
||
real(kind=8) :: dvydx, dvydy, dvydz
|
||
real(kind=8) :: dvzdx, dvzdy, dvzdz
|
||
|
||
integer, save :: nt
|
||
!$ integer :: omp_get_thread_num
|
||
|
||
real(kind=8), dimension(3) :: ga, dh
|
||
real(kind=8), dimension(nn) :: x, y
|
||
#if NDIMS == 3
|
||
real(kind=8), dimension(nn) :: z
|
||
#else /* NDIMS == 3 */
|
||
real(kind=8), dimension( 1) :: z
|
||
#endif /* NDIMS == 3 */
|
||
real(kind=8), dimension(:,:,:) , pointer, save :: db
|
||
real(kind=8), dimension(:,:,:,:,:), pointer, save :: tmp
|
||
!$omp threadprivate(first, nt, db, tmp)
|
||
|
||
character(len=*), parameter :: loc = 'SOURCES:update_sources()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
nt = 0
|
||
!$ nt = omp_get_thread_num()
|
||
if (first) then
|
||
i = nn**NDIMS
|
||
j = 10 * i
|
||
|
||
call resize_workspace(j, status)
|
||
if (status /= 0) then
|
||
call print_message(loc, "Could not resize the workspace!")
|
||
go to 100
|
||
end if
|
||
|
||
#if NDIMS == 3
|
||
db(1:nn,1:nn,1:nn) => work( 1:i,nt)
|
||
tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j,nt)
|
||
#else /* NDIMS == 3 */
|
||
db(1:nn,1:nn,1: 1) => work( 1:i,nt)
|
||
tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j,nt)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
|
||
! proceed only if the gravitational term is enabled
|
||
!
|
||
if (gravity_enabled) then
|
||
|
||
! prepare block coordinates
|
||
!
|
||
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
||
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
||
#if NDIMS == 3
|
||
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! iterate over all positions in the YZ plane
|
||
!
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
! get gravitational acceleration components
|
||
!
|
||
call gravitational_acceleration(t, dt, x(i), y(j), z(k), ga(1:3))
|
||
|
||
! calculate the gravitational source terms
|
||
!
|
||
gx = pdata%q(idn,i,j,k) * ga(1)
|
||
gy = pdata%q(idn,i,j,k) * ga(2)
|
||
#if NDIMS == 3
|
||
gz = pdata%q(idn,i,j,k) * ga(3)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! add source terms to momentum equations
|
||
!
|
||
pdata%du(imx,i,j,k) = pdata%du(imx,i,j,k) + gx
|
||
pdata%du(imy,i,j,k) = pdata%du(imy,i,j,k) + gy
|
||
#if NDIMS == 3
|
||
pdata%du(imz,i,j,k) = pdata%du(imz,i,j,k) + gz
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! add source terms to total energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
#if NDIMS == 2
|
||
pdata%du(ien,i,j,k) = pdata%du(ien,i,j,k) &
|
||
+ gx * pdata%q(ivx,i,j,k) &
|
||
+ gy * pdata%q(ivy,i,j,k)
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
pdata%du(ien,i,j,k) = pdata%du(ien,i,j,k) &
|
||
+ gx * pdata%q(ivx,i,j,k) &
|
||
+ gy * pdata%q(ivy,i,j,k) &
|
||
+ gz * pdata%q(ivz,i,j,k)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
end if
|
||
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
end if ! gravity enabled
|
||
|
||
! proceed only if the viscosity coefficient is not zero
|
||
!
|
||
if (viscosity > 0.0d+00) then
|
||
|
||
if (work_in_use(nt)) &
|
||
call print_message(loc, "Workspace is being used right now! " // &
|
||
"Corruptions can occur!")
|
||
|
||
work_in_use(nt) = .true.
|
||
|
||
! prepare coordinate increments
|
||
!
|
||
dh(1) = adx(pdata%meta%level)
|
||
dh(2) = ady(pdata%meta%level)
|
||
dh(3) = adz(pdata%meta%level)
|
||
|
||
! calculate the velocity Jacobian
|
||
!
|
||
call gradient(dh(:), pdata%q(ivx,:,:,:), tmp(inx,inx:inz,:,:,:))
|
||
call gradient(dh(:), pdata%q(ivy,:,:,:), tmp(iny,inx:inz,:,:,:))
|
||
call gradient(dh(:), pdata%q(ivz,:,:,:), tmp(inz,inx:inz,:,:,:))
|
||
|
||
! iterate over all cells
|
||
!
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
! prepare the νρ factor
|
||
!
|
||
gc = viscosity * pdata%q(idn,i,j,k)
|
||
fc = 2.0d+00 * gc
|
||
|
||
! get the velocity Jacobian elements
|
||
!
|
||
dvxdx = tmp(inx,inx,i,j,k)
|
||
dvxdy = tmp(inx,iny,i,j,k)
|
||
dvxdz = tmp(inx,inz,i,j,k)
|
||
dvydx = tmp(iny,inx,i,j,k)
|
||
dvydy = tmp(iny,iny,i,j,k)
|
||
dvydz = tmp(iny,inz,i,j,k)
|
||
dvzdx = tmp(inz,inx,i,j,k)
|
||
dvzdy = tmp(inz,iny,i,j,k)
|
||
dvzdz = tmp(inz,inz,i,j,k)
|
||
divv = (dvxdx + dvydy + dvzdz) / 3.0d+00
|
||
|
||
! calculate elements of the viscous stress tensor
|
||
!
|
||
tmp(inx,inx,i,j,k) = fc * (dvxdx - divv)
|
||
tmp(iny,iny,i,j,k) = fc * (dvydy - divv)
|
||
tmp(inz,inz,i,j,k) = fc * (dvzdz - divv)
|
||
tmp(inx,iny,i,j,k) = gc * (dvxdy + dvydx)
|
||
tmp(inx,inz,i,j,k) = gc * (dvxdz + dvzdx)
|
||
tmp(iny,inz,i,j,k) = gc * (dvydz + dvzdy)
|
||
tmp(iny,inx,i,j,k) = tmp(inx,iny,i,j,k)
|
||
tmp(inz,inx,i,j,k) = tmp(inx,inz,i,j,k)
|
||
tmp(inz,iny,i,j,k) = tmp(iny,inz,i,j,k)
|
||
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! calculate the divergence of the first tensor row
|
||
!
|
||
call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))
|
||
|
||
! add viscous source terms to the X momentum equation
|
||
!
|
||
pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) + db(:,:,:)
|
||
|
||
! calculate the divergence of the second tensor row
|
||
!
|
||
call divergence(dh(:), tmp(iny,inx:inz,:,:,:), db(:,:,:))
|
||
|
||
! add viscous source terms to the Y momentum equation
|
||
!
|
||
pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) + db(:,:,:)
|
||
|
||
! calculate the divergence of the third tensor row
|
||
!
|
||
call divergence(dh(:), tmp(inz,inx:inz,:,:,:), db(:,:,:))
|
||
|
||
! add viscous source terms to the Z momentum equation
|
||
!
|
||
pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) + db(:,:,:)
|
||
|
||
! add viscous source term to total energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! iterate over all cells
|
||
!
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
! calculate scalar product of v and viscous stress tensor τ
|
||
!
|
||
gx = pdata%q(ivx,i,j,k) * tmp(inx,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(inx,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(inx,inz,i,j,k)
|
||
gy = pdata%q(ivx,i,j,k) * tmp(iny,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(iny,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(iny,inz,i,j,k)
|
||
gz = pdata%q(ivx,i,j,k) * tmp(inz,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(inz,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(inz,inz,i,j,k)
|
||
|
||
! update (v.τ), use the first row of the tensor tmp
|
||
!
|
||
tmp(inx,inx,i,j,k) = gx
|
||
tmp(inx,iny,i,j,k) = gy
|
||
tmp(inx,inz,i,j,k) = gz
|
||
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! calculate the divergence of (v.τ)
|
||
!
|
||
call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))
|
||
|
||
! update the energy increment
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + db(:,:,:)
|
||
|
||
end if ! ien > 0
|
||
|
||
work_in_use(nt) = .false.
|
||
|
||
end if ! viscosity is not zero
|
||
|
||
!=== add magnetic field related source terms ===
|
||
!
|
||
if (magnetized) then
|
||
|
||
if (work_in_use(nt)) &
|
||
call print_message(loc, "Workspace is being used right now! " // &
|
||
"Corruptions can occur!")
|
||
|
||
work_in_use(nt) = .true.
|
||
|
||
! prepare coordinate increments
|
||
!
|
||
dh(1) = adx(pdata%meta%level)
|
||
dh(2) = ady(pdata%meta%level)
|
||
dh(3) = adz(pdata%meta%level)
|
||
|
||
! add the EGLM-MHD and KEPES source terms
|
||
!
|
||
if (glm_type > 0) then
|
||
|
||
! calculate the magnetic field divergence
|
||
!
|
||
call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))
|
||
|
||
! update the momentum component increments, i.e.
|
||
! d/dt (ρv) + ∇.F = - (∇.B)B
|
||
!
|
||
pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) &
|
||
- db(:,:,:) * pdata%q(ibx,:,:,:)
|
||
pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) &
|
||
- db(:,:,:) * pdata%q(iby,:,:,:)
|
||
pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) &
|
||
- db(:,:,:) * pdata%q(ibz,:,:,:)
|
||
|
||
! add the HEGLM-MHD and KEPES source terms
|
||
!
|
||
if (glm_type > 1) then
|
||
|
||
! update magnetic field component increments, i.e.
|
||
! d/dt B + ∇.F = - (∇.B)v
|
||
!
|
||
pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) &
|
||
- db(:,:,:) * pdata%q(ivx,:,:,:)
|
||
pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) &
|
||
- db(:,:,:) * pdata%q(ivy,:,:,:)
|
||
pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) &
|
||
- db(:,:,:) * pdata%q(ivz,:,:,:)
|
||
|
||
! update the energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! calculate scalar product of velocity and magnetic field
|
||
!
|
||
tmp(inx,inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) &
|
||
* pdata%q(ibx:ibz,:,:,:), 1)
|
||
|
||
! add the divergence potential source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = - (∇.B) (v.B)
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
|
||
- db(:,:,:) * tmp(inx,inx,:,:,:)
|
||
|
||
end if ! ien > 0
|
||
|
||
end if ! glm_type > 1
|
||
|
||
if (ibp > 0) then
|
||
if (glm_type == 3) then
|
||
|
||
! calculate the divergence of velocity
|
||
!
|
||
call divergence(dh(:), pdata%q(ivx:ivz,:,:,:), db(:,:,:))
|
||
|
||
! add the divergence potential source term to the energy equation, i.e.
|
||
! d/dt ψ + ∇.F = ½ψ(∇.v)
|
||
!
|
||
pdata%du(ibp,:,:,:) = pdata%du(ibp,:,:,:) &
|
||
+ 5.0d-01 * pdata%q(ibp,:,:,:) * db(:,:,:)
|
||
|
||
else if (ien > 0) then
|
||
|
||
! calculate the gradient of the 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.∇)ψ
|
||
!
|
||
pdata%du(ibp,:,:,:) = pdata%du(ibp,:,:,:) - db(:,:,:)
|
||
|
||
! add the divergence potential source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = - B.(∇ψ)
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
|
||
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
|
||
|
||
end if ! glm == 3
|
||
end if ! ibp
|
||
end if ! glm_type > 0
|
||
|
||
! if anomalous resistivity is enabled
|
||
!
|
||
if (anomalous > 0.0d+00) then
|
||
|
||
! calculate current density (J = ∇xB)
|
||
!
|
||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inx:inz,inx,:,:,:))
|
||
|
||
! calculate the normalized absolute value of current density (|J|/Jcrit)
|
||
!
|
||
tmp(inx,iny,:,:,:) = sqrt(sum(tmp(inx:inz,inx,:,:,:)**2, 1)) / jcrit
|
||
|
||
! calculate the local resistivity [ηu + ηa (|J|/Jcrit - 1) H(|J|/Jcrit)]
|
||
!
|
||
tmp(iny,iny,:,:,:) = resistivity + &
|
||
anomalous * max(0.0d+00, (tmp(inx,iny,:,:,:) - 1.0d+00))
|
||
|
||
! multiply the current density vector by the local resistivity (ηJ)
|
||
!
|
||
tmp(inx,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(inx,inx,:,:,:)
|
||
tmp(iny,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(iny,inx,:,:,:)
|
||
tmp(inz,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(inz,inx,:,:,:)
|
||
|
||
! calculate the curl of (ηJ)
|
||
!
|
||
call curl(dh(:), tmp(inx:inz,inz,:,:,:), tmp(inx:inz,iny,:,:,:))
|
||
|
||
! update magnetic field component increments
|
||
!
|
||
pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) - tmp(inx,iny,:,:,:)
|
||
pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) - tmp(iny,iny,:,:,:)
|
||
pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) - tmp(inz,iny,:,:,:)
|
||
|
||
! update energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! calculate the vector product Bx(η ∇xB)
|
||
!
|
||
tmp(inx,iny,:,:,:) = pdata%q(iby,:,:,:) * tmp(inz,inz,:,:,:) &
|
||
- pdata%q(ibz,:,:,:) * tmp(iny,inz,:,:,:)
|
||
tmp(iny,iny,:,:,:) = pdata%q(ibz,:,:,:) * tmp(inx,inz,:,:,:) &
|
||
- pdata%q(ibx,:,:,:) * tmp(inz,inz,:,:,:)
|
||
tmp(inz,iny,:,:,:) = pdata%q(ibx,:,:,:) * tmp(iny,inz,:,:,:) &
|
||
- pdata%q(iby,:,:,:) * tmp(inx,inz,:,:,:)
|
||
|
||
! calculate the divergence ∇.[Bx(η ∇xB)]
|
||
!
|
||
call divergence(dh(:), tmp(inx:inz,iny,:,:,:), db(:,:,:))
|
||
|
||
! add the second resistive source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = η J²
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + db(:,:,:)
|
||
|
||
end if ! energy equation present
|
||
|
||
else if (resistivity > 0.0d+00) then
|
||
|
||
! calculate the Laplace operator of B, i.e. Δ(B)
|
||
!
|
||
call laplace(dh(:), pdata%q(ibx,:,:,:), tmp(inx,inx,:,:,:))
|
||
call laplace(dh(:), pdata%q(iby,:,:,:), tmp(inx,iny,:,:,:))
|
||
call laplace(dh(:), pdata%q(ibz,:,:,:), tmp(inx,inz,:,:,:))
|
||
|
||
! multiply by the resistivity coefficient
|
||
!
|
||
tmp(iny,inx:inz,:,:,:) = resistivity * tmp(inx,inx:inz,:,:,:)
|
||
|
||
! update magnetic field component increments
|
||
!
|
||
pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) + tmp(iny,inx,:,:,:)
|
||
pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) + tmp(iny,iny,:,:,:)
|
||
pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) + tmp(iny,inz,:,:,:)
|
||
|
||
! update energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! add the first resistive source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = η B.[Δ(B)]
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
|
||
+ (pdata%q(ibx,:,:,:) * tmp(iny,inx,:,:,:) &
|
||
+ pdata%q(iby,:,:,:) * tmp(iny,iny,:,:,:) &
|
||
+ pdata%q(ibz,:,:,:) * tmp(iny,inz,:,:,:))
|
||
|
||
! calculate current density J = ∇xB
|
||
!
|
||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inz,inx:inz,:,:,:))
|
||
|
||
! calculate J²
|
||
!
|
||
db(:,:,:) = tmp(inz,inx,:,:,:)**2 &
|
||
+ tmp(inz,iny,:,:,:)**2 &
|
||
+ tmp(inz,inz,:,:,:)**2
|
||
|
||
! add the second resistive source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = η J²
|
||
!
|
||
pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + resistivity * db(:,:,:)
|
||
|
||
end if ! energy equation present
|
||
|
||
end if ! resistivity is not zero
|
||
|
||
work_in_use(nt) = .false.
|
||
|
||
end if ! magnetized
|
||
|
||
! add extra source terms
|
||
!
|
||
if (associated(update_extra_sources)) &
|
||
call update_extra_sources(t, dt, pdata)
|
||
|
||
100 continue
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine update_sources
|
||
|
||
!===============================================================================
|
||
!
|
||
end module sources
|