dtopo file behavior if start time <= t0, and initial time step

[rjleveque]

If a dtopo file is specified with the dtopo t0 in the header less than or equal to t0 for the job run, then the change to topo may be applied before qinit is called, so that the initial ocean is flat and never changes if there are no subsequent dtopo changes specified (particularly if mt = 1 so only one dtopo time is specified). This may or may not be what is desired, but for instantaneous earthquake models it probably is not desired, instead you might want to see the initial surface displacement matching the specified dtopo. At any rate this should be clarified in the documentation.

There might also be a bug in this, since in some tests I did the initial surface was not entirely flat but had some small perturbations near the coast on coarser grid levels (but not on the finest level patches). This needs more exploration.

At any rate, for earthquake models it is best to specify t0 in the dtopo file to be larger than the t0 set in setrun.py for the simulation, e.g. to 1 second as we often do for "instantaneous uplift" models. And then setting e.g.

      rundata.clawdata.dt_initial = 0.5 

      rundata.dtopo_data.dt_max_dtopo = 0.5

insures that it takes small enough time steps initially on the coarsest grid to capture the initial deformation before it starts taking time steps governed by rundata.clawdata.cfl_desired, which can be several hundred seconds if the coarsest level is one or two degrees.

Also, if dt_max_dtopo <= dt_initial then it should be using this time step initially, one would think, but in fact it uses dt_initial for the first time step regardless of dt_max_dtopo and this should be fixed.

[rjleveque]

Update on this issue:

If an instantaneous deformation is specified at time t0 or (before) and if you set

rundata.qinit_data.variable_eta_init = True

then I think everything works ok and the full dtopo deformation appears in the initial conditions that are printed out at time t0.

However, if you set variable_eta_init = False (or don't set it at all in setrun.py, since the default is False), then it sometimes does not capture the full deformation on all grid levels (e.g. in some tests it seems to be captured on levels 1 and 2 but on level 3 the ocean is flat).

In retrospect this is not surprising, and it is always best to set variable_eta_init = True and use set_eta_init whenever you want an instantaneous deformation at time 0 (and also more generally for problems with coastal uplift or subsidence, especially for problems where you have fine grids near the coast in the region that experiences vertical land motion).

If the vertical motion dz is incorporated in the initial topography, then you do not want these cells to be filled with water up to the original sea_level (which is the GeoClaw default in qinit.f90) since the ocean near the coast subsides or rises along with the land, the the initial water level eta in these cells should be set to sea_level + dz where dz varies with location. This is exactly why variable_eta_init was introduced, although originally we were thinking of coastal regions where that do not get refined until after the earthquake happened (but before the tsunami arrives), in which case new coastal cells again have to be filled to sea_level as modified by the vertical motion dz as specified by the dtopo file. This is what the default version of $CLAW/geoclaw/src/2d/shallow/set_eta_init.f90 does.

Also note that if an instantaneous deformation is specified at time t0 (only one time frame in the dtopo file) then it should be ok to not set rundata.dtopo_data.dt_max_dtopo. The initial data should already include the deformation in the topo and in eta (if variable_eta_init = True) and so it should be possible to take any initial time step allowed by the CFL conditions.