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 2015-12-13 11:43:51 -02:00
commit 600e93239b
1 changed files with 134 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

233
src/interpolations.F90
View File

@ -184,6 +184,7 @@ module interpolations
if (verbose .and. ng < 2) & if (verbose .and. ng < 2) &
call print_warning("interpolations:initialize_interpolation" & call print_warning("interpolations:initialize_interpolation" &
, "Increase the number of ghost cells (at least 2).") , "Increase the number of ghost cells (at least 2).")
eps = max(1.0d-12, eps)
case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z") case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z")
name_rec = "5th order WENO-Z (Borges et al. 2008)" name_rec = "5th order WENO-Z (Borges et al. 2008)"
reconstruct_states => reconstruct_weno5z reconstruct_states => reconstruct_weno5z
@ -484,12 +485,12 @@ module interpolations
! !
! calculate the left- and right-side interface interpolations ! calculate the left- and right-side interface interpolations
! !
do i = 1, n do i = 2, n - 1
! calculate left and right indices ! calculate left and right indices
! !
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
! calculate left and right side derivatives ! calculate left and right side derivatives
! !
@ -505,11 +506,14 @@ module interpolations
fl(i ) = f(i) + df fl(i ) = f(i) + df
fr(im1) = f(i) - df fr(im1) = f(i) - df
end do ! i = 1, n end do ! i = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = f(1) i = n - 1
fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n) fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
@ -563,12 +567,12 @@ module interpolations
! iterate along the vector ! iterate along the vector
! !
do i = 1, n do i = 2, n - 1
! prepare neighbour indices ! prepare neighbour indices
! !
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
! calculate the left and right derivatives ! calculate the left and right derivatives
! !
@ -625,12 +629,15 @@ module interpolations
! !
fr(im1) = (wp * fp + wm * fm) / ww fr(im1) = (wp * fp + wm * fm) / ww
end do ! i = 1, n end do ! i = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = f (1) i = n - 1
fr(n) = fl(n) fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -686,12 +693,12 @@ module interpolations
! iterate over positions and interpolate states ! iterate over positions and interpolate states
! !
do i = 1, n do i = 2, n - 1
! prepare neighbour indices ! prepare neighbour indices
! !
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
! prepare left and right differences ! prepare left and right differences
! !
@ -709,11 +716,7 @@ module interpolations
! calculate values at i + ½ ! calculate values at i + ½
! !
if (dfr == 0.0d+00) then if (abs(dfr) > eps) then
fl(i) = f(i)
else
! calculate the slope ratio (eq. 2.8 in [1]) ! calculate the slope ratio (eq. 2.8 in [1])
! !
@ -735,15 +738,15 @@ module interpolations
! !
fl(i) = f(i) + dfr * (xl * f1 + xi * f2) fl(i) = f(i) + dfr * (xl * f1 + xi * f2)
else
fl(i) = f(i)
end if end if
! calculate values at i - ½ ! calculate values at i - ½
! !
if (dfl == 0.0d+00) then if (abs(dfl) > eps) then
fr(im1) = f(i)
else
! calculate the slope ratio (eq. 2.8 in [1]) ! calculate the slope ratio (eq. 2.8 in [1])
! !
@ -765,14 +768,21 @@ module interpolations
! !
fr(im1) = f(i) - dfl * (xl * f1 + xi * f2) fr(im1) = f(i) - dfl * (xl * f1 + xi * f2)
else
fr(im1) = f(i)
end if end if
end do ! i = 1, n end do ! i = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = f (1) i = n - 1
fr(n) = fl(n) fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -863,13 +873,13 @@ module interpolations
! iterate along the vector ! iterate along the vector
! !
do i = 1, n do i = 2, n - 1
! prepare neighbour indices ! prepare neighbour indices
! !
im1 = max(1, i - 1)
im2 = max(1, i - 2) im2 = max(1, i - 2)
ip1 = min(n, i + 1) im1 = i - 1
ip1 = i + 1
ip2 = min(n, i + 2) ip2 = min(n, i + 2)
! calculate βₖ (eqs. 9-11 in [1]) ! calculate βₖ (eqs. 9-11 in [1])
@ -928,12 +938,15 @@ module interpolations
! !
fr(im1) = (wl * ql + wr * qr) + wc * qc fr(im1) = (wl * ql + wr * qr) + wc * qc
end do ! i = 1, n end do ! i = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = fr(1) i = n - 1
fr(n) = fl(n) fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -2489,11 +2502,11 @@ module interpolations
! obtain the face values using high order interpolation ! obtain the face values using high order interpolation
! !
do i = 1, n do i = 2, n - 1
im2 = max(1, i - 2) im2 = max(1, i - 2)
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
ip2 = min(n, i + 2) ip2 = min(n, i + 2)
fl(i) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & fl(i) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) &
@ -2503,14 +2516,14 @@ module interpolations
- (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) &
/ 6.0d+01 / 6.0d+01
end do ! i = 1, n end do ! i = 2, n - 1
! apply monotonicity preserving limiting ! apply monotonicity preserving limiting
! !
do i = 1, n do i = 2, n - 1
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
if (dfm(i) * dfp(i) >= 0.0d+00) then if (dfm(i) * dfp(i) >= 0.0d+00) then
sigma = kappa sigma = kappa
@ -2580,12 +2593,15 @@ module interpolations
! !
fr(im1) = fr(i) fr(im1) = fr(i)
end do end do ! n = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = fr(1) i = n - 1
fr(n) = fl(n) fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -2662,10 +2678,10 @@ module interpolations
! prepare the tridiagonal system coefficients for the interior ! prepare the tridiagonal system coefficients for the interior
! !
do i = 1, n do i = ng, n - ng + 1
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
a(i,1) = 3.0d-01 a(i,1) = 3.0d-01
b(i,1) = 6.0d-01 b(i,1) = 6.0d-01
@ -2678,15 +2694,36 @@ module interpolations
r(i,1) = (f(im1) + 1.9d+01 * f(i ) + 1.0d+01 * f(ip1)) / 3.0d+01 r(i,1) = (f(im1) + 1.9d+01 * f(i ) + 1.0d+01 * f(ip1)) / 3.0d+01
r(i,2) = (f(ip1) + 1.9d+01 * f(i ) + 1.0d+01 * f(im1)) / 3.0d+01 r(i,2) = (f(ip1) + 1.9d+01 * f(i ) + 1.0d+01 * f(im1)) / 3.0d+01
end do ! i = 1, n end do ! i = ng, n - ng + 1
! interpolate ghost zones using explicit method (left-side reconstruction) ! interpolate ghost zones using explicit method (left-side reconstruction)
! !
do i = 1, ng do i = 2, ng
im2 = max(1, i - 2) im2 = max(1, i - 2)
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
ip2 = i + 2
a(i,1) = 0.0d+00
b(i,1) = 1.0d+00
c(i,1) = 0.0d+00
r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) &
- (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) &
/ 6.0d+01
end do ! i = 2, ng
a(1,1) = 0.0d+00
b(1,1) = 1.0d+00
c(1,1) = 0.0d+00
r(1,1) = 0.5d+00 * (f(1) + f(2))
do i = n - ng, n - 1
im2 = i - 2
im1 = i - 1
ip1 = i + 1
ip2 = min(n, i + 2) ip2 = min(n, i + 2)
a(i,1) = 0.0d+00 a(i,1) = 0.0d+00
@ -2697,32 +2734,40 @@ module interpolations
- (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) &
/ 6.0d+01 / 6.0d+01
end do ! i = 1, ng end do ! i = n - ng, n - 1
a(n,1) = 0.0d+00
do i = n - ng, n b(n,1) = 1.0d+00
c(n,1) = 0.0d+00
im2 = max(1, i - 2) r(n,1) = f(n)
im1 = max(1, i - 1)
ip1 = min(n, i + 1)
ip2 = min(n, i + 2)
a(i,1) = 0.0d+00
b(i,1) = 1.0d+00
c(i,1) = 0.0d+00
r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) &
- (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) &
/ 6.0d+01
end do ! i = n - ng, n
! interpolate ghost zones using explicit method (right-side reconstruction) ! interpolate ghost zones using explicit method (right-side reconstruction)
! !
do i = 1, ng + 1 do i = 2, ng + 1
im2 = max(1, i - 2) im2 = max(1, i - 2)
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
ip2 = i + 2
a(i,2) = 0.0d+00
b(i,2) = 1.0d+00
c(i,2) = 0.0d+00
r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) &
- (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) &
/ 6.0d+01
end do ! i = 2, ng + 1
a(1,2) = 0.0d+00
b(1,2) = 1.0d+00
c(1,2) = 0.0d+00
r(1,2) = f(1)
do i = n - ng + 1, n - 1
im2 = i - 2
im1 = i - 1
ip1 = i + 1
ip2 = min(n, i + 2) ip2 = min(n, i + 2)
a(i,2) = 0.0d+00 a(i,2) = 0.0d+00
@ -2733,24 +2778,11 @@ module interpolations
- (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) &
/ 6.0d+01 / 6.0d+01
end do ! i = 1, ng + 1 end do ! i = n - ng + 1, n - 1
a(n,2) = 0.0d+00
do i = n - ng + 1, n b(n,2) = 1.0d+00
c(n,2) = 0.0d+00
im2 = max(1, i - 2) r(n,2) = 0.5d+00 * (f(n-1) + f(n))
im1 = max(1, i - 1)
ip1 = min(n, i + 1)
ip2 = min(n, i + 2)
a(i,2) = 0.0d+00
b(i,2) = 1.0d+00
c(i,2) = 0.0d+00
r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) &
- (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) &
/ 6.0d+01
end do ! i = n - ng + 1, n
! solve the tridiagonal system of equations for the left-side interpolation ! solve the tridiagonal system of equations for the left-side interpolation
! !
@ -2758,10 +2790,10 @@ module interpolations
! apply the monotonicity preserving limiting ! apply the monotonicity preserving limiting
! !
do i = 1, n do i = 2, n - 1
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
if (dfm(i) * dfp(i) >= 0.0d+00) then if (dfm(i) * dfp(i) >= 0.0d+00) then
sigma = kappa sigma = kappa
@ -2798,7 +2830,7 @@ module interpolations
end if end if
end do end do ! i = 2, n - 1
! solve the tridiagonal system of equations for the right-side interpolation ! solve the tridiagonal system of equations for the right-side interpolation
! !
@ -2806,10 +2838,10 @@ module interpolations
! apply the monotonicity preserving limiting ! apply the monotonicity preserving limiting
! !
do i = 1, n do i = 2, n - 1
im1 = max(1, i - 1) im1 = i - 1
ip1 = min(n, i + 1) ip1 = i + 1
if (dfm(i) * dfp(i) >= 0.0d+00) then if (dfm(i) * dfp(i) >= 0.0d+00) then
sigma = kappa sigma = kappa
@ -2851,12 +2883,15 @@ module interpolations
! !
fr(im1) = fr(i) fr(im1) = fr(i)
end do end do ! i = 2, n - 1
! update the interpolation of the first and last points ! update the interpolation of the first and last points
! !
fl(1) = fr(1) i = n - 1
fr(n) = fl(n) fl(1) = 0.5d+00 * (f(1) + f(2))
fr(i) = 0.5d+00 * (f(i) + f(n))
fl(n) = f(n)
fr(n) = f(n)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !