diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 8ede69d5ff..a5f3589bd1 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -18,6 +18,8 @@ module MOM_remapping use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 use MOM_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs +use Recon1d_PCM, only : PCM + implicit none ; private !> Container for remapping parameters @@ -1624,6 +1626,7 @@ logical function remapping_unit_tests(verbose) type(testing) :: test ! Unit testing convenience functions integer :: om4 character(len=4) :: om4_tag + type(PCM) :: PCM call test%set( verbose=verbose ) ! Sets the verbosity flag in test @@ -2257,6 +2260,9 @@ logical function remapping_unit_tests(verbose) 3, (/0.,0.,0./), (/0.,0.,0./), & 3, (/0.,0.,0./), (/0.,0.,0./) ) + if (verbose) write(test%stdout,*) '- - - - - - - - - - Recon1d PCM tests - - - - - - - - -' + call test%test( PCM%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test') + remapping_unit_tests = test%summarize('remapping_unit_tests') end function remapping_unit_tests diff --git a/src/ALE/Recon1d_PCM.F90 b/src/ALE/Recon1d_PCM.F90 new file mode 100644 index 0000000000..21b1cadafc --- /dev/null +++ b/src/ALE/Recon1d_PCM.F90 @@ -0,0 +1,163 @@ +!> The equation of state using the Jackett and McDougall fits to the UNESCO EOS +module Recon1d_PCM + +! This file is part of MOM6. See LICENSE.md for the license. + +use Recon1d_type, only : Recon1d, testing + +implicit none ; private + +public PCM + +!> The EOS_base implementation of the UNESCO equation of state +type, extends (Recon1d) :: PCM + + real, allocatable :: h(:) !< Grid spacing (thickness) [typically H] + +contains + !> Implementation of the PCM constructor + procedure :: construct => construct_PCM + !> Implementation of the PCM reconstruction + procedure :: reconstruct => PCM_reconstruct + !> Implementation of function returning the PCM edge values + procedure :: lr_edge => PCM_lr_edge + !> Implementation of the PCM average over an interval [A] + procedure :: average => PCM_average + !> Implementation of finding the PCM position of a value + procedure :: x_for_f => PCM_x_for_f + !> Implementation of unit tests for the PCM reconstruction + procedure :: unit_tests => PCM_unit_tests + +end type PCM + +contains + +!> Initialize a 1D PCM reconstruction for n cells +subroutine construct_PCM(this, n) + class(PCM), intent(out) :: this !< This reconstruction + integer, intent(in) :: n !< Number of cells in this column + + this%n = n + + allocate( this%u_mean(n) ) + allocate( this%h(n) ) + + if (this%do_legacy) then + this%degree = 1 + allocate( this%poly_coef(n,1) ) + endif + +end subroutine construct_PCM + +!> Calculate a 1D PCM reconstructions based on h(:) and u(:) +subroutine PCM_reconstruct(this, h, u) + class(PCM), intent(inout) :: this !< This reconstruction + real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H] + real, intent(in) :: u(*) !< Cell mean values [A] + ! Local variables + integer :: k + + do k = 1, this%n + this%u_mean(k) = u(k) + this%h(k) = h(k) + enddo + + if (this%do_legacy) then + do k = 1, this%n + this%poly_coef(k,1) = u(k) + enddo + endif + +end subroutine PCM_reconstruct + +!> Returns the left/right edge values in cell k of a 1D PCM reconstruction +subroutine PCM_lr_edge(this, k, u_left, u_right) + class(PCM), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(out) :: u_left !< Left edge value [A] + real, intent(out) :: u_right !< Right edge value [A] + + u_left = this%u_mean(k) + u_right = this%u_mean(k) + +end subroutine PCM_lr_edge + +!> Position at which 1D PCM reconstruction has a value f in cell k [0..1] +real function PCM_x_for_f(this, k, f) + class(PCM), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: f !< Value of reconstruction to solve for [A] + + PCM_x_for_f = 0.5 ! For PCM, every value between x=0 and x=1 has the same value + ! but x=0.5 is the average position + if ( f < this%u_mean(k) ) then + PCM_x_for_f = 0. + elseif ( f > this%u_mean(k) ) then + PCM_x_for_f = 1. + endif + +end function PCM_x_for_f + +!> Average between xa and xb for cell k of a 1D PCM reconstruction [A] +real function PCM_average(this, k, xa, xb) + class(PCM), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: xa !< Start of averageing interval on element (0 to 1) + real, intent(in) :: xb !< End of averageing interval on element (0 to 1) + + PCM_average = this%u_mean(k) + +end function PCM_average + +!> Runs PCM reconstruction unit tests and returns True for any fails, False otherwise +logical function PCM_unit_tests(this, verbose, stdout, stderr) + class(PCM), intent(inout) :: this !< This reconstruction + logical, intent(in) :: verbose !< True, if verbose + integer, intent(in) :: stdout !< I/O channel for stdout + integer, intent(in) :: stderr !< I/O channel for stderr + ! Local variables + real, allocatable :: ul(:), ur(:), um(:) ! test values [A] + type(testing) :: test ! convenience functions + integer :: k + + call test%set( verbose=verbose ) ! Sets the verbosity flag in test + + call this%construct(3) + call test%test( this%n /= 3, 'Setting number of levels') + allocate( um(3), ul(3), ur(3) ) + + call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) ) + call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values') + do k = 1, 3 + um(k) = this%cell_mean(k) + enddo + call test%real_arr(3, um, (/1.,3.,5./), 'Return cell mean') + + do k = 1, 3 + call this%lr_edge(k, ul(k), ur(k)) + enddo + call test%real_arr(3, ul, (/1.,3.,5./), 'Return left edge') + call test%real_arr(3, ur, (/1.,3.,5./), 'Return right edge') + + do k = 1, 3 + um(k) = this%average(k, 0.5, 0.75) + enddo + call test%real_arr(3, um, (/1.,3.,5./), 'Return interval average') + + do k = 1, 3 + ul(k) = this%x_for_f(k, real(2*k-2)) + um(k) = this%x_for_f(k, real(2*k-1)) + ur(k) = this%x_for_f(k, real(2*k-0)) + enddo + call test%real_arr(3, ul, (/0.,0.,0./), 'Return position of f<') + call test%real_arr(3, um, (/0.5,0.5,0.5/), 'Return position of f=') + call test%real_arr(3, ur, (/1.,1.,1./), 'Return position of f>') + + PCM_unit_tests = test%summarize('PCM_unit_tests') + +end function PCM_unit_tests + +!> \namespace recon1d_pcm +!! + +end module Recon1d_PCM diff --git a/src/ALE/Recon1d_type.F90 b/src/ALE/Recon1d_type.F90 new file mode 100644 index 0000000000..6d91dc1c72 --- /dev/null +++ b/src/ALE/Recon1d_type.F90 @@ -0,0 +1,349 @@ +!> A generic type for vertical 1D reconstructions +module Recon1d_type + +! This file is part of MOM6. See LICENSE.md for the license. + +implicit none ; private + +public Recon1d +public testing + +!> The base class for implementations of 1D reconsrtuctions +type, abstract :: Recon1d + + integer :: n !< Number of cells in column + real, allocatable, dimension(:) :: u_mean !< Cell mean [A] + + logical :: do_legacy = .true. !< Whether to use legacy constructs or not + + ! Legacy representation + integer :: degree !< Degree of polynomial used in legacy representation + real, allocatable, dimension(:,:) :: poly_coef !< Polynomial coefficients in legacy representation + +contains + + ! The following functions/subroutines are deferred and must be provided specifically by each scheme + + !> Deferred implementation constructor + procedure(i_construct), deferred :: construct + !> Deferred implementation of reconstruction function + procedure(i_reconstruct), deferred :: reconstruct + !> Deferred implementation of function returning the edge values + procedure(i_lr_edge), deferred :: lr_edge + !> Deferred implementation of the average over an interval + procedure(i_average), deferred :: average + !> Deferred implementation of the finding position of a value + procedure(i_x_for_f), deferred :: x_for_f + !> Deferred implementation of unit tests for the reconstruction + procedure(i_unit_tests), deferred :: unit_tests + + ! The following functions/subroutines are shared across all reconstructions and provided by this module + !> Cell mean of cell k [A] + procedure :: cell_mean => a_cell_mean + +end type Recon1d + +!> Class to assist in unit tests +type :: testing + private + !> True if any fail has been encountered since instantiation of "testing" + logical :: state = .false. + !> Count of tests checked + integer :: num_tests_checked = 0 + !> Count of tests failed + integer :: num_tests_failed = 0 + !> If true, be verbose and write results to stdout. Default True. + logical :: verbose = .true. + !> Error channel + integer :: stderr = 0 + !> Standard output channel + integer :: stdout = 6 + !> If true, stop instantly + logical :: stop_instantly = .false. + !> Record instances that fail + integer :: ifailed(100) = 0. + !> Record label of first instance that failed + character(len=:), allocatable :: label_first_fail + + contains + procedure :: test => test !< Update the testing state + procedure :: set => set !< Set attributes + procedure :: outcome => outcome !< Return current outcome + procedure :: summarize => summarize !< Summarize testing state +! procedure :: real_scalar => real_scalar !< Compare two reals + procedure :: real_arr => real_arr !< Compare array of reals + procedure :: int_arr => int_arr !< Compare array of integers +end type + +interface + + !> Initialize a 1D reconstruction for n cells + subroutine i_construct(this, n) + import :: Recon1d + class(Recon1d), intent(out) :: this !< This reconstruction + integer, intent(in) :: n !< Number of cells in this column + end subroutine i_construct + + !> Calculate a 1D reconstructions based on h(:) and u(:) + subroutine i_reconstruct(this, h, u) + import :: Recon1d + class(Recon1d), intent(inout) :: this !< This reconstruction + real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H] + real, intent(in) :: u(*) !< Cell mean values [A] + end subroutine i_reconstruct + + !> Returns the left/right edge values in cell k of a 1D reconstruction + subroutine i_lr_edge(this, k, u_left, u_right) + import :: Recon1d + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(out) :: u_left !< Left edge value [A] + real, intent(out) :: u_right !< Right edge value [A] + end subroutine i_lr_edge + + !> Average between xa and xb for cell k of a 1D reconstruction [A] + !! + !! It is assumed that 0<=xa<=1, 0<=xb<=1, and xa<=xb + real function i_average(this, k, xa, xb) + import :: Recon1d + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: xa !< Start of averageing interval on element (0 to 1) + real, intent(in) :: xb !< End of averageing interval on element (0 to 1) + end function i_average + + !> Position at which 1D reconstruction has a value f [0..1] + !! + !! It is assumed the reconstruction is not decreasing with cell index. + real function i_x_for_f(this, k, f) + import :: Recon1d + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: f !< Value of reconstruction to solve for [A] + end function i_x_for_f + + !> Runs reconstruction unit tests and returns True for any fails, False otherwise + !! + !! Assumes single process/thread context + logical function i_unit_tests(this, verbose, stdout, stderr) + import :: Recon1d + class(Recon1d), intent(inout) :: this !< This reconstruction + logical, intent(in) :: verbose !< True, if verbose + integer, intent(in) :: stdout !< I/O channel for stdout + integer, intent(in) :: stderr !< I/O channel for stderr + end function i_unit_tests + +end interface + +contains + +!> Mean value for cell k [A] +real function a_cell_mean(this, k) + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + + a_cell_mean = this%u_mean(k) + +end function a_cell_mean + +! ========================================================================================= +! The following provide the function for the testing_type helper class +! ========================================================================================= + +!> Update the state with "test" +subroutine test(this, state, label) + class(testing), intent(inout) :: this !< This testing class + logical, intent(in) :: state !< True to indicate a fail, false otherwise + character(len=*), intent(in) :: label !< Message + + this%num_tests_checked = this%num_tests_checked + 1 + if (state) then + this%state = .true. + this%num_tests_failed = this%num_tests_failed + 1 + this%ifailed( this%num_tests_failed ) = this%num_tests_checked + if (this%num_tests_failed == 1) this%label_first_fail = label + endif + if (this%stop_instantly .and. this%state) stop 1 +end subroutine test + +!> Set attributes +subroutine set(this, verbose, stdout, stderr, stop_instantly) + class(testing), intent(inout) :: this !< This testing class + logical, optional, intent(in) :: verbose !< True or false setting to assign to verbosity + integer, optional, intent(in) :: stdout !< The stdout channel to use + integer, optional, intent(in) :: stderr !< The stderr channel to use + logical, optional, intent(in) :: stop_instantly !< If true, stop immediately on error detection + + if (present(verbose)) then + this%verbose = verbose + endif + if (present(stdout)) then + this%stdout = stdout + endif + if (present(stderr)) then + this%stderr = stderr + endif + if (present(stop_instantly)) then + this%stop_instantly = stop_instantly + endif +end subroutine set + +!> Returns state +logical function outcome(this) + class(testing), intent(inout) :: this !< This testing class + outcome = this%state +end function outcome + +!> Summarize results +logical function summarize(this, label) + class(testing), intent(inout) :: this !< This testing class + character(len=*), intent(in) :: label !< Message + integer :: i + + if (this%state) then + write(this%stdout,'(a," : ",a,", ",i4," failed of ",i4," tested")') & + red('FAIL'), trim(label), this%num_tests_failed, this%num_tests_checked + write(this%stdout,'(a,100i4)') 'Failed tests:',(this%ifailed(i),i=1,this%num_tests_failed) + write(this%stdout,'(a,a)') 'First failed test: ',trim(this%label_first_fail) + write(this%stderr,'(a,100i4)') 'Failed tests:',(this%ifailed(i),i=1,this%num_tests_failed) + write(this%stderr,'(a,a)') 'First failed test: ',trim(this%label_first_fail) + write(this%stderr,'(a," : ",a)') trim(label),'FAILED' + else + write(this%stdout,'(a," : ",a,", all ",i4," tests passed")') & + green('Pass'), trim(label), this%num_tests_checked + endif + summarize = this%state +end function summarize + +!> Pads string with the color codes for red/reset +function red(string) result(newstr) + character(len=*), intent(in) :: string !< Input string + character(len=:), allocatable :: newstr !< The modified output string + newstr = achar(27)//'[31m'//trim(string)//achar(27)//'[0m' +end function red + +!> Pads string with the color codes for green/reset +function green(string) result(newstr) + character(len=*), intent(in) :: string !< Input string + character(len=:), allocatable :: newstr !< The modified output string + newstr = achar(27)//'[32m'//trim(string)//achar(27)//'[0m' +end function green + +!! !> Compare u_test to u_true, report, and return true if a difference larger than tol is measured +!! !! +!! !! If in verbose mode, display results to stdout +!! !! If a difference is measured, display results to stdout and stderr +!! subroutine real_scalar(this, u_test, u_true, label) +!! class(testing), intent(inout) :: this !< This testing class +!! real, intent(in) :: u_test !< Values to test [A] +!! real, intent(in) :: u_true !< Values to test against (correct answer) [A] +!! character(len=*), intent(in) :: label !< Message +!! ! Local variables +!! logical :: this_test +!! real :: err +!! +!! this_test = .false. +!! +!! ! Scan for any mismatch between u_test and u_true +!! err = u_test - u_true +!! if (abs(err) > 0.) this_test = .true. +!! +!! ! If either being verbose, or an error was measured then display results +!! if (this_test .or. this%verbose) then +!! if (this_test) then +!! write(this%stdout,'(3(a,1pe24.16,1x),a)') 'Calculated value =',u_test,'correct value =',u_true,'error =',err,label +!! write(this%stderr,'(3(a,1pe24.16,1x),a)') 'Calculated value =',u_test,'correct value =',u_true,'error =',err,label +!! else +!! write(this%stdout,'(2(a,1pe24.16,1x),a)') 'Calculated value =',u_test,'correct value =',u_true,label +!! endif +!! endif +!! +!! call this%test( this_test, label ) ! Updates state and counters in this +!! end subroutine real_scalar + +!> Compare u_test to u_true, report, and return true if a difference larger than tol is measured +!! +!! If in verbose mode, display results to stdout +!! If a difference is measured, display results to stdout and stderr +subroutine real_arr(this, n, u_test, u_true, label, tol) + class(testing), intent(inout) :: this !< This testing class + integer, intent(in) :: n !< Number of cells in u + real, dimension(n), intent(in) :: u_test !< Values to test [A] + real, dimension(n), intent(in) :: u_true !< Values to test against (correct answer) [A] + character(len=*), intent(in) :: label !< Message + real, optional, intent(in) :: tol !< The tolerance for differences between u and u_true [A] + ! Local variables + integer :: k + logical :: this_test + real :: tolerance, err + + tolerance = 0.0 + if (present(tol)) tolerance = tol + this_test = .false. + + ! Scan for any mismatch between u_test and u_true + do k = 1, n + if (abs(u_test(k) - u_true(k)) > tolerance) this_test = .true. + enddo + + ! If either being verbose, or an error was measured then display results + if (this_test .or. this%verbose) then + write(this%stdout,'(a4,2a24,1x,a)') 'k','Calculated value','Correct value',label + if (this_test) write(this%stderr,'(a4,2a24,1x,a)') 'k','Calculated value','Correct value',label + do k = 1, n + err = u_test(k) - u_true(k) + if (abs(err) > tolerance) then + write(this%stdout,'(i4,1p2e24.16,a,1pe24.16,a)') k, u_test(k), u_true(k), & + ' err=', err, red(' <--- WRONG') + write(this%stderr,'(i4,1p2e24.16,a,1pe24.16,a)') k, u_test(k), u_true(k), & + ' err=', err, ' <--- WRONG' + else + write(this%stdout,'(i4,1p2e24.16)') k, u_test(k), u_true(k) + endif + enddo + endif + + call this%test( this_test, label ) ! Updates state and counters in this +end subroutine real_arr + +!> Compare i_test to i_true and report and return true if a difference is found +!! +!! If in verbose mode, display results to stdout +!! If a difference is measured, display results to stdout and stderr +subroutine int_arr(this, n, i_test, i_true, label) + class(testing), intent(inout) :: this !< This testing class + integer, intent(in) :: n !< Number of cells in u + integer, dimension(n), intent(in) :: i_test !< Values to test [A] + integer, dimension(n), intent(in) :: i_true !< Values to test against (correct answer) [A] + character(len=*), intent(in) :: label !< Message + ! Local variables + integer :: k + logical :: this_test + + this_test = .false. + + ! Scan for any mismatch between u_test and u_true + do k = 1, n + if (i_test(k) .ne. i_true(k)) this_test = .true. + enddo + + if (this%verbose) then + write(this%stdout,'(a12," : calculated =",30i3)') label, i_test + write(this%stdout,'(12x," correct =",30i3)') i_true + if (this_test) write(this%stdout,'(3x,a,8x,"error =",30i3)') red('FAIL --->'), i_test(:) - i_true(:) + endif + if (this_test) then + write(this%stderr,'(a12," : calculated =",30i3)') label, i_test + write(this%stderr,'(12x," correct =",30i3)') i_true + write(this%stderr,'(" FAIL ---> error =",30i3)') i_test(:) - i_true(:) + endif + + call this%test( this_test, label ) ! Updates state and counters in this +end subroutine int_arr + +!> \namespace recon1d_type +!! +!! \section section_recon1d_type Generic vertical reconstruction type +!! + +end module Recon1d_type