Inaccurate documentation for iqinit = 4?

[piyueh]

On GeoClaw's website (https://www.clawpack.org/setrun_geoclaw.html#setrun-qinit), it says

iqinit : integer
...
... 4 = Perturbation to surface level

Maybe I misunderstand this description, but it seems the source code is simply taking the values as the water surface level eta, rather than the perturbation to eta. For example, in qinit_module.f90:

geoclaw/src/2d/shallow/qinit_module.f90

Lines 183 to 184 in 3431f65

	! if iqinit = 4, the z column corresponds to the definition of the
	! surface elevation eta. The depth is then set as q(i,j,1)=max(eta-b,0)

And also from here:

geoclaw/src/2d/shallow/qinit_module.f90

Lines 168 to 170 in 3431f65

	else if (qinit_type == 4) then
	q(1,i,j) = max(dq-aux(1,i,j),0.d0)
	endif

So I guess the documentation on GeoClaw's website is not accurate? (Actually, on another page here, the description of iqinit=4 in the last paragraph matches the source code.)

[mandli]

You are correct that this is simply "perturbing" a state at rest but GeoClaw currently has no other way to initialize the sea-surface without additional code modifications so there really is no alternative anyway EXCEPT when an instantaneous earthquake is specified. In that case I think it really does perturb the surface like you were thinking although I am not entirely certain, it was not really designed to be used simultaneously like that.

[rjleveque]

@piyueh: Yes, I think you're right that there are some inconsistencies in the documentation, and probably unintended behavior in the code in some cases. Thanks for pointing this out.

This ability to read a qinit file dates from early versions of GeoClaw where we always assumed sea_level == 0 and so any non-zero eta could also be viewed as the perturbation to 0. I see that it needs some updating now that users are not only allowed to set a different constant sea_level but also a variable initial eta using the set_eta_init function.

I think the lines you flagged in

geoclaw/src/2d/shallow/qinit_module.f90

Lines 168 to 170 in 3431f65

	else if (qinit_type == 4) then
	q(1,i,j) = max(dq-aux(1,i,j),0.d0)
	endif

should probably be replaced by something like:

     else if (qinit_type == 4) then 
        eta = aux(1,i,j) + q(1,i,j)   # the surface elevation that was set in qinit.f90
        new_eta = eta + dq   # add the perturbation
        q(1,i,j) = max(new_eta - aux(1,i,j), 0.d0)   # the new depth
    endif

Does that seem right?

We should think a bit about what the desired behavior is on land. If q(1,i,j) = 0 is set in qinit.f90 then eta is equal to the topo value aux(1,i,j). If this gets perturbed by a dq that is positive then the new water depth will be equal to dq. That's probably the sensible thing? But note that it could have the effect that some low areas onshore get turned into lakes that are separated from the offshore water by higher ground in between. Also trying to use this in conjunction with a force_dry_init array might not work as expected. I'm not sure what would be expected, since that's not a use case I envisioned, but we should probably add a warning to the documentation at least.

[piyueh]

Thanks for your replies! I was just trying to make sure I understand it correctly.

Another thing is that I just realized the iqinit on GeoClaw's website is not available anymore in the python object QinitData. I guess it might be renamed to qinit_type at some point in the past?

[mandli]

Yeah, trying to make things a bit more explicitly named.

[rjleveque]

Looking at the earlier lines in qinit_module.f90:

geoclaw/src/2d/shallow/qinit_module.f90

Lines 164 to 170 in 3431f65

	if (qinit_type < 4) then
	if (aux(1,i,j) <= sea_level) then
	q(qinit_type,i,j) = q(qinit_type,i,j) + dq
	endif
	else if (qinit_type == 4) then
	q(1,i,j) = max(dq-aux(1,i,j),0.d0)
	endif

I see that for qinit_type < 4 (i.e. when specifying a perturbation to h, hu, or hv) the perturbation is only apply to cells if the topo is below sea level, presumably meant to be only the wet cells. But now that var_eta_init is allowed this won't necessarily have the expected behavior. I think instead we could just check if q(1,i,j) > 0, i.e. if qinit.f90 initialized the cell as wet.

For consistency we could also do the same thing for qinit_type == 4, e.g. combine this with my earlier suggestion and rewrite this as:

                        if (q(1,i,j) > 0.d0) then
                            if (qinit_type < 4) then
                                q(qinit_type,i,j) = q(qinit_type,i,j) + dq
                            else if (qinit_type == 4) then
                                eta = aux(1,i,j) + q(1,i,j) + dq
                                q(1,i,j) = max(eta - aux(1,i,j), 0.d0)
                            endif
                        endif

And note that in this case any cell forced to be dry in qinit would remain dry after this perturbation is applied.

[rjleveque]

Oops -- thinking more about this, I think that we actually wanted qinit_type == 4 to mean the value of eta is specified, not a perturbation as in the other fields. If we make the change I suggested above, handling dq as a perturbation to eta but only where q(1,i,j) > 0 from qinit, then I think qinit_type==4 would give exactly the same behavior as qinit_type==1, i.e. you can get the same surface perturbation by instead perturbing the depth in wet cells.

So the whole point of introducing qinit_type==4 was to allow the user to specify a desired surface elevation everywhere directly, ignoring what was set in qinit (at least in cells covered by the qinit array).

So maybe the code is right and only the documentation needs tweaking.

[mandli]

That makes sense to me given what we generally do here. The subroutine name is a bit misleading but I think it is probably still better just to clarify the documentation.