From eb118745a198dd852df4fec257de6a6af6457867 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 08:03:08 -0200 Subject: [PATCH 1/6] ALGEBRA: Add tridiagonal linear system solver. Signed-off-by: Grzegorz Kowal --- src/algebra.F90 | 75 +++++++++++++++++++++++++++++++++++++++++++++++++ src/makefile | 2 +- 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/src/algebra.F90 b/src/algebra.F90 index 6086965..9086059 100644 --- a/src/algebra.F90 +++ b/src/algebra.F90 @@ -42,6 +42,7 @@ module algebra public :: quadratic, quadratic_normalized public :: cubic, cubic_normalized public :: quartic + public :: tridiag !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -1080,6 +1081,80 @@ module algebra !------------------------------------------------------------------------------- ! end subroutine sort +! +!=============================================================================== +! +! subroutine TRIDIAG: +! ------------------ +! +! Subroutine solves the system of tridiagonal linear equations. +! +! 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; +! 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 +! +!=============================================================================== +! + subroutine tridiag(n, l, d, u, r, x) + +! import external procedures +! + use error, only : print_error + +! 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(out) :: x + +! local variables +! + integer :: i, j + real(kind=8) :: t + real(kind=8), dimension(n) :: g +! +!------------------------------------------------------------------------------- +! +! decomposition and forward substitution +! + t = d(1) + x(1) = r(1) / t + do i = 2, n + j = i - 1 + g(i) = u(j) / t + t = d(i) - l(i) * g(i) + if (t == 0.0d+00) then + call print_error("algebra::tridiag", "solution failed!") + stop + end if + x(i) = (r(i) - l(i) * x(j)) / t + end do + +! backsubstitution +! + do i = n - 1, 1, -1 + j = i + 1 + x(i) = x(i) - g(j) * x(j) + end do + +!------------------------------------------------------------------------------- +! + end subroutine tridiag !=============================================================================== ! diff --git a/src/makefile b/src/makefile index dc313a9..7bc4d6c 100644 --- a/src/makefile +++ b/src/makefile @@ -120,7 +120,7 @@ clean-all: clean-bak clean-data clean-exec clean-logs clean-modules \ #------------------------------------------------------------------------------- -algebra.o : algebra.F90 constants.o +algebra.o : algebra.F90 constants.o error.o blocks.o : blocks.F90 error.o timers.o boundaries.o : boundaries.F90 blocks.o coordinates.o equations.o error.o \ interpolations.o mpitools.o timers.o From a00dc32c71c4897d90a21eb76830b8ab79daee54 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 08:21:39 -0200 Subject: [PATCH 2/6] INTERPOLATIONS: Use tridiag() solver in CRMP5 reconstruction. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 34 ++++++++-------------------------- src/makefile | 2 +- 2 files changed, 9 insertions(+), 27 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 777075f..a7fd155 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -2650,6 +2650,10 @@ module interpolations ! subroutine reconstruct_crmp5(n, h, f, fl, fr) +! include external procedures +! + use algebra , only : tridiag + ! local variables are not implicit by default ! implicit none @@ -2666,12 +2670,12 @@ module interpolations integer :: i, im1, ip1, im2, ip2 real(kind=8) :: df, ds, dc0, dc4, dm1, dp1, dml, dmr real(kind=8) :: flc, fmd, fmp, fmn, fmx, ful - real(kind=8) :: sigma, bt + real(kind=8) :: sigma ! local arrays for derivatives ! real(kind=8), dimension(n) :: dfm, dfp - real(kind=8), dimension(n) :: u, g + real(kind=8), dimension(n) :: u real(kind=8), dimension(n,2) :: a, b, c, r ! !------------------------------------------------------------------------------- @@ -2780,18 +2784,7 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - bt = b(1,1) - u(1) = r(1,1) / bt - do i = 2, n - im1 = i - 1 - g(i) = c(im1,1) / bt - bt = b(i,1) - a(i,1) * g(i) - u(i) = (r(i,1) - a(i,1) * u(im1)) / bt - end do - do i = n - 1, 1, -1 - ip1 = i + 1 - u(i) = u(i) - g(ip1) * u(ip1) - end do + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) ! apply the monotonicity preserving limiting ! @@ -2839,18 +2832,7 @@ module interpolations ! solve the tridiagonal system of equations for the right-side interpolation ! - bt = b(n,2) - u(n) = r(n,2) / bt - do i = n - 1, 1, -1 - ip1 = i + 1 - g(i) = a(ip1,2) / bt - bt = b(i,2) - c(i,2) * g(i) - u(i) = (r(i,2) - c(i,2) * u(ip1)) / bt - end do - do i = 2, n - im1 = i - 1 - u(i) = u(i) - g(im1) * u(im1) - end do + call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n)) ! apply the monotonicity preserving limiting ! diff --git a/src/makefile b/src/makefile index 7bc4d6c..fb1e425 100644 --- a/src/makefile +++ b/src/makefile @@ -138,7 +138,7 @@ evolution.o : evolution.F90 blocks.o boundaries.o coordinates.o \ domains.o : domains.F90 blocks.o boundaries.o coordinates.o parameters.o integrals.o : integrals.F90 blocks.o coordinates.o equations.o error.o \ evolution.o mpitools.o parameters.o timers.o -interpolations.o : interpolations.F90 blocks.o coordinates.o error.o \ +interpolations.o : interpolations.F90 algebra.o blocks.o coordinates.o error.o \ parameters.o timers.o io.o : io.F90 blocks.o coordinates.o equations.o error.o \ evolution.o mesh.o mpitools.o random.o refinement.o timers.o From 52f520f2156918400661b77e3c27ab3c168e4647 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 09:14:48 -0200 Subject: [PATCH 3/6] INTERPOLATIONS: Use tridiag() solver in CRWENO5Z reconstruction. Also fix the matrix coefficient calculation. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 76 ++++++++++++++++++------------------------ 1 file changed, 33 insertions(+), 43 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index a7fd155..58c8643 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -1315,6 +1315,10 @@ module interpolations ! subroutine reconstruct_crweno5z(n, h, f, fl, fr) +! include external procedures +! + use algebra , only : tridiag + ! local variables are not implicit by default ! implicit none @@ -1337,7 +1341,7 @@ module interpolations ! real(kind=8), dimension(n) :: dfm, dfp, df2 real(kind=8), dimension(n) :: al, ac, ar - real(kind=8), dimension(n) :: u, g + real(kind=8), dimension(n) :: u real(kind=8), dimension(n,2) :: a, b, c, r ! smoothness indicator coefficients @@ -1353,7 +1357,7 @@ module interpolations ! implicit method coefficients ! - real(kind=8), parameter :: dq = 2.5d-01 + real(kind=8), parameter :: dq = 5.0d-01 ! interpolation coefficients ! @@ -1431,20 +1435,20 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,1) = 2.0d+00 * wl + 5.0d-01 * wc - c(i,1) = 5.0d-01 * wr + a(i,1) = 2.0d+00 * wl + wc + b(i,1) = wl + 2.0d+00 * (wc + wr) + c(i,1) = wr ! prepare right hand side of tridiagonal equation ! - r(i,1) = ( 2.0d+00 * wl * f(im1) & - + (1.0d+01 * wl + 5.0d+00 * wc + wr) * f(i ) & - + ( wc + 5.0d+00 * wr) * f(ip1)) * dq + r(i,1) = (wl * f(im1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(ip1)) * dq ! calculate weights ! - wl = cr * al(i) + wl = cl * ar(i) wc = cc * ac(i) - wr = cl * ar(i) + wr = cr * al(i) ww = (wl + wr) + wc wl = wl / ww wr = wr / ww @@ -1452,14 +1456,14 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,2) = 5.0d-01 * wl - c(i,2) = 5.0d-01 * wc + 2.0d+00 * wr + a(i,2) = wr + b(i,2) = wl + 2.0d+00 * (wc + wr) + c(i,2) = 2.0d+00 * wl + wc ! prepare right hand side of tridiagonal equation ! - r(i,2) = ((5.0d+00 * wl + wc ) * f(im1) & - + ( wl + 5.0d+00 * wc + 1.0d+01 * wr) * f(i ) & - + 2.0d+00 * wr * f(ip1)) * dq + r(i,2) = (wl * f(ip1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(im1)) * dq end do ! i = 1, n @@ -1497,6 +1501,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -1536,6 +1541,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -1575,6 +1581,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -1614,6 +1621,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -1621,37 +1629,19 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - ib = 1 - ie = n - tt = 1.0d+00 - u(ib) = r(ib,1) - do i = ib + 1, ie - im1 = i - 1 - g(i) = c(im1,1) / tt - tt = 1.0d+00 - a(i,1) * g(i) - u(i) = (r(i,1) - a(i,1) * u(im1)) / tt - end do - do i = ie - 1, ib, -1 - ip1 = i + 1 - u(i) = u(i) - g(ip1) * u(ip1) - end do - fl(ib:ie) = u(ib:ie) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) -! solve the tridiagonal system of equations for the right-side interpolation +! substitute the left-side values ! - tt = 1.0d+00 - u(ie) = r(ie,2) - do i = ie - 1, ib, -1 - ip1 = i + 1 - g(i) = a(ip1,2) / tt - tt = 1.0d+00 - c(i,2) * g(i) - u(i) = (r(i,2) - c(i,2) * u(ip1)) / tt - end do - do i = ib + 1, ie - im1 = i - 1 - u(i) = u(i) - g(im1) * u(im1) - end do - fr(ib:ie-1) = u(ib+1:ie) + fl(1:n ) = u(1:n) + +! 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)) + +! substitute the right-side values +! + fr(1:n-1) = u(2:n) ! update the interpolation of the first and last points ! From e3d91d48f1feb65242495ce63bc5e62d8705b38f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 09:23:13 -0200 Subject: [PATCH 4/6] INTERPOLATIONS: Use tridiag() solver in CRWENO5YC reconstruction. Also fix the matrix coefficient calculation. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 76 ++++++++++++++++++------------------------ 1 file changed, 33 insertions(+), 43 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 58c8643..44c9c20 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -1688,6 +1688,10 @@ module interpolations ! subroutine reconstruct_crweno5yc(n, h, f, fl, fr) +! include external procedures +! + use algebra , only : tridiag + ! local variables are not implicit by default ! implicit none @@ -1710,7 +1714,7 @@ module interpolations ! real(kind=8), dimension(n) :: dfm, dfp, df2 real(kind=8), dimension(n) :: al, ac, ar - real(kind=8), dimension(n) :: u, g + real(kind=8), dimension(n) :: u real(kind=8), dimension(n,2) :: a, b, c, r ! smoothness indicator coefficients @@ -1726,7 +1730,7 @@ module interpolations ! implicit method coefficients ! - real(kind=8), parameter :: dq = 2.5d-01 + real(kind=8), parameter :: dq = 5.0d-01 ! interpolation coefficients ! @@ -1807,20 +1811,20 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,1) = 2.0d+00 * wl + 5.0d-01 * wc - c(i,1) = 5.0d-01 * wr + a(i,1) = 2.0d+00 * wl + wc + b(i,1) = wl + 2.0d+00 * (wc + wr) + c(i,1) = wr ! prepare right hand side of tridiagonal equation ! - r(i,1) = ( 2.0d+00 * wl * f(im1) & - + (1.0d+01 * wl + 5.0d+00 * wc + wr) * f(i ) & - + ( wc + 5.0d+00 * wr) * f(ip1)) * dq + r(i,1) = (wl * f(im1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(ip1)) * dq ! calculate weights ! - wl = cr * al(i) + wl = cl * ar(i) wc = cc * ac(i) - wr = cl * ar(i) + wr = cr * al(i) ww = (wl + wr) + wc wl = wl / ww wr = wr / ww @@ -1828,14 +1832,14 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,2) = 5.0d-01 * wl - c(i,2) = 5.0d-01 * wc + 2.0d+00 * wr + a(i,2) = wr + b(i,2) = wl + 2.0d+00 * (wc + wr) + c(i,2) = 2.0d+00 * wl + wc ! prepare right hand side of tridiagonal equation ! - r(i,2) = ((5.0d+00 * wl + wc ) * f(im1) & - + ( wl + 5.0d+00 * wc + 1.0d+01 * wr) * f(i ) & - + 2.0d+00 * wr * f(ip1)) * dq + r(i,2) = (wl * f(ip1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(im1)) * dq end do ! i = 1, n @@ -1873,6 +1877,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -1912,6 +1917,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -1951,6 +1957,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -1990,6 +1997,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -1997,37 +2005,19 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - ib = 1 - ie = n - tt = 1.0d+00 - u(ib) = r(ib,1) - do i = ib + 1, ie - im1 = i - 1 - g(i) = c(im1,1) / tt - tt = 1.0d+00 - a(i,1) * g(i) - u(i) = (r(i,1) - a(i,1) * u(im1)) / tt - end do - do i = ie - 1, ib, -1 - ip1 = i + 1 - u(i) = u(i) - g(ip1) * u(ip1) - end do - fl(ib:ie) = u(ib:ie) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) -! solve the tridiagonal system of equations for the right-side interpolation +! substitute the left-side values ! - tt = 1.0d+00 - u(ie) = r(ie,2) - do i = ie - 1, ib, -1 - ip1 = i + 1 - g(i) = a(ip1,2) / tt - tt = 1.0d+00 - c(i,2) * g(i) - u(i) = (r(i,2) - c(i,2) * u(ip1)) / tt - end do - do i = ib + 1, ie - im1 = i - 1 - u(i) = u(i) - g(im1) * u(im1) - end do - fr(ib:ie-1) = u(ib+1:ie) + fl(1:n ) = u(1:n) + +! 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)) + +! substitute the right-side values +! + fr(1:n-1) = u(2:n) ! update the interpolation of the first and last points ! From c382b714ad50e19e4e38795816c9b4c2510e2278 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 09:53:34 -0200 Subject: [PATCH 5/6] INTERPOLATIONS: Use tridiag() solver in CRWENO5NS reconstruction. Also fix the matrix coefficient calculation. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 98 +++++++++++++++++++----------------------- 1 file changed, 44 insertions(+), 54 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 44c9c20..4cac6f3 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -2064,6 +2064,10 @@ module interpolations ! subroutine reconstruct_crweno5ns(n, h, f, fl, fr) +! include external procedures +! + use algebra , only : tridiag + ! local variables are not implicit by default ! implicit none @@ -2077,7 +2081,7 @@ module interpolations ! local variables ! - integer :: i, im1, ip1, im2, ip2, ib, ie + integer :: i, im1, ip1, im2, ip2 real(kind=8) :: bl, bc, br, tt real(kind=8) :: wl, wc, wr, ww real(kind=8) :: df, lq, l3, zt @@ -2087,7 +2091,7 @@ module interpolations ! real(kind=8), dimension(n) :: dfm, dfp, df2 real(kind=8), dimension(n,2) :: al, ac, ar - real(kind=8), dimension(n) :: u, g + real(kind=8), dimension(n) :: u real(kind=8), dimension(n,2) :: a, b, c, r ! the free parameter for smoothness indicators (see eq. 3.6 in [3]) @@ -2103,7 +2107,7 @@ module interpolations ! implicit method coefficients ! - real(kind=8), parameter :: dq = 2.5d-01 + real(kind=8), parameter :: dq = 5.0d-01 ! 3rd order interpolation coefficients for three stencils ! @@ -2203,20 +2207,20 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,1) = 2.0d+00 * wl + 5.0d-01 * wc - c(i,1) = 5.0d-01 * wr + a(i,1) = 2.0d+00 * wl + wc + b(i,1) = wl + 2.0d+00 * (wc + wr) + c(i,1) = wr ! prepare right hand side of tridiagonal equation ! - r(i,1) = ( 2.0d+00 * wl * f(im1) & - + (1.0d+01 * wl + 5.0d+00 * wc + wr) * f(i ) & - + ( wc + 5.0d+00 * wr) * f(ip1)) * dq + r(i,1) = (wl * f(im1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(ip1)) * dq ! calculate weights ! - wl = cr * al(i,2) + wl = cl * ar(i,2) wc = cc * ac(i,2) - wr = cl * ar(i,2) + wr = cr * al(i,2) ww = (wl + wr) + wc wl = wl / ww wr = wr / ww @@ -2224,14 +2228,14 @@ module interpolations ! calculate tridiagonal matrix coefficients ! - a(i,2) = 5.0d-01 * wl - c(i,2) = 5.0d-01 * wc + 2.0d+00 * wr + a(i,2) = wr + b(i,2) = wl + 2.0d+00 * (wc + wr) + c(i,2) = 2.0d+00 * wl + wc ! prepare right hand side of tridiagonal equation ! - r(i,2) = ((5.0d+00 * wl + wc ) * f(im1) & - + ( wl + 5.0d+00 * wc + 1.0d+01 * wr) * f(i ) & - + 2.0d+00 * wr * f(ip1)) * dq + r(i,2) = (wl * f(ip1) + (5.0d+00 * (wl + wc) + wr) * f(i ) & + + (wc + 5.0d+00 * wr) * f(im1)) * dq end do ! i = 1, n @@ -2269,6 +2273,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -2308,6 +2313,7 @@ module interpolations ! prepare coefficients of the tridiagonal system ! a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 c(i,1) = 0.0d+00 r(i,1) = fl(i) @@ -2326,9 +2332,9 @@ module interpolations ! normalize weights ! - wl = dr * al(i,2) + wl = dl * ar(i,2) wc = dc * ac(i,2) - wr = dl * ar(i,2) + wr = dr * al(i,2) ww = (wl + wr) + wc wl = wl / ww wr = wr / ww @@ -2336,17 +2342,18 @@ module interpolations ! calculate the interpolations of the right state ! - ql = a31 * f(i ) + a32 * f(im1) + a33 * f(im2) + ql = a11 * f(ip2) + a12 * f(ip1) + a13 * f(i ) qc = a21 * f(ip1) + a22 * f(i ) + a23 * f(im1) - qr = a11 * f(ip2) + a12 * f(ip1) + a13 * f(i ) + qr = a31 * f(i ) + a32 * f(im1) + a33 * f(im2) ! calculate the right state ! - fr(i) = (wr * qr + wl * ql) + wc * qc + fr(i) = (wl * ql + wr * qr) + wc * qc ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -2365,9 +2372,9 @@ module interpolations ! normalize weights ! - wl = dr * al(i,2) + wl = dl * ar(i,2) wc = dc * ac(i,2) - wr = dl * ar(i,2) + wr = dr * al(i,2) ww = (wl + wr) + wc wl = wl / ww wr = wr / ww @@ -2375,17 +2382,18 @@ module interpolations ! calculate the interpolations of the right state ! - ql = a31 * f(i ) + a32 * f(im1) + a33 * f(im2) + ql = a11 * f(ip2) + a12 * f(ip1) + a13 * f(i ) qc = a21 * f(ip1) + a22 * f(i ) + a23 * f(im1) - qr = a11 * f(ip2) + a12 * f(ip1) + a13 * f(i ) + qr = a31 * f(i ) + a32 * f(im1) + a33 * f(im2) ! calculate the right state ! - fr(i) = (wr * qr + wl * ql) + wc * qc + fr(i) = (wl * ql + wr * qr) + wc * qc ! prepare coefficients of the tridiagonal system ! a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 c(i,2) = 0.0d+00 r(i,2) = fr(i) @@ -2393,37 +2401,19 @@ module interpolations ! solve the tridiagonal system of equations for the left-side interpolation ! - ib = 1 - ie = n - tt = 1.0d+00 - u(ib) = r(ib,1) - do i = ib + 1, ie - im1 = i - 1 - g(i) = c(im1,1) / tt - tt = 1.0d+00 - a(i,1) * g(i) - u(i) = (r(i,1) - a(i,1) * u(im1)) / tt - end do - do i = ie - 1, ib, -1 - ip1 = i + 1 - u(i) = u(i) - g(ip1) * u(ip1) - end do - fl(ib:ie) = u(ib:ie) + call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n)) -! solve the tridiagonal system of equations for the right-side interpolation +! substitute the left-side values ! - tt = 1.0d+00 - u(ie) = r(ie,2) - do i = ie - 1, ib, -1 - ip1 = i + 1 - g(i) = a(ip1,2) / tt - tt = 1.0d+00 - c(i,2) * g(i) - u(i) = (r(i,2) - c(i,2) * u(ip1)) / tt - end do - do i = ib + 1, ie - im1 = i - 1 - u(i) = u(i) - g(im1) * u(im1) - end do - fr(ib:ie-1) = u(ib+1:ie) + fl(1:n ) = u(1:n) + +! 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)) + +! substitute the right-side values +! + fr(1:n-1) = u(2:n) ! update the interpolation of the first and last points ! From 224cc37b531ec213e36a2bcebff6dddf5aa01b17 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 8 Dec 2015 09:54:35 -0200 Subject: [PATCH 6/6] INTERPOLATIONS: Remove unused variables. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 4cac6f3..d74511b 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -1332,7 +1332,7 @@ module interpolations ! local variables ! - integer :: i, im1, ip1, im2, ip2, ib, ie + integer :: i, im1, ip1, im2, ip2 real(kind=8) :: bl, bc, br, tt real(kind=8) :: wl, wc, wr, ww real(kind=8) :: ql, qc, qr @@ -1705,7 +1705,7 @@ module interpolations ! local variables ! - integer :: i, im1, ip1, im2, ip2, ib, ie + integer :: i, im1, ip1, im2, ip2 real(kind=8) :: bl, bc, br, tt real(kind=8) :: wl, wc, wr, ww real(kind=8) :: ql, qc, qr