Merge branch 'master' into reconnection
This commit is contained in:
commit
adc354f448
@ -1165,8 +1165,8 @@ module equations
|
||||
!
|
||||
write(sid,'(i15)') id
|
||||
write(snc,'(i15)') nc
|
||||
write(msg,'("Unphysical cells in block ID:",a," (",a," counted)")') &
|
||||
trim(sid), trim(snc)
|
||||
write(msg,'("Unphysical cells in block ID:",a," (",a," counted).")') &
|
||||
trim(adjustl(sid)), trim(adjustl(snc))
|
||||
call print_warning(loc, trim(msg))
|
||||
|
||||
! allocate temporary vectors for unphysical states
|
||||
@ -3522,7 +3522,7 @@ module equations
|
||||
en = u(ien,i) + u(idn,i)
|
||||
dn = u(idn,i)
|
||||
|
||||
! find the exact W using an Newton-Ralphson interative method
|
||||
! find the exact W using the Newton-Ralphson interative method
|
||||
!
|
||||
call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
|
||||
|
||||
@ -3873,124 +3873,135 @@ module equations
|
||||
call start_timer(imp)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! check if the state is physical; this can save some time if unphysical state
|
||||
! is considered
|
||||
!
|
||||
wl = sqrt(mm + dn * dn)
|
||||
if (en > wl) then
|
||||
|
||||
! prepare the initial brackets
|
||||
!
|
||||
wl = sqrt(mm + dn * dn) + gammaxi * pmin
|
||||
wu = en + pmin
|
||||
wl = wl + gammaxi * pmin
|
||||
wu = en + pmin
|
||||
|
||||
! calculate the value of function for the lower bracket
|
||||
!
|
||||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||||
|
||||
! the lower bracket gives negative function, so there is chance it bounds
|
||||
! the root
|
||||
!
|
||||
if (fl < 0.0d+00) then
|
||||
if (fl < 0.0d+00) then
|
||||
|
||||
! make sure that the upper bracket is larger than the lower one and
|
||||
! the function has positive value
|
||||
!
|
||||
if (wu <= wl) wu = 2.0d+00 * wl
|
||||
if (wu <= wl) wu = 2.0d+00 * wl
|
||||
|
||||
! check if the brackets bound the root region, if not proceed until
|
||||
! opposite function signs are found for the brackets
|
||||
!
|
||||
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 = 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
|
||||
|
||||
! the upper bracket was found, so proceed with determining the root
|
||||
!
|
||||
if (it > 0 .and. fu >= 0.0d+00) then
|
||||
if (it > 0 .and. fu >= 0.0d+00) then
|
||||
|
||||
! 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
|
||||
!
|
||||
info = .true.
|
||||
keep = .true.
|
||||
it = nrmax
|
||||
cn = nrext
|
||||
info = .true.
|
||||
keep = .true.
|
||||
it = nrmax
|
||||
cn = nrext
|
||||
|
||||
! iterate using the Newton-Raphson method in order to find a root w of the
|
||||
! function
|
||||
!
|
||||
do while(keep)
|
||||
do while(keep)
|
||||
|
||||
! 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
|
||||
!
|
||||
if (f > fl .and. f < 0.0d+00) then
|
||||
wl = w
|
||||
fl = f
|
||||
end if
|
||||
if (f < fu .and. f > 0.0d+00) then
|
||||
wu = w
|
||||
fu = f
|
||||
end if
|
||||
if (f > fl .and. f < 0.0d+00) then
|
||||
wl = w
|
||||
fl = f
|
||||
end if
|
||||
if (f < fu .and. f > 0.0d+00) then
|
||||
wu = w
|
||||
fu = f
|
||||
end if
|
||||
|
||||
! calculate the increment dW, update the solution, and estimate the error
|
||||
!
|
||||
dw = f / df
|
||||
w = w - dw
|
||||
err = abs(dw / w)
|
||||
dw = f / df
|
||||
w = w - dw
|
||||
err = abs(dw / w)
|
||||
|
||||
! check the convergence, if the convergence is not reached, iterate until
|
||||
! the maximum number of iteration is reached
|
||||
!
|
||||
if (err < tol) then
|
||||
keep = cn > 0
|
||||
cn = cn - 1
|
||||
else
|
||||
keep = it > 0
|
||||
end if
|
||||
if (err < tol) then
|
||||
keep = cn > 0
|
||||
cn = cn - 1
|
||||
else
|
||||
keep = it > 0
|
||||
end if
|
||||
|
||||
! if new W leaves the brackets, use the bisection method to estimate
|
||||
! the new guess
|
||||
!
|
||||
if (w < wl .or. w > wu) then
|
||||
w = 0.5d+00 * (wl + wu)
|
||||
end if
|
||||
if (w < wl .or. w > wu) then
|
||||
w = 0.5d+00 * (wl + wu)
|
||||
end if
|
||||
|
||||
it = it - 1
|
||||
end do ! NR iterations
|
||||
it = it - 1
|
||||
end do ! NR iterations
|
||||
|
||||
! print information about failed convergence or unphysical variables
|
||||
!
|
||||
if (err >= tol) then
|
||||
call print_warning(loc, "Convergence not reached!")
|
||||
write(*,"(a,1x,1e24.16e3)") "Error: ", err
|
||||
end if
|
||||
if (err >= tol) then
|
||||
call print_warning(loc, "Convergence not reached!")
|
||||
write(*,"(a,1x,1e24.16e3)") "Error: ", err
|
||||
end if
|
||||
|
||||
! calculate |V|² from W
|
||||
!
|
||||
vv = mm / (w * w)
|
||||
vv = mm / (w * w)
|
||||
|
||||
else ! the upper brack not found
|
||||
call print_warning(loc, "Could not find the upper bracket!")
|
||||
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
|
||||
|
||||
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!")
|
||||
else ! the state is unphysical
|
||||
call print_warning(loc, "The state is not physical!")
|
||||
info = .false.
|
||||
end if
|
||||
|
||||
@ -5363,7 +5374,8 @@ module equations
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use algebra , only : cubic_normalized, quartic
|
||||
use algebra, only : cubic_normalized, quartic
|
||||
use error , only : print_warning
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -5387,6 +5399,11 @@ module equations
|
||||
!
|
||||
real(kind=8), dimension(5) :: a
|
||||
real(kind=8), dimension(4) :: x
|
||||
|
||||
! local parameters
|
||||
!
|
||||
character(len=*), parameter :: loc = &
|
||||
'EQUATIONS::nr_initial_brackets_srmhd_adi()'
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -5450,29 +5467,19 @@ module equations
|
||||
keep = .true.
|
||||
it = nrmax
|
||||
wl = wu
|
||||
|
||||
do while(keep)
|
||||
|
||||
call nr_pressure_srmhd_adi_1d(mm, bb, mb, ec, dn, wl, f, df)
|
||||
|
||||
dw = f / df
|
||||
wl = wl - dw
|
||||
|
||||
err = abs(dw / wl)
|
||||
it = it - 1
|
||||
keep = (err > tol) .and. it > 0
|
||||
|
||||
end do
|
||||
|
||||
if (it <= 0) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)") "ERROR in" &
|
||||
, "EQUATIONS::nr_initial_brackets_srmhd_adi()"
|
||||
write(*,"(a)" ) "Could not find the lower limit for the enthalpy!"
|
||||
call print_warning(loc, &
|
||||
"Iterative solver failed to find the lower bracket!")
|
||||
write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = " &
|
||||
, dn, mm, mb, bb, en
|
||||
info = .false.
|
||||
return
|
||||
end if
|
||||
|
||||
end if ! nr > 0
|
||||
@ -5481,71 +5488,53 @@ module equations
|
||||
!
|
||||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wl, fl)
|
||||
|
||||
if (fl > 0.0d+00) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)") "ERROR in" &
|
||||
, "EQUATIONS::nr_initial_brackets_srmhd_adi()"
|
||||
write(*,"(a)" ) "Lower limit positive!"
|
||||
write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = " &
|
||||
, dn, mm, mb, bb, en, wl
|
||||
info = .false.
|
||||
return
|
||||
end if
|
||||
if (fl <= 0.0d+00) then
|
||||
|
||||
! make sure that the upper limit is larger than the lower one
|
||||
!
|
||||
keep = wl >= wu
|
||||
it = nrmax
|
||||
do while(keep)
|
||||
wu = 2.0d+00 * wu
|
||||
it = it - 1
|
||||
keep = (wl >= wu) .and. it > 0
|
||||
end do
|
||||
if (it <= 0) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)") "ERROR in" &
|
||||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||||
write(*,"(a)" ) "Could not find the upper limit for enthalpy!"
|
||||
write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = " &
|
||||
, dn, mm, mb, bb, en, wl
|
||||
info = .false.
|
||||
return
|
||||
end if
|
||||
if (wu <= wl) wu = 2.0d+00 * wl
|
||||
|
||||
! check if the brackets bound the root region, if not proceed until
|
||||
! opposite function signs are found for the brackets
|
||||
!
|
||||
call nr_function_srmhd_adi_1d(mm, bb, mb, 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_srmhd_adi_1d(mm, bb, mb, 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_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
|
||||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||||
end do
|
||||
|
||||
it = it - 1
|
||||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||||
|
||||
end do
|
||||
if (it <= 0) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)") "ERROR in" &
|
||||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||||
write(*,"(a)" ) "No initial brackets found!"
|
||||
write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = " &
|
||||
, dn, mm, mb, bb, en
|
||||
info = .false.
|
||||
return
|
||||
end if
|
||||
! the upper bracket was found, so proceed with determining the root
|
||||
!
|
||||
if (it > 0 .and. fu >= 0.0d+00) then
|
||||
|
||||
! estimate the enthalpy value close to the root
|
||||
!
|
||||
wc = wl - fl * (wu - wl) / (fu - fl)
|
||||
wc = wl - fl * (wu - wl) / (fu - fl)
|
||||
|
||||
! we have good brackets and guess, so good to go
|
||||
!
|
||||
info = .true.
|
||||
|
||||
else ! the upper brack not found
|
||||
call print_warning(loc, "Could not find the upper bracket!")
|
||||
write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = " &
|
||||
, dn, mm, mb, bb, en
|
||||
info = .false.
|
||||
|
||||
end if
|
||||
|
||||
else ! the root cannot be found, since it is below the lower bracket
|
||||
call print_warning(loc, "Positive function for lower bracket!")
|
||||
write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = " &
|
||||
, dn, mm, mb, bb, en, wl
|
||||
info = .false.
|
||||
end if
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the initial brackets
|
||||
@ -5587,6 +5576,10 @@ module equations
|
||||
!
|
||||
subroutine nr_iterate_srmhd_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
|
||||
!
|
||||
implicit none
|
||||
@ -5604,6 +5597,10 @@ module equations
|
||||
real(kind=8) :: wl, wu, fl, fu
|
||||
real(kind=8) :: f , df, dw
|
||||
real(kind=8) :: err
|
||||
|
||||
! local parameters
|
||||
!
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srmhd_adi_1dw()'
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -5616,104 +5613,80 @@ module equations
|
||||
! find the initial brackets and estimate the initial enthalpy
|
||||
!
|
||||
call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn &
|
||||
, wl, wu, w, fl, fu, info)
|
||||
, wl, wu, w, fl, fu, info)
|
||||
|
||||
! if the brackets could not be found, return the lower bracket as the solution
|
||||
! continue if brackets found
|
||||
!
|
||||
if (.not. info) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||||
write(*,"(a,1x)" ) "The solution lays in unphysical regime."
|
||||
write(*,"(a,1x,1e24.16e3)") "Using the lower bracket as solution: ", wl
|
||||
if (info) then
|
||||
|
||||
! use the lower bracket, since it guarantees the positive pressure
|
||||
! initialize iteration parameters
|
||||
!
|
||||
w = wl
|
||||
info = .true.
|
||||
keep = .true.
|
||||
it = nrmax
|
||||
cn = nrext
|
||||
|
||||
! iterate using the Newton-Raphson method in order to find a root w of the
|
||||
! function
|
||||
!
|
||||
do while(keep)
|
||||
|
||||
! calculate F(W) and dF(W)/dW
|
||||
!
|
||||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
|
||||
|
||||
! update brackets
|
||||
!
|
||||
if (f > fl .and. f < 0.0d+00) then
|
||||
wl = w
|
||||
fl = f
|
||||
end if
|
||||
if (f < fu .and. f > 0.0d+00) then
|
||||
wu = w
|
||||
fu = f
|
||||
end if
|
||||
|
||||
! calculate the increment dW, update the solution, and estimate the error
|
||||
!
|
||||
dw = f / df
|
||||
w = w - dw
|
||||
err = abs(dw / w)
|
||||
|
||||
! check the convergence, if the convergence is not reached, iterate until
|
||||
! the maximum number of iteration is reached
|
||||
!
|
||||
if (err < tol) then
|
||||
keep = cn > 0
|
||||
cn = cn - 1
|
||||
else
|
||||
keep = it > 0
|
||||
end if
|
||||
|
||||
! if new W lays out of the brackets, use the bisection method to estimate
|
||||
! the new guess
|
||||
!
|
||||
if (w < wl .or. w > wu) then
|
||||
w = 0.5d+00 * (wl + wu)
|
||||
end if
|
||||
|
||||
! decrease the number of remaining iterations
|
||||
!
|
||||
it = it - 1
|
||||
|
||||
end do ! NR iterations
|
||||
|
||||
! let know the user if the convergence failed
|
||||
!
|
||||
if (err >= tol) then
|
||||
call print_warning(loc, "Convergence not reached!")
|
||||
write(*,"(a,1x,1e24.16e3)") "Error: ", err
|
||||
end if
|
||||
|
||||
! calculate |V|² from W
|
||||
!
|
||||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||||
|
||||
info = .true.
|
||||
return
|
||||
|
||||
end if
|
||||
|
||||
! initialize iteration parameters
|
||||
!
|
||||
info = .true.
|
||||
keep = .true.
|
||||
it = nrmax
|
||||
cn = nrext
|
||||
|
||||
! iterate using the Newton-Raphson method in order to find a root w of the
|
||||
! function
|
||||
!
|
||||
do while(keep)
|
||||
|
||||
! calculate F(W) and dF(W)/dW
|
||||
!
|
||||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
|
||||
|
||||
! update brackets
|
||||
!
|
||||
if (f > fl .and. f < 0.0d+00) then
|
||||
wl = w
|
||||
fl = f
|
||||
end if
|
||||
if (f < fu .and. f > 0.0d+00) then
|
||||
wu = w
|
||||
fu = f
|
||||
end if
|
||||
|
||||
! calculate the increment dW
|
||||
!
|
||||
dw = f / df
|
||||
|
||||
! update the solution
|
||||
!
|
||||
w = w - dw
|
||||
|
||||
! calculate the error
|
||||
!
|
||||
err = abs(dw / w)
|
||||
|
||||
! check the convergence, if the convergence is not reached, iterate until
|
||||
! the maximum number of iteration is reached
|
||||
!
|
||||
if (err < tol) then
|
||||
keep = cn > 0
|
||||
cn = cn - 1
|
||||
else
|
||||
keep = it > 0
|
||||
end if
|
||||
|
||||
! if new W lays out of the brackets, use the bisection method to estimate
|
||||
! the new guess
|
||||
!
|
||||
if (w < wl .or. w > wu) then
|
||||
w = 0.5d+00 * (wl + wu)
|
||||
end if
|
||||
|
||||
! decrease the number of remaining iterations
|
||||
!
|
||||
it = it - 1
|
||||
|
||||
end do ! NR iterations
|
||||
|
||||
! calculate |V|² from W
|
||||
!
|
||||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||||
|
||||
! let know the user if the convergence failed
|
||||
!
|
||||
if (err >= tol) then
|
||||
write(*,*)
|
||||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||||
write(*,"(a,1x,1e24.16e3)") "Convergence not reached: ", err
|
||||
end if
|
||||
end if ! correct brackets
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for variable solver
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user