 Okay, let's go ahead and get started. Again, welcome everybody to the webinar for November 2023. My name is Greg Tucker. I'm the director of CSDMS and I am absolutely thrilled to introduce Daniel Shapiro from the Polar Science Center. Daniel is a researcher there. He has a PhD in applied mathematics from University of Washington in 2017. And among other things, he is the lead developer of the software package, which is a tool for modeling the flow of glaciers and ice sheets and solving inverse problems. And that's what we're gonna hear about today. So without further ado, I will turn the floor over to you, Daniel. Okay, thanks. So as Greg said, I'll be talking in part about the software package called IcePak that I develop. So IcePak is in turn built on another library called FireDrake for finding element analysis. And I'll talk about that some too. It's possible that the bit I say about FireDrake is gonna be more interesting to this crowd than what I say about IcePak. So I was supposed to give this talk in April, which was good timing because the birth of our first child was due around early in May. And then this little Joker thought that it would be funny to be born six weeks early, which was fine. But it also meant I had to delay this talk. So thanks to Greg and Lynn for being flexible and to the nurses of UW Northwest Hospital for keeping him alive. And now he has grown to be, fat as can be 90th percent of for height and weight. So premature birth aside, everything's totally great. And yeah, you can see the personality already coming out there. So yeah, that's what I've really been working on lately, not a lot of math, mostly just this creature. Okay, so the motivation for pretty much all of this is that at least in glaciology, traditionally for a long time, doing any numerical modeling has usually been the province of experts in numerical modeling as such. And what I was trying to set out to do with all of this is make a tool that a broader audience could use, whether or not they were experts. So just a brief overview, I'm gonna talk about the finite element method and fire drake. I'll give a little demo of it. I'll give a brief overview of the physics of glaciers going into some detail, but not a ton. I'll talk about IcePak, some of the design decisions that I made along the way. And then I'll talk about something that I've been working on lately that I'm pretty excited about. Okay, so the finite element method, very briefly. Simulating Glacier flow means solving some big nasty partial differential equations numerically, but a solution of a partial differential equation can have infinitely many degrees of freedom. And the finite element method is a way to approximate something that has infinitely many degrees of freedom with only finitely many. So we express the solution U of the differential equation we're looking for as a finite linear combination of a set of basis functions, phi i, that we choose in some judicious way. And then we need to come up with some way to choose the coefficients. And the finite element method is one way of doing that. There are also collocation methods and spectral methods and all sorts of other things, the details of which mostly won't matter for the rest of this. So this is really embarrassing that I see that Wolfgang Geingerth is in the room because I used deal two for a while and I don't use it any longer. So the kernel of developing the finite element method is, okay, all right, well, at least we have a sense of humor about this. So I started writing all of this by hand in C++. This is a couple of language standards ago, but very often the complaint about this is that it can be error prone to write. There can be a lot of boilerplate and most especially it can be very difficult to customize by the end user. And that was kind of one of the big things that I was trying to address when writing Icepack. I will say, so I wrote a ton of code completely by hand in Fortran, wrote my own sparse matrix data structures and such like, and then I started using, I started trying out deal two, whoops. And I was able to get everything rewritten in about six months. And I would say that deal two has the most fantastic documentation of any software project I've ever used. And yeah, I was able to use it and actually contribute some code to the library in six months and Wolfgang and Timo and all the others who work on it have been really amazing about it. So briefly, if you're not familiar with it, some pseudo code for what you might do to assemble the stiffness matrix and the finite element method involves a really gigantic nest of loop. So for every cell of a mesh, for every quadrature point in that cell and for every pair of degrees of freedom or of basis functions that are non-zero over that cell, you assemble some big thing that involves a lot of complicated tensor algebra and pack it into a larger matrix. And here is some code that I used during my thesis that would actually do all of that. So if you know how the finite element method works, then this loop is pretty comprehensible. Again, it roughly corresponds to that pseudo code. I'm gonna make this a little bigger. So again, for every cell of a mesh and for every quadrature point in that cell and for every pair of degrees of freedom in that cell. So this, with a lot of work, I managed to get down to roughly one laptop screen of code, but there are a lot of things like this and it became, it was easy enough for me to develop with practice, but I had sort of set a higher bar for myself, which is that I wanted PhD students to be able to easily modify some of the guts of what goes in here. So just again, a couple of helpful graphics. This is from the CommSol website. This is a triangular mesh of a heat sink. So the principle describing, say the temperature field throughout this entire heat sink would require infinitely many numbers to do. So instead, we discretize the problem by breaking up the geometry and do a set of small pieces or triangles and define a set of basis functions. In this case, this is showing two basis functions defined on two of the vertices of that triangular mesh. These two happen to not overlap, but they would overlap if we were looking at the basis function defined, say on this vertex instead. And finally, by writing functions as a linear combination of those basis functions, we can do things like solve the heat equation. So this shows in colors a heat field on that. Okay, so where FireDrake kicks off is, instead of writing a lot of the low-level C or C++ code, could we instead use a high-level symbolic description language to describe the problem that we wanna solve and have someone write a compiler, basically, which will automatically generate all the low-level C and C++ code to fill matrices and vectors. The result of this originally in 2004 was called the Unified Form Language and a software project called Phoenix. I am using an offshoot of Phoenix called FireDrake, but there have been others. So DuneFem also uses the Unified Form Language. There's another software project called DeVito, which is mostly aimed at finding difference as opposed to finding element methods, but it too uses this kind of higher-level symbolic description language and code generation. And most of what they're interested in with DeVito is seismic wave propagation and seismic conversion problems. So Phoenix, for example, they have kind of a heavy bias towards problems in mathematical biology. So one of their, I think the thing they like to brag about is doing simulations of the heart and the aorta and many arteries, which is super cool. And FireDrake has been oriented more at doing geophysics problems. So just as a brief demonstration, I've written down here the variational form of the diffusion equation for a field U, which states that for every test function V, this integral relation holds. And without belaboring a ton of the details, what I wanted to write down here is the code that expresses this variational form. And what's so nice about this is that the code looks a lot like the math. So we replace dot products that we might write in math notation with a call to a function called inner and we replace some upside down triangles with the word grad, but more or less there's a one-to-one correspondence between one expression and the other. And that's really nice. It saves you a lot of typing. And on top of that, I find like it eliminates a lot of the low level fiddling with the details of what goes on a lot of inner loops. So I wanted to give a brief demonstration of this. So on an example that's a little more sophisticated than just the diffusion equation. So all of this is using FireDrake. And what I wanted to do was compute the drainage area, which is a quantity that you often want to use in the stream power incision model, which describes the evolution of landscapes under brainwash and other sources of erosion. So I took a digital elevation model from the state of Oregon for this region called the Sullivan Creek. And I picked that because one of the papers on hillslope diffusion uses this exact region. So this is a contour map of the height field of the area we're looking at. And what I wanted to do was express the upslope drainage area as a differential equation rather than using one of the flow routing models that people commonly use. So a lot of this is, and I can go back and explain as much of this as people want to see. But really the important thing I want people to get out of this is how closely the code that you write corresponds to the math. So the idealized model that we're trying to solve is that the divergence of this area field that we're looking for times the slope is equal to one. Now we have to turn that into a variational form. So we multiply by some test function phi and integrate. And here the things I want to draw your attention to are this first expression, the dot product of the upslope area times the slope with the gradient of the test function that corresponds exactly to this bit of code here. And some boundary terms are again here. And what's nice about these is it functions in a lot of ways like a computer algebra system. So we can define like store subsets of these expressions and variables that we define piecemeal rather than have to define the entire statement of the problem all at once. I did a bit of extra things to add a small bit of diffusion to regularize things using the symmetric interior penalty discontinuous Galerkin method, which sure is a bit of a mouthful, but the end result is a fairly nice plot, which is maybe a little mushier than what you get out of a flow routing method, but still pretty decent. And the plot should show up in just a sec. Here we go. So these darker areas are where you have the greatest amount of upslope drainage area that are flowing into there, which is one of the source terms in the string power incision model. So again, really the thing that's so nice about these is how much all of the math that you write down corresponds exactly to the code and how succinct it is to express what problem it is that you wanna solve without having to worry about a lot. Well, without having to manage a lot of the low level details yourself, which has both good parts and bad, as I'll allude to. So again, the code looks a lot like the math. One of the other features that's really nice about this is that it makes doing automatic differentiation very simple and it's baked into the language itself. So if you have a nonlinear problem and you wanna have the linearization of it derived for you automatically, you don't have to do any of that yourself. This is possible in other languages using there's operator overloading techniques that you can do in C++ that will also do the same thing. So it's not like this is impossible to have otherwise but it is really nice to have it built in. And on top of that, it's very easy to modify, which was one of my big criteria. Now there are drawbacks to this. The first is a loss of low level control. So when I was writing everything in C++, I could do a lot of hand optimization of the loops I was writing in order to make sure that they ran as fast as possible. There is a decent amount of that that you lose when you go to using domain specific languages, which can be good if the people writing the compiler are really focused on that and I find often they are or it can be bad if there are low level optimizations that you wanna do but cannot. Which is a thing that probably I'm just gonna go ahead and guess that Wolfgang and I are gonna argue about that at the end of this talk. The second con and this is really a big one that I have no answer for is that you have a much more complex software stack. So with deal two, well I don't know how much of this is a function of their wonderful documentation but at least it's pretty easy to get your head around most of the data structures that you really need to work with. Whereas with FireDrake, I don't understand the low level guts of it. There's basically a compiler built in there and it's unlikely that I ever will. So this comes with trade-offs and as I hope I'll convince you some of them are worth making and a lot of that is gonna be application dependent. So FireDrake, which is the software package I was showing you just a second ago. It started out as an experiment on Phoenix. David Hamm wanted to see how much of the internals he could rewrite to make more of the code generated as opposed to defined beforehand. And at a certain point, they realized that it was just going to become its own project. It wasn't going to get merged back into the mainline Phoenix. So it became a friendly fork and the people who developed FireDrake and Phoenix all worked together on maintaining the outer level symbolic description language but they effectively, you could think of it as having two different compilers for the same program language. So some of the really big advantages of FireDrake over many other tools in the same kind of space are that it has all the weird hipster elements. So if you wanna do, for example, our gear is finding elements for plate bending, it has that. One of the other kind of cool things that people have been doing with this are new and clever and energy and extra fee conserving ways to solve the shallow water equations. So if that's interesting to you, then you can ask me about it. The other thing that has been a really big deal, especially for me is tensor product meshes. So when you have a 2D footprint mesh that you wanna lift up into 3D and solve problems on this layered domain, it is really nice for that and they have that baked in from the start. And there are a lot of problems in the geophysical sciences that are naturally defined on tensor product meshes. So that to me is really just a kind of killer out for FireDrake. Okay, so before I talk about IcePak, which is a software package for ice flow modeling, I should probably talk about Glacier physics a little bit. So I always love starting with this quote, this is from JD Forbes, who I guess at the time they would have called him a natural philosopher and it's from the transactions of the Royal Society of London. And he was one of the first people to express the idea that Glacier's flow by viscous deformation. And he kind of came to this conclusion from some observations that he was making in the Swiss Alps. And one of the things that he realized that he wrote down in his notebooks is that if that's true, then it means that you're using the same mathematics to describe the flow of glaciers as to describe the flow of lava. So totally different kinds of things, but it's all really the same math underlying. And he says there's something pleasing to the imagination and the unexpected analogies presented by a torrent of fiery lava in the icy stream of a glacier. And I love that quote, so I always throw that in wherever I can. So ice flow is like a viscous fluid, but what's unusual about it is that it's not a Newtonian fluid, it's sheer thinning. So if you double the applied stress, then you more than double the resulting strain rate response. And this is part of what makes solving the momentum balance equations of glaciers difficult is that you now have a nonlinear differential equation. So I love showing a couple observations of this just because I think it really nicely illustrates the point. So this is the time lapse from Landsat imagery going back to the 70s of Malaspina Glacier in Alaska. So Malaspina is a Piedmont glacier. So all of the ice is flowing down from this very narrow channel to the top and it spreads out into this very wide plain. There are actually forests that are growing on top of the ice at the terminus there, which is just wild to me. So that's kind of a long time series. This next one is from some more recent imagery. I think this is just from 2017 to 2021, but that does a nice job showing the viscous flow out of that channel and into the... Okay, one more, I swear to God, almost done. This is the Conta Glacier in Alaska and this is some time lapse imagery from in situ photographs that I think does a nice job showing off like, oh yeah, it really does kind of look like a fluid if you take a lot of pictures over a really long time. Okay, fun's over, time for math. So there are really two big parts of the physics of glacier flow, which is just mass and momentum balance. So the mass balance part is kind of the easier of the two to describe. So first, how does the thickness of the glacier change? So the first is flux divergence. So if there's more ice flowing into some control volume that you put on the ground and there is flowing out, then the thickness will tend to increase and the sources are accumulation. So that can come from snowfall or meteoric ice that slowly compacts, snowfalls that can have a density anywhere between around 50 to 350 kilograms per meter cubed and under the weight of more and more snow, it'll densify to about 900 kilograms per meter cubed. But you can also have accumulation on floating ice shelves that comes from seawater that freezes underneath it. So that's common in Antarctica. So like the Filtrenorani ice shelf a significant fraction of the entire thicknesses frozen seawater and melting. So for say a mountain glacier, that's all from the sun at the surface for floating ice shelves in Antarctica and for tidewater glaciers in Greenland that can also come from ocean melt. So you can have accumulation and melt both at the top and bottom depending on what region you're looking at. So this, at least if we're only thinking about the thickness is just a linear hyperbolic equation which is relatively easier to deal with the real bad boy in the room is looking at momentum balance. So what I'm describing here is a simplified glacier flow model where we've taken the full viscous stress balance. We've assumed that it's a thin film flow and done a bit of perturbation theory to remove a couple of terms. And also that most of the flow is by horizontal extension rather than vertical shear. So I'm just going to describe this one model. There are loads of different perturbative models. I'm just sticking to this one for brevity. So breaking down all the terms here, first is the balance of membrane stresses. So M is the membrane stress which is basically a little ranked two tensor that you get that is the result of what happens when you apply this thin film approximation to the full Stokes equations. This is balanced by drag at the glacier base. So if you have a glacier that's flowing on land, it's experiencing resistance as it flows over bedrock or sediments. And this kind of ties into some of the geology and hydrology of the region. And finally gravity, which is that ice flows downhill. The kind of real, the operative term here is just the gradient of the surface slope S. But that, you know, just these two equations isn't a complete description. I haven't told you really what the membrane stress tensor is. So what I've written down here is the famous Glenn flow law which states that the rate of strain or the symmetrized gradient of the ice velocity U is proportional to some power of the membrane stress tensor. So roughly the membrane stress tensor raised to the power M where N is about three. And there is a team N equals four. And we get into playful fights with each other of the N equals three and the N equals four people. Myself, I have no opinion because I just write the modeling tools. I only care that there is a fight and that I can help both sides. It's sort of like, you know, when there's a gold rush, self-pick axis, that's me. The next is the sliding law. So how do we relate the basal shear stress or the basal drag to the basal velocity? And once again, the simple model that people have liked to use is some kind of power law so that the basal velocity is proportional to the drag to power M where M is again greater than one. So there's some kind of shear thinning behavior happening at the base. This is actually, if you chase all the papers down, this is really a fascinating story that played out largely in the 80s and 90s. So for a while, people thought that the basal sliding exponent M had to be the same as the flow law exponent N because sliding occurred largely through processes of internal deformation within the ice. And then people started getting observations of glaciers in Antarctica, like the Cyple Coast ice streams, where they were flowing really fast, like upwards of a kilometer a year, and yet the surface slope was almost zero. And the reason that this was able to happen is because they learned all sorts of new things about subglacial geology, that you could have nearly plastic flow. So where this exponent M gets very, very large or even in some models, infinity, and then you have to do strange things. So really the upshot of all of this is that glacier flow is quite nonlinear because it's a shear thinning material. And the ways that it interacts with the substrate is also nonlinear as well. And just some more observations. This is from the NASA Scientific Visualization Studio. This is a log plot of the velocity of all of Antarctica. So you'll notice these purple regions, those are the really fast ones. This is the Filchner and Ronnie ice shelves, the Larson Sea ice shelf, the Ross ice shelf, and Pine Island and Fwates. They've been in the news a lot lately. So all the fastest flowing parts are where the ice is floating because there's no drag from the bed. And so they're essentially unconfined other than by internal resistant stresses. And so there are a lot of ice streams that feed these ice shelves sort of like tributaries feeding a larger river. And this, if you've read any of the news articles about, and I hate this term, the doomsday glacier, this is the weights, this is the doomsday glacier, and people call it that because they think that it could destabilize all the way back here to the ice streams that feed the Ross ice shelf. So one of the best ways that we have of validating these nonlinear flow laws of saying that we think the strain rate is proportional to the stress to some power is by looking at these areas of floating ice shelves and basically doing regression on the measurements there. Because when there's no basal drag then that kind of removes a lot of the complications and you can test these assumptions. Okay, so there are a couple of missing pieces. So I said that the strain rate is proportional to the stress to some power but what's the proportionality constant in there and what does it depend on? So well, first of all, there are some laboratory measurements that can tell you roughly what the temperature dependence is and it's an Arrhenius law. But there are loads of other factors that come in here too. So the first is impurities, which I just find really amusing. The atmosphere of the last glacial period was a lot dustier. It didn't rain as much. There wasn't as much precipitation so you get much less scavenging out of all the dust that gets kicked up and for example, the Sahara Desert. So much more of that dust reaches the ice sheets and mechanically it can cause like an enhancement factor of this flow law by a factor of two or three, which is really, really big for something that you definitely would not expect just going from first principles. And the final one is crystal orientation. So I've written this as sort of an isotropic kind of relation where the strain rate is proportional to the stress to some power but there can actually be a directional dependence due to the fact that eventually as you have enough sustained strain along the same direction over really, really long periods of time, you can get preferential alignment of the ice crystals that make them easier to deform in some directions as opposed to others. And coming up with good models for this is an area of active research. Okay, so now what's missing in the sliding law? The big ones are geology. So is the ice flowing over hard crystalline granite or is it flowing over glacial till or sediments and really kind of the holy grail problem is basal hydrology. So how much meltwater from the surface is making it to the bed? How much meltwater is being generated at the bed itself? Where is that meltwater routed afterwards? And how does it affect sliding? This is a very challenging problem. My understanding from talking to people who really study this is that the current forefront of research is getting a model of basal hydrology to couple to a model of glacier flow in such a way that they don't explode. So just keeping the thing from not blowing up is current research. And it's especially hard because we can get very few observations of basal hydrology. The final one that I think is, this is basically the holy grail problem of glaciology is coming up with some good mathematical model for iceberg calving. So what makes calving go? What are the relevant fields? What is the frequency and how large are the calving events? So I'll show you a little movie. This is a sped up time lapse image from Helheim Glacier in East Greenland. And you can see this block of ice breaking off and tipping over, which is kind of an important part of the whole process. So what's interesting is it's really the same physical process, but in Greenland, the events are frequent and small. And in Antarctica, the events are infrequent. So maybe every couple of decades, but gigantic. So periodically, you know, for example, you might've seen a few years ago the Larson Sea ice shelf calved off an iceberg the size of Manhattan Island. So just roughly some big questions that people are asking in the field. The main one that I'm kind of attacking here is, can you simulate how big the ice sheets are gonna be in 100 or 200 years under a changing climate? How much sea level rise can we expect? The one that I kind of wish we paid more attention to is better managing glaciers as a water resource. So for example, in Washington state where I live, there's a very large proportion of the agricultural land that's irrigated by water that comes from glacial melt in the summer. And glaciers cover only 2% of the land area of Washington state. So how are we going to deal with this in 50 or 100 years? And a lot of the people here at UW work on paleoclimate and paleoglaciology. So how can we use our understanding of ice flow to get a better idea of what we're seeing in the ice core record? Okay, so where do I come into this? So I started writing ice pack in 2013 and it was a C++ library and it was pretty much just for me. I got a postdoc to turn this into usable code and I was told to my face that no glaciologist is going to learn C++ which was a little harsh but empirically I found that to be true. And I know this differs by subfield within the earth sciences so I'm speaking mostly just for glaciology here. I was basically told in as many words that it had to have a Python interface. So I started taking all this deal two code that I'd written and grown quite attached to and then I started running into, you know what do I do to make the C++ memory model comprehensible to the Python and proliferation of shared pointers occurred and at some point I thought maybe I should just rewrite the whole thing using these fancy new tools with the main specific languages. So to fire Drake's credit I was able to rewrite a whole lot of thesis code in about six weeks over a winter vacation. Which I thought was pretty impressive just for them. So it wasn't long before a couple of people started using it and some of them even managed to do so without mentioning to me that they were trying. So they were able to do it without getting any help from me which I thought was a pretty good sign that I might have something on my hands that could be more useful to the community. And I got an NSF grant to develop this and I became a PI at the lab here. So just a few kind of observations that sort of informed a lot of my thinking when I was trying to develop this. So I got a degree in applied math and I of course had taken a whole lot of classes and numerical methods. And largely the things that we were taught were how to solve certain common families and differential equations really fast and how to solve them on big parallel computers. But a lot of the problems that I encountered in glaciology and in looking at other areas of the geosciences is they're only in 2D. So the need for really big HPC resources isn't quite there. Now I know this isn't true for all fields. If you're doing atmospheric dynamics or ocean dynamics you really do have honest to God 3D problems and you need the big iron. But many of the problems that we're interested in are 2D but they're just weird. So when I show mathematicians the full stream power incision model usually they're very, very weirded out by it. They don't know how to classify it. It's extremely confusing to them and they look at the dynamics of what it does and usually by then they're jaws on the floor. So this kind of, when you view this as sort of the normal then a lot of what I had been taught in an applied math degree wasn't preparing me for the reality of what I would actually be doing which is finding ways to discretize relatively small but incredibly odd problems. At the last part, which is more of a human factors thing is that a lot of the people in my grad school cohort who were doing earth sciences rather than applied math they really hated doing numerical modeling. They found it frustrating. And we ended up having to deal with an emotional problem of people don't like doing this. They don't find it fun. How do we deliver something to them that makes you feel like there's raw power sluicing through your blood? So they hated doing it and they had other options. So of course I have to show that I do in fact get outside sometimes this is from the Juno ice field in Southeast Alaska and people would rather be doing fieldwork than messing around with a lot of computer stuff. So when you factor in, people have other options they don't necessarily have to use these big sophisticated flow models they can do other things. So if you want people to use them then you're gonna have to make some changes. Okay, so the kind of operating hypothesis I was working off of was that a lot of students have a harder time translating very odd or unusual physics into code in the first place than they do making the code run fast enough. So a lot of what I learned in my degree was kind of under this assumption that performance and speed is the rate limiting factor for domain scientists getting work done. And that is true a lot of the time but it also is a completely different set of difficulties that they're encountering the remainder of the time. And it's those difficulties that I was kind of interested in addressing. Okay, so now we get into weird awful things that scientists really don't like talking about because they're hard to quantify and we have to talk about people. So what are, if you take an average domain scientists in field X and ask them what is their mental model of the physics or the dynamics of the thing they're interested in? How can you design simulation tools that best correspond to the mental model that they have to how they think of it in their head? And this isn't, I don't wanna knock people doing performance optimization or high performance computing too hard. But one of the, you can measure the performance of code. It actually is quite difficult but measuring something like this is basically impossible. We're talking about squishy human factors here. So it's sort of an adjustment to start thinking about this. And I think I can at least answer this question for glaciology, maybe not for other sub-disciplines but by and large, if you ask a glaciologist to describe roughly what's going on with glacier flow or how they think about it in their head, they'd say there's, so there's one equation which is the momentum balance that sets the velocity. It has no time derivatives in it because they're basically minuscule. So it sets the velocity instantaneously as a function of the thickness and other forcing data. And next is a prognostic equation which is the mass balance equation that updates the thickness and time. So my idea is could we write a software tool where the application programming interface for it more or less matched exactly the domain scientists mental model for how it worked. And I've written a bit of pseudocode here but as I'll show you in a second and I should have been running this demo earlier. I'm just going to start it now and it should be done in a second. Right, okay. So could you design a software tool where the application programming interface matches closely as possible people's mental models of it. And this is roughly what I was trying to do. And the other thing that was really important to me is that people using it should have fairly fine grain control of what happens in the internals of the simulation. So the kind of abstraction that I present to them is that I give them a solver for the mass of a metam balance equations and it's up to them to put it in a loop that does the time stepping which is sort of more responsibility on them but it also means that they're able to intervene in this loop and debug it or put in whatever other analysis or visualization bits they want in the middle of it. However they see fit. And this is a design choice that has trade offs but many of the other tools that are aimed at doing similar things the only interface that you're presented is solve the entire couple of equations of glacier flow for this length of time with this many time steps and you don't get to do anything in the middle of the simulation you have a fairly coarse grain control which is good if the entire thing works really well and you trust that it's doing internally what you want but many people find that this feels like a black box so for me it was really important that you should have a fairly fine grain level of control over it. So there's a second part of this which I had alluded to earlier which is not only do I want to have something that matches well people's metal model of what the physics is that they have fairly fine grain control but also there should be a pretty wide flexibility in your ability as an end user to make changes to what's going on. The reason being that we don't understand all of glacier physics and the example I like to use here is some glacial hydrology. So for example, a lot of users are going to want to mess around with the sliding law. Maybe they just want to change the exponent which is pretty simple to do from outside you change one floating point number to another but what if instead we wanted to couple the sliding law to a model of some glacial hydrology well now we have to add a whole other field into the physics that has its own dynamics as well and the couples into them amount of balance. So designing software that is amenable to that kind of degree of modifiability is quite hard. I have a very simple and stupid solution for this which you can do in Python through the power of keyword arguments. I feel like I'm posing this great grand problem and then I have the absolute silliest and dumbest solution for it. So roughly speaking, all of the models that I've developed an ice pack the main kind of function you want to call is this diagnostic solve. And really what that's doing internally is it's summing up many terms. So there's say a viscous power dissipation, gravitational power, frictional power and so forth making a sum of all of those many components and solving a differential equation or minimizing some functional if you're familiar with variational calculus and the modifiability comes in because every single part of those terms those are also attributes of this model that the user can modify. So I've kind of roughly sketched out here what the viscous power dissipation does or looks like but the user can substitute in their own choices for how to compute each term of that model. So for example, if they want to change how the viscosity is computed to include say a fabric model, they can do it. If they wanted to give their own function for how you compute the frictional power dissipation they can do the same thing. So again, just to summarize that the specification of what physics you want to solve essentially boils down to a sum of many terms and I've given people access to you can change individual terms of those or more than one. So you can replace any one of them. And the reason that you can do that which again is really silly is that all of the arguments are passed by keyword and every single term of that model can see every field that is passed into the entire solve. So this is a bit of a code smell. It means that for example, the friction law can see fields that are related to the flow law. So that's not ideal from a software engineering point of view but it was a trade-off that I elected to make because it means that people can very easily replace any of the terms of that model. And I think that's been pretty nice. So, and I can show you some, a couple of papers where people have actually done that before. The other thing I wanted to mention is I did a decent amount of reading of the literature on human computer interaction and I'll bring up the Wikipedia page for it because it summarizes it really nicely. This comes from a paper by Green and Peter in 1996. They were looking at visual programming languages but it applies to a lot of other things. And they laid out these concepts for what makes it easy for people to grasp. In their case of visual programming language but I think it applies to a lot of other software. So there were a couple of the kind of criteria that they point out that I think are really helpful for thinking about. So first is how closely does the notation or does the application programming interface correspond to the problem? And I'd spoken about that already in terms of giving people an interface where what they can do is solve one step of the prognostic model or the diagnostic model. Progressive evaluation. So how easy is it to obtain feedback on an incomplete solution? Again, this is one of those grandiose sounding criteria that has a very silly solution which is, all of this is written in an interpreted programming language. And I showed you that I like using Jupyter notebooks for this but you can also just do all of it in an IPython console. And for example, if you're not comfortable with writing a whole loop that'll do the iterated diagnostic and prognostic solves, you can do a single step of it and visualize the results and see if they look sane or not and then do another step. You can manually walk through every time step of the simulation quite easily. And finally, what are the minimum and maximum levels of abstraction that are exposed and can details be encapsulated? So this is a thing that, well, so I'm on one level presenting an interface for you just solve the diagnostic model but it's a big nonlinear PDE and can users also go in and say substitute different solver options for how that nonlinear equation gets solved. So I've presented that interface I'm not gonna show it in this talk but if need be you do, you are able to go down to a lower level. And for example, if a user wants to say use a multi-grid solver to get a faster convergence rate, then they can do so. But this is kind of one of those other things where there are trade-offs and no one really has a perfect answer here. So for example, earlier I'd mentioned that when you use domain-specific languages like FireDrake you give up some of the low-level control. You're sort of handing off a lot of the work to the compiler. Okay, so I wanted to have a demonstration of this. So this is from the documentation for Icepack. What I did is I created this sort of silly geometry where ice is flowing in from this larger semi-circle and it's flowing out from the higher radius one. And this is defining the initial mesh and I had to create some silly input data. So this is a completely goofy initial thickness that I made as a giant nasty algebraic expression. And what I'm gonna do is I'm gonna spin this model up to steady state and we'll see the thickness propagate through the end of the domain and it's still running. Oh, here we go. Oh, one of the things I get the most compliments on seems kind of funny is that I militantly use progress bars in all time-dependent simulations. So people know, do I need to make a whole new pot of coffee or not before this thing finishes running? So to actually look at the details here, again, I had alluded to the fact that you're given fairly fine-grained control over what goes into a simulation. And I think this is on display here. So for every step and the number of time steps, we first call this prognostic solve procedure, which takes in a time step, the current thickness, the current velocity, the accumulation rate, and the inflow thickness values, and it gives you a new value of the ice thickness H, which you can then substitute into the diagnostic solve procedure that takes in the current guess for the velocity, the current thickness, and the fluidity, which is one of the factors that goes into the flow law and it gives you a new value of the thickness. So again, this like, I tried to design this in such a way that the mental model of it that most glaciologists have corresponds as closely as possible to all of the code. And once again, the fact that we're passing all of these arguments in by keyword means that it's relatively easy to replace some of them internally if need be. And, okay, and so this is the final value of the thickness where it's propagated to the end of the domain. And the stream plot of the final velocity, this might take a second to render the movie. I'll show you a movie in a sec, but just for now we'll move on. Okay, so I just wanted to briefly touch on some things that I'm doing lately that I'm pretty excited about. So IcePak includes solvers for the mass and energy balance equations as well. We have, I mentioned earlier that there's more than one momentum balance equation that people like to use. I have all the 2D ones. I have one simplified 3D momentum balance equation and I'm working on the full-stokes equations. I had to do actual math for that. That was hard. And data assimilation, including assimilating sparse data, which took a lot of work. That was by a student, Ruben Hill, who's an Imperial, which I'm also very excited about. It's been used in a couple of publications. So I'll just bring up this one because I'm very pleased about it. This was by a PhD student, Chris Mielly. And he was able to get more or less all of this done with very little help for me. And he's now a postdoc here. And I had a student tell me that it was ridiculously easy to work with. That felt good. And so Ben Hills is a student here. He mostly works on radio glaciology, but he wanted to take some of the 2D and 3D models and simplify them to work in just an X or an XZ geometry. So he was able to make some contributions to the package that he's been using it a little bit as well. So I set out to do this whole thing with the idea of I really wanted to make it so PhD students can use it and they have. So at least in that metric, I think I've largely set out to do what I wanted to. So as I mentioned, there are some aspects of the physics of glacier flow that are not entirely understood. So the one I wanna focus on here is what do we do about the fact that we're really confronting a problem where the geometry has its own dynamics as well. So what do we do about iceberg calving? And if you're familiar with a lot of computational physics, you know, whenever there are dynamic geometries, your life is really, really awful. So there's very dramatic differences in the calving regime between Greenland and Antarctica. So it's sort of a holy grail problem of how do you come up with a model that can describe both of them? We don't really know what physical fields are relevant. So we're going to need a model that's very easy to modify. And yes, dynamic geometries means weeping numerical analysts. So again, I have this hypothesis that's more about the meta science and the science itself, which is that one of the rate limiting factors in us coming up with better calving models is that it's just really awful to simulate. It's very hard if you have some idea for how calving works to actually test it in silico. So the kind of conventional orthodox approach or phase field models or level set methods. And again, I'll point out the Wikipedia page. They're really cool. I've written one before. I think they're fascinating and they're worth knowing about. But the idea is that to describe the dynamic geometry, you can do so implicitly by adding a new spatial field. I'm not going to talk about this a whole lot because I have an even crazier idea, which is can we make the ice flow models work even in the limit of zero thickness? So if we think of solving the momentum balance equation as one big nonlinear equation, we always have to regularize that problem. The reason being that when the thickness goes to zero, this problem becomes degenerate. We get an equation of the form zero equals zero. And it's singular when the gradient of view goes to zero. So there's a term in this big fancy F that goes to infinity. So that's really nasty. And the worst part is they're multiplied by each other. So you get, if you try and do this in floating point of risk, your life is going to be a living hell. So the idea that I had is using what is called the dual form of the momentum balance equation, where rather than eliminate the membrane stress and basal stress variables from the entire problem, we make them unknowns. If you're familiar with, say, the mixed form of the Poisson equation, where you solve for both the field and it's gradient at the same time, this is a similar procedure, or if you know about the Helinger-Reissner principle of elasticity, it's the same idea. Now, the nice part about this dual form is that it changes the character of all the nonlinearities. And what it means is that we can still solve this in the limit of the ice thickness going to zero, really honestly, actually going to zero. So this is a simulation I did of the Larson-C ice shelf in the Antarctic Peninsula. Before the talk started, Greg was giving me a hard time for doing live demos. And of all the things that I didn't expect to break the YouTube video, which is not a live demo at all, is not what I expected. So let's see if we can make this go. Okay, so it starts out advancing slowly. There's a calving event and then it continues advancing again. So again, the only way to do this before would be to oppose some minimum ice thickness, which introduces this awful mask conservation error into the whole problem, and it is really not ideal. So this is kind of new stuff that's hot off the presses and what I'm hoping I've delivered is a nice and new and easy way to do iceberg calving. And just because hopefully it is finally finished. Okay, so this is the promised movie of that spin-up experiment. So I started with a very silly initial ice thickness and propagate it until it reaches the end of the domain by doing that iterated mass moment balance solves. So this is just taken straight from the documentation for ice pack. And if you want to see more of those, this is the website and we have a whole bunch of tutorials, I once again need to credit DL2 because this idea of having a whole bunch of tutorials with a graph that links them as something I blatantly stole from the DL2 documentation. And that's all I've got. Thank you so much, Daniel. You make me simultaneously worried and hopeful. I think worry just when I hear that the mathematicians jaw is falling on the floor with the stream power equation, like what must happen when we actually talk about getting more realistic with physics and the fact that the whole, I think all of geomorphology is probably a dynamic geometry problem. But it's often- You actually work on the stream power incision model. Have you ever shown it to a mathematician and did their jaws actually fall on the floor? Oh, they grumble. I've seen a lot of grumpiness come out of this. Okay, so I'm not wrong there. I don't think so. But it's so cool to see this. So, questions for Daniel? You can either just use the raise hand feature or use the chat if you prefer. And I'll keep an eye on the chat as we go. Looks like there's a question from Mark, go ahead. Hey, Daniel, thanks. I was really cool to learn about your work. I have two really trivial questions and this needs a little embarrassing, but one is that I heard that fire drake is hard to install and configure. Is that still the case? And then secondly, what did you use to make your slides? Oh, yeah, it's still hard to install. You know, I'm not gonna bullshit you there. Everything is a pain in the ass. Are you to use a Mac? Yeah, okay, I'm so sorry. A lot of that is Mac OS' fault. Okay. Do you use Docker? Have you used it at all? Yeah. And then it's fine. Okay. All right. If you go with Docker, then you can use it through that. Yeah. It's, you know, Docker is annoying, but it can be a lifesaver and in this case, it is. Oh, the slides. I'll send, here, I'll send a link to the thing. Cool, thanks. I know it's trivial, but it was neat to watch. Oh, okay. Hey, you know, at least the graphic design is on point even if the math and software engineering is not. Wolfgang. So that was fun to hear then. I really enjoyed it. And for all your self-consciousness about not using deal two, I think you're doing exactly the right thing, right? There's enough space for all of these tools. And... Big background has entered the chat and I was like, oh God, he's gonna eviscerate me. No, no, no, that's neither my intention nor my character. No, the way I really think about it is that, you know, there's different tools for different purposes, right? And I think you describe this, you know, well, that there are trade-offs to all of these sort of things, right? And some of these are technical. Like, do you really need to solve 3D problems that are expensive or is it a 2D problem that you can solve at Python? And the other part is human. And I think that's how I felt about it for a very long time, right? That at the end of the day, we don't write software for computers, we write software for humans to use, right? I mean, that's really what it boils down to. And if a community has a hard time learning certain tools, perhaps because they have a very different background, right? I mean, a lot of glaciologists are probably geologists by training. And do not come out of sort of the PDE of PyPath community, right? And so- Our reform climbers skier, dirtbags. Yeah, so there you go. They just don't have the background to really deal with, you know, the details. And that's totally fine, right? And so then I think what you need to do is choose the right tool and that's exactly what you did. Yeah. Yeah, and I mean, my understanding, so aspect is built on deal too, right? That's right, yes. Right, like my, I guess my impression, and you can correct me if I'm wrong, is that the people who do metal convection are like generally more comfortable with C++ and doing really big parallel things. I think that's, let's say more, so within a glaciologist probably. I think the people who are fascinated by things that happen 2000 kilometers down are perhaps not the reformed dirtbags of skiers, right? They, I mean, there's just not very much connection to what you can see on the surface, right? And so I think these are more people, you know, you think about sort of what would have people published 30 years ago. There's a lot of benchmark studies, right? There's sort of the things that applied mathematicians like to do. So I think there's a substantial overlap. Are they all expert D++ programmers now absolutely not, right? But I think it's perhaps a tool that works better in that community. And then also the problem is 3D, right? And so that limits what you can do, yeah. But I have absolutely no complaints about what you do and certainly not the character inclination to rate you. I know, I'm just giving you a hard time. So we're almost out of time, but I have to ask, I mean, I'll preface my question by saying that one of the things I really love about this is that you're trying to solve, I think a really important problem, which is, you know, we need models, models are half of science, like we have theory and models and we have observations and we're constantly trying to put those two things together to demonstrate our understanding of the natural world. And when it comes to numerical models, there's kind of a dynamic tension between, it takes a lot of work to create a numerical code, whatever it is. And so you put a tremendous amount of investment in this object that represents an idea about how nature works. But of course it's wrong. It's always wrong until you really are really done. And so it needs to be iteratively changed and updated. And I think you've addressed that really nicely here by saying, okay, the stuff that we really feel confident about, we're gonna like really build in, but I'm gonna give you the flexibility to explore the things that we don't understand yet. And I guess the question I have about this is a practical one, given the conversation that you guys just had, a lot of geoscientists don't have a lot of training in applied math and the standard track in math that science students tend to take, takes a really long time just to get to your very first basic PDE. Is there another way we should be approaching quantitative education for scientists that gets people to the things they really need to know faster than how we do it now? Oh, this sounds like a three-tier conversation. I don't know and I also don't, it's a thing I think about, but I also don't feel like I have the wherewithal to make significant waves in that direction. I think there's a lot that I, like if you ask me, are there things you would change? Yeah, plenty. I'll give you the list. For example, well, okay, so I'm gonna teach the Earth Sciences Numerical Methods course here in the spring. And I mean, you know, I can say this because I haven't applied math degree, but there are many respects in which I find courses on, courses on numerical PDEs to be really unsatisfying because it feels like you're doing a lot of work to solve very elementary problems and to solve them really fast, which is of course cool and useful if you have really big problems where speed is the rate limiting factor. So at least one of the things that I'm trying to do with this is have it be like focused first on applications and then we'll go back and do what theory you need to understand like what are good tools for each particular application with the understanding that some of that is going to be superficial and I'm going to have to end up sending them off to read a book at some point, but to like focus more on, first of all, I often find like one of the things that a lot of students have struggled with the hardest is just even understanding how do you derive the variational form of the problem they care about in the first place? So even just being able to specify the problem which really requires practice on how do I come up with the right specification for a ton of different problems? But it requires a lot of wide exposure rather than focusing on a couple of different benchmarks. So that's a very different style of how like how would you teach the class if it were up to you to teach it? No, I haven't taught it yet. So maybe come back in a couple of months and see what my student evaluations are. They might all say that it was total garbage in a waste of their time. But I think at least on some level it requires a change in thinking and the hard part at least, again, I'm coming at this from having done a degree in applied math is that I'd like to see more, I think the geoscientists shouldn't change and I think the applied mathematicians should. It's the short answer, which is a little, well, you know, maybe not the answer they wanna hear but I would love to see more engagement from people who do applied math in the actual domain areas because I think that's really what it's going to take. Or maybe that's more just because I feel more comfortable criticizing the community that I feel like I come from than the one that I've come into. But, you know, I mean, in principle, that's like the joy of doing applied math is that you get to play in everyone else's backyard. So play more, you know. Awesome, thank you. So I think at this point we should release everybody but Wolfgang has a comment so I won't start off the Zoom quite yet. I mean, I'm not gonna keep anybody for longer than necessary but I do want to add to what Dan just said. I think one of the big problems that I find with teaching people computational methods is that they show up for these courses and generally can't program. And I mean, I see this in mathematics because our students have to take two credit hours of programming in the undergraduate program. Well, that's just never going to be enough for anything they're going to do in life, right? And so when they show up in computational modeling courses, you have to teach them the computational modeling and the software and physics, right? And so that's just too much and it doesn't work well. I think that these things, go ahead. Do you teach them get and version control at the same time? Yeah, and all of this, right? And so one of the things that I am working on changing in my department is that we require substantially more programming sort of as a foundation upon which we can build. Because I think that people are going to show up and, well, okay, so you teach them the bound on your form but at least they now have to implement it, right? Whereas right now, you teach them the bound on your form and then you need to tell them that if you run a for loop or a while loop at the condition, the variables in the while loop change in the middle. It's not like the while loop instantly, you know, you jump out, but you have to first get back to the while before you terminate it, that sort of stuff, right? And I think that's not very helpful. And I think that most applied disciplines are no better than applied mathematics at teaching programming as sort of a foundational skill. Yeah. Anyway, I'm not gonna keep anybody any longer. Thanks. It just does indeed sound like a three beer conversation, but for now we're gonna cut it off. Thank you again so much, Daniel, for that fantastic talk and thanks everybody for joining us. Yeah, thanks for having me.