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

USER_PROBLEMS: Add reconnection rate terms through YZ boundaries.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-04-04 08:52:10 -03:00
parent 58075613e6
commit 3def71f6ed
1 changed files with 97 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
sources/user_problem.F90
View File

@ -1531,7 +1531,8 @@ module user_problem
subroutine user_statistics(time) subroutine user_statistics(time)
use blocks , only : block_meta, block_data, list_data use blocks , only : block_meta, block_data, list_data
use coordinates, only : nn => bcells, nb, ne use coordinates, only : nn => bcells, nb, ne, nbm, nep, periodic
use coordinates, only : xmin, xmax
use coordinates, only : adx, ady, adz, advol, ay, yarea use coordinates, only : adx, ady, adz, advol, ay, yarea
use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
use helpers , only : print_message, flush_and_sync use helpers , only : print_message, flush_and_sync
@ -1554,7 +1555,7 @@ module user_problem
!$ integer :: omp_get_thread_num !$ integer :: omp_get_thread_num
integer :: ni, nl, nu, nm, np, status integer :: ni, nl, nu, nm, np, status
real(kind=8) :: dvol, dxz real(kind=8) :: dvol, dyz, dxz
real(kind=8), dimension(3) :: dh real(kind=8), dimension(3) :: dh
real(kind=8), dimension(nn) :: yc real(kind=8), dimension(nn) :: yc
real(kind=8), dimension(32) :: rterms real(kind=8), dimension(32) :: rterms
@ -1604,6 +1605,7 @@ module user_problem
dvol = advol(pmeta%level) dvol = advol(pmeta%level)
dyz = ady(pmeta%level) * adz(pmeta%level)
dxz = adx(pmeta%level) * adz(pmeta%level) dxz = adx(pmeta%level) * adz(pmeta%level)
yc(:) = pmeta%ymin + ay(pmeta%level,:) yc(:) = pmeta%ymin + ay(pmeta%level,:)
@ -1641,6 +1643,98 @@ module user_problem
+ sum(pdata%u(ibx:ibz,nb:ne,nl:nu, : )**2) * dvol + sum(pdata%u(ibx:ibz,nb:ne,nl:nu, : )**2) * dvol
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
if (.not. periodic(1)) then
if (pmeta%xmin <= xmin) then
! advection of magnetic energy across the lower X-boundary
!
rterms(14) = rterms(14) &
#if NDIMS == 3
+ sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne)**2,1) &
* pdata%q(ivx ,nbm:nb,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
+ sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu, : )**2,1) &
* pdata%q(ivx ,nbm:nb,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
rterms(17) = rterms(17) &
#if NDIMS == 3
- sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu,nb:ne) &
* pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne),1) &
* pdata%q(ibx ,nbm:nb,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
- sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu, : ) &
* pdata%q(ibx:ibz,nbm:nb,nl:nu, : ),1) &
* pdata%q(ibx ,nbm:nb,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
if (resistive) then
! diffusion of magnetic energy across the lower X-boundary
!
rterms(20) = rterms(20) &
#if NDIMS == 3
- sum(cur(2,nbm:nb,nl:nu,nb:ne) &
* pdata%q(ibz,nbm:nb,nl:nu,nb:ne) &
- cur(3,nbm:nb,nl:nu,nb:ne) &
* pdata%q(iby,nbm:nb,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
- sum(cur(2,nbm:nb,nl:nu, : ) &
* pdata%q(ibz,nbm:nb,nl:nu, : ) &
- cur(3,nbm:nb,nl:nu, : ) &
* pdata%q(iby,nbm:nb,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
end if ! resistivity
end if
if (pmeta%xmax >= xmax) then
! advection of magnetic energy across the upper X-boundary
!
rterms(14) = rterms(14) &
#if NDIMS == 3
- sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne)**2,1) &
* pdata%q(ivx ,ne:nep,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
- sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu, : )**2,1) &
* pdata%q(ivx ,ne:nep,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
rterms(17) = rterms(17) &
#if NDIMS == 3
+ sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu,nb:ne) &
* pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne),1) &
* pdata%q(ibx ,ne:nep,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
+ sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu, : ) &
* pdata%q(ibx:ibz,ne:nep,nl:nu, : ),1) &
* pdata%q(ibx ,ne:nep,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
if (resistive) then
! diffusion of magnetic energy across the upper X-boundary
!
rterms(20) = rterms(20) &
#if NDIMS == 3
+ sum(cur(2,ne:nep,nl:nu,nb:ne) &
* pdata%q(ibz,ne:nep,nl:nu,nb:ne) &
- cur(3,ne:nep,nl:nu,nb:ne) &
* pdata%q(iby,ne:nep,nl:nu,nb:ne)) * dyz
#else /* NDIMS == 3 */
+ sum(cur(2,ne:nep,nl:nu, : ) &
* pdata%q(ibz,ne:nep,nl:nu, : ) &
- cur(3,ne:nep,nl:nu, : ) &
* pdata%q(iby,ne:nep,nl:nu, : )) * dyz
#endif /* NDIMS == 3 */
end if ! resistivity
end if
end if
if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then
ni = nl ni = nl
nm = ni - 1 nm = ni - 1
@ -1696,7 +1790,7 @@ module user_problem
+ sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2) * dxz + sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! advection of magnetic energy ! advection of magnetic energy across the lower Y-boundary
! !
rterms(15) = rterms(15) & rterms(15) = rterms(15) &
#if NDIMS == 3 #if NDIMS == 3