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

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2014-11-30 18:51:19 -02:00
commit 0293bd06fd
3 changed files with 2201 additions and 2055 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3907
src/boundaries.F90
View File

File diff suppressed because it is too large Load Diff

82
src/mpitools.F90
View File

@ -49,12 +49,16 @@ module mpitools
! MPI global variables
!
integer(kind=4), save :: comm3d
integer(kind=4), save :: nproc, nprocs, npmax
integer(kind=4), save :: nproc, nprocs, npmax, npairs
integer(kind=4), save, dimension(3) :: pdims, pcoords, pparity
integer(kind=4), save, dimension(3,2) :: pneighs
logical , save, dimension(3) :: periodic
logical , save :: master = .true.
! allocatable array for processor pairs
!
integer(kind=4), dimension(:,:), allocatable, save :: pairs
! by default everything is public
!
public
@ -92,7 +96,11 @@ module mpitools
! local variables
!
#ifdef MPI
integer :: iret
integer :: mprocs, i, j, l, n, iret
! allocatable array for processors order
!
integer(kind=4), dimension(:), allocatable :: procs
#endif /* MPI */
!
!-------------------------------------------------------------------------------
@ -170,6 +178,72 @@ module mpitools
!
npmax = nprocs - 1
! roung up the number of processors to even number
!
mprocs = nprocs + mod(nprocs, 2)
! calculate the number of processor pairs for data exchange
!
npairs = nprocs * npmax
! allocate space for all processor pairs
!
allocate(pairs(npairs, 2))
! allocate space for the processor order
!
allocate(procs(mprocs))
! fill the processor order array
!
procs(:) = (/(l, l = 0, mprocs - 1)/)
! generate processor pairs
!
n = 0
! iterate over turns
!
do l = 1, mprocs - 1
! generate pairs for a given turn
!
do i = 1, mprocs / 2
! calculate the pair for the current processor
!
j = mprocs - i + 1
! continue, if the process number is correct (for odd nprocs case)
!
if (procs(i) < nprocs .and. procs(j) < nprocs) then
! increase the pair number
!
n = n + 1
! substitute the processor numbers for the current pair
!
pairs(n,1:2) = (/ procs(i), procs(j) /)
end if ! max(procs(i), procs(j)) < nprocs
end do ! i = 1, mprocs / 2
! shift elements in the processor order array
!
procs(2:mprocs) = cshift(procs(2:mprocs), -1)
end do ! l = 1, mprocs - 1
! fill out the remaining pairs (swapped)
!
pairs(npairs/2+1:npairs,1:2) = pairs(1:npairs/2,2:1:-1)
! allocate space for the processor order
!
deallocate(procs)
! store the MPI pool handles
!
comm3d = mpi_comm_world
@ -232,6 +306,10 @@ module mpitools
stop
end if
! deallocate space used for processor pairs
!
if (allocated(pairs)) deallocate(pairs)
! stop time accounting for the MPI initialization
!
call stop_timer(imi)

267
src/sources.F90
View File

@ -238,7 +238,7 @@ module sources
!
use blocks , only : block_data
use coordinates , only : im, jm, km
use coordinates , only : ax, ay, az, adx, ady, adz, adxi, adyi, adzi
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
use equations , only : ibx, iby, ibz, ibp
@ -256,10 +256,12 @@ module sources
! local variables
!
integer :: i , j , k
integer :: im1, jm1, km1
integer :: ip1, jp1, kp1
real(kind=8) :: r2, r3, gc, gx, gy, gz
real(kind=8) :: dxi2, dyi2, dzi2, dvx, dvy, dvz, dbx, dby, dbz, dvv
real(kind=8) :: fc, gc
real(kind=8) :: r2, r3, gx, gy, gz
real(kind=8) :: dbx, dby, dbz
real(kind=8) :: dvxdx, dvxdy, dvxdz, divv
real(kind=8) :: dvydx, dvydy, dvydz
real(kind=8) :: dvzdx, dvzdy, dvzdz
! local arrays
!
@ -267,8 +269,8 @@ module sources
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,im,jm,km) :: jc
real(kind=8), dimension(im,jm,km) :: db
real(kind=8), dimension(3,3,im,jm,km) :: tmp
!
!-------------------------------------------------------------------------------
!
@ -351,120 +353,129 @@ module sources
! prepare coordinate increments
!
dxi2 = viscosity * adxi(pdata%meta%level)**2
dyi2 = viscosity * adyi(pdata%meta%level)**2
#if NDIMS == 3
dzi2 = viscosity * adzi(pdata%meta%level)**2
#endif /* NDIMS == 3 */
dh(1) = adx(pdata%meta%level)
dh(2) = ady(pdata%meta%level)
dh(3) = adz(pdata%meta%level)
! iterate over all positions in the YZ plane
! 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))
! iterate over all cells
!
do k = 1, km
#if NDIMS == 3
km1 = max( 1, k - 1)
kp1 = min(km, k + 1)
#endif /* NDIMS == 3 */
do j = 1, jm
jm1 = max( 1, j - 1)
jp1 = min(jm, j + 1)
do i = 1, jm
im1 = max( 1, i - 1)
ip1 = min(im, i + 1)
! calculate second order derivatives of Vx
! prepare the νρ factor
!
dvx = (pdata%q(ivx,ip1,j,k) + pdata%q(ivx,im1,j,k)) &
- 2.0d+00 * pdata%q(ivx,i,j,k)
dvy = (pdata%q(ivx,i,jp1,k) + pdata%q(ivx,i,jm1,k)) &
- 2.0d+00 * pdata%q(ivx,i,j,k)
#if NDIMS == 3
dvz = (pdata%q(ivx,i,j,kp1) + pdata%q(ivx,i,j,km1)) &
- 2.0d+00 * pdata%q(ivx,i,j,k)
#endif /* NDIMS == 3 */
gc = viscosity * pdata%q(idn,i,j,k)
fc = 2.0d+00 * gc
! calculate the source term for Vx
! get the velocity Jacobian elements
!
#if NDIMS == 2
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy)
#endif /* NDIMS == 2 */
#if NDIMS == 3
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy + dzi2 * dvz)
#endif /* NDIMS == 3 */
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
! add viscous source terms to X-momentum equation
! calculate elements of the viscous stress tensor
!
du(imx,i,j,k) = du(imx,i,j,k) + dvv
! add viscous source term to total energy equation
!
if (ien > 0) then
du(ien,i,j,k) = du(ien,i,j,k) + pdata%q(ivx,i,j,k) * dvv
end if
! calculate second order derivatives of Vy
!
dvx = (pdata%q(ivy,ip1,j,k) + pdata%q(ivy,im1,j,k)) &
- 2.0d+00 * pdata%q(ivy,i,j,k)
dvy = (pdata%q(ivy,i,jp1,k) + pdata%q(ivy,i,jm1,k)) &
- 2.0d+00 * pdata%q(ivy,i,j,k)
#if NDIMS == 3
dvz = (pdata%q(ivy,i,j,kp1) + pdata%q(ivy,i,j,km1)) &
- 2.0d+00 * pdata%q(ivy,i,j,k)
#endif /* NDIMS == 3 */
! calculate the source term for Vy
!
#if NDIMS == 2
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy)
#endif /* NDIMS == 2 */
#if NDIMS == 3
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy + dzi2 * dvz)
#endif /* NDIMS == 3 */
! add viscous source terms to Y-momentum equation
!
du(imy,i,j,k) = du(imy,i,j,k) + dvv
! add viscous source term to total energy equation
!
if (ien > 0) then
du(ien,i,j,k) = du(ien,i,j,k) + pdata%q(ivy,i,j,k) * dvv
end if
! calculate second order derivatives of Vz
!
dvx = (pdata%q(ivz,ip1,j,k) + pdata%q(ivz,im1,j,k)) &
- 2.0d+00 * pdata%q(ivz,i,j,k)
dvy = (pdata%q(ivz,i,jp1,k) + pdata%q(ivz,i,jm1,k)) &
- 2.0d+00 * pdata%q(ivz,i,j,k)
#if NDIMS == 3
dvz = (pdata%q(ivz,i,j,kp1) + pdata%q(ivz,i,j,km1)) &
- 2.0d+00 * pdata%q(ivz,i,j,k)
#endif /* NDIMS == 3 */
! calculate the source term for Vz
!
#if NDIMS == 2
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy)
#endif /* NDIMS == 2 */
#if NDIMS == 3
dvv = pdata%q(idn,i,j,k) * (dxi2 * dvx + dyi2 * dvy + dzi2 * dvz)
#endif /* NDIMS == 3 */
! add viscous source terms to Y-momentum equation
!
du(imz,i,j,k) = du(imz,i,j,k) + dvv
! add viscous source term to total energy equation
!
if (ien > 0) then
du(ien,i,j,k) = du(ien,i,j,k) + pdata%q(ivz,i,j,k) * dvv
end if
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, im
end do ! j = 1, jm
end do ! k = 1, km
! 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))
! 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)
! 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))
! 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)
! 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))
! 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)
! add viscous source term to total energy equation
!
if (ien > 0) then
! iterate over all cells
!
do k = 1, km
do j = 1, jm
do i = 1, jm
! 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, im
end do ! j = 1, jm
end do ! k = 1, km
! 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))
! 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)
end if ! ien > 0
end if ! viscosity is not zero
!=== add magnetic field related source terms ===
@ -498,13 +509,13 @@ module sources
! calculate the gradient of divergence potential
!
call gradient(dh(:), pdata%q(ibp,:,:,:), jc(inx:inz,:,:,:))
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,:,:,:) * jc(inx:inz,:,:,:), 1)
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
end if ! ien > 0
@ -525,13 +536,13 @@ module sources
! calculate scalar product of velocity and magnetic field
!
jc(inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) &
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)
!
du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * jc(inx,:,:,:)
du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * tmp(inx,inx,:,:,:)
end if ! ien > 0
@ -545,15 +556,26 @@ module sources
! calculate the Laplace operator of B, i.e. Δ(B)
!
call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:))
call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:))
call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:))
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))
! 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)
! update magnetic field component increments
!
du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:)
du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:)
du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:)
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)
! update energy equation
!
@ -562,20 +584,21 @@ module sources
! add the first resistive source term to the energy equation, i.e.
! d/dt E + ∇.F = η B.[Δ(B)]
!
du(ien,:,:,:) = du(ien,:,:,:) &
+ resistivity * (pdata%q(ibx,:,:,:) * jc(inx,:,:,:) &
+ pdata%q(iby,:,:,:) * jc(iny,:,:,:) &
+ pdata%q(ibz,:,:,:) * jc(inz,:,:,:))
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))
! calculate current density J = ∇xB
!
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:))
call curl(dh(:), pdata%q(ibx:ibz,1:im,1:jm,1:km) &
, tmp(inz,inx:inz,1:im,1:jm,1:km))
! add the second resistive source term to the energy equation, i.e.
! d/dt E + ∇.F = η J²
!
du(ien,:,:,:) = du(ien,:,:,:) &
+ resistivity * sum(jc(:,:,:,:) * jc(:,:,:,:), 1)
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) &
+ resistivity * sum(tmp(inz,inx:inz,1:im,1:jm,1:km)**2, 2)
end if ! energy equation present