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 2020-09-29 12:05:16 -03:00
commit 163edd0025
2 changed files with 40 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
sources/algebra.F90
View File

@ -1093,83 +1093,55 @@ module algebra
! subroutine TRIDIAG:
! ------------------
!
! Subroutine solves the system of tridiagonal linear equations.
! Subroutine solves the system of tridiagonal linear equations using
! Thomas algorithm.
!
! Arguments:
!
! n - the size of the matrix;
! l(:) - the vector of the lower off-diagonal values;
! d(:) - the vector of the diagonal values;
! u(:) - the vector of the upper off-diagonal values;
! r(:) - the vector of the right side values;
! a(:) - the vector of the lower off-diagonal values;
! b(:) - the vector of the diagonal values;
! c(:) - the vector of the upper off-diagonal values;
! d(:) - the vector of the right side values;
! x(:) - the solution vector;
!
! References:
!
! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992
! [1] https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
!
!===============================================================================
!
subroutine tridiag(n, l, d, u, r, x, iret)
subroutine tridiag(n, a, b, c, d, x)
! import external procedures
!
use iso_fortran_env, only : error_unit
! local variables are not implicit by default
!
implicit none
! input/output arguments
!
integer , intent(in) :: n
real(kind=8), dimension(n), intent(in) :: l, d, u, r
real(kind=8), dimension(n), intent(in) :: a, b, c, d
real(kind=8), dimension(n), intent(out) :: x
integer , intent(out) :: iret
! local variables
!
integer :: i, j
integer :: i
real(kind=8) :: t
real(kind=8), dimension(n) :: g
! local parameters
!
character(len=*), parameter :: loc = 'ALGEBRA::tridiag()'
real(kind=8), dimension(n) :: cp, dp
!
!-------------------------------------------------------------------------------
!
! decomposition and forward substitution
!
t = d(1)
x(1) = r(1) / t
cp(1) = c(1) / b(1)
dp(1) = d(1) / b(1)
do i = 2, n
j = i - 1
g(i) = u(j) / t
t = d(i) - l(i) * g(i)
if (abs(t) > 0.0d+00) then
x(i) = (r(i) - l(i) * x(j)) / t
else
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Could not find solution!"
iret = 1
return
end if
t = b(i) - a(i) * cp(i-1)
cp(i) = c(i) / t
dp(i) = (d(i) - a(i) * dp(i-1)) / t
end do
! backsubstitution
! backward substitution
!
x(n) = dp(n)
do i = n - 1, 1, -1
j = i + 1
x(i) = x(i) - g(j) * x(j)
x(i) = dp(i) - cp(i) * x(i+1)
end do
! set return value to success
!
iret = 0
!-------------------------------------------------------------------------------
!
end subroutine tridiag

45
sources/interpolations.F90
View File

@ -2631,7 +2631,6 @@ module interpolations
! local variables
!
integer :: n, i, im1, ip1, im2, ip2
integer :: iret
real(kind=8) :: bl, bc, br, tt
real(kind=8) :: wl, wc, wr, ww
real(kind=8) :: ql, qc, qr
@ -2948,7 +2947,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
! substitute the left-side values
!
@ -2956,7 +2955,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
! substitute the right-side values
!
@ -3027,7 +3026,6 @@ module interpolations
! local variables
!
integer :: n, i, im1, ip1, im2, ip2
integer :: iret
real(kind=8) :: bl, bc, br, tt
real(kind=8) :: wl, wc, wr, ww
real(kind=8) :: ql, qc, qr
@ -3347,7 +3345,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
! substitute the left-side values
!
@ -3355,7 +3353,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
! substitute the right-side values
!
@ -3426,7 +3424,6 @@ module interpolations
! local variables
!
integer :: n, i, im1, ip1, im2, ip2
integer :: iret
real(kind=8) :: bl, bc, br
real(kind=8) :: wl, wc, wr, ww
real(kind=8) :: df, lq, l3, zt
@ -3766,7 +3763,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
! substitute the left-side values
!
@ -3774,7 +3771,7 @@ module interpolations
! solve the tridiagonal system of equations for the left-side interpolation
!
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
! substitute the right-side values
!
@ -4649,7 +4646,7 @@ module interpolations
! local variables
!
integer :: n, i, iret
integer :: n, i
! local arrays for derivatives
!
@ -4718,7 +4715,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -4755,7 +4752,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -4829,7 +4826,7 @@ module interpolations
! local variables
!
integer :: n, i, iret
integer :: n, i
! local arrays for derivatives
!
@ -4899,7 +4896,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -4936,7 +4933,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -5003,7 +5000,7 @@ module interpolations
! local variables
!
integer :: n, i, iret
integer :: n, i
! local arrays for derivatives
!
@ -5078,7 +5075,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -5117,7 +5114,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -5185,7 +5182,7 @@ module interpolations
! local variables
!
integer :: n, i, iret
integer :: n, i
! local arrays for derivatives
!
@ -5261,7 +5258,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -5300,7 +5297,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiting
!
@ -5353,7 +5350,7 @@ module interpolations
real(kind=8), dimension(:), intent(in) :: fc
real(kind=8), dimension(:), intent(out) :: fl, fr
integer :: n, i, iret
integer :: n, i
real(kind=8), dimension(size(fc)) :: fi
real(kind=8), dimension(size(fc)) :: a, b, c
@ -5417,7 +5414,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiter
!
@ -5454,7 +5451,7 @@ module interpolations
! solve the tridiagonal system of equations
!
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
! apply the monotonicity preserving limiter
!