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

SOURCES: Use bcells instead of im, jm, km.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2019-02-05 11:17:59 -02:00
parent e63d3cff0d
commit 969afab57f
1 changed files with 76 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

151
sources/sources.F90
View File

@ -268,7 +268,7 @@ module sources
! include external variables
!
use blocks , only : block_data
use coordinates , only : im, jm, km
use coordinates , only : nb => bcells
use coordinates , only : ax, ay, az, adx, ady, adz
use equations , only : nv, inx, iny, inz
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
@ -283,13 +283,13 @@ module sources
! subroutine arguments
!
type(block_data), pointer , intent(inout) :: pdata
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(nv,im,jm,km), intent(inout) :: du
type(block_data), pointer , intent(inout) :: pdata
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(:,:,:,:), intent(inout) :: du
! local variables
!
integer :: i , j , k
integer :: i, j, k = 1
real(kind=8) :: fc, gc
real(kind=8) :: gx, gy, gz
real(kind=8) :: dbx, dby, dbz
@ -300,11 +300,17 @@ module sources
! local arrays
!
real(kind=8), dimension(3) :: ga, dh
real(kind=8), dimension(im) :: x
real(kind=8), dimension(jm) :: y
real(kind=8), dimension(km) :: z
real(kind=8), dimension(im,jm,km) :: db
real(kind=8), dimension(3,3,im,jm,km) :: tmp
real(kind=8), dimension(nb) :: x
real(kind=8), dimension(nb) :: y
#if NDIMS == 3
real(kind=8), dimension(nb) :: z
real(kind=8), dimension(nb,nb,nb) :: db
real(kind=8), dimension(3,3,nb,nb,nb) :: tmp
#else /* NDIMS == 3 */
real(kind=8), dimension( 1) :: z
real(kind=8), dimension(nb,nb, 1) :: db
real(kind=8), dimension(3,3,nb,nb, 1) :: tmp
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
@ -320,17 +326,19 @@ module sources
! prepare block coordinates
!
x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km)
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */
! iterate over all positions in the YZ plane
!
do k = 1, km
do j = 1, jm
do i = 1, im
#if NDIMS == 3
do k = 1, nb
#endif /* NDIMS == 3 */
do j = 1, nb
do i = 1, nb
! get gravitational acceleration components
!
@ -368,9 +376,11 @@ module sources
end if
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
end do ! i = 1, nb
end do ! j = 1, nb
#if NDIMS == 3
end do ! k = 1, nb
#endif /* NDIMS == 3 */
end if ! gravity enabled
@ -386,18 +396,17 @@ module sources
! calculate the velocity Jacobian
!
call gradient(dh(:), pdata%q(ivx,1:im,1:jm,1:km) &
, tmp(inx,inx:inz,1:im,1:jm,1:km))
call gradient(dh(:), pdata%q(ivy,1:im,1:jm,1:km) &
, tmp(iny,inx:inz,1:im,1:jm,1:km))
call gradient(dh(:), pdata%q(ivz,1:im,1:jm,1:km) &
, tmp(inz,inx:inz,1:im,1:jm,1:km))
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
!
do k = 1, km
do j = 1, jm
do i = 1, im
#if NDIMS == 3
do k = 1, nb
#endif /* NDIMS == 3 */
do j = 1, nb
do i = 1, nb
! prepare the νρ factor
!
@ -429,36 +438,35 @@ module sources
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, im
end do ! j = 1, jm
end do ! k = 1, km
end do ! i = 1, nb
end do ! j = 1, nb
#if NDIMS == 3
end do ! k = 1, nb
#endif /* NDIMS == 3 */
! calculate the divergence of the first tensor row
!
call divergence(dh(:), tmp(inx,inx:inz,1:im,1:jm,1:km) &
, db(1:im,1:jm,1:km))
call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))
! add viscous source terms to the X momentum equation
!
du(imx,1:im,1:jm,1:km) = du(imx,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
du(imx,:,:,:) = du(imx,:,:,:) + db(:,:,:)
! calculate the divergence of the second tensor row
!
call divergence(dh(:), tmp(iny,inx:inz,1:im,1:jm,1:km) &
, db(1:im,1:jm,1:km))
call divergence(dh(:), tmp(iny,inx:inz,:,:,:), db(:,:,:))
! add viscous source terms to the Y momentum equation
!
du(imy,1:im,1:jm,1:km) = du(imy,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
du(imy,:,:,:) = du(imy,:,:,:) + db(:,:,:)
! calculate the divergence of the third tensor row
!
call divergence(dh(:), tmp(inz,inx:inz,1:im,1:jm,1:km) &
, db(1:im,1:jm,1:km))
call divergence(dh(:), tmp(inz,inx:inz,:,:,:), db(:,:,:))
! add viscous source terms to the Z momentum equation
!
du(imz,1:im,1:jm,1:km) = du(imz,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
du(imz,:,:,:) = du(imz,:,:,:) + db(:,:,:)
! add viscous source term to total energy equation
!
@ -466,9 +474,11 @@ module sources
! iterate over all cells
!
do k = 1, km
do j = 1, jm
do i = 1, im
#if NDIMS == 3
do k = 1, nb
#endif /* NDIMS == 3 */
do j = 1, nb
do i = 1, nb
! calculate scalar product of v and viscous stress tensor τ
!
@ -488,18 +498,19 @@ module sources
tmp(inx,iny,i,j,k) = gy
tmp(inx,inz,i,j,k) = gz
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
end do ! i = 1, nb
end do ! j = 1, nb
#if NDIMS == 3
end do ! k = 1, nb
#endif /* NDIMS == 3 */
! calculate the divergence of (v.τ)
!
call divergence(dh(:), tmp(inx,inx:inz,1:im,1:jm,1:km) &
, db(1:im,1:jm,1:km))
call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))
! update the energy increment
!
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
du(ien,:,:,:) = du(ien,:,:,:) + db(:,:,:)
end if ! ien > 0
@ -583,26 +594,19 @@ module sources
! calculate the Laplace operator of B, i.e. Δ(B)
!
call laplace(dh(:), pdata%q(ibx,1:im,1:jm,1:km) &
, tmp(inx,inx,1:im,1:jm,1:km))
call laplace(dh(:), pdata%q(iby,1:im,1:jm,1:km) &
, tmp(inx,iny,1:im,1:jm,1:km))
call laplace(dh(:), pdata%q(ibz,1:im,1:jm,1:km) &
, tmp(inx,inz,1:im,1:jm,1:km))
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,1:im,1:jm,1:km) = &
resistivity * tmp(inx,inx:inz,1:im,1:jm,1:km)
tmp(iny,inx:inz,:,:,:) = resistivity * tmp(inx,inx:inz,:,:,:)
! update magnetic field component increments
!
du(ibx,1:im,1:jm,1:km) = du(ibx,1:im,1:jm,1:km) &
+ tmp(iny,inx,1:im,1:jm,1:km)
du(iby,1:im,1:jm,1:km) = du(iby,1:im,1:jm,1:km) &
+ tmp(iny,iny,1:im,1:jm,1:km)
du(ibz,1:im,1:jm,1:km) = du(ibz,1:im,1:jm,1:km) &
+ tmp(iny,inz,1:im,1:jm,1:km)
du(ibx,:,:,:) = du(ibx,:,:,:) + tmp(iny,inx,:,:,:)
du(iby,:,:,:) = du(iby,:,:,:) + tmp(iny,iny,:,:,:)
du(ibz,:,:,:) = du(ibz,:,:,:) + tmp(iny,inz,:,:,:)
! update energy equation
!
@ -611,27 +615,24 @@ module sources
! add the first resistive source term to the energy equation, i.e.
! d/dt E + ∇.F = η B.[Δ(B)]
!
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) &
+ (pdata%q(ibx,1:im,1:jm,1:km) * tmp(iny,inx,1:im,1:jm,1:km) &
+ pdata%q(iby,1:im,1:jm,1:km) * tmp(iny,iny,1:im,1:jm,1:km) &
+ pdata%q(ibz,1:im,1:jm,1:km) * tmp(iny,inz,1:im,1:jm,1:km))
du(ien,:,:,:) = 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,1:im,1:jm,1:km) &
, tmp(inz,inx:inz,1:im,1:jm,1:km))
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inz,inx:inz,:,:,:))
! calculate J²
!
db(1:im,1:jm,1:km) = tmp(inz,inx,1:im,1:jm,1:km)**2 &
+ tmp(inz,iny,1:im,1:jm,1:km)**2 &
+ tmp(inz,inz,1:im,1:jm,1:km)**2
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²
!
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) &
+ resistivity * db(1:im,1:jm,1:km)
du(ien,:,:,:) = du(ien,:,:,:) + resistivity * db(:,:,:)
end if ! energy equation present
@ -641,7 +642,7 @@ module sources
! add user defined source terms
!
call update_sources_user(pdata, t, dt, du(1:nv,1:im,1:jm,1:km))
call update_sources_user(pdata, t, dt, du(:,:,:,:))
#ifdef PROFILE
! stop accounting time for source terms