diff --git a/lessons/python/yet_another_oop.ipynb b/lessons/python/yet_another_oop.ipynb
index 470ef49..6ed6db2 100644
--- a/lessons/python/yet_another_oop.ipynb
+++ b/lessons/python/yet_another_oop.ipynb
@@ -507,7 +507,9 @@
"id": "4d4b08d0-7027-4fae-909e-4fe8d43c2113",
"metadata": {},
"source": [
- "Just kidding - here's a template, but feel free to draw your own owl :)"
+ "Just kidding - here's a template, but feel free to draw your own owl :)\n",
+ "\n",
+ "Note: where I've written pass, you will probably want to return something or save a new class attribute (self.whatever), which means you can remove the pass statement."
]
},
{
@@ -535,19 +537,27 @@
" self.temperature = initial_temperature\n",
"\n",
" # We probably want to calculate dz --> try using np.gradient\n",
+ " # Note: I didn't realize this when I originally wrote the class, but we won't actually need dz\n",
+ " # Want to know why? Check out the documentation for np.gradient -> https://numpy.org/doc/stable/reference/generated/numpy.gradient.html\n",
+ " # (What does the second argument to np.gradient do?)\n",
"\n",
" # We'll need attributes for rho, c, and k\n",
" # You could also store info about boundary conditions at this point, if you want\n",
+ " self.rho = 917\n",
+ " self.c = 2093\n",
+ " self.k = 2.1\n",
+ "\n",
+ " # Let's store boundary conditions too\n",
"\n",
" # Maybe we should go ahead and calculate diffusivity right away?\n",
"\n",
" # Let's keep track of the elapsed time\n",
"\n",
- " def calc_diffusivity(self):\n",
+ " def calc_diffusivity(self) -> float:\n",
" \"\"\"From above, kappa = k / (rho * c).\"\"\"\n",
" pass\n",
"\n",
- " def calc_heat_flux(self):\n",
+ " def calc_heat_flux(self) -> np.ndarray:\n",
" \"\"\"The heat flux is -kappa * dT / dz.\"\"\"\n",
"\n",
" # How should we calculate the difference in temperature with depth? (hint: see dz, above)\n",
@@ -622,13 +632,19 @@
"source": [
"model = SimpleGlacier(z, T0)\n",
"\n",
- "dt = 60 * 60 * 24\n",
+ "dt = 60 * 60 * 24 * 20\n",
"\n",
- "for i in range(365):\n",
+ "for i in range(30000):\n",
" model.run_one_step(dt)\n",
- " print(f\"{model.time_elapsed / 31556926} years.\")\n",
"\n",
- "model.plot()"
+ " if i % 1000 == 0:\n",
+ " print(f\"{model.time_elapsed / 31556926} years.\")\n",
+ "\n",
+ "plt.plot(model.temperature, model.z)\n",
+ "plt.xlabel(r\"Temperature ($^\\circ$ C)\")\n",
+ "plt.ylabel(\"Depth (m)\")\n",
+ "plt.gca().invert_yaxis() # This forces matplotlib to invert the y-axis\n",
+ "plt.show()"
]
},
{
@@ -636,17 +652,145 @@
"id": "7ae3c5ca-a9eb-4481-961f-1cc0d07911a5",
"metadata": {},
"source": [
- "Once you've done that, we could think about how to add a better boundary condition at the bed. A typical geothermal heat flux in Greenland is $0.04$ W m$^{-2}$ - but don't forget to also divide by $(\\rho * c)$ to keep the units intact! Maybe we should also let the user pass their own geothermal heat flux to our model. At this point, I bet you can figure out how to do that."
+ "Once you've done that, we could think about how to add a better boundary condition at the bed. Let's add a geothermal heat flux of $0.01$ W m$^{-2}$. Maybe we should also let the user pass their own geothermal heat flux to our model. At this point, I bet you can figure out how to do that."
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "2ebbb16c-f568-4358-992f-4bfe38817afc",
+ "id": "dccba120-feb5-4ae6-b726-99146df0700f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Initialize your model with a new bottom boundary condition"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4cc607e6-10bb-4dd5-9994-b3076be2a093",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dt = 60 * 60 * 24 * 20\n",
+ "\n",
+ "for i in range(30000):\n",
+ " model.run_one_step(dt)\n",
+ "\n",
+ " if i % 1000 == 0:\n",
+ " print(f\"{model.time_elapsed / 31556926} years.\")\n",
+ "\n",
+ "plt.plot(model.temperature, model.z)\n",
+ "plt.xlabel(r\"Temperature ($^\\circ$ C)\")\n",
+ "plt.ylabel(\"Depth (m)\")\n",
+ "plt.gca().invert_yaxis() # This forces matplotlib to invert the y-axis\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "038ac06d-9307-4717-987a-5dece0a6b682",
"metadata": {},
+ "source": [
+ "### Spoilers below!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "038be991-e91b-488f-b423-f9ad2d833d8b",
+ "metadata": {},
+ "source": [
+ "This is one example of how you could choose to write this class."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a74ed31d-6be6-40cb-a018-89830b53b923",
+ "metadata": {
+ "jupyter": {
+ "source_hidden": true
+ }
+ },
"outputs": [],
"source": [
- "# Plot results here!"
+ "class SimpleGlacier:\n",
+ " \"\"\"Models the temperature profile with a glacier.\n",
+ "\n",
+ " This model is based off of:\n",
+ " The Physics of Glaciers (Cuffey and Paterson, 2010).\n",
+ " Lecture notes from Andy Aschwanden (McCarthy school, summer 2012).\n",
+ "\n",
+ " Attributes:\n",
+ " z: an array of z-coordinates\n",
+ " temperature: an array representing the temperature profile with depth\n",
+ " \"\"\"\n",
+ "\n",
+ " def __init__(self, z: np.ndarray, initial_temperature: np.ndarray):\n",
+ " \"\"\"Initialize the model with arrays of z-coordinates and initial temperature profile.\"\"\"\n",
+ " self.z = z\n",
+ " self.temperature = initial_temperature\n",
+ "\n",
+ " # We probably want to calculate dz --> try using np.gradient\n",
+ " self.dz = np.gradient(self.z)\n",
+ " # Note: I didn't realize this when I originally wrote the class, but we won't need dz\n",
+ " # Want to know why? Check out the documentation for np.gradient -> https://numpy.org/doc/stable/reference/generated/numpy.gradient.html\n",
+ "\n",
+ " # We'll need attributes for rho, c, and k\n",
+ " # You could also store info about boundary conditions at this point, if you want\n",
+ " self.rho = 917\n",
+ " self.c = 2093\n",
+ " self.k = 2.1\n",
+ "\n",
+ " # Let's store boundary conditions too\n",
+ " self.surface_temperature = -10.0\n",
+ " self.basal_heat_flux = 0.0\n",
+ "\n",
+ " # Maybe we should go ahead and calculate diffusivity right away?\n",
+ " self.kappa = self.calc_diffusivity()\n",
+ "\n",
+ " # Let's keep track of the elapsed time\n",
+ " self.time_elapsed = 0.0\n",
+ "\n",
+ " def calc_diffusivity(self) -> float:\n",
+ " \"\"\"From above, kappa = k / (rho * c).\"\"\"\n",
+ " return self.k / (self.rho * self.c)\n",
+ "\n",
+ " def calc_heat_flux(self) -> np.ndarray:\n",
+ " \"\"\"The heat flux is -kappa * dT / dz.\"\"\"\n",
+ "\n",
+ " # How should we calculate the difference in temperature with depth? (hint: see dz, above)\n",
+ " temperature_gradient = np.gradient(self.temperature, self.z)\n",
+ "\n",
+ " # Are dT and dz the same size? Are they the same size as z?\n",
+ " assert temperature_gradient.shape == self.dz.shape == self.z.shape\n",
+ "\n",
+ " # Don't forget to apply boundary conditions! The heat flux at the bed should be zero, for now.\n",
+ " temperature_gradient[-1] = self.basal_heat_flux\n",
+ "\n",
+ " return -self.kappa * temperature_gradient\n",
+ "\n",
+ " def calc_divergence(self) -> np.ndarray:\n",
+ " \"\"\"In 1D, divergence is just the derivative. yay!\"\"\"\n",
+ " heat_flux = self.calc_heat_flux()\n",
+ "\n",
+ " flux_divergence = np.gradient(heat_flux, self.z)\n",
+ "\n",
+ " return flux_divergence\n",
+ "\n",
+ " def run_one_step(self, dt: float):\n",
+ " \"\"\"Advance the model by one step of size dt.\"\"\"\n",
+ " flux_divergence = self.calc_divergence()\n",
+ "\n",
+ " # updated temperature = current temperature - the divergence * dt\n",
+ " updated_temperature = self.temperature - flux_divergence * dt\n",
+ "\n",
+ " # Don't forget to apply boundary conditions! The temperature at the surface is equal to a fixed value.\n",
+ " updated_temperature[0] = self.surface_temperature\n",
+ "\n",
+ " self.temperature = updated_temperature\n",
+ "\n",
+ " self.time_elapsed += dt"
]
},
{