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

Add documentation vor v1.6.1 to tutorials #95

Merged
merged 5 commits into from
Aug 16, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
396 changes: 396 additions & 0 deletions notebooks/10minutes/numeric_solvers.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,396 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "003a7438-9867-41cf-b221-ad5e54e26d10",
"metadata": {},
"source": [
"# 10 minutes to... understand difference between ice dynamic models"
]
},
{
"cell_type": "markdown",
"id": "3dc919c7-f3f2-4193-b22c-270d4054cd09",
"metadata": {},
"source": [
"In version 1.6, OGGM changed the default numeric solver to the **Semi-Implicit** model. In this notebook, we explore the main differences compared to the old default, the **Flux-Based** model."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d6b6d9cb-8320-4219-af12-763eace92772",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import time\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"from oggm import cfg, utils, workflow, graphics, tasks\n",
"from oggm.core.flowline import FluxBasedModel, SemiImplicitModel"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e2186deb-7f89-403e-98c1-5b93eabd4f5a",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Initialize OGGM and set up the default run parameters\n",
"cfg.initialize(logging_level='WARNING')\n",
"\n",
"# Define our test glacier (Baltoro)\n",
"rgi_ids = ['RGI60-14.06794']\n",
"\n",
"# load elevation band representation\n",
"cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_dynamic_solvers_elevation_bands', reset=True)\n",
"base_url_eb = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/elev_bands/W5E5/'\n",
"gdir_eb = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_base_url=base_url_eb)[0]\n",
"\n",
"# load centerline representation\n",
"cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_dynamic_solvers_centerliens', reset=True)\n",
"base_url_cl = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/centerlines/W5E5/'\n",
"gdir_cl = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_base_url=base_url_cl)[0]"
]
},
{
"cell_type": "markdown",
"id": "234d6445-db4c-4738-a060-47de49e2f90c",
"metadata": {},
"source": [
"## Flux-Based model is more flexible, but unstable"
]
},
{
"cell_type": "markdown",
"id": "f2d798e7-21f2-4ad8-9322-fcb53f13620a",
"metadata": {},
"source": [
"The big advantage of the Flux-Based model is that it works for all flowline representations (multiple flowlines and different bed shapes). See also the tutorial [\"elevation band\" and \"centerline\" flowlines](elevation_bands_vs_centerlines.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9815e1f4-b41f-46cb-b170-9996d7d3ad99",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# run Flux-Based with centerlines\n",
"tasks.run_random_climate(gdir_cl,\n",
" evolution_model=FluxBasedModel,\n",
" nyears=300,\n",
" y0=2000,\n",
" seed=0,\n",
" store_fl_diagnostics=True,\n",
" output_filesuffix='_flux_based')\n",
"\n",
"# plot result\n",
"with xr.open_dataset(gdir_cl.get_filepath('model_diagnostics', filesuffix='_flux_based')) as ds:\n",
" ds_trap = ds.load()\n",
"ds_trap.volume_m3.plot();"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "278f2d7a-bc0d-4b3f-99fb-ffdb5bea0231",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# run Flux-Based with elevation bands\n",
"start_time = time.time() # time it for later comparision\n",
"tasks.run_random_climate(gdir_eb,\n",
" evolution_model=FluxBasedModel,\n",
" nyears=300,\n",
" y0=2000,\n",
" seed=0,\n",
" store_fl_diagnostics=True,\n",
" output_filesuffix='_flux_based')\n",
"flux_based_time = time.time() - start_time\n",
"\n",
"# plot result\n",
"with xr.open_dataset(gdir_eb.get_filepath('model_diagnostics', filesuffix='_flux_based')) as ds:\n",
" ds_flux_eb = ds.load()\n",
"ds_flux_eb.volume_m3.plot();"
]
},
{
"cell_type": "markdown",
"id": "c8987c1b-3c40-4339-a96f-6fcfdb8e9961",
"metadata": {},
"source": [
"Whereas the Semi-Impicit model only works for single trapezoidal flowlines (elevation bands)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "316c1c8f-8b78-4a04-ad21-2e15775b3e71",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# run Semi-Implicit with centerlines raise an Error\n",
"tasks.run_random_climate(gdir_cl,\n",
" evolution_model=SemiImplicitModel,\n",
" y0=2000,\n",
" seed=0,\n",
" store_fl_diagnostics=True,\n",
" output_filesuffix='_semi_implicit')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "053aa403-1659-4c5c-bbb9-baa3edba4995",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# run Semi-Implicit with elevation bands\n",
"start_time = time.time() # time it for later comparision\n",
"tasks.run_random_climate(gdir_eb,\n",
" evolution_model=SemiImplicitModel,\n",
" nyears=300,\n",
" y0=2000,\n",
" seed=0,\n",
" store_fl_diagnostics=True,\n",
" output_filesuffix='_semi_implicit')\n",
"semi_implicit_time = time.time() - start_time\n",
"\n",
"# plot result\n",
"with xr.open_dataset(gdir_eb.get_filepath('model_diagnostics', filesuffix='_semi_implicit')) as ds:\n",
" ds_impl_eb = ds.load()\n",
"\n",
"ds_impl_eb.volume_m3.plot(label='SemiImplicitModel', lw=4)\n",
"ds_flux_eb.volume_m3.plot(label='FluxBasedModel')\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"id": "5b20b896-bbaa-4e1e-9a10-26e20ef7515f",
"metadata": {},
"source": [
"You see that for the elevation band flowlines, both produce similar results. The differences arise from numeric instabilities in the Flux-Based model (see next paragraph). You can redo the experiment with a glacier where these instabilities are not that severe (e.g. RGI60-11.00897 Hintereisferner) and you will see both models produce the same result."
]
},
{
"cell_type": "markdown",
"id": "7130bacd-9442-4303-afdf-f498949645f0",
"metadata": {
"tags": []
},
"source": [
"## Semi-Implicit model is faster and more stable, but less flexible"
]
},
{
"cell_type": "markdown",
"id": "90a2680e-324b-4830-9249-cbf094b518ec",
"metadata": {},
"source": [
"Even the Semi-Implicit model is not as flexible as the Flux-Based one, we see it is faster when comparing the computing time:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b247bdf0-7ba6-4fa5-b2e7-7f20a5cec333",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"print(f'Semi-Implicit time needed: {semi_implicit_time:.1f} s')\n",
"print(f'Flux-Based time needed: {flux_based_time:.1f} s')"
]
},
{
"cell_type": "markdown",
"id": "fea7afda-0f06-4947-8565-78821ba157e7",
"metadata": {},
"source": [
"For a single glacier, this speed-up is probably not that important, but when thinking about regional to global simulations it can save you a lot of time.\n",
"\n",
"One reason for the speed-up is that the Semi-Implicit model is numerically more stable and can take larger time steps without producing instabilities:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "308b247a-2da4-4184-9109-bbf7fb6829e8",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# open flowline diagnostics\n",
"f_impl = gdir_eb.get_filepath('fl_diagnostics', filesuffix='_semi_implicit')\n",
"f_flux = gdir_eb.get_filepath('fl_diagnostics', filesuffix='_flux_based')\n",
"with xr.open_dataset(f_impl, group=f'fl_0') as ds:\n",
" ds_fl_impl = ds.load()\n",
"with xr.open_dataset(f_flux, group=f'fl_0') as ds:\n",
" ds_fl_flux = ds.load()\n",
" \n",
"# compare velocities along flowline\n",
"year = 100\n",
"ds_fl_impl.sel(time=year).ice_velocity_myr.plot(label='SemiImplicitModel')\n",
"ds_fl_flux.sel(time=year).ice_velocity_myr.plot(label='FluxBasedModel')\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"id": "12675cac-31f1-4579-88fb-65f0d384c33b",
"metadata": {},
"source": [
"The instabilities are visible for the FluxBasedModel at around 30 km distance along the flowline. They can lead to very large velocities which reduce the maximum possible step size due to the cfl-criterion (see also in the [documentation](https://docs.oggm.org/en/latest/faq.html#ice-velocities-in-oggm-are-sometimes-noisy-or-unrealistic-how-so)).\n",
"\n",
"The computing speed up and, even more, the increased stability are the reasons why we switched to the SemiImplicitModel in OGGM v1.6.\n",
"\n",
"However, if you want to set the FluxBasedModel as the default you can do so with:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "60e58f51-2078-4673-90df-1dd06bb1567f",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cfg.PARAMS['evolution_model'] = 'FluxBased' # default is 'SemiImplicit'"
]
},
{
"cell_type": "markdown",
"id": "2ad59e71-865a-4598-a8e4-63767974b28e",
"metadata": {},
"source": [
"## Have 5 minutes more? The bed shape of the downstream line"
]
},
{
"cell_type": "markdown",
"id": "0b2447c9-f876-4c85-bc4d-0bb9355e6e43",
"metadata": {},
"source": [
"This paragraph deals with the downstream line, the initially ice-free part of the glacier. You can see it below as the red line connecting the end of the outline with the left border of the figure:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b19e38fd-cdff-47a7-9a27-eba764340090",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"graphics.plot_centerlines(gdir_cl,\n",
" use_flowlines=True,\n",
" add_downstream=True)"
]
},
{
"cell_type": "markdown",
"id": "b88982a2-1d05-4d77-9fc4-9ae246ac6f69",
"metadata": {},
"source": [
"In OGGM before v1.6, with the FluxBasedModel, the shape of this downstream line was defined by fitting a parabola to the valley walls. However, for the SemiImplicitModel we had to change the shape to a trapezoidal, even a parabola approximates a mountain valley arguably better. We checked the influence of this change on advancing glaciers and found negligibly small differences in the volume on a regional scale. There might be some differences in the area.\n",
"\n",
"By default, we use a trapezoidal bed shape for the downstream line:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2af788f7-6c92-47df-aa94-c7e6317a0b15",
"metadata": {},
"outputs": [],
"source": [
"fl_trap = gdir_eb.read_pickle('model_flowlines')\n",
"fl_trap[-1].is_trapezoid[fl_trap[-1].thick == 0]"
]
},
{
"cell_type": "markdown",
"id": "3af0bbfe-2530-4037-8354-a34924c8551b",
"metadata": {},
"source": [
"But if for any reason you decided to use the FluxBasedModel you also can switch back to a parabolic downstream line using `cfg.PARAMS['downstream_line_shape'] = 'parabola'`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a18eed87-d4b3-4b60-8f3c-0d6dfd1d1ec6",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# change the downstream line shape\n",
"cfg.PARAMS['downstream_line_shape'] = 'parabola' # default is 'trapezoidal'\n",
"\n",
"# IMPORTANT: need to call init_present_time_glacier to take effect\n",
"tasks.init_present_time_glacier(gdir_eb)\n",
"\n",
"fl_trap = gdir_eb.read_pickle('model_flowlines')\n",
"fl_trap[-1].is_trapezoid[fl_trap[-1].thick == 0]"
]
},
{
"cell_type": "markdown",
"id": "861c7f1e-4a04-4cd8-b90e-d8f694dc18e5",
"metadata": {},
"source": [
"# What's next?"
]
},
{
"cell_type": "markdown",
"id": "2f70547c-1646-4ce6-9951-b8bb34a2b30d",
"metadata": {},
"source": [
"- return to the [OGGM documentation](https://docs.oggm.org)\n",
"- back to the [table of contents](../welcome.ipynb)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading