forked from NOAA-GFDL/MOM6
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Refactor of vertical reconstruction with new schemes
This refactor is motivated by i) adding a set of new reconstructions schemes that are ready-to-use-now for remapping, but ii) will be an essential part of a new grid generator that will use some of the new methods available as a result of the refactor. The "refactor" actually duplicates many existing schemes, leaving the original code in place to provide a period of transition - just in case there's some major reason to revert to the old style of code. Six new schemes are added: C_PLM_CW, C_PLM_CWK, C_MPLM_CWK, C_EMPLM_CWK, C_PPM_CWK, and C_EPPM_CWK. Eight existing schemes were ported to the new class (these are the most widely used ones), along with their non-monotonic behaviors or ill-formed expressions. Others can be ported later if needed, although one result of the refactor was the discvovery that none of the existing schemes worked as desired. The new form uses a class (Recon1d) to define a uniform API for all the reconstruction schemes. The methods of the class include: - init() that creates work arrays and stores parameters such as `h_neglect` - reconstruct() that causes the reconstruction parameters to be calculated - average() that returns the average between two points in a reconstruction (used in the remap_src_to_sub_grid() functions) - unit_tests() method - check_reconstruction() method that checks tha a reconstruction has the properties that it should (e.g. reconstruction values are monotonic) - f() that evaluates the reconstruction at a point (this will be used in the new grid-generator) - dfdx() that evaluates the derivative of the reconstruction at a point (this will be used in the new grid-generator) An advantage of the new class approach is that the check_reconstruction() method can be customized to each reconstruction, which allows us to document which properties (e.g. monotonicity) are satisfied and which are not. The base class also provides a testing type for convenience. In time, this could be moved out to a stand-alone module, so long as it does not become cumbersome. Currently this type is hidden to avoid expanding the network of dependencies, particularly to the `framework/` directory which is not needed by the class. Performance - the new class-based schemes are all faster than their counterparts in the original code style - example numbers can be found in the build-test-perfmon job (expand the "Display timing results") - e.g. timings for PPM_HYBGEN (8.68) and C_PPM_HYBGEN (7.98) show a marginal speed up due to code style, and C_PPM_CWK (6.91) has a clear speed up due to a new algorithm - I imagine we will see further speed ups if we work on arrays rather than columns which can be implemented as methods later Documentation - the new module APIs are documented - a new page specific to vertical reconstruction added that tabulates the schemes (see https://mom6-doctest.readthedocs.io/en/latest/api/generated/pages/Vertical_Reconstruction.html)
- Loading branch information
Showing
26 changed files
with
6,423 additions
and
954 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,148 @@ | ||
| !> Piecewise Linear Method 1D reconstruction in index space and boundary extrapolation | ||
| !! | ||
| !! This implementation of PLM follows Colella and Woodward, 1984 \cite colella1984, except for assuming | ||
| !! uniform resolution so that the method is independent of grid spacing. The cell-wise reconstructions | ||
| !! are limited so that the edge values (which are also the extrema in a cell) are bounded by the neighbors. | ||
| !! The slope of the first and last cells are set so that the first interior edge values match the interior | ||
| !! cell (i.e. extrapolates from the interior). | ||
| module Recon1d_EMPLM_CWK | ||
|
|
||
| ! This file is part of MOM6. See LICENSE.md for the license. | ||
|
|
||
| use Recon1d_type, only : testing | ||
| use Recon1d_MPLM_CWK, only : MPLM_CWK | ||
|
|
||
| implicit none ; private | ||
|
|
||
| public EMPLM_CWK, testing | ||
|
|
||
| !> PLM reconstruction following Colella and Woodward, 1984 | ||
| !! | ||
| !! Implemented by extending recon1d_mplm_cwk. | ||
| !! | ||
| !! The source for the methods ultimately used by this class are: | ||
| !! - init() -> recon1d_mplm_cwk -> recon1d_plm_cw.init() | ||
| !! - reconstruct() *locally defined | ||
| !! - average() -> recon1d_mplm_cwk -> recon1d_plm_cw.average() | ||
| !! - f() -> recon1d_mplm_cwk -> recon1d_plm_cw.f() | ||
| !! - dfdx() -> recon1d_mplm_cwk -> recon1d_plm_cw.dfdx() | ||
| !! - check_reconstruction() -> recon1d_mplm_cwk.check_reconstruction() | ||
| !! - unit_tests() *locally defined | ||
| !! - destroy() -> recon1d_mplm_cwk -> recon1d_plm_cw.destroy() | ||
| !! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid() | ||
| !! - init_parent() -> init() | ||
| !! - reconstruct_parent() -> recon1d_mplm_cwk.reconstruct() | ||
| type, extends (MPLM_CWK) :: EMPLM_CWK | ||
|
|
||
| contains | ||
| !> Implementation of the EMPLM_CWK reconstruction | ||
| procedure :: reconstruct => reconstruct | ||
| !> Implementation of unit tests for the EMPLM_CWK reconstruction | ||
| procedure :: unit_tests => unit_tests | ||
|
|
||
| end type EMPLM_CWK | ||
|
|
||
| contains | ||
|
|
||
| !> Calculate a 1D PLM reconstructions based on h(:) and u(:) | ||
| subroutine reconstruct(this, h, u) | ||
| class(EMPLM_CWK), intent(inout) :: this !< This reconstruction | ||
| real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H] | ||
| real, intent(in) :: u(*) !< Cell mean values [A] | ||
| ! Local variables | ||
| real :: slp ! The PLM slopes (difference across cell) [A] | ||
| real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as | ||
| ! differences across the cell [A] | ||
| real :: u_min, u_max ! Minimum and maximum value across cell [A] | ||
| real :: u_l, u_r, u_c ! Left, right, and center values [A] | ||
| real :: u_e(this%n+1) ! Average of edge values [A] | ||
| integer :: k, n | ||
|
|
||
| n = this%n | ||
|
|
||
| call this%reconstruct_parent(h, u) | ||
|
|
||
| this%ur(1) = this%ul(2) | ||
| this%ul(1) = u(1) + ( u(1) - this%ur(1) ) | ||
|
|
||
| this%ul(n) = this%ur(n-1) | ||
| this%ur(n) = u(n) + ( u(n) - this%ul(n) ) | ||
|
|
||
| end subroutine reconstruct | ||
|
|
||
| !> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise | ||
| logical function unit_tests(this, verbose, stdout, stderr) | ||
| class(EMPLM_CWK), 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] | ||
| real, allocatable :: ull(:), urr(:) ! test values [A] | ||
| type(testing) :: test ! convenience functions | ||
| integer :: k | ||
|
|
||
| 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 this%init(3) | ||
| call test%test( this%n /= 3, 'Setting number of levels') | ||
| allocate( um(3), ul(3), ur(3), ull(3), urr(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 | ||
| ul(k) = this%f(k, 0.) | ||
| um(k) = this%f(k, 0.5) | ||
| ur(k) = this%f(k, 1.) | ||
| enddo | ||
| call test%real_arr(3, ul, (/0.,2.,4./), 'Evaluation on left edge') | ||
| call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center') | ||
| call test%real_arr(3, ur, (/2.,4.,6./), 'Evaluation on right edge') | ||
|
|
||
| do k = 1, 3 | ||
| ul(k) = this%dfdx(k, 0.) | ||
| um(k) = this%dfdx(k, 0.5) | ||
| ur(k) = this%dfdx(k, 1.) | ||
| enddo | ||
| call test%real_arr(3, ul, (/2.,2.,2./), 'dfdx on left edge') | ||
| call test%real_arr(3, um, (/2.,2.,2./), 'dfdx in center') | ||
| call test%real_arr(3, ur, (/2.,2.,2./), 'dfdx on right edge') | ||
|
|
||
| do k = 1, 3 | ||
| um(k) = this%average(k, 0.25, 0.75) ! Average from x=0.25 to 0.75 in each cell | ||
| enddo | ||
| call test%real_arr(3, um, (/1.,3.,5./), 'Return interval average') | ||
|
|
||
| call this%destroy() | ||
| deallocate( um, ul, ur, ull, urr ) | ||
|
|
||
| allocate( um(4), ul(4), ur(4) ) | ||
| call this%init(4) | ||
|
|
||
| ! These values lead to non-monotonic reconstuctions which are | ||
| ! valid for transport problems but not always appropriate for | ||
| ! remapping to arbitrary resolution grids. | ||
| ! The O(h^2) slopes are -, 2, 2, - and the limited | ||
| ! slopes are 0, 1, 1, 0 so the everywhere the reconstructions | ||
| ! are bounded by neighbors but ur(2) and ul(3) are out-of-order. | ||
| call this%reconstruct( (/1.,1.,1.,1./), (/0.,3.,4.,7./) ) | ||
| do k = 1, 4 | ||
| ul(k) = this%f(k, 0.) | ||
| ur(k) = this%f(k, 1.) | ||
| enddo | ||
| call test%real_arr(4, ul, (/-2.5,2.5,3.5,4.5/), 'Evaluation on left edge') | ||
| call test%real_arr(4, ur, (/2.5,3.5,4.5,9.5/), 'Evaluation on right edge') | ||
|
|
||
| deallocate( um, ul, ur ) | ||
|
|
||
| unit_tests = test%summarize('EMPLM_CWK:unit_tests') | ||
|
|
||
| end function unit_tests | ||
|
|
||
| !> \namespace recon1d_emplm_cwk | ||
| !! | ||
|
|
||
| end module Recon1d_EMPLM_CWK |
Oops, something went wrong.