!!****************************************************************************** !! !! 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 !! !! 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 . !! !!****************************************************************************** !! !! 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