Skip to content

Commit

Permalink
Refactor of vertical reconstruction with six new schemes
Browse files Browse the repository at this point in the history
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 member methods available only
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.

What changed
------------
- removed remapping_attic.F90
- added 16 modules name Recon1d_*.F90, which extend a base Recon1d_type.F90
- MOM_remapping.F90
  - rewrote remapping_unit_tests to be tidier
  - added tests for
    1) invariance when remapping to/from the same grid;
    2) preservation of a uniform value (many of the old schemes fail);
    3) internal consistency such as monotonic ordering of values at
       edges, no internal extrema, etc., but only for the new Recon1d_* schemes.
- added _Vertical_Reconstruction.dox

Six new schemes are added: C_PLM_CW, C_PLM_CWK, C_MPLM_CWK, C_EMPLM_CWK,
C_PPM_CWK, and C_EPPM_CWK.

Ten existing schemes were ported to the new class, 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 discovery 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() checks that a reconstruction has the
  properties that it should (e.g. reconstruction values are monotonic)
  Theses properties are defined by each scheme, although many are shared
  between schemes.
- 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)

Testing
-------
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.
These tests use seeded random numbers which has proven useful for
debugging/development, but are not necessarily useful in the CI. Nevertheless,
I've made it part of remapping_unit_tests() so there is a template for
debugging/testing new schemes when we add them.

The base class also provides a testing type for convenience. In time, this
could (should) 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 using anything in 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  of this commit
  (expand the "Display timing results");
- e.g. timings for PPM_HYBGEN and C_PPM_HYBGEN show a marginal (~10%) speed
  up due to code style, and C_PPM_CWK has a more substantial speed up (~20%)
  due to a lighter algorithm;
- I imagine we will see further speed ups if we worked on arrays rather than
  columns, which is an extension that could be implemented later.

Documentation
-------------
- the new module APIs are documented;
- a new page, `_Vertical_Reconstruction.dox`, details the vertical
  reconstruction schemes, tabulates the values of `REMPAPPING_SCHEME`,
  and summarize some properties.
  • Loading branch information
adcroft committed Oct 16, 2024
1 parent 80d8b5f commit 061857b
Show file tree
Hide file tree
Showing 25 changed files with 6,410 additions and 951 deletions.
6 changes: 3 additions & 3 deletions .testing/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -269,16 +269,16 @@ $(BUILD)/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests
# Build executables
.NOTPARALLEL:$(foreach e,$(UNIT_EXECS),$(BUILD)/unit/$(e))
$(BUILD)/unit/test_%: $(BUILD)/unit/Makefile FORCE
cd $(@D) && $(TIME) $(MAKE) $(@F) -j
cd $(@D) && $(TIME) $(MAKE) $(@F)
$(BUILD)/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90)

.NOTPARALLEL:$(foreach e,$(TIMING_EXECS),$(BUILD)/timing/$(e))
$(BUILD)/timing/time_%: $(BUILD)/timing/Makefile FORCE
cd $(@D) && $(TIME) $(MAKE) $(@F) -j
cd $(@D) && $(TIME) $(MAKE) $(@F)
$(BUILD)/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90)

$(BUILD)/%/MOM6: $(BUILD)/%/Makefile FORCE
cd $(@D) && $(TIME) $(MAKE) $(@F) -j
cd $(@D) && $(TIME) $(MAKE) $(@F)

# Target codebase should use its own build system
$(BUILD)/target/MOM6: $(BUILD)/target FORCE | $(TARGET_CODEBASE)
Expand Down
28 changes: 22 additions & 6 deletions config_src/drivers/timing_tests/time_MOM_remapping.F90
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,27 @@ program time_MOM_remapping
implicit none

type(remapping_CS) :: CS
integer, parameter :: nk=75, nij=20*20, nits=10, nsamp=100, nschemes = 2
character(len=10) :: scheme_labels(nschemes)
integer, parameter :: nk=75, nij=20*20, nits=10, nsamp=100, nschemes = 19
character(len=16) :: scheme_labels(nschemes) = [ character(len=16) :: &
'PCM', &
'C_PCM', &
'PLM', &
'C_MPLM_WA', &
'C_EMPLM_WA', &
'C_PLM_HYBGEN', &
'C_PLM_CW', &
'C_PLM_CWK', &
'C_MPLM_WA_POLY', &
'C_EMPLM_WA_POLY', &
'C_MPLM_CWK', &
'PPM_CW', &
'PPM_HYBGEN', &
'C_PPM_H4_2018', &
'C_PPM_H4_2019', &
'C_PPM_HYBGEN', &
'C_PPM_CW', &
'C_PPM_CWK', &
'C_EPPM_CWK' ]
real, dimension(nschemes) :: timings ! Time for nits of nij calls for each scheme [s]
real, dimension(nschemes) :: tmean ! Mean time for a call [s]
real, dimension(nschemes) :: tstd ! Standard deviation of time for a call [s]
Expand All @@ -30,9 +49,6 @@ program time_MOM_remapping
seed(:) = 102030405
call random_seed(put=seed)

scheme_labels(1) = 'PCM'
scheme_labels(2) = 'PLM'

! Set up some test data (note: using k,i indexing rather than i,k)
allocate( u0(nk,nij), h0(nk,nij), u1(nk,nij), h1(nk,nij) )
call random_number(u0) ! In range 0-1
Expand All @@ -59,7 +75,7 @@ program time_MOM_remapping
do isamp = 1, nsamp
! Time reconstruction + remapping
do ischeme = 1, nschemes
call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)))
call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)), nk=nk)
call cpu_time(start)
do iter = 1, nits ! Make many passes to reduce sampling error
do ij = 1, nij ! Calling nij times to make similar to cost in MOM_ALE()
Expand Down
14 changes: 13 additions & 1 deletion config_src/drivers/unit_tests/test_MOM_remapping.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,18 @@ program test_MOM_remapping

use MOM_remapping, only : remapping_unit_tests

if (remapping_unit_tests(.true.)) stop 1
integer :: n !< Number of arguments, or tests
character(len=12) :: cmd_ln_arg !< Command line argument (if any)

n = command_argument_count()

if (n==1) then
call get_command_argument(1, cmd_ln_arg)
read(cmd_ln_arg,*) n
else
n = 3000 ! Fallback value if no argument provided
endif

if (remapping_unit_tests(.true., num_comp_samp=n)) stop 1

end program test_MOM_remapping
1 change: 1 addition & 0 deletions docs/discrete_space.rst
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,5 @@ algorithm.
api/generated/pages/Discrete_Coriolis
api/generated/pages/Discrete_PG
api/generated/pages/Energetic_Consistency
api/generated/pages/Vertical_Reconstruction
api/generated/pages/Discrete_OBC
14 changes: 14 additions & 0 deletions docs/zotero.bib
Original file line number Diff line number Diff line change
Expand Up @@ -2946,3 +2946,17 @@ @article{Young1994
pages={1812--1826},
year={1994}
}

@article{van_leer_1977,
title = {Towards the ultimate conservative difference scheme. {IV}. {A} new approach to numerical convection},
volume = {23},
issn = {0021-9991},
doi = {10.1016/0021-9991(77)90095-X},
number = {3},
journal = {Journal of Computational Physics},
author = {Van Leer, Bram},
month = mar,
year = {1977},
pages = {276--299},
}

1,060 changes: 780 additions & 280 deletions src/ALE/MOM_remapping.F90

Large diffs are not rendered by default.

148 changes: 148 additions & 0 deletions src/ALE/Recon1d_EMPLM_CWK.F90
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
Loading

0 comments on commit 061857b

Please sign in to comment.