Ponded water/alternative solutions

[billsacks]

NOTE: This PR was originally created in the temporary CTSM repository by @martynpclark, Sep 6, 2017, where it was NCAR/clm-ctsm#12. I am moving over most of the comments from there: both my own and those of @martynpclark .

Original comment from @martynpclark Sep 6, 2017

NOTE: This PR has not been compiled/tested. It is included to illustrate concepts for alternative numerical solutions, for discussion at the CTSM meeting on 7 Sep 2017.

[billsacks]

Comment from @billsacks Sep 7, 2017

@martynpclark @dlawrenncar @swensosc @barlage @bandre-ucar @mvertens

Here is the pdf I referenced in today's meeting:

Modern.Fortran.in.Practice.-.section.5.3.pdf

Specifically, see the "Type­-Bound Procedures" implementation.

There is a more involved example that follows a similar pattern in the book "Scientific Software Design: The Object-Oriented Way" - see section 6.2.1. That book is available in the NCAR library (http://opencat.library.ucar.edu/cgi-bin/koha/opac-detail.pl?biblionumber=58255&query_desc=kw%2Cwrdl%3A%20scientific%20software%20design). (Disclaimer: I have only given that a quick skim, so may be missing something.)

After giving this more thought, while I like Martyn's proposed solution, I think I'd prefer something like the object-oriented one illustrated in the two above references. (Or the solution that @bandre-ucar proposed, which feels like it would end up looking somewhat similar in the end.) A big advantage is that it is more generally applicable: you don't need to have a containing subroutine that sets all of the input parameters. But a lot of my feeling is based on the problems I have found with using internal subroutines like this, where you access variables from the enclosing routine:

1. It's hard to trace the data flow in the containing routine. If I have something like this simple example:

   foo = 17
   call my_internal_routine()
   bar = 2 * foo

   where my_internal_routine accesses foo: If I go to modify this routine, a common thing I do is try to trace the data flow through the routine, by searching for where foo is set and used. It's very easy to miss the fact that foo is accessed by my_internal_routine here, and thus it is not safe to move the setting of foo to after that call.

   The same issue would apply to the example in this PR.

2. The above problem gets even worse if any variables are set by the internal routine. And because these variables aren't passed as arguments, you can't enforce an intent(in) expectation.

3. It's impossible to unit test internal routines in isolation.

[billsacks]

Comment from @billsacks Sep 7, 2017

@martynpclark I think you said that drainMax is supposed to be inputRate (or vice versa); am I remembering right? (So that I can show some alternatives correctly.)

[billsacks]

Comment from @martynpclark Sep 7, 2017

@billsacks -- yes, inputRate should be drainMax -- sorry about that.

[billsacks]

Comment from @martynpclark Sep 7, 2017

inputRate should be drainMax

[billsacks]

Comment from @billsacks Sep 7, 2017

Is this a general pattern - that, to get the final flux with implicit Euler, we do a final function call like this? If so, can this final call be incorporated into the implicitEuler routine itself?

[billsacks]

Comment from @martynpclark Sep 7, 2017

Yes. The final flux should be a function of the state at the end of the time step.

[billsacks]

Comment from @billsacks Sep 11, 2017

Addressed in NCAR/clm-ctsm#14 in NCAR/clm-ctsm@315c14b

Update: This is now #192 , commit b872d21 . The comment referenced in that commit is the one right above: https://github.com/ESCOMP/ctsm/pull/191/files#r159001381

[billsacks]

Comment from @billsacks Sep 7, 2017

In today's meeting, we discussed the fact that this implementation will have performance problems - and particularly will not allow vectorization - because the routine is called separately for each point, rather than operating on all points. We discussed refactoring this so that the flux routine operates on arrays rather than scalars.

However, I'm not sure how to fit that in to the integration routine, which has an outer loop over iterations, with a convergence criteria for each point. Does anyone know how this is done elsewhere in CESM (or other models) - so that you can have vectorizable code in the context of a solver like this?

[billsacks]

Comment from @martynpclark Sep 8, 2017

@billsacks -- there are iterative solutions in other parts of CLM (e.g., see line 754 in CanopyFluxesMod.F90. Are there problems with vectorization in these parts of the code also?

[billsacks]

Comment from @billsacks Sep 8, 2017

@martynpclark Thanks for that pointer. That example handles this well, from a vectorization perspective, by passing a temporary filter to each science routine. Each time through the iteration loop, the filter is recreated via this code at the end of the iteration loop:

         ! Test for convergence

         itlef = itlef+1
         if (itlef > itmin) then
            do f = 1, fn
               p = filterp(f)
               dele(p) = abs(efe(p)-efeb(p))
               efeb(p) = efe(p)
               det(p)  = max(del(p),del2(p))
            end do
            fnold = fn
            fn = 0
            do f = 1, fnold
               p = filterp(f)
               if (.not. (det(p) < dtmin .and. dele(p) < dlemin)) then
                  fn = fn + 1
                  filterp(fn) = p
               end if
            end do
         end if

That solution makes sense to me in a case like CanopyFluxes, where there is a lot of work that happens in the iteration loop. In the example in this PR, the overhead of recreating the filter every time through the iteration loop may not pay for itself, since there isn't much work done in the science routine. That said, it may be worth going with the filter solution anyway.

Another alternative (simpler, and good enough or better in terms of performance when there isn't much work done in the flux calculation routine) would be: Have a logical array passed to the flux calculation routine that says whether we have converged at each point. Then the loop in the science routine would look like:

  do fn = 1, num_f
     n = filter_f(fn)
     if (not_converged(fn)) then
        ! do calculations
     end if
  end do

My understanding is that loops with conditionals can sometimes be vectorized, though maybe not efficiently. (I think the machine may compute the result for all elements, then throw out the ones where the conditional was false.)

Before settling on a solution, we should probably do some timing tests of a few different approaches:

1. Scalar approach like you have here
2. Array approach with filter
3. Array approach with conditional
4. Others?

[billsacks]

Comment from @billsacks Sep 8, 2017

@martynpclark I was going to spend part of today illustrating the object-oriented approach I discussed yesterday. Would it be more useful if I:

1. First helped rewrite your current implementation to use arrays, then showed a side-by-side comparison with an object-oriented approach also using arrays?

or

2. Showed an object-oriented approach using scalars - then we can treat the conversion to array-based as a separate step?

I think you can see the differences between the internal subroutine approach and OO approach about equally in scalar or array form, so I'm happy to do this in either order if one seems more useful to you.

[billsacks]

Comment from @martynpclark Sep 8, 2017

@billsacks -- thanks. Just FYSA, the example I presented was deliberately simple, which was done in order to make it easy to illustrate the implicit Euler approach. In practice most implementations will be much more complex (i.e., a large amount of code in the flux+derivative subroutine). We still may have some cases where the function call in the implicit Euler is trivial.

[billsacks]

Comment from @billsacks Sep 8, 2017

In practice most implementations will be much more complex (i.e., a large amount of code in the flux+derivative subroutine)

Okay, then going with a filter-based approach like the one in CanopyFluxes is probably a good way to go, though I still think we should do some timings before settling on this completely.

[billsacks]

Comment from @martynpclark Sep 8, 2017

@billsacks -- we would never actually use the implicit Euler in this case, since we are implementing operator splitting and the analytical solution is available. This will become more clear once we move to solutions for the complete state equation.

[billsacks]

Comment from @martynpclark Sep 8, 2017

Agreed that we should look at the timings.

[billsacks]

Comment from @martynpclark Sep 8, 2017

@billsacks -- in terms of the path forward:
1. First helped rewrite your current implementation to use arrays, then showed a side-by-side comparison with an object-oriented approach also using arrays?
or
2 Showed an object-oriented approach using scalars - then we can treat the conversion to array-based as a separate step?

Either way is fine, whatever is easiest for you. I've read the book chapter that you sent and I think I understand the approach.

[billsacks]

Comment from @martynpclark Sep 8, 2017

@billsacks -- your code block

  do fn = 1, num_f
     n = filter_f(fn)
     if (not_converged(fn)) then
        ! do calculations
     end if
  end do

seems much easier to understand. Is there a reason why the other (longer) code block is preferable?

[billsacks]

Comment from @billsacks Sep 8, 2017

@martynpclark I should first point out that the key part of the longer code block is just the second half:

            fnold = fn
            fn = 0
            do f = 1, fnold
               p = filterp(f)
               if (.not. (det(p) < dtmin .and. dele(p) < dlemin)) then
                  fn = fn + 1
                  filterp(fn) = p
               end if
            end do

(The other part is needed in either case, I think.)

This rebuilding of the filter replaces the not_converged logical.

The reason this is preferable (in principle) is the same reason that we have filters at all in CLM: It's more efficient in general - and especially in a context of vectorization - to loop over some subset of indices rather than looping over all indices (or a larger subset of indices) with a conditional inside the loop.

That said, I don't think anyone has done timing tests of the actual impact of this with recent machines / compilers.

[billsacks]

I'm closing this PR that was really mainly opened for discussion purposes.

[billsacks]

I have deleted this branch. Its entire history is contained in the pondedWater/alternativeSolutions_oo branch, which currently resides on my fork. The head of the pondedWater/alternativeSolutions branch was commit 0041a92.