Skip to content

Commit

Permalink
Wetting and Drying Master commit
Browse files Browse the repository at this point in the history 
tmp commit for example

Adds Gaussian wave test case for Delaware

Adds an initial test case of a gravity wave entrying the Delaware Bay

Wetting and Drying Master commit

Adds overflow to wetting / drying test suite stub

Working example of wave going onshore

Adds passive tracer (tracer1) to onshore wave case

Enhances strength of depression wave (stable)

Fixes bug in onshow wave case, drying causes break

Adds wetting and drying (non/and) working tests

Adds check to RK4 for wetting and drying

Draft working case for RK4 (mostly)

minor bug left that condition for minimum can
be violated by a little bit-- implies that
hard-line threshold needs to have an eps to fix

Working test case for RK4 for drying

Strange that works best for thiner layer, which implies
that more damping is applied earlier in the simulation (stronger
damping response to prevent layer thickness from getting too small).

Outputs minimum thickness as a diagnostic

Improved wetting / drying version

Serious wave reflection because drying prevention
is not coupled back to the velocity field via a tendency.

Could enforce this by directly applying the result to the velocity
field too... (or use some type of damped restoring).

Working code, more realistic view of drying

Problem: still don't have clean constraint on early warning for drying

Algorithm best approach known so far for RK4

Works with drying case

Fixed drying condition check

Verifies no drying has occurred until after
timestep completes.

Changes to ensure stability under strong drying

Note, new config_zero_drying_velocity can be used to
estimate an appropriate config_drying_min_cell_height value.

Clean up removing print debugs

Make flow restriction go to zero faster

This prevents threshold from being exceeded in standard use.

Adds new drying slope test case

Adds and fixes up test cases

Note, potential issue with layerEdgeThickness centered vs
upwinding approaches

Adds 10km drying_slope cases

Note, need to verify that the 10km cases are comparable with the 20km
in the future...

Adds in flag for selection of thickness flux

Note, code is Not fully vetted

Includes path for test suite output

Setup for split explicit barotropic drying dev

Configured split explicit test cases for debug output

Defaults debugging to off for split explicit

Fixes indendation

f, fixed indentation

Initial/broken split explicit draft (error check)

We have a test for failure of the algorithm and draft code.
This is an initial commit as a coding "sketch".

Quasi-working draft of split-explicit wetting/drying

Uses simple first-order upwinding approach with flux
limiting of outflowing velocity

f, updates comment

f, cleans up comments

Working version for barotropic fluxes

Fixes bug when using RK4 for compatability with SE

Need to select appropriate layer thickness for RK4 or split
explicit algorithms

Fixes merge bug

Moves wave test case to RK4 for wetting/drying

Fixes wetting and drying parallelization bug

Previously, normalVelocity and normalTransportVelocity
were not updated with a halo exchange following computation
to prevent drying.

Modifies vel tendency to be zero for drying case

This ensures that for RK4 no spurious velocity forcing can
occur that allows a drying cell velocity to increase without
bound.

Zero out velocity for drying cells

Example of draining Delaware Bay case

Memory leak was observed, will try to setup
working inundation case now

Configuration update Delaware Bay inundation

Note, spurious gravity waves need additional examination

Whitespace cleanup

Adds two new resolutions to Delaware case

Also updates 3km resolution parameters

Fixes to Delaware workflows from cross-comparisons

Updates high res Delaware del4 for better scaling

Tweak to timestep

Adds higher resolution Delaware cases

Draft spreading mound test case

Thacker test case of analytical solution of quadratic mount movement on
flat surface, includes nonlinearity but not Coriolis.  Preliminary, not
compltetely general, implementation.

Added missing forward file (squash with previous)

Temp commits for merge (revisit)
  • Loading branch information
@pwolfram
pwolfram committed Feb 21, 2019
1 parent e898586 commit d242282
Show file tree
Hide file tree
Showing 167 changed files with 7,701 additions and 101 deletions.
2 changes: 1 addition & 1 deletion Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ gfortran:
"CFLAGS_OPT = -O3 -m64" \
"CXXFLAGS_OPT = -O3 -m64" \
"LDFLAGS_OPT = -O3 -m64" \
"FFLAGS_DEBUG = -g -m64 -ffree-line-length-none -fconvert=big-endian -ffree-form -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow" \
"FFLAGS_DEBUG = -g -m64 -ffree-line-length-none -fconvert=big-endian -ffree-form -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow -Wall" \
"CFLAGS_DEBUG = -g -m64" \
"CXXFLAGS_DEBUG = -O3 -m64" \
"LDFLAGS_DEBUG = -g -m64" \
Expand Down
43 changes: 43 additions & 0 deletions src/core_ocean/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -916,6 +916,40 @@
possible_values="any positive real, typically 1.0e-3"
/>
</nml_record>
<nml_record name="wetting_drying" mode="init;forward">
<nml_option name="config_use_wetting_drying" type="logical" default_value=".false." units="unitless"
description="If true, use wetting and drying algorithm to allow for dry cells to config_min_cell_height."
possible_values=".true. or .false."
/>
<nml_option name="config_prevent_drying" type="logical" default_value=".true." units="unitless"
description="If true, prevent cells from drying past config_min_cell_height."
possible_values=".true. or .false."
/>
<nml_option name="config_drying_min_cell_height" type="real" default_value="1.0e-2" units="m"
description="Minimum allowable cell height under drying. Cell to be kept wet to at least this thickness."
possible_values="any positive real, typically 1.0e-3"
/>
<nml_option name="config_zero_drying_velocity" type="logical" default_value=".false." units="unitless"
description="If true, just zero out velocity that is contributing to drying for cell that is drying. This option can be used to estimate acceptable minimum thicknesses for a run."
possible_values=".true. or .false."
/>
<nml_option name="config_debug_split_explicit" type="logical" default_value=".false." units="unitless"
description="If true, output ssh from split explicit (1 proc only)."
possible_values=".true. or .false."
/>
<nml_option name="config_debug_split_explicit_filename" type="character" default_value="debugSplitExplicit.dat" units="unitless"
description="File name for output for debugging split explicit."
possible_values="Any acceptable filename pattern."
/>
<nml_option name="config_verify_not_dry" type="logical" default_value=".true." units="unitless"
description="If true, verify that cells are at least config_min_cell_height thick."
possible_values=".true. or .false."
/>
<nml_option name="config_thickness_flux_type" type="character" default_value="upwind" units="unitless"
description="If 'upwind', use upwind to evaluate the edge-value for layerThickness, i.e., layerThicknessEdge. The standard MPAS-O approach is 'centered'."
possible_values="'upwind', 'centered'"
/>
</nml_record>
<nml_record name="ocean_constants">
<nml_option name="config_density0" type="real" default_value="1026.0" units="kg m^{-3}"
description="Density used as a coefficient of the pressure gradient terms, $\rho_0$. This is a constant due to the Boussinesq approximation."
Expand Down Expand Up @@ -1309,6 +1343,7 @@
<var name="fEdge"/>
<var name="fVertex"/>
<var name="fCell"/>
<var name="bottomDepth"/>
</stream>

<stream name="output_init"
Expand Down Expand Up @@ -2313,6 +2348,14 @@
description="horizontal velocity used to transport mass and tracers"
packages="forwardMode;analysisMode"
/>
<var name="barotropicWettingVelocity" type="real" dimensions="nEdges Time" units="m s^{-1}"
description="Barotropic velocity to prevent drying of column."
packages="forwardMode"
/>
<var name="wettingVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
description="Velocity to prevent drying of cell."
packages="forwardMode"
/>
<var name="vertAleTransportTop" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^{-1}"
description="vertical transport through the layer interface at the top of the cell"
packages="forwardMode;analysisMode"
Expand Down
25 changes: 25 additions & 0 deletions src/core_ocean/mode_forward/mpas_ocn_forward_mode.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ module ocn_forward_mode
use mpas_log
use mpas_decomp
use mpas_log
use mpas_io_units

use ocn_analysis_driver
use ocn_init_routines
Expand Down Expand Up @@ -79,6 +80,7 @@ module ocn_forward_mode
use ocn_time_varying_forcing

use ocn_constants
use ocn_wetting_drying

implicit none
private
Expand Down Expand Up @@ -127,9 +129,11 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
logical, pointer :: config_do_restart, config_read_nearest_restart, config_filter_btr_mode, config_conduct_tests
character (len=StrKIND), pointer :: config_vert_coord_movement, config_pressure_gradient_type
character (len=StrKIND), pointer :: config_surface_salinity_monthly_restoring_compute_interval
character (len=StrKIND), pointer :: config_debug_split_explicit_filename

real (kind=RKIND), pointer :: config_maxMeshDensity
logical, pointer :: config_use_surface_salinity_monthly_restoring
logical, pointer :: config_debug_split_explicit
ierr = 0

!
Expand All @@ -153,6 +157,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
config_use_surface_salinity_monthly_restoring)
call mpas_pool_get_config(domain % configs, 'config_surface_salinity_monthly_restoring_compute_interval', &
config_surface_salinity_monthly_restoring_compute_interval)
call mpas_pool_get_config(domain % configs, 'config_debug_split_explicit', config_debug_split_explicit)
call mpas_pool_get_config(domain % configs, 'config_debug_split_explicit_filename', config_debug_split_explicit_filename)
!
! Read input data for model
!
Expand Down Expand Up @@ -378,6 +384,14 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
call mpas_dmpar_exch_group_add_field(domain, 'subcycleFields', 'normalBarotropicVelocitySubcycle')
call mpas_dmpar_exch_group_build_reusable_buffers(domain, 'subcycleFields')

if (config_debug_split_explicit) then
! open up temp file
call mpas_new_unit(debugSplitExplicitUnit)
! note, presumes we are running on a single core here...
open(unit=debugSplitExplicitUnit, file=trim(config_debug_split_explicit_filename), &
action='write', status='unknown')
end if

end function ocn_forward_mode_init!}}}

!***********************************************************************
Expand Down Expand Up @@ -708,6 +722,7 @@ end function ocn_forward_mode_run!}}}
function ocn_forward_mode_finalize(domain) result(iErr)!{{{

type (domain_type), intent(inout) :: domain
logical, pointer :: config_debug_split_explicit

integer :: ierr

Expand All @@ -718,6 +733,16 @@ function ocn_forward_mode_finalize(domain) result(iErr)!{{{
call mpas_destroy_clock(domain % clock, ierr)

call mpas_decomp_destroy_decomp_list(domain % decompositions)
call mpas_pool_get_config(domain % configs, 'config_debug_split_explicit', config_debug_split_explicit)

! assumes running on a single core
if (config_debug_split_explicit) then
! close file
flush(debugSplitExplicitUnit)
close(debugSplitExplicitUnit)
call mpas_release_unit(debugSplitExplicitUnit)
end if


end function ocn_forward_mode_finalize!}}}

Expand Down
84 changes: 80 additions & 4 deletions src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ module ocn_time_integration_rk4
use ocn_equation_of_state
use ocn_vmix
use ocn_time_average_coupled
use ocn_wetting_drying

use ocn_effective_density_in_land_ice

Expand Down Expand Up @@ -128,6 +129,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
logical, pointer :: config_disable_vel_all_tend
logical, pointer :: config_disable_tr_all_tend
real (kind=RKIND), pointer :: config_mom_del4
real (kind=RKIND), pointer :: config_drying_min_cell_height
logical, pointer :: config_use_wetting_drying, config_verify_not_dry, config_prevent_drying, config_zero_drying_velocity
character (len=StrKIND), pointer :: config_land_ice_flux_mode

! State indices
Expand Down Expand Up @@ -177,7 +180,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
type (field2DReal), pointer :: normalizedRelativeVorticityEdgeField, divergenceField, relativeVorticityField

! State/Tend Field Pointers
type (field2DReal), pointer :: normalVelocityField, layerThicknessField
type (field2DReal), pointer :: normalVelocityField, normalTransportVelocityField, wettingVelocityField
type (field3DReal), pointer :: tracersGroupField

! Tracer Group Iteartion
Expand All @@ -197,6 +200,11 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call mpas_pool_get_config(domain % configs, 'config_disable_vel_all_tend', config_disable_vel_all_tend)
call mpas_pool_get_config(domain % configs, 'config_disable_thick_all_tend', config_disable_thick_all_tend)
call mpas_pool_get_config(domain % configs, 'config_disable_tr_all_tend', config_disable_tr_all_tend)
call mpas_pool_get_config(domain % configs, 'config_use_wetting_drying', config_use_wetting_drying)
call mpas_pool_get_config(domain % configs, 'config_prevent_drying', config_prevent_drying)
call mpas_pool_get_config(domain % configs, 'config_verify_not_dry', config_verify_not_dry)
call mpas_pool_get_config(domain % configs, 'config_drying_min_cell_height', config_drying_min_cell_height)
call mpas_pool_get_config(domain % configs, 'config_zero_drying_velocity', config_zero_drying_velocity)

!
! Initialize time_levs(2) with state at current time
Expand Down Expand Up @@ -445,6 +453,29 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
endif
! require that cells don't dry out
if (config_use_wetting_drying) then
call mpas_timer_start("RK4-prevent drying")

! compute wetting velocity to prevent drying of cell (sets up start of next iterate to not dry)
if (config_prevent_drying) then
block => domain % blocklist
do while (associated(block))
call ocn_prevent_drying_rk4(block, dt, rk_substep_weights(rk_step), config_zero_drying_velocity, err)
block => block % next
end do
! exchange fields for parallelization, may not be needed-- this should be checked for efficiency
call mpas_pool_get_field(diagnosticsPool, 'normalTransportVelocity', normalTransportVelocityField)
call mpas_pool_get_field(statePool, 'normalVelocity', normalVelocityField, 1)
call mpas_dmpar_exch_halo_field(normalVelocityField)
call mpas_dmpar_exch_halo_field(normalTransportVelocityField)
call mpas_pool_get_field(diagnosticsPool, 'wettingVelocity', wettingVelocityField)
call mpas_dmpar_exch_halo_field(wettingVelocityField)
end if

call mpas_timer_stop("RK4-prevent drying")
call mpas_threading_barrier()
end if

! Compute tendencies for velocity, thickness, and tracers.
! In RK4 notation, we are computing the right hand side f(t,y),
Expand All @@ -456,6 +487,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call ocn_time_integrator_rk4_compute_vel_tends( block, dt, rk_substep_weights(rk_step), err )

call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err )

block => block % next
end do

Expand Down Expand Up @@ -495,6 +527,24 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call mpas_timer_stop("RK4 tracer prognostic halo update")
call mpas_threading_barrier()

!! ensure drying momentum tendency is zero
!if (config_use_wetting_drying) then
! call mpas_timer_start("RK4-prevent drying momentum tendency")

! ! compute wetting velocity to prevent drying of cell (sets up start of next iterate to not dry)
! if (config_prevent_drying) then
! block => domain % blocklist
! do while (associated(block))
! call ocn_modify_tendencies_drying_rk4(block, dt, rk_substep_weights(rk_step), config_zero_drying_velocity, err)
! block => block % next
! end do
! end if

! call mpas_timer_stop("RK4-prevent drying momentum tendency")
! call mpas_threading_barrier()
!end if


! Compute next substep state for velocity, thickness, and tracers.
! In RK4 notation, we are computing y_n + a_j k_j.

Expand Down Expand Up @@ -533,6 +583,23 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
! END RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! verify that cells are not dry at conclusion of time step
if (config_use_wetting_drying) then
call mpas_timer_start("RK4- check wet cells")

! ensure existing layerThickness is valid
if (config_verify_not_dry) then
block => domain % blocklist
do while (associated(block))
call ocn_wetting_drying_verify(block, config_drying_min_cell_height, err)
block => block % next
end do
end if

call mpas_timer_stop("RK4- check wet cells")
call mpas_threading_barrier()
end if

call mpas_timer_stop("RK4-main loop")
call mpas_threading_barrier()

Expand Down Expand Up @@ -776,7 +843,7 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(block, dt, rkWeight, err)!{
endif
call mpas_threading_barrier()
call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1)
call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt)
call mpas_threading_barrier()
end subroutine ocn_time_integrator_rk4_compute_vel_tends!}}}
Expand Down Expand Up @@ -977,7 +1044,7 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{
endif
call mpas_threading_barrier()
call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1)
call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt)
call mpas_threading_barrier()
if (associated(highFreqThicknessProvis)) then
Expand Down Expand Up @@ -1025,6 +1092,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{
real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, layerThicknessProvis, layerThicknessTend
real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceCur, lowFreqDivergenceProvis, lowFreqDivergenceTend
real (kind=RKIND), dimension(:, :), pointer :: normalTransportVelocity, normalGMBolusVelocity
real (kind=RKIND), dimension(:, :), pointer :: wettingVelocity
real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupCur, tracersGroupProvis, tracersGroupTend
Expand Down Expand Up @@ -1075,13 +1143,15 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{
call mpas_pool_get_array(diagnosticsPool, 'normalTransportVelocity', normalTransportVelocity)
call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', normalGMBolusVelocity)
call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity)
call mpas_threading_barrier()
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdges
do k = 1, maxLevelEdgeTop(iEdge)
normalVelocityProvis(k, iEdge) = normalVelocityCur(k, iEdge) + rkWeight &
* normalVelocityTend(k, iEdge)
normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge)*(1.0_RKIND - wettingVelocity(k, iEdge))
end do
end do
!$omp end do
Expand All @@ -1092,6 +1162,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkWeight, err)!{
do k = 1, maxLevelCell(iCell)
layerThicknessProvis(k, iCell) = layerThicknessCur(k, iCell) + rkWeight &
* layerThicknessTend(k, iCell)
!call mpas_log_write('layer tend = $r', realArgs=(/layerThicknessTend(k,iCell)/))
end do
end do
!$omp end do
Expand Down Expand Up @@ -1191,13 +1262,14 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
integer, pointer :: nCells, nEdges
integer :: iCell, iEdge, k
type (mpas_pool_type), pointer :: statePool, tendPool, meshPool
type (mpas_pool_type), pointer :: statePool, tendPool, meshPool, diagnosticsPool
type (mpas_pool_type), pointer :: tracersPool, tracersTendPool
real (kind=RKIND), dimension(:, :), pointer :: normalVelocityNew, normalVelocityTend
real (kind=RKIND), dimension(:, :), pointer :: layerThicknessNew, layerThicknessTend
real (kind=RKIND), dimension(:, :), pointer :: highFreqThicknessNew, highFreqThicknessTend
real (kind=RKIND), dimension(:, :), pointer :: lowFreqDivergenceNew, lowFreqDivergenceTend
real (kind=RKIND), dimension(:, :), pointer :: wettingVelocity
real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew, tracersGroupTend
Expand All @@ -1216,6 +1288,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
call mpas_pool_get_subpool(block % structs, 'state', statePool)
call mpas_pool_get_subpool(block % structs, 'tend', tendPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
call mpas_pool_get_subpool(tendPool, 'tracersTend', tracersTendPool)
Expand All @@ -1235,20 +1308,23 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
call mpas_pool_get_array(tendPool, 'highFreqThickness', highFreqThicknessTend)
call mpas_pool_get_array(tendPool, 'lowFreqDivergence', lowFreqDivergenceTend)
call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity)
call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
!$omp do schedule(runtime) private(k)
do iCell = 1, nCells
do k = 1, maxLevelCell(iCell)
layerThicknessNew(k, iCell) = layerThicknessNew(k, iCell) + rkWeight * layerThicknessTend(k, iCell)
!call mpas_log_write('update tend = $r', realArgs=(/layerThicknessTend(k,iCell)/))
end do
end do
!$omp end do
!$omp do schedule(runtime)
do iEdge = 1, nEdges
normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) + rkWeight * normalVelocityTend(:, iEdge)
normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge)*(1.0_RKIND - wettingVelocity(:, iEdge))
end do
!$omp end do
Expand Down
Loading

0 comments on commit d242282

Please sign in to comment.