Skip to content
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

Ponded water/alternative solutions - object-oriented #192

Closed
wants to merge 11 commits into from
Closed

Ponded water/alternative solutions - object-oriented #192

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
3932256
alternative numerical solutions for surface runoff of ponded water
Sep 5, 2017
3ce6202
change scope of numerical options
Sep 5, 2017
0041a92
alternative numerical solutions for drainage of ponded water
Sep 6, 2017
dd3ec25
Fixes to get Martyn's code to compile
billsacks Sep 8, 2017
98da3bf
Reimplement Martyn's idea in an object-oriented way
billsacks Sep 8, 2017
cba171d
Rename flux_type to computeFlux_type
billsacks Sep 11, 2017
87406e2
Rename flux_inst to computeFlux_inst for consistency
billsacks Sep 11, 2017
b872d21
Move final flux calculation into implicitEuler
billsacks Sep 11, 2017
ea15490
Move different drainPond calculation options to their own methods
billsacks Sep 12, 2017
e4d56cd
Revert "Move different drainPond calculation options to their own met…
billsacks Sep 14, 2017
9133ffd
Merge branch 'unified_land_model' into pondedWater/alternativeSolutio…
billsacks Sep 14, 2017
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
111 changes: 105 additions & 6 deletions src/biogeophys/SoilHydrologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,10 @@ module SoilHydrologyMod
use LandunitType , only : lun
use ColumnType , only : col
use PatchType , only : patch

use computeFluxMod ! template for the flux module
use implicitEulerMod , only : implicitEuler ! implicit Euler solver for a single state

!
! !PUBLIC TYPES:
implicit none
Expand All @@ -54,9 +58,22 @@ module SoilHydrologyMod
private :: QflxH2osfcSurf ! Compute qflx_h2osfc_surf
private :: QflxH2osfcDrain ! Compute qflx_h2osfc_drain

type, extends(flux_type) :: drainPond_type
real(r8) :: smoothScale
real(r8) :: drainMax
contains
procedure :: getFlux => drainPondFlux
end type drainPond_type

!-----------------------------------------------------------------------
real(r8), private :: baseflow_scalar = 1.e-2_r8

! !ALTERNATIVE NUMERICAL SOLUTIONS
integer, parameter :: ixExplicitEuler=1001 ! constrained Explicit Euler solution with operator splitting (original)
integer, parameter :: ixImplicitEuler=1002 ! implicit Euler solution with operator splitting
integer, parameter :: ixAnalytical=1003 ! analytical solution with operator splitting
integer :: ixSolution=ixExplicitEuler

Copy link
Member Author

Choose a reason for hiding this comment

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

Comment from @martynpclark Sep 13, 2017

We should talk about how this is implemented. In principle it is possible to have multiple different numerical solutions for every single flux. This would provide scope to mix-and-match different numerical solutions in different parts of the code (sort of the way that CLM is developed now, although CLM does not have multiple numerical options for each flux). This would be incredibly cumbersome and likely confusing. At the other end of the spectrum we use the same numerical options everywhere in the code. This could be too restrictive.

I initially included these alternative options as (1) an example of how we can increase the accuracy of the operator splitting constrained explicit Euler solution without any increase in cost, expecting that there may be cases where people desire to retain operator splitting; and (2) begin the process of separating the physics from the numerical solution.

In my opinion the alternative numerical options for a single flux (as implemented here) may not be a representative use case.

Copy link
Member Author

Choose a reason for hiding this comment

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

Comment from @billsacks Sep 14, 2017

I'd also like to discuss this point further. And I'd be interested to see a more representative use case at some point.

One other concern I have is: If each flux has multiple possible parameterization options, each with multiple numerical solution options, the code is going to get very complex. (Much of the work I was doing earlier this week was to try to play around with ways to reduce that complexity, but there's no way to remove it entirely.)

Originally I had imagined that, for the most part, the choice of numerical solution method would just appear at high level in the code, not down at the level of individual fluxes. That would reduce this problem of a combinatorial explosion of complexity. I think that's what you're saying, too.

character(len=*), parameter, private :: sourcefile = &
__FILE__

Expand Down Expand Up @@ -478,6 +495,9 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, &
real(r8) :: dtime ! land model time step (sec)
real(r8) :: frac_infclust ! fraction of submerged area that is connected
real(r8) :: k_wet ! linear reservoir coefficient for h2osfc
real(r8) :: dRate ! effective linear reservoir coefficient for h2osfc (account for fraction of submerged area that is connected)
real(r8) :: activePond ! ponded water above threshold = h2osfc - h2osfc_thresh
real(r8) :: h2osfc1 ! h2osfc at the end of the time step, given qflx_h2osfc_surf only

character(len=*), parameter :: subname = 'QflxH2osfcSurf'
!-----------------------------------------------------------------------
Expand All @@ -503,11 +523,36 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, &

! limit runoff to value of storage above S(pc)
if(h2osfc(c) > h2osfc_thresh(c) .and. h2osfcflag/=0) then

! spatially variable k_wet
k_wet=1.0e-4_r8 * sin((rpi/180.) * topo_slope(c))
qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c))
dRate=k_wet*frac_infclust

! ponding above threshold
activePond = h2osfc(c) - h2osfc_thresh(c)

! switch between different numerical solutions
select case(ixSolution)

! constrained Explicit Euler solution with operator splitting (original)
case(ixExplicitEuler)
qflx_h2osfc_surf(c)=min(dRate*activePond,(h2osfc(c) - h2osfc_thresh(c))/dtime)

! implicit Euler solution with operator splitting
case(ixImplicitEuler)
qflx_h2osfc_surf(c) = dRate*activePond/(1._r8 + dRate*dtime)

! analytical solution with operator splitting
case(ixAnalytical)
h2osfc1 = activePond*exp(-dRate*dtime) + h2osfc_thresh(c)
qflx_h2osfc_surf(c) = (h2osfc1 - h2osfc(c))/dtime

! check case identfied
case default
call endrun(subname // ':: the numerical solution must be specified!')

end select ! (switch between different numerical solutions)

qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime)
else
qflx_h2osfc_surf(c)= 0._r8
endif
Expand Down Expand Up @@ -547,8 +592,16 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
logical , intent(inout) :: truncate_h2osfc_to_zero( bounds%begc: ) ! whether h2osfc should be truncated to 0 to correct for roundoff errors, in order to maintain bit-for-bit the same answers as the old code
!
! !LOCAL VARIABLES:
integer :: fc, c
real(r8) :: dtime ! land model time step (sec)
type(drainPond_type) :: drainPond_inst
integer :: fc, c
real(r8) :: dtime ! land model time step (sec)
real(r8) :: h2osfc1 ! h2osfc at the end of the time step, given qflx_h2osfc_drain only
real(r8) :: drainMax ! maximum drainage rate of ponded water
real(r8) :: const ! constant in analytical integral
real(r8) :: uFunc ! analytical integral
real(r8), parameter :: smoothScale=0.05_r8 ! smoothing scale
integer :: err ! error code
character(len=128) :: message ! error message

character(len=*), parameter :: subname = 'QflxH2osfcDrain'
!-----------------------------------------------------------------------
Expand All @@ -568,17 +621,63 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
qflx_h2osfc_drain(c) = h2osfc(c)/dtime
truncate_h2osfc_to_zero(c) = .true.
else
qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax(c),h2osfc(c)/dtime)

! define maximum drainage
! NOTE: check fraction
drainMax = frac_h2osfc(c)*qinmax(c)

! switch between different numerical solutions
select case(ixSolution)

! constrained Explicit Euler solution with operator splitting (original)
case(ixExplicitEuler)
qflx_h2osfc_drain(c)=min(drainMax,h2osfc(c)/dtime)

! implicit Euler solution with operator splitting
case(ixImplicitEuler)
! The next two lines look like extra overhead, but in a real implementation
! (1) smoothScale could probably be local to the drainPond type; (2)
! drainMax would probably also be local to the drainPond type, and if not,
! it would be an array that is just set once outside a loop, rather than
! being set for each loop iteration.
drainPond_inst%smoothScale = smoothScale
drainPond_inst%drainMax = drainMax
call implicitEuler(drainPond_inst, dtime, h2osfc(c), h2osfc1, err, message)
if(err/=0) call endrun(subname // ':: '//trim(message))
call drainPond_inst%getFlux(h2osfc1, qflx_h2osfc_drain(c))

! analytical solution with operator splitting
case(ixAnalytical)
const = 1._r8 - 1._r8/(1._r8 - exp(-h2osfc(c)/smoothScale))
uFunc = -1._r8/(const*exp(drainMax*dtime/smoothScale) - 1._r8)
h2osfc1 = -log(1._r8 - uFunc)*smoothScale
qflx_h2osfc_drain(c)= (h2osfc1 - h2osfc(c))/dtime
end select

if(h2osfcflag==0) then
! ensure no h2osfc
qflx_h2osfc_drain(c)= max(0._r8,h2osfc(c)/dtime)
end if
truncate_h2osfc_to_zero(c) = .false.
end if
end do

end subroutine QflxH2osfcDrain

subroutine drainPondFlux(this, x, f, dfdx)
! dummy variables
class(drainPond_type), intent(in) :: this
real(r8), intent(in) :: x ! storage
real(r8), intent(out) :: f ! drainage flux
real(r8), intent(out),optional :: dfdx ! derivative
! local variables
real(r8) :: arg ! temporary argument
! procedure starts here
arg = exp(-x/this%smoothScale)
f = -this%drainMax*(1._r8 - arg)
if(present(dfdx)) dfdx = -this%drainMax*arg/this%smoothScale
end subroutine drainPondFlux


!-----------------------------------------------------------------------
subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
Expand Down
25 changes: 25 additions & 0 deletions src/utils/computeFluxMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
module computeFluxMod

! provides interface for the flux modules

use shr_kind_mod , only : r8 => shr_kind_r8
implicit none

type, abstract :: flux_type
Copy link
Member Author

Choose a reason for hiding this comment

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

Comment from @martynpclark Sep 8, 2017

This is pedantic, I know, though can we call this computeFlux_type (or something similar)? The reason is that it is easy to confuse flux_type with WaterFlux_type and EnergyFlux_type, and folks could assume that the type only includes fluxes (whereas, in fact, it will include all variables required to compute fluxes).

Copy link
Member Author

Choose a reason for hiding this comment

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

Comment from @billsacks Sep 11, 2017

Fixed - thanks. Not pedantic at all - this is a very good point.

contains
procedure(getFlux_interface), deferred :: getFlux
end type flux_type

abstract interface
subroutine getFlux_interface(this,x,f,dfdx)
use shr_kind_mod , only : r8 => shr_kind_r8
import :: flux_type
implicit none
class(flux_type), intent(in) :: this
real(r8), intent(in) :: x ! state
real(r8), intent(out) :: f ! f(x)
real(r8), intent(out), optional :: dfdx ! derivative
end subroutine getFlux_interface
end interface

end module computeFluxMod
49 changes: 49 additions & 0 deletions src/utils/implicitEulerMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
module implicitEulerMod

use shr_kind_mod , only : r8 => shr_kind_r8
use computeFluxMod , only : flux_type
implicit none

contains

subroutine implicitEuler(flux_inst, dt, xInit, xNew, err, message)
! dummy variables
class(flux_type), intent(in) :: flux_inst
real(r8), intent(in) :: dt ! time step
real(r8), intent(in) :: xInit ! initial state
real(r8), intent(out) :: xNew ! updated state
integer,intent(out) :: err ! error code
character(*),intent(out) :: message ! error message
! local variables
integer :: iter ! iteration index
integer,parameter :: maxiter=20 ! maximum number of iterations
real(r8), parameter :: xTol=1.e-8_r8 ! convergence tolerance
real(r8) :: flux ! flux
real(r8) :: dfdx ! derivative
real(r8) :: xRes ! residual error
real(r8) :: delX ! state update
! initialiuze routine
err=0; message='implicitEuler/'

! initialize
xNew = xInit

! iterate
do iter=1,maxiter
! get the flux and derivative
call flux_inst%getFlux(xNew, flux, dfdx)
! get the residual and the update
xRes = xNew - (xInit + dt*flux)
delX = -xRes/(1._r8 - dt*dfdx)
! update and check for convergence
xNew = xNew + delX
if(abs(delX) < xTol) exit
! check for non-convergence
if(iter==maxiter)then
message=trim(message)//'the implicit Euler solution did not converge!'
err=20; return
endif
end do
end subroutine implicitEuler

end module implicitEulerMod