Change Timing for Gauge Updating

[mandli]

This PR moves the update_gauge call to after a step is taken. It also puts a condition on an initial time abs(t0 - time) < TOLERANCE so that the initial time is still recorded at the initial time.

[mandli]

This seems to work although there are two questions:

• Does the impact the min_time_output functionality for gauges?
• Should we not record the gauge secondarily after the first step?

My guess is no to both, but it would be good to have a double check.

This should be merged before #642.

[rjleveque]

This modification can make a significant change in the gauge values printed, e.g. for the chile2010 example comparing this to the master branch version commit 271bea3 that it is based on, there are changes in the 2nd significant digit in many values, particularly in hu and hv. E.g. master gives these lines as part of gauge32412.txt:

   02  0.7967691E+04  0.4425462E+04  0.5695225E-02 -0.7138274E-02 -0.4718117E-04
   03  0.8159613E+04  0.4382718E+04  0.7401617E-02 -0.8644394E-02 -0.6137179E-04
   03  0.8197998E+04  0.4382718E+04  0.8354446E-02 -0.1245157E-01 -0.7849258E-04

while the modified code gives

   02  0.7967691E+04  0.4425462E+04  0.5695225E-02 -0.6724048E-02 -0.4665426E-04
   02  0.8159613E+04  0.4425461E+04  0.1331391E-01 -0.1412212E-01 -0.9894767E-04
   03  0.8197998E+04  0.4382718E+04  0.8354446E-02 -0.1245157E-01 -0.7849258E-04

The third line is identical, but there are significant differences in the first line, and in the middle line shown, hu changes by 85% from 0.74e-2 to 1.33e-2 and even the level is different.

I think I know why, and I'm remembering now why we (including @mjberger) previously decided it wasn't so easy to move the update_gauges call to after the call to stepgrid in advanc. The original code calls it after ghost cells have been set on the patch but before the time step is taken, and so interpolating to gauge locations near the edge of patch works fine. The new version does this interpolation after the interior values on the grid have been updated by a time step, but the ghost cells are still those from the previous time step, so the interpolation in space is using data at two different times near the edge of the patch. In the example above, at the first time the gauge was presumably near the edge of a patch while at the third time shown it was interior to the patch and so interpolation gave the same thing in both codes.

I think when we discussed this before, we determined it would be a lot more work to update the ghost cells after the step in order for this interpolation to work properly and was not worth doing every time step just so that the gauges would be output at the final time of the simulation. An alternative would be to do something special at the end of the simulation to fill ghost cells and output gauges properly, but that also didn't seem worth doing since the gauge value is rarely actually needed at the final time. If it is, the user should just increase tfinal a bit to capture the gauge at the desired time properly.

But we could reconsider this. I notice that the routine valout has a call to bound to set boundary values for each patch if we are writing out in a binary format rather than ascii, because we do a binary dump that includes the ghost cells and so we want these to be at the same time as the interior points, since the resulting ghost cell values might actually be used (e.g. if the user wants to plot vorticity, which can then be computed by applying central differencing to all the interior grid cells to compute $v_x$ and $u_y$ properly). So we could do a final call to update_gauges after then final time frame is output, if we changed valout to always set the ghost cell values even for ascii output. (And assuming the user asks for the full output at the final time -- in some cases we only want fgmax and/or fgout output and set clawdata.num_output_times = 0 to suppress all valout output.)

The reason the middle line in the example above is so different is that in the new code the gauge value is output immediately after taking a time step, while in the original code there is a regridding that takes place in between the end of one time step and printing the gauge at the beginning of the next, and a recomputation of mbestsrc determining which grid is the best source of data for this gauge. One could debate which value is more accurate, the value just computed on level 2 or the value interpolated to the new level 3 grid that gets introduced around this gauge, but at any rate they will generally be very different values, as seen here.

[mandli]

So I decided to go through the exercise of adding time to setaux and to my surprise, it did not fix the problem. I am currently trying to trace down where exactly we would need to add code in bound and likely filrecur to fill in the ghost cells with aux values, but this seems to be part of the problem. Just to be clear, this was working with the switch in recording time as in b4step the ghost cell aux values were also being set.

I will continue to get a better handle again as to exactly what is going on but there are two questions that this raises:

• Do we want to put time into setaux? This was not nearly as hard as I thought but I was not able to guarantee that whatever time is there is correct.
• There are a number of places where we copy aux values from the old grids (mostly for topography). This was modified to handle moving topography, but not general time-dependent aux values. While I do not think the problem here is directly impacted by this, we may need to reconsider this point.
• This may impact how @kbarnhart implements temporal aux arrays.
• Do we want to change how dtopo is done so that it can be handled directly within setaux rather than having an entirely separate apparatus?

In any case, I think this PR can be closed, but we should reopen a new issue to maintain the discussion and make PRs from there.