EQUATIONS: Rewrite nr_iterate_srhd_adi_1dw().
Make sure all cases are considered. Also, in the case of profiling, the subroutine counts timing correctly. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
c7a4159af7
commit
fe117327d1
@ -3835,6 +3835,10 @@ module equations
|
|||||||
!
|
!
|
||||||
subroutine nr_iterate_srhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)
|
subroutine nr_iterate_srhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)
|
||||||
|
|
||||||
|
! include external procedures
|
||||||
|
!
|
||||||
|
use error, only : print_warning
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
implicit none
|
implicit none
|
||||||
@ -3852,6 +3856,10 @@ module equations
|
|||||||
real(kind=8) :: wl, wu, fl, fu
|
real(kind=8) :: wl, wu, fl, fu
|
||||||
real(kind=8) :: f, df, dw
|
real(kind=8) :: f, df, dw
|
||||||
real(kind=8) :: err
|
real(kind=8) :: err
|
||||||
|
|
||||||
|
! local parameters
|
||||||
|
!
|
||||||
|
character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srhd_adi_1dw()'
|
||||||
!
|
!
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
@ -3866,132 +3874,120 @@ module equations
|
|||||||
wl = sqrt(mm + dn * dn) + gammaxi * pmin
|
wl = sqrt(mm + dn * dn) + gammaxi * pmin
|
||||||
wu = en + pmin
|
wu = en + pmin
|
||||||
|
|
||||||
! make sure that the upper bracket is larger than the lower one
|
! calculate the value of function for the lower bracket
|
||||||
!
|
!
|
||||||
keep = wl >= wu
|
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||||||
it = nrmax
|
|
||||||
do while(keep)
|
! the lower bracket gives negative function, so there is chance it bounds
|
||||||
wu = 2.0d+00 * wu
|
! the root
|
||||||
it = it - 1
|
!
|
||||||
keep = (wl >= wu) .and. it > 0
|
if (fl < 0.0d+00) then
|
||||||
end do
|
|
||||||
if (it <= 0) then
|
! make sure that the upper bracket is larger than the lower one and
|
||||||
write(*,*)
|
! the function has positive value
|
||||||
write(*,"(a,1x,a)") "ERROR in" &
|
!
|
||||||
, "EQUATIONS::nr_iterate_srhd_adi_1dw()"
|
if (wu <= wl) wu = 2.0d+00 * wl
|
||||||
write(*,"(a)" ) "Could not find the upper limit for enthalpy!"
|
|
||||||
info = .false.
|
|
||||||
return
|
|
||||||
end if
|
|
||||||
|
|
||||||
! check if the brackets bound the root region, if not proceed until
|
! check if the brackets bound the root region, if not proceed until
|
||||||
! opposite function signs are found for the brackets
|
! opposite function signs are found for the brackets
|
||||||
!
|
!
|
||||||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
|
||||||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
|
||||||
|
|
||||||
keep = (fl * fu > 0.0d+00)
|
|
||||||
it = nrmax
|
|
||||||
|
|
||||||
do while (keep)
|
|
||||||
|
|
||||||
wl = wu
|
|
||||||
fl = fu
|
|
||||||
wu = 2.0d+00 * wu
|
|
||||||
|
|
||||||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||||||
|
it = nrmax
|
||||||
|
keep = fl * fu > 0.0d+00
|
||||||
|
do while (keep)
|
||||||
|
it = it - 1
|
||||||
|
wl = wu
|
||||||
|
fl = fu
|
||||||
|
wu = 2.0d+00 * wu
|
||||||
|
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||||||
|
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||||||
|
end do
|
||||||
|
|
||||||
it = it - 1
|
! the upper bracket was found, so proceed with determining the root
|
||||||
|
!
|
||||||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
if (it > 0 .and. fu >= 0.0d+00) then
|
||||||
|
|
||||||
end do
|
|
||||||
if (it <= 0) then
|
|
||||||
write(*,*)
|
|
||||||
write(*,"(a,1x,a)") "ERROR in" &
|
|
||||||
, "EQUATIONS::nr_iterate_srhd_adi_1dw()"
|
|
||||||
write(*,"(a)" ) "No initial brackets found!"
|
|
||||||
info = .false.
|
|
||||||
return
|
|
||||||
end if
|
|
||||||
|
|
||||||
! estimate the value of enthalpy close to the root and corresponding v²
|
! estimate the value of enthalpy close to the root and corresponding v²
|
||||||
!
|
!
|
||||||
w = wl - fl * (wu - wl) / (fu - fl)
|
w = wl - fl * (wu - wl) / (fu - fl)
|
||||||
|
|
||||||
! initialize iteration parameters
|
! initialize iteration parameters
|
||||||
!
|
!
|
||||||
info = .true.
|
info = .true.
|
||||||
keep = .true.
|
keep = .true.
|
||||||
it = nrmax
|
it = nrmax
|
||||||
cn = nrext
|
cn = nrext
|
||||||
|
|
||||||
! iterate using the Newton-Raphson method in order to find a root w of the
|
! iterate using the Newton-Raphson method in order to find a root w of the
|
||||||
! function
|
! function
|
||||||
!
|
!
|
||||||
do while(keep)
|
do while(keep)
|
||||||
|
|
||||||
! calculate F(W) and dF(W)/dW
|
! calculate F(W) and dF(W)/dW
|
||||||
!
|
!
|
||||||
call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
|
call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
|
||||||
|
|
||||||
! update brackets
|
! update brackets
|
||||||
!
|
!
|
||||||
if (f > fl .and. f < 0.0d+00) then
|
if (f > fl .and. f < 0.0d+00) then
|
||||||
wl = w
|
wl = w
|
||||||
fl = f
|
fl = f
|
||||||
end if
|
end if
|
||||||
if (f < fu .and. f > 0.0d+00) then
|
if (f < fu .and. f > 0.0d+00) then
|
||||||
wu = w
|
wu = w
|
||||||
fu = f
|
fu = f
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! calculate the increment dW
|
! calculate the increment dW, update the solution, and estimate the error
|
||||||
!
|
!
|
||||||
dw = f / df
|
dw = f / df
|
||||||
|
w = w - dw
|
||||||
! update the solution
|
err = abs(dw / w)
|
||||||
!
|
|
||||||
w = w - dw
|
|
||||||
|
|
||||||
! calculate the error
|
|
||||||
!
|
|
||||||
err = abs(dw / w)
|
|
||||||
|
|
||||||
! check the convergence, if the convergence is not reached, iterate until
|
! check the convergence, if the convergence is not reached, iterate until
|
||||||
! the maximum number of iteration is reached
|
! the maximum number of iteration is reached
|
||||||
!
|
!
|
||||||
if (err < tol) then
|
if (err < tol) then
|
||||||
keep = cn > 0
|
keep = cn > 0
|
||||||
cn = cn - 1
|
cn = cn - 1
|
||||||
else
|
else
|
||||||
keep = it > 0
|
keep = it > 0
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! if new W lays out of the brackets, use the bisection method to estimate
|
! if new W leaves the brackets, use the bisection method to estimate
|
||||||
! the new guess
|
! the new guess
|
||||||
!
|
!
|
||||||
if (w < wl .or. w > wu) then
|
if (w < wl .or. w > wu) then
|
||||||
w = 0.5d+00 * (wl + wu)
|
w = 0.5d+00 * (wl + wu)
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! decrease the number of remaining iterations
|
it = it - 1
|
||||||
!
|
end do ! NR iterations
|
||||||
it = it - 1
|
|
||||||
|
|
||||||
end do ! NR iterations
|
|
||||||
|
|
||||||
! calculate |V|² from W
|
|
||||||
!
|
|
||||||
vv = mm / (w * w)
|
|
||||||
|
|
||||||
! print information about failed convergence or unphysical variables
|
! print information about failed convergence or unphysical variables
|
||||||
!
|
!
|
||||||
if (err >= tol) then
|
if (err >= tol) then
|
||||||
write(*,*)
|
call print_warning(loc, "Convergence not reached!")
|
||||||
write(*,"(a,1x,a)" ) "WARNING in" &
|
write(*,"(a,1x,1e24.16e3)") "Error: ", err
|
||||||
, "EQUATIONS::nr_iterate_srhd_adi_1dw()"
|
end if
|
||||||
write(*,"(a,1x,1e24.16e3)") "Convergence not reached: ", err
|
|
||||||
|
! calculate |V|² from W
|
||||||
|
!
|
||||||
|
vv = mm / (w * w)
|
||||||
|
|
||||||
|
else ! the upper brack not found
|
||||||
|
call print_warning(loc, "Could not find the upper bracket!")
|
||||||
|
info = .false.
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
else if (fl == 0.0d+00) then ! the lower bracket is a root, so return it
|
||||||
|
w = wl
|
||||||
|
info = .true.
|
||||||
|
|
||||||
|
else ! the root cannot be found, since it is below the lower bracket
|
||||||
|
call print_warning(loc, "Positive function for lower bracket!")
|
||||||
|
info = .false.
|
||||||
end if
|
end if
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user