Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lateral drag / Shape factors #429

Merged
merged 14 commits into from Mar 27, 2018
Merged

Lateral drag / Shape factors #429

merged 14 commits into from Mar 27, 2018

Conversation

@phigre
Copy link
Contributor

@phigre phigre commented Mar 8, 2018

  • Closes #423
  • Passes Inversion Tests for runs without shape factors
  • Tests for inversion with shape factors added/passed
  • Tests for shape-factor interpolation from Adhikari /Nye added and passed
  • Fully documented, including whats-new.rst for all changes (remove if this change is small enough to be undocumented)
    not yet to a sufficient amount

For first review (@fmaussion)

Philipp Gregor added 2 commits Mar 8, 2018
…nservation inversion.

Shape factors are used as given in Adhikari 2012.
A new cfg-param "use_shape_factor" is introduced.
Needs further tests and eventually adaption in forward model.
Changed one comment/TODO in inversion routine.
Copy link
Member

@fmaussion fmaussion left a comment

Thanks so much! This already looks very good. I have a few minor comments only.

For the rest (and for future reference) it would be super nice to have some plots pasted in this PR:

  • curves for various zetas and for both methods
  • maybe some application test at HEF?

This we could talk bout it or I'll do it in a subsequent PR

oggm/cfg.py Outdated
@@ -274,6 +274,7 @@ def initialize(file=None):
PARAMS['invert_with_sliding'] = cp.as_bool('invert_with_sliding')
_k = 'optimize_inversion_params'
PARAMS[_k] = cp.as_bool(_k)
PARAMS['use_shape_factor'] = cp.as_bool('use_shape_factor')
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

I think we should rename it to use_shape_factor_for_inversion. It is likely that we will want to separate between inversion and dynamical model's factors.

For testing between the factors it might be useful to have this option as a str: 'Huss' (or the original paper name), "Adhikari" or empty for no shape factor

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Yes, I already also thought about an option like this, but just did not spend more time on it. Will do so.

according to Huss and Farinotti (2012)
Discouraged to use, as this does not make any difference
between parabolic and rectangular beds
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

I'm not sure to understand: wouldn't it be better to return a factor of 1 when the bed is rectangular in this case?

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

I could not really find any information on the assumed bed form in the paper. Therefore I did not go into too much detail yet. But you're right, I will adapt it this way.

Adhikari (2012)
TODO: should we expand this here to also include
the factors suggested for sliding?
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

Not yet I would say

factors_parabolic = np.array([0.251, 0.448, 0.653, 0.748,
0.803, 0.839, 0.917])
f_parab = interp1d(tabular_zetas, factors_parabolic,
fill_value='extrapolate')
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

This entire block could live at the module scope instead

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Agreed

is_rectangular = is_rectangular.astype(bool)

# TODO: could check for division by 0, but at the moment
# this is covered by interpolation and clip, resulting in a factor of 1
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

It's just a problem if it throws numerical warnings. If you want you can put this statement in a catch as done here

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

In case is_rectangular is of type integer, these integers (0/1) are directly interpreted as indices, which is not what is intended here. In case of floats I have to admit that I did not check.

-------
"""
a0s_shaped = a0s / (shape_factor ** 3)
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

you can safely overwrite a0s here

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

But then I cannot exploit the precomputed a0s in subsequent loops

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Sorry, my mistake. Did not realize I was inside the function.

@@ -161,6 +273,11 @@ def mass_conservation_inversion(gdir, glen_a=cfg.A, fs=0., write=True,
fd = 2. / (cfg.N+2) * glen_a
a3 = fs / fd

# Shape factor params
sf_func = shape_factor_adhikari
sf_tol = 1e-2 # TODO: better as params in cfg?
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

This is too specific it's OK to leave it here

# Start iteration for shape factor with guess of 1
i = 0
sf = np.ones(slope.shape)
sf_diff = np.ones(slope.shape)
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

i = 0 and sf_diff could be defined in the if block below

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Absolutely right


shape_factors = np.ones(widths.shape)

# TODO: higher order interpolation? (e.g. via InterpolatedUnivariateSpline)
Copy link
Member

@fmaussion fmaussion Mar 8, 2018

Maybe? Depends on how the curve looks like

Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Already looks quite ok, i would say?

grafik

I can provide the notebook to this figure on demand

@fmaussion
Copy link
Member

@fmaussion fmaussion commented Mar 8, 2018

Ping @alexjarosch for info


# TODO: how to handle zeta > 10? at the moment extrapolation
# Table 1 from Adhikari (2012) and corresponding interpolation functions
_ADHIKARI_TABLE_ZETAS = np.array([0.5, 1, 2, 3, 4, 5, 10])
Copy link
Contributor Author

@phigre phigre Mar 8, 2018

Move this to top of file?

Copy link
Member

@fmaussion fmaussion Mar 20, 2018

I'm not sure there is a pep rule for this, but yes I think it's slightly better

Philipp Gregor added 2 commits Mar 19, 2018
…ion of Nye-shape factor.

Updated test to match change in shape factor calculation
oggm/params.cfg Outdated
@@ -189,6 +189,9 @@ optimize_inversion_params = False
# If false, tell OGGM which should be used
inversion_glen_a = 2.4e-24
inversion_fs = 0.
# Do you want to use shape factors to account for lateral drag?
# Allowed is empty, "Adhikari", "Nye" equal to "Adhikari" or "Huss"
Copy link
Member

@fmaussion fmaussion Mar 20, 2018

"Nye" (equivalent to "Adhikari")

@@ -161,6 +287,16 @@ def mass_conservation_inversion(gdir, glen_a=cfg.A, fs=0., write=True,
fd = 2. / (cfg.N+2) * glen_a
a3 = fs / fd

# Shape factor params
sf_func = None
if cfg.PARAMS['use_shape_factor_for_inversion'] == 'Adhikari' or \
Copy link
Member

@fmaussion fmaussion Mar 20, 2018

create a new variable, like use_shape_factor = cfg.PARAMS['use_shape_factor_for_inversion'] and use it throughout the function to shorten the code a little bit

@fmaussion
Copy link
Member

@fmaussion fmaussion commented Mar 20, 2018

@phigre should we merge this as is and add the tests at a later stage or are you still working on this?

Philipp Gregor added 2 commits Mar 20, 2018
…shape factor in FluxbasedModel.step

Currently, the shape factor methods of the inversion module are used for the calculation of shape factors.
…he shape factors for sections with height = 0 are no longer np.nan but are set to 1. (so shape factor is actually not influencing this section)
@fmaussion
Copy link
Member

@fmaussion fmaussion commented Mar 20, 2018

Ah, and don't forget to acknowledge your work in whats-new.rst !

sf = sf_func(fl.widths_m, thick, fl.is_rectangular)
# TODO: avoid unnecessary new array initialization by adding
# sf_stag to self._stags
sf_stag = np.zeros(thick_stag.shape)
Copy link
Contributor Author

@phigre phigre Mar 20, 2018

Staggered shape factors need to be computed to match with thick_stag and u_stag.
I decided to use shape factors from the unstaggered grid and combine upstream and downstream shape factors to obtain staggered values.
Directly computing shape factors on the staggered grid would make an additional width_stag necessary and the question arises what is_rectangular_stag would look like.

@@ -21,6 +21,7 @@
from oggm import entity_task
import oggm.core.massbalance as mbmods
from oggm.core.centerlines import Centerline, line_order
from oggm.core.inversion import shape_factor_adhikari, shape_factor_huss
Copy link
Contributor Author

@phigre phigre Mar 20, 2018

I referenced the shape factor functions from the inversion module. IMHO, having them in a separate dedicated module and referencing this module from inversion and flowline would make more sense in regards of code structure.

Copy link
Member

@fmaussion fmaussion Mar 20, 2018

Yes, probably. It could be added to utils for now (and the test moved to test_utils for now, even if the utils module is definitely overcrowded and could need a split-up one day. What do you think?

Copy link
Contributor Author

@phigre phigre Mar 20, 2018

e.g. move to oggm.utils?

Philipp Gregor added 5 commits Mar 20, 2018
Introduced use_sf in flowline and inversion to shorten lines and avoid multiple calls to cfg.PARAMS[use_shape_...']
Using a proxy variable `use_sf` instead of `cfg.PARAMS['use_shape_...']`. This variable is set to None, if the acquired key is not found in PARAMS.
@phigre
Copy link
Contributor Author

@phigre phigre commented Mar 21, 2018

There are still some issues like the storage of sf_stag in self._stag and how to handle not set fl.is_rectangular (yes, that is possible), which I didn't want to decide on my own.

Otherwise ready to merge now?
A simple test for creating a synthetic glacier with shape factors in the FluxBasedModel and trying to recover the bed with mass-conservation inversion and the adequate shape factors for inversion worked out, so at least we ca assume that either both are right or broken ;)

…s. Moved shape factor function choice to init of class. Added tests with shape factors for inversion and FluxBasedModel doing a roundtrip.
Copy link
Contributor Author

@phigre phigre left a comment

Stiil a lot of testing necessary and some easy performance possible, but getting better

d = np.zeros(nx-1)
e = np.zeros(nx)
self._stags.append((a, b, c, d, e))
d = np.ones(nx+1) # shape factor default is 1
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

could set to 1.0 in case shape factors are not wanted to make the array-array mutliplication an array-scalar mutliplication afterwards. May give some perfomance advantages (especially if numpy could identify multiplication by 1 as skippable)

rhogh = (RHO*G*slope_stag*sf_stag)**N
u_stag = (thick_stag**(N+1)) * self._fd * rhogh + \
rhogh = (RHO*G*slope_stag)**N
u_stag = (thick_stag**(N+1)) * self._fd * rhogh * sf_stag**N + \
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

shape factor as we use it only makes sense for the non-sliding part

@@ -1377,6 +1377,96 @@ def test_inversion_parabolic(self):
plt.plot(bed_shape_gl[:-3])
plt.show()

def test_inversion_parabolic_sf_adhikari(self):
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

Tests for each case would probably be much better with one single callable function and write wrappers, which are actually called for testing and set the shape factor mode before calling the function.

@@ -1377,6 +1377,96 @@ def test_inversion_parabolic(self):
plt.plot(bed_shape_gl[:-3])
plt.show()

def test_inversion_parabolic_sf_adhikari(self):
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

Much better to define one single function and then write wrappers for different shape factor modes, which are actually called for testing and do nothing else than setting shape factor cfg and call the function doing the work. Would avoid a lot of ugly copied code

d = np.zeros(nx-1)
e = np.zeros(nx)
self._stags.append((a, b, c, d, e))
d = np.ones(nx+1) # shape factor default is 1
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

Set to 1.0 as a scalar if scalararray is faster than arrayarray in numpy. Would be cool, if numpy could skip multiplication by 1 automatically

oggm/utils.py Outdated
@@ -3161,10 +3161,10 @@ def shape_factor_huss(widths, heights, is_rectangular):

# Ensure bool (for masking)
is_rectangular = is_rectangular.astype(bool)
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

Avoid typecast necessary? Did not look into the implementation and cannot say anything about perfomance impacts

oggm/utils.py Outdated
@@ -3161,10 +3161,10 @@ def shape_factor_huss(widths, heights, is_rectangular):

# Ensure bool (for masking)
is_rectangular = is_rectangular.astype(bool)

shape_factors = np.ones(widths.shape)
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

always repeating initialization could be substituted by storing shape factor array in the flowmodel and passing the reference to the shape factor functions

oggm/utils.py Outdated
shape_factors = np.ones(widths.shape)
shape_factors[~is_rectangular] = \
(widths / (2 * heights + widths))[~is_rectangular]
shape_factors[heights <= 0.] = 1.
Copy link
Contributor Author

@phigre phigre Mar 22, 2018

Need to avoid nans for areas without ice

@fmaussion
Copy link
Member

@fmaussion fmaussion commented Mar 27, 2018

This is a brilliant piece of work, thanks a lot! I will merge this now, we can always go back to some details later.

Experience shows that the global applications will require some tuning, but I'm eager to test this out. Note that this probably should become the default some day.

@fmaussion fmaussion merged commit 5f98d9a into OGGM:master Mar 27, 2018
3 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants