Skip to content

Add t17d and t20d diagnostics #1001

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
adcroft wants to merge 1 commit into
base: dev/gfdl
Choose a base branch
from
Open

Add t17d and t20d diagnostics #1001

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions src/ALE/Recon1d_PCM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ module Recon1d_PCM
!! average() *locally defined
!! f() *locally defined
!! dfdx() *locally defined
!! - x() *locally defined
!! check_reconstruction() *locally defined
!! unit_tests() *locally defined
!! destroy() *locally defined
Expand All @@ -36,6 +37,8 @@ module Recon1d_PCM
procedure :: f => f
!> Implementation of the derivative of the PCM reconstruction at a point [A]
procedure :: dfdx => dfdx
!> Implementation of solver for x: f(x)=t
procedure :: x => x
!> Implementation of deallocation for PCM
procedure :: destroy => destroy
!> Implementation of check reconstruction for the PCM reconstruction
Expand Down Expand Up @@ -105,6 +108,24 @@ real function dfdx(this, k, x)

end function dfdx

!> Solver for x: f(x)=t
real function x(this, k, t)
class(PCM), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(in) :: t !< Value to solve for [A]
real :: slp ! Difference across cell [A]

slp = this%u_mean(min(k+1,this%n)) - this%u_mean(max(k-1,1))
if ( abs(slp) > 0. ) slp = sign(1., slp)
x = 0.5 ! Fall back if t==u_mean
! if t>u_mean & slp=1 then x=1
! if t<u_mean & slp=1 then x=0
! if t>u_mean & slp=-1 then x=0
! if t<u_mean & slp=-1 then x=1
! if slp=0 then x=0.5
if ( abs(t - this%u_mean(k)) > 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
end function x

!> Average between xa and xb for cell k of a 1D PCM reconstruction [A]
real function average(this, k, xa, xb)
class(PCM), intent(in) :: this !< This reconstruction
Expand Down Expand Up @@ -181,6 +202,16 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(3, um, (/0.,0.,0./), 'dfdx in center')
call test%real_arr(3, ur, (/0.,0.,0./), 'dfdx on right edge')

call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0')
call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5')
call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1')
call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5')
call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1')
call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0')
call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5')
call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1')

do k = 1, 3
um(k) = this%average(k, 0.5, 0.75)
enddo
Expand Down
38 changes: 38 additions & 0 deletions src/ALE/Recon1d_PLM_CW.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ module Recon1d_PLM_CW
!! - average() *locally defined
!! - f() *locally defined
!! - dfdx() *locally defined
!! - x() *locally defined
!! - check_reconstruction() *locally defined
!! - unit_tests() *locally defined
!! - destroy() *locally defined
Expand All @@ -45,6 +46,8 @@ module Recon1d_PLM_CW
procedure :: f => f
!> Implementation of the derivative of the PLM_CW reconstruction at a point [A]
procedure :: dfdx => dfdx
!> Implementation of solver for x: f(x)=t
procedure :: x => x
!> Implementation of deallocation for PLM_CW
procedure :: destroy => destroy
!> Implementation of check reconstruction for the PLM_CW reconstruction
Expand Down Expand Up @@ -194,6 +197,31 @@ real function dfdx(this, k, x)

end function dfdx

!> Solver for x such that f(x)=t
real function x(this, k, t)
class(PLM_CW), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(in) :: t !< Value to solve for [A]
real :: slp ! Difference across cell [A]

slp = this%ur(k) - this%ul(k)
if ( abs(slp) > 0. ) then
x = ( t - this%ul(k) ) / slp
x = max( 0., min( x, 1. ) )
else
slp = this%ul(min(k+1,this%n)) - this%ur(max(k-1,1))
if ( abs(slp) > 0. ) slp = sign(1., slp)
x = 0.5 ! Fall back if t==u_mean
! if t>u_mean & slp=1 then x=1
! if t<u_mean & slp=1 then x=0
! if t>u_mean & slp=-1 then x=0
! if t<u_mean & slp=-1 then x=1
! if t=u_mean then slp=0
! if slp=0 then x=0.5
if ( abs(t - this%u_mean(k)) > 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
endif
end function x

!> Average between xa and xb for cell k of a 1D PLM reconstruction [A]
real function average(this, k, xa, xb)
class(PLM_CW), intent(in) :: this !< This reconstruction
Expand Down Expand Up @@ -334,6 +362,16 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center')
call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge')

call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0')
call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5')
call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1')
call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5')
call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1')
call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0')
call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5')
call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1')

do k = 1, 3
um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
enddo
Expand Down
1 change: 1 addition & 0 deletions src/ALE/Recon1d_PLM_CWK.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ module Recon1d_PLM_CWK
!! - average() -> recon1d_plm_cw.average()
!! - f() -> recon1d_plm_cw.f()
!! - dfdx() -> recon1d_plm_cw.dfdx()
!! - x() -> recon1d_plm_cw.x()
!! - check_reconstruction() -> recon1d_plm_cw.check_reconstruction()
!! - unit_tests() -> recon1d_plm_cw.unit_tests()
!! - destroy() -> recon1d_plm_cw.destroy()
Expand Down
11 changes: 11 additions & 0 deletions src/ALE/Recon1d_PLM_WLS.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ module Recon1d_PLM_WLS
!! - average() *locally defined
!! - f() *locally defined
!! - dfdx() *locally defined
!! - x() -> recon1d_type.x()
!! - check_reconstruction() *locally defined
!! - unit_tests() *locally defined
!! - destroy() *locally defined
Expand Down Expand Up @@ -371,6 +372,16 @@ logical function unit_tests(this, verbose, stdout, stderr)
enddo
call test%real_arr(3, um, (/1.25,3.25,5.25/), "Return interval average")

call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0')
call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5')
call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1')
call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5')
call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1')
call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0')
call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5')
call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1')

call this%destroy()
deallocate( um, ul, ur, ull, urr )

Expand Down
11 changes: 11 additions & 0 deletions src/ALE/Recon1d_PLM_hybgen.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ module Recon1d_PLM_hybgen
!! - average() *locally defined
!! - f() *locally defined
!! - dfdx() *locally defined
!! - x() -> recon1d_plm_cw.x()
!! - check_reconstruction() *locally defined
!! - unit_tests() *locally defined
!! - destroy() *locally defined
Expand Down Expand Up @@ -358,6 +359,16 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center')
call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge')

call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0')
call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5')
call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1')
call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5')
call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1')
call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0')
call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5')
call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1')

do k = 1, 3
um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
enddo
Expand Down
14 changes: 13 additions & 1 deletion src/ALE/Recon1d_PPM_CW.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module Recon1d_PPM_CW
!! - average() *locally defined
!! - f() *locally defined
!! - dfdx() *locally defined
!! - x() *locally defined
!! - check_reconstruction() *locally defined
!! - unit_tests() *locally defined
!! - destroy() *locally defined
Expand All @@ -49,6 +50,8 @@ module Recon1d_PPM_CW
procedure :: f => f
!> Implementation of the derivative of the PPM_CW reconstruction at a point [A]
procedure :: dfdx => dfdx
!> Implementation of solver for x: f(x)=t
! procedure :: x => x
!> Implementation of deallocation for PPM_CW
procedure :: destroy => destroy
!> Implementation of check reconstruction for the PPM_CW reconstruction
Expand Down Expand Up @@ -152,7 +155,7 @@ subroutine reconstruct(this, h, u)
this%ur(n) = u(n) ! PCM
this%ul(n) = u(n) ! PCM

do K = 2, n ! K=2 is interface between cells 1 and 2
do K = 2, n-1 ! K=2 is interface between cells 1 and 2
u0 = u(k-1)
u1 = u(k)
u2 = u(k+1)
Expand Down Expand Up @@ -333,6 +336,7 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%set( stdout=stdout ) ! Sets the stdout channel in test
call test%set( stderr=stderr ) ! Sets the stderr channel in test
call test%set( verbose=verbose ) ! Sets the verbosity flag in test
call test%set( stop_instantly=.true. )

if (verbose) write(stdout,'(a)') 'PPM_CW:unit_tests testing with linear fn'

Expand Down Expand Up @@ -366,6 +370,10 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(5, um, (/0.,3.,3.,3.,0./), 'dfdx in center')
call test%real_arr(5, ur, (/0.,3.,3.,3.,0./), 'dfdx on right edge')

call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
call test%real_scalar( this%x(2,4.), 0.5, 'f-1(2,4)=0.5')
call test%real_scalar( this%x(2,5.5), 1., 'f-1(2,5.5)=1')

do k = 1, 5
um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
enddo
Expand All @@ -392,6 +400,10 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center')
call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge')

call test%real_scalar( this%x(3,12.), 0., 'f-1(3,12)=0')
call test%real_scalar( this%x(3,18.75), 0.5, 'f-1(3,18.75)=0.5', robits=1)
call test%real_scalar( this%x(3,27.), 1., 'f-1(3,27)=1')

! x = 3 i i=0 at origin
! f(x) = x^2 / 3 = 3 i^2
! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
Expand Down
51 changes: 50 additions & 1 deletion src/ALE/Recon1d_PPM_CWK.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ module Recon1d_PPM_CWK
!! - average() *locally defined
!! - f() *locally defined
!! - dfdx() *locally defined
!! - x() *locally defined
!! - check_reconstruction() *locally defined
!! - unit_tests() *locally defined
!! - destroy() *locally defined
Expand All @@ -50,6 +51,8 @@ module Recon1d_PPM_CWK
procedure :: f => f
!> Implementation of the derivative of the PPM_CWK reconstruction at a point [A]
procedure :: dfdx => dfdx
!> Implementation of solver for x: f(x)=t
procedure :: x => x
!> Implementation of deallocation for PPM_CWK
procedure :: destroy => destroy
!> Implementation of check reconstruction for the PPM_CWK reconstruction
Expand Down Expand Up @@ -137,7 +140,7 @@ subroutine reconstruct(this, h, u)
this%ur(n) = u(n) ! PCM
this%ul(n) = u(n) ! PCM

do K = 2, n ! K=2 is interface between cells 1 and 2
do K = 2, n-1 ! K=2 is interface between cells 1 and 2
u0 = u(k-1)
u1 = u(k)
u2 = u(k+1)
Expand Down Expand Up @@ -215,6 +218,48 @@ real function dfdx(this, k, x)

end function dfdx

!> Solver for x: f(x)=t
real function x(this, k, t)
class(PPM_CWK), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(in) :: t !< Value to solve for [A]
real :: xl, xr, xo ! Left/right bounds and guess [nondim]
real :: fl, fr ! Left right values [A]
real :: slp ! Difference across cell or derivative wrt nondim x [A]
real :: f_at_x ! Value at current x [A]
integer :: iter

slp = this%ur(k) - this%ul(k)
if ( abs(slp) > 0. ) then
xl = 0.
xr = 1.
xo = ( t - this%ul(k) ) / slp ! First guess by regula falsi
f_at_x = this%f(k, xo)
do iter = 1,5
Copy link
Member

@Hallberg-NOAA Hallberg-NOAA Dec 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a comment explaining why this code is using Newton's method to solve for the root of a quadratic equation rather than just using the closed-form solution. (Could it be that the profile is quadratic in T and S, but with a nonlinear equation of state the function might be nonlinear for density, for example?)

slp = this%dfdx(k, xo)
x = xo - ( f_at_x - t ) / slp ! Newton-Raphson step
if ( x < xl ) x = 0.5 * ( xl + xo ) ! Replace with bi-section
if ( x > xr ) x = 0.5 * ( xr + xo ) ! Replace with bi-section
f_at_x = this%f(k, x)
if ( abs(f_at_x - t) <= 0. .or. abs(x - xo) < this%x_tolerance ) return
if ( f_at_x < t ) xl = x ! Replace left bound
if ( f_at_x > t ) xr = x ! Replace right bound
xo = x
enddo
else
x = 0.5 ! Fall back if t==u_mean
slp = this%ul(min(k+1,this%n)) - this%ur(max(k-1,1))
if ( abs(slp) > 0. ) slp = sign(1., slp)
! if t>u_mean & slp=1 then x=1
! if t<u_mean & slp=1 then x=0
! if t>u_mean & slp=-1 then x=0
! if t<u_mean & slp=-1 then x=1
! if t=u_mean then slp=0
! if slp=0 then x=0.5
x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
endif
end function x

!> Average between xa and xb for cell k of a 1D PPM reconstruction [A]
real function average(this, k, xa, xb)
class(PPM_CWK), intent(in) :: this !< This reconstruction
Expand Down Expand Up @@ -373,6 +418,10 @@ logical function unit_tests(this, verbose, stdout, stderr)
call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center')
call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge')

call test%real_scalar( this%x(3,12.), 0., 'f-1(3,12)=0')
call test%real_scalar( this%x(3,18.75), 0.5, 'f-1(3,18.75)=0.5', robits=1)
call test%real_scalar( this%x(3,27.), 1., 'f-1(3,27)=1')

! x = 3 i i=0 at origin
! f(x) = x^2 / 3 = 3 i^2
! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
Expand Down
2 changes: 1 addition & 1 deletion src/ALE/Recon1d_PPM_hybgen.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ subroutine reconstruct(this, h, u)
this%ur(n) = u(n) ! PCM
this%ul(n) = u(n) ! PCM

do K = 2, n ! K=2 is interface between cells 1 and 2
do K = 2, n-1 ! K=2 is interface between cells 1 and 2
u0 = u(k-1)
u1 = u(k)
u2 = u(k+1)
Expand Down
Loading