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 answer key to yet another OOP lesson #143

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
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
164 changes: 154 additions & 10 deletions lessons/python/yet_another_oop.ipynb
Original file line number Diff line number Diff line change
@@ -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 <tt>pass</tt>, you will probably want to <tt>return</tt> something or save a new class attribute (<tt>self.whatever</tt>), which means you can remove the <tt>pass</tt> 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,31 +632,165 @@
"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()"
]
},
{
"cell_type": "markdown",
"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"
]
},
{