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 2019-03-07 11:14:05 -03:00
commit 99d5cc5069
2 changed files with 217 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
sources/interpolations.F90
View File

@ -292,6 +292,11 @@ module interpolations
interfaces => interfaces_dir
reconstruct_states => reconstruct_mp7
nghosts = max(nghosts, 4)
case ("mp9", "MP9")
name_rec = "9th order Monotonicity Preserving"
interfaces => interfaces_dir
reconstruct_states => reconstruct_mp9
nghosts = max(nghosts, 6)
case ("crmp5", "CRMP5")
name_rec = "5th order Compact Monotonicity Preserving"
interfaces => interfaces_dir
@ -3988,6 +3993,146 @@ module interpolations
!
!===============================================================================
!
! subroutine RECONSTRUCT_MP9:
! --------------------------
!
! Subroutine reconstructs the interface states using the ninth order
! Monotonicity Preserving (MP) method.
!
! Arguments are described in subroutine reconstruct().
!
! References:
!
! [1] Suresh, A. & Huynh, H. T.,
! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
! Time Stepping"
! Journal on Computational Physics,
! 1997, vol. 136, pp. 83-99,
! http://dx.doi.org/10.1006/jcph.1997.5745
! [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
! "A 5th order monotonicity-preserving upwind compact difference
! scheme",
! Science China Physics, Mechanics and Astronomy,
! Volume 54, Issue 3, pp. 511-522,
! http://dx.doi.org/10.1007/s11433-010-4220-x
!
!===============================================================================
!
subroutine reconstruct_mp9(h, fc, fl, fr)
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8) , intent(in) :: h
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
! local variables
!
integer :: n, i
! local arrays for derivatives
!
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: u
! local parameters
!
real(kind=8), dimension(9), parameter :: &
ce9 = (/ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03, &
1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 /) / 2.520d+03
real(kind=8), dimension(7), parameter :: &
ce7 = (/-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02, &
-3.800d+01, 4.000d+00 /) / 4.200d+02
real(kind=8), dimension(5), parameter :: &
ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 /) / 6.0d+01
real(kind=8), dimension(3), parameter :: &
ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
!
!-------------------------------------------------------------------------------
!
! get the input vector length
!
n = size(fc)
!! === left-side interpolation ===
!!
! reconstruct the interface state using the 9th order interpolation
!
do i = 5, n - 4
u(i) = sum(ce9(:) * fc(i-4:i+4))
end do
! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
u( 1) = sum(ce2(:) * fc( 1: 2))
u( 2) = sum(ce3(:) * fc( 1: 3))
u( 3) = sum(ce5(:) * fc( 1: 5))
u( 4) = sum(ce7(:) * fc( 1: 7))
u(n-3) = sum(ce7(:) * fc(n-6: n))
u(n-2) = sum(ce5(:) * fc(n-4: n))
u(n-1) = sum(ce3(:) * fc(n-2: n))
u(n ) = fc(n )
! apply the monotonicity preserving limiting
!
call mp_limiting(fc(:), u(:))
! copy the interpolation to the respective vector
!
fl(1:n) = u(1:n)
!! === right-side interpolation ===
!!
! invert the cell-centered value vector
!
fi(1:n) = fc(n:1:-1)
! reconstruct the interface state using the 9th order interpolation
!
do i = 5, n - 4
u(i) = sum(ce9(:) * fi(i-4:i+4))
end do
! interpolate the interface state of the ghost zones using the interpolations
! of lower orders
!
u( 1) = sum(ce2(:) * fi( 1: 2))
u( 2) = sum(ce3(:) * fi( 1: 3))
u( 3) = sum(ce5(:) * fi( 1: 5))
u( 4) = sum(ce7(:) * fi( 1: 7))
u(n-3) = sum(ce7(:) * fi(n-6: n))
u(n-2) = sum(ce5(:) * fi(n-4: n))
u(n-1) = sum(ce3(:) * fi(n-2: n))
u(n ) = fi(n )
! apply the monotonicity preserving limiting
!
call mp_limiting(fi(:), u(:))
! copy the interpolation to the respective vector
!
fr(1:n-1) = u(n-1:1:-1)
! update the interpolation of the first and last points
!
i = n - 1
fl(1) = 0.5d+00 * (fc(1) + fc(2))
fr(i) = 0.5d+00 * (fc(i) + fc(n))
fl(n) = fc(n)
fr(n) = fc(n)
!-------------------------------------------------------------------------------
!
end subroutine reconstruct_mp9
!
!===============================================================================
!
! subroutine RECONSTRUCT_CRMP5:
! ----------------------------
!

118
sources/mesh.F90
View File

@ -47,6 +47,10 @@ module mesh
!
integer, save :: funit = 11
! allocatable array for prolongation
!
real(kind=8), dimension(:,:,:,:), allocatable :: up
! by default everything is private
!
private
@ -86,8 +90,10 @@ module mesh
! import external procedures and variables
!
use coordinates, only : toplev
use mpitools , only : master, nprocs
use coordinates , only : ni => ncells, ng => nghosts, toplev
use equations , only : nv
use iso_fortran_env, only : error_unit
use mpitools , only : master, nprocs
! local variables are not implicit by default
!
@ -102,11 +108,15 @@ module mesh
! local variables
!
character(len=64) :: fn
integer(kind=4) :: i, j, k, l, n
integer :: l, n
! local arrays
!
integer(kind=4), dimension(3) :: bm, rm, dm
integer, dimension(3) :: pm = 1
! local parameters
!
character(len=*), parameter :: loc = 'MESH::initialize_mesh()'
!
!-------------------------------------------------------------------------------
!
@ -170,6 +180,19 @@ module mesh
end if ! master
! prepare dimensions
!
pm(1:NDIMS) = 2 * (ni + ng)
! allocate array for the prolongation array
!
allocate(up(nv, pm(1), pm(2), pm(3)), stat = status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Prolongation array could not be allocated!"
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
@ -197,7 +220,8 @@ module mesh
! import external procedures and variables
!
use mpitools, only : master
use iso_fortran_env, only : error_unit
use mpitools , only : master
! local variables are not implicit by default
!
@ -206,6 +230,10 @@ module mesh
! subroutine arguments
!
integer, intent(out) :: status
! local parameters
!
character(len=*), parameter :: loc = 'MESH::finalize_mesh()'
!
!-------------------------------------------------------------------------------
!
@ -221,13 +249,17 @@ module mesh
! only master posses a handler to the mesh statistics file
!
if (master) then
if (master) close(funit)
! close the statistics file
! deallocate prolongation array
!
close(funit)
end if ! master
if (allocated(up)) then
deallocate(up, stat = status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Prolongation array could not be deallocated!"
end if
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
@ -1062,17 +1094,18 @@ module mesh
!
! Arguments:
!
! pmeta - the input meta block;
! pmeta - the input meta block;
! status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
subroutine prolong_block(pmeta)
subroutine prolong_block(pmeta, status)
! import external procedures and variables
!
use blocks , only : block_meta, block_data, nchildren
use coordinates , only : ni => ncells, nn => bcells
use coordinates , only : ng => nghosts, nh => nghosts_half
use coordinates , only : nh => nghosts_half
use coordinates , only : nb, ne
use equations , only : nv, idn, ien
use interpolations , only : limiter_prol
@ -1085,6 +1118,7 @@ module mesh
! input arguments
!
type(block_meta), pointer, intent(inout) :: pmeta
integer , intent(out) :: status
! local variables
!
@ -1100,13 +1134,8 @@ module mesh
! local arrays
!
integer , dimension(3) :: dm = 1
real(kind=8), dimension(NDIMS) :: du
! local allocatable arrays
!
real(kind=8), dimension(:,:,:,:), allocatable :: u
! local parameters
!
character(len=*), parameter :: loc = 'MESH::prolong_block()'
@ -1119,18 +1148,14 @@ module mesh
call start_timer(imp)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! assign the pdata pointer
!
pdata => pmeta%data
! prepare dimensions
!
dm(1:NDIMS) = 2 * (ni + ng)
! allocate array for the expanded arrays
!
allocate(u(nv, dm(1), dm(2), dm(3)))
! prepare indices
!
il = nb - nh
@ -1186,6 +1211,7 @@ module mesh
, "Positive variable is not positive for block ID:" &
, pmeta%id, " at ( ", i, j, k, " )"
du(:) = 0.0d+00
status = 1
end if
end if
@ -1194,11 +1220,10 @@ module mesh
#if NDIMS == 2
du1 = du(1) + du(2)
du2 = du(1) - du(2)
u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2
u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du1
up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2
up(p,ip,jp,kc) = pdata%u(p,i,j,k) + du1
#endif /* NDIMS == 2 */
#if NDIMS == 3
@ -1206,14 +1231,14 @@ module mesh
du2 = du(1) - du(2) - du(3)
du3 = du(1) - du(2) + du(3)
du4 = du(1) + du(2) - du(3)
u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3
u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4
u(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4
u(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3
u(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2
u(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1
up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3
up(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4
up(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4
up(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3
up(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2
up(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1
#endif /* NDIMS == 3 */
end do
end do
@ -1255,18 +1280,14 @@ module mesh
! copy data to the current child
!
#if NDIMS == 2
pchild%data%u(1:nv,:,:,:) = u(1:nv,il:iu,jl:ju, : )
pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju, : )
#endif /* NDIMS == 2 */
#if NDIMS == 3
pchild%data%u(1:nv,:,:,:) = u(1:nv,il:iu,jl:ju,kl:ku)
pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju,kl:ku)
#endif /* NDIMS == 3 */
end do ! nchildren
! deallocate local arrays
!
if (allocated(u)) deallocate(u)
#ifdef PROFILE
! stop accounting time for the block prolongation
!
@ -2322,7 +2343,12 @@ module mesh
! perform the data prolongation
!
call prolong_block(pparent)
call prolong_block(pparent, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot prolong current data block!"
go to 100
end if
! remove the data block associated with the new parent
!