Draft: Small stepgrid change

[kbarnhart]

Another change based on D-Claw development: In D-Claw we have one variable in q that should not be set to zero when q(1)<drytol.

This suggestion is intended to retain geoclaw functionality, but not require a D-Claw specific routine that differs only by this line.

@rjleveque @mandli @dlgeorge

[mandli]

This is probably true for any field that is momentum related I suppose. Most of the other fields, unless they somehow have velocity issues, would not necessarily want to zero it out either. In the multilayer code we do the following:

      forall(i=1:mitot, j=1:mjtot, k=1:num_layers,
     &       q(3*(k-1)+1,i,j) / rho(k) < dry_tolerance(k))
        q(3*(k-1)+1,i,j) = max(q(3*(k-1)+1,i,j), 0.d0)
        q(3*(k-1)+2,i,j) = 0.d0
        q(3*(k-1)+3,i,j) = 0.d0
      end forall

Given the momentum equations are interspersed with the other fields, this also required another stepgrid.f.

In any case, modulo the bouss code, I think being explicit about this is probably fine. The alternative would be to call something that would do only this custom bit. I suppose this might be a recurring task that may be easier to implement that way, but that may be too complex a solution for this.

[kbarnhart]

@mandli -

I think it be worth having stepgrid call a custom bit here would make sense. If Geoclaw, D-Claw, multi-layer all have custom stepgrid that differs only in these lines, reducing duplication seems worth it to me.

e.g., for dclaw it would look like this:

subroutine fixdry(q, nvar, mitot, mjtot)
      use geoclaw_module, only dry_tolerance
      use digclaw_module, only: i_h,i_hu,i_hv,i_hm,i_pb,i_hchi

      ! Input parameters
      integer, intent(in) :: nvar,mitot,mjtot

      ! Output
      double precision, intent(inout) :: q(nvar,mitot,mjtot)

      forall(i=1:my, j=1:mx, q(i_h,i,j) < dry_tolerance)
        q(i_h,i,j) = max(q(i_h,i,j),0.d0)
        q(i_hu,i,j) = 0.d0
        q(i_hv,i,j) = 0.d0
        q(i_hm,i,j) = 0.d0
        q(i_pb,i,j) = 0.d0
        q(i_hchi,i,j) = 0.d0
        ! do not reset i_bdif
      end forall
end subroutine fixdry

And for geoclaw it would look like:

subroutine fixdry(q, nvar, mitot, mjtot)
      use geoclaw_module, only dry_tolerance

      ! Input parameters
      integer, intent(in) :: nvar,mitot,mjtot

      ! Output
      double precision, intent(inout) :: q(nvar,mitot,mjtot)

    forall(i=1:mitot, j=1:mjtot, q(1,i,j) < dry_tolerance)
        q(1,i,j) = max(q(1,i,j),0.d0)
        q(2:nvar,i,j) = 0.d0
      end forall

end subroutine fixdry

I can implement this in Geoclaw (new function, updates to the Makefile) if you think it would make sense. (suggestions for the subroutine name are welcome).

[mandli]

@rjleveque and I chatted about this and were recalling that this type of operation is done a couple of places so having a fixdry subroutine that could be called from stepgrid and b4step would be fine. My only thought is to change the signature so that we have:

 ! =====================================================
 subroutine fixdry(meqn, mbc, mx, my, q, maux, aux)
 ! =====================================================
 !
 ! Reset cells whose depth is below dry_tolerance to a consistent dry state:
 ! zero negative depths to zero and zero the momenta.  Operates over the full
 ! patch, ghost cells included.
 !
 ! This routine is an override point.  Applications and GeoClaw variants may
 ! shadow it by placing their own fixdry.f90 in the application SOURCES list
 ! (see the GeoClaw docs on replacing library routines).  Two constraints on
 ! any replacement:
 !   * It must stay a plain external subroutine, not a module procedure
 !   * It must be thread safe.  Both call sites run inside an !$OMP PARALLEL DO
 !     over grids

     use geoclaw_module, only: dry_tolerance

     implicit none

     ! Arguments
     integer, intent(in) :: meqn, mbc, mx, my, maux
     real(kind=8), intent(inout) :: q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc)
     real(kind=8), intent(in) :: aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc)

     ! Locals
     integer :: i, j

     do j = 1-mbc, my+mbc
         do i = 1-mbc, mx+mbc
             if (q(1,i,j) < dry_tolerance) then
                 q(1,i,j) = max(q(1,i,j), 0.d0)
                 q(2,i,j) = 0.d0
                 q(3,i,j) = 0.d0
             end if
         end do
     end do

 end subroutine fixdry

The big difference is the addition of aux in case that is needed.

A couple of things to perhaps note:

• Since this is a call per-patch it should carry very little overhead (< 1%). I think using do loops may be fine here and avoids a possible logical array creation (although I have never quite understood when this does and does not happen).
• There are other places where this happens but may not be so easy to plugin. I had claude find some places where it may be worth a thought:
  1. src/2d/shallow/upbnd.f:120-122 — val(:,icrse,jcrse) = 0.d0 zeroes every component including h, contradicting the comment above it ("Reset small h to zeros") and skipping the max(h,0) every other site uses. For D-Claw this destroys i_bdif, the one component Draft: Small stepgrid change #733 exists to preserve. multilayer/upbnd.f:101-115 is the same code and additionally imports the scalar geoclaw_module::dry_tolerance and tests val(1,...) = ρ₁h₁ against it — off by ρ₁.
  2. src/2d/shallow/update.f90:168-170 — fine→coarse averaging writes only iadd(1..3). Components ≥4 on a coarse cell keep pre-update values forever (stale, not zeroed — worse to debug).
  3. src/2d/shallow/coarsen.f:109-111 — writes only valbgc(1..3); components ≥4 are fed uninitialized to the Richardson estimator via prepbigstep.f:84. Only reachable when flag_richardson = .true., which is False in every in-tree setrun.py.
  4. filpatch.f90:368 (do n = 2, nvar) and filval.f90:267 (do ivar = 2,nvar) — treat all of 2:nvar as h-scaled momentum: divide by h, velocity-bound against neighbours, remultiply by h_fine. Correct for hm/hchi, meaningless for a pressure like pb.
  5. prepregstep.f:24 and prepbigstep.f:88 — pass mitot/mjtot (resp. nx/ny against a mi2tot-sized array) as mx/my to b4step2, so it addresses out of bounds. fixdry inherits this verbatim; it is where a -fcheck=bounds build will trip. Unreachable today (flag_richardson = False everywhere in tree).
  6. amr2.f90 setprob ordering — setprob() runs before set_geo/set_storm/set_multilayer on fresh start (:521 vs :524/:531/:532) but after on restart (:486 vs :470/:471), so the existing procedure-pointer hooks cannot reliably be assigned by an app. A no-op setprob_final() stub in COMMON_SOURCES, called at the end of both branches, fixes it in ~8 lines.

[kbarnhart]

@mandli @rjleveque - I've now updated this with Kyle's suggested function, put it into stepgrid, updated the Geoclaw Makefile, and tested it against both the Geoclaw and D-Claw bowl-slosh example.

I didn't think it made sense for me to replace it in other locations in Geoclaw, but if one of you wants to do that on this branch, you should have push access.

Let me know if there are any other improvements to make or alternative approaches to implementation to discuss.

[kbarnhart]

Working on some larger simulations and this (or something else I did) has seemed to make things slower.

This is mostly to say that I think a timing test on full probably makes sense to do for this one.

I can plan to do that (sometime this week at the latest, I hope) and report back.

[kbarnhart]

OK, probably was something else. On a small (but not super small) problem, I get the following timing.

259.063 s CPU time on current master, 259.480 s with this change.

Without this change:

============================== Timing Data ==============================

Integration Time (stepgrid + BC + overhead)
Level           Wall Time (seconds)    CPU Time (seconds)   Total Cell Updates
  1                    11.991                 22.318            0.463E+08
  2                    61.271                206.457            0.283E+09
total                  73.262                228.775            0.329E+09

All levels:
stepgrid               72.460                227.093    
BC/ghost cells          0.480                  1.349
Regridding             13.876                 20.281  
Output (valout)         1.337                  0.721  

Total time:            97.658                259.063  
Using  4 thread(s)

Note: The CPU times are summed over all threads.
      Total time includes more than the subroutines listed above
Note: timings are also recorded for each output step
      in the file timing.csv.

With this change

============================== Timing Data ==============================

Integration Time (stepgrid + BC + overhead)
Level           Wall Time (seconds)    CPU Time (seconds)   Total Cell Updates
  1                    12.009                 22.359            0.463E+08
  2                    61.310                206.685            0.283E+09
total                  73.319                229.044            0.329E+09

All levels:
stepgrid               72.512                227.349    
BC/ghost cells          0.484                  1.361
Regridding             13.901                 20.309  
Output (valout)         1.743                  0.725  

Total time:            98.357                259.480  
Using  4 thread(s)

Note: The CPU times are summed over all threads.
      Total time includes more than the subroutines listed above
Note: timings are also recorded for each output step
      in the file timing.csv.