High performance climate modelling: mimetic finite differences, and beyond?
Dubos, Thomas
Climate models simulate atmospheric flows interacting with many physical processes. Because they address long time scales, from centuries to millennia, they need to be efficient, but not at the expense of certain desirable properties, especially conservation of total mass and energy. Most of my talk will explain the design principles behind DYNAMICO, a highly scalable unstructured-mesh energy-conserving finite volume/mimetic finite difference atmospheric flow solver and potential successor of LMD-Z, a structured-mesh (longitude-latitude) solver currently operational as part of IPSL-CM, the Earth System Model developed by Institut Pierre Simon Laplace (IPSL). Specifically, the design exploits the variational structure of the equations of motion and their Hamiltonian formulation, so that the conservation of energy requires only that the discrete grad and div operators be compatible, i.e. that a discrete integration by parts formula holds. I will finish my talk by sketching how the desirable properties of DYNAMICO may be obtained with a different approach based on mixed finite elements (FEM). Indeed while DYNAMICO is very fast and scalable, it is low-order and higher-order accuracy may be desirable. While FEM methods can provide higher-order accuracy, they are computationally more expensive. They offer a viable path only if the performance gap compared to finite differences is not too large. The aim of the CEMRACS project A-HA is to evaluate how wide this gap may be, and whether it can be narrowed by using a recently proposed duality-based approach to assemble the various matrices involved in a FEM method.
 The work I will present is the result of several years of work, including a former PhD student Marine Torx, Sarvesh Dubey, who is now in the US, Yanmer Dessoif from LSE, another lab from the Institut Pierre Simon la Place Federation, and Eva Gellos who will be here in two weeks. And I will talk about climate modeling with an ion high performance. I will explain how we designed an atmospheric flow solver based on mimetic finite differences. And I will sketch ideas that we have to possibly go beyond these finite differences with similar ideas but based on finite elements. So first I will give some background on climate modeling and atmospheric modeling before going into the details of the model that we have been developing the last few years and finish with these finite element ideas. So I think last week you had a presentation by my colleague Fabrice Voitus who is working at Meteor France and doing weather forecasting. So forecasting is really you model the atmospheric flow but it is an initial value problem. You try, you do your best to know the flow tomorrow given an initial condition today. So it's a short term in our time scales. It's a very short term prediction. And I think Fabrice has emphasized the very tight time constraints they have to do this exercise. So if you do the calculation typically they have one hour to give their forecast which means the model must run about 100 times faster than reality to be useful. So for climate we are not looking at the next few days. We are interested in several decades or several centuries. And it's not an initial value problem. It's more a boundary value problem where we look for the statistical equilibrium of the system given some external conditions like solar forcing or the composition of the atmosphere things like this. So we look for the statistics of the model but we do run an initial condition problem, initial value problem for a long time. But the initial value is not an important daytime in this exercise. And because we are interested in these long time scales we need an even faster model. So for modern climate if you have a patient, a PhD student who waits for one month until the simulation gets out of the machine you need to run a thousand times faster than reality. So modern climate is very interesting and important but you may also be interested in paleoclimate. So typically even so-called modern paleoclimate it's how climate changes over a few thousands of years under say changing insulation conditions. And even if you're very patient if you want to have 10,000 years of simulation you need an even faster model. So we will aim at high performance to reach these kind of time scales but we don't want to do this at the expense of the quality of the simulation. So I will explain what we focus on later on. And also in order to reach these high throughput rates we can't really achieve the kind of resolution that you will have in a weather forecasting model which resolves a few kilometers now. In a typical climate model today, sorry, it will be more in maybe 100 kilometers maybe now it's a little bit below that depending on which group is running its model. It used to be 500 kilometers 30, 40 years ago. That means also that the fluid solver is only a small part of the model because very important things happen at subgrid scale. So what we really call a climate model now it's actually what better described as an earth system model. So it's a modeling platform which couples the atmosphere, the ocean, sea ice, possibly biogeochemistry, plankton in the ocean, possibly usually you need for climate continental surfaces that's subsurface hydrology and vegetation. All of these different physics coupled together. So what I will talk about is only the atmospheric component and in fact what I talk about is just the fluid solver inside this atmospheric component. And 90% of an atmospheric model is actually devoted to everything that's subgrid which I won't talk at all about. And this kind of model can be used for terrestrial climate but also for planetary climates like Mars, Venus which are also very interesting scientific objects. Okay, so now more to the specifics. So I'll talk now about the flow solver. So the atmosphere, it's just, you know, the loads of physics are valid for the atmosphere as well. So it's a compressible flow and it obeys Navier-Stokes equations for compressible flow. I will talk only about the inviscid part of the equations for the talk. One first difference compared to standard CFD is that we work in a rotating frame. So we work in a frame attached to the Earth which is rotating. And this implies that in addition to the usual pressure gradient and say gravity, so that would be if the Earth was not rotating, you would have here the gradient of the gravitational potential of the planet. So because of rotation, we have this Coriolis term here. So r is the velocity of the planet itself, the velocity field, a solid body rotation. And the curl of r is twice the rotation rate of the planet. So that's the Coriolis term, 2 omega cross velocity. That's the fluid velocity. And in addition to this Coriolis term, we have a centrifugal term which also derives from a potential. So you can put it together with the gravitational potential to form what's called the geopotential. So even if the Earth were perfectly spherical, and even if V had a spherical symmetry, because of that centrifugal term, the isosurfaces of geopotential would be slightly elliptical. And in fact, V itself isn't perfectly spherical because the Earth is slightly flat. Okay, so that's what we want to solve. And we also have the mass budget here and some thermodynamic equations. So s is entropy. And for what I'd be interested in, this heating term will be zero. So that's the second principle of thermodynamics. Okay. But we don't solve it, the atmosphere is not just a random standard flow. It's a flow in a very specific regime. And to understand that, it's good to have some scales in mind. So there are some velocity scales given by the speed of sound and wind. And some time scales. So the buoyancy oscillation is the oscillation of a piece of air. That imagine that you have some piece of air here. And the air above is slightly warmer. And if I bring this up, this cold air will be in a warmer environment and be heavier and it will sink back. And this creates oscillations at a certain frequency, which is this one. And we have another time scale, which is given by the rotation rate of the planet. So the velocity together with the compressibility of the gas measured by the speed of sound gives rise to a very important scale, which is the high scale, which is essentially the thickness of the atmosphere. So if you go above by a number of these scales, essentially you're in outer space. You can construct this other scale, which is given by the speed of sound on the Coriolis force, which is about 1,000 kilometers on Earth. And it's essentially the size of the cyclones and anti-cyclones. So this is from Jupiter, actually. But that would be the cyclones and anti-cyclones that you find typically at the mid-latitudes. So on Earth, there are some small numbers. Like the Mach number is quite small, the centrifugal acceleration compared to the gravity of the Earth is small, the atmosphere is quite shallow, a few tens of kilometers compared to thousands of kilometers. So there's also a scale separation between the scale height and the Rosby radius. So that means that typically a climate model will capture all these large scales, but in current conditions it will not resolve these smaller scales here, which are widely separated. So this means that there will be room for approximations in order to simplify a little bit the equations we're solving. And the most important one probably is the hydrostatic approximation. So if I take the vertical momentum budget, typically this term here, the vertical acceleration, is negligible at the scales at which a climate model is operating. So you can do some standard estimates based on physics. That's the physical way of doing estimates. And you end up figuring out that this acceleration term is negligible compared to the other ones if the horizontal scale of your flow is large enough compared to the scale height. So for a climate model, this term is small and we typically will neglect vertical acceleration, which has profound consequences on the structure of the equations we saw. But for weather forecasting, they don't do that. Also because the atmosphere is a very thin shell, we don't work in, say, Cartesian coordinating in outer space. We use coordinates that are fitted to the atmosphere. And especially what has been used for a while is to use coordinates which follow the geopotential surfaces. So that means that if you stay on this surface over which phi is a constant, you don't feel gravity because you're on an outer surface of the total potential, the geopotential. And gravity by definition, the vertical by definition is the direction along the gradient of phi. So by using this kind of coordinates, you separate the gravity, which is the dominant force, from the other forces. That also means we need to choose some coordinate system in the horizontal, which cannot be done everywhere on a sphere. So either you use, say, lat-long coordinates, which have a singularity at the pole, or you use more this kind of coordinates, but you need to connect together several patches, like as is shown here. That's the so-called cube sphere. OK, so typically in the next slides, I will be working in curvilinear coordinates, but it doesn't mean that everything becomes more complicated, especially for transport. If you, so when you change from Cartesian to curvilinear coordinates, and if you work with the right quantities, you don't even see that you're working in curvilinear coordinates. So first, your velocity has to be, I'm not sure you see that very well, your velocity is, so the motion of the fluid has to be described by the so-called contravariant velocity component, which are just the Lagrangian derivative of a fluid parcel in your coordinate system. And then you have to take into account the fact that the volume is not given by just dx dy, the volume is given by dx dy dz, but because you have this curvilinear coordinate, there is a Jacobian here. But if you take, if you look at the so-called pseudo density, which is the density of the fluid times the Jacobian, which is essentially the amount of mass per unit psi coordinates, then your transport equations look exactly the same in curvilinear coordinates as in Cartesian coordinates. You don't see very well here. So this would be lambda dot. So lambda is longitude, lambda dot, phi is latitude, phi dot, and r is the distance from the center of the planet. So here I've just taken as an example, latitude, longitude, radius, so three coordinates. So for transport, this looks exactly like Udxs plus, so the three components of Cartesian velocities times the three components of the Cartesian gradient. So if you do things right, you don't really see the coordinate system. So where you do see the coordinate system is in the dynamics. In the momentum equation, here you have to take into account this metric tensor. The metric tensor tells you the distances between neighboring points, and only the dynamics really need to know about that. So you have nice Christopher symbols all over the place, and the rest doesn't change too much. That also means that once you've done this, once you've changed to these curvilinear coordinates, you don't really have to be, well, you can decide how accurate you want to be on these metric terms. And typically what is done in atmospheric models is that because these geopotential surfaces are very close to perfect spheres, we just pretend they are perfect spheres. So in the end, the metric factor we have here are just the metric factors of the spherical coordinates. So you don't even see that the earth is not completely flat. When you forget about this idea that you need to incorporate together gravitational and centrifugal potentials. So that's called the spherical geoid approximation. You can also take into account the fact that the atmospheric is shallow to further simplify the metric terms. In the end, so the real equations of motion are up there, the compressible Euler equation in 3D space. And you can simplify the metric like this as I explained, spherical geoid. You can further simplify by pretending basically that any air parcel is at a fixed distance from the center of the earth. And you can also do this hydrostatic approximation I talked about. So in the end, climate models actually do solve these equations here, which look a lot like standard Euler equations, but they are slightly different. And an important consequence of the hydrostatic approximation is that because it introduces a constraint in the system, because the balance, so the vertical momentum balance, which was an equation used to predict vertical velocity, is now changed into purely balanced equations. That means the system has fewer degrees of freedom when you do the hydrostatic approximation. And this suppresses acoustic waves, which is good because we are not, we don't care about acoustic waves. And from a numerical point of view, they cause additional difficulties because they are fast oscillations. So it's easier to treat these systems without these fast oscillations rather than the full system that Fabrice has been talking about. So that's basically what all climate models do regarding the flow solver. So what do we do to solve it numerically? So what we did before is we had a formulation of these equations in lat-long coordinates, which is nice because it has a Cartesian structure. So you can do centered finite differences and have easily second order numerics with some conservation properties. And for a while, we've had finite volume transport. So as I explained, the transport part doesn't really see anything special. So you can use standard methods like finite volumes or your preferred transport method. So we've had positive definite finite volume transport for a while, basically van Lier schemes, and some more ancient finite difference schemes, but which had some interesting conservation properties, which I will develop now. The only problem with this is that we have this singularity at the poles, which is a feature of the spherical topology. And so from a computational point of view, you can't really distribute the work. So we need to put filters at the pole for stability to remove the fastest oscillations there to have a decent Courant stability condition. So that means we need to apply some Fourier transforms back and forth there, which prevents parallelization along the zonal direction. So you can't really distribute the work in parallel as much as you would like with this kind of model. Essentially, you distribute along the latitude but not longitude. So it creates a bottleneck in scalability for higher resolutions. So the solution, in a way, is well known. You want to go away from this singular mesh and have a more uniform mesh. Then what you lose is the regularity. So essentially, you have a mesh where every cell is different from its neighbors, and it's harder to get second order or even first order accuracy there. So what you get easily is scalability, because then it's just straightforward to distribute the work. What you need to work a little bit for is to preserve this consistency, so properties, essentially halting discrete conservation laws at the discrete level. And also, that's another thing that we wanted. We didn't want to be locked in forever with this hydrostatic primitive equation set. So we wanted to have an approach which would be adaptable to different sets of equations of motion. So especially, we will conserve energy. So just what I call discrete conservation law is the idea that after discretization, you have some quantity that is exactly preserved in the numerical model. So for instance, it's pretty easy to conserve total mass. If you use a finite-volume method, mass is going from one cell to the neighboring cell. And so if you do your communications properly, everything that leaves one cell goes into another cell. And since we don't have real boundaries or we have no flux boundary conditions, total mass is exactly preserved in the discrete model. So that's pretty easy to do this for, say, total mass, total mass of water, total entropy, because you can cast them as flux form. d dt plus divergence of something is zero. What's more difficult is energy, because energy is typically a complicated expression of the variables that you're involving in time. You have internal energy, which is very nonlinear, kinetic energy, which is quadratic, et cetera. So what I will discuss is essentially how to conserve energy at the discrete level in the model. So first, maybe we want a good reason to make that effort, maybe, because if it's not a worthy effort, we're not going to do it. So there are physical reasons to do it. If we do climate, climate modeling is essentially looking at how energy flows in the system. The atmosphere receives energy from the sun, and then this creates motion, and this motion will transport energy here and there from the equator to the pole, for instance. And this will be sent back to space by radiation. So we want to have a good energy budget in the model. So that's a physical motivation to have an energy-conserving numeric. But also from a purely numerical point of view, I think it's interesting to stress that if you conserve energy, which is a convex function of the prognostic variables, then you gain some stability guarantees. So that's energy, right? Kinetic internal potential. And it's typically will prognose rho, u, and s, and this is all convex. So you have a convex functional, and it's conserved. So if you are close to a minimum, you should stay close to that minimum. The thing is, this thing doesn't really have an interesting minimum, but if you have other conserved quantities like total mass or total entropy, you can form what is called an absurdo energy, which now has an interesting minimum. So if you preserve all these three quantities, you can demonstrate that any isothermal state of rest is stable. So physically it's stable, but numerically it will be stable. So it gives you some stability guarantee to your model, which means that conservation, real conservation laws will create limits to what the model can do. Even if it makes errors, these errors will remain within certain physically realizable states. So if now we're convinced we should conserve energy, we can reflect on where this conservation of energy comes from. And it comes from the fact that our equations of motion have a variational structure. So they come from a least action principle. And if you have equations of motion deriving from a least action principle with Lagrangian being time independent, then you conserve energy. So to see this, I rewrote these equations of motion in Cartesian vector form. And you see that if I write Lagrangian, which is kinetic energy plus this Coriolis term, plus internal and potential energy, then I can write the equations exactly in this form. So dL dx dot is just the acceleration, the only term here, which depends on x dot is this one. So if I derive L with respect to x dot, I get x dot and so on. So this internal energy will produce your pressure term, potential energy will produce your gravity term and so on. So you can show that all the equations that I've shown derive from an least action principle. And you can rewrite this in curvilinear coordinates too, but it's just a rewriting of the same thing. So we have a good start. We know that all our equations of motion can be written from a least action principle. So we know why they have the conservation of energy. So that's more a physical interesting thing. Numerically, you will not discretize the least action principle. Some people do that, but it's difficult. Another approach is to use the Hamiltonian formulation. So if you have in classical mechanics, if you have a system which derives from a least action principle, you can derive a Hamiltonian formulation for it. So there is a very elegant algebraic way to do it. I will use a much lower level presentation. So what does it mean that we have a Hamiltonian formulation? So the Hamiltonian formulation says that instead of looking at individual evolution equations for each variable of the system, we will derive a single equation which is valid for any quantity depending on the current state of the system. So if we do have those equations of motion, the kind of which I've already shown, and we have some expression F, which depends on velocity, density, etc., typically some integral, then we know how to compute this time derivative. We just do the derivation under the integral and we can compute everything. So what that means is in order to know how F evolves in time, we only need to know how velocity, density, and entropy evolve in time, and how F depends on all these quantities. That's just a chain rule. There's nothing here but the chain rule. So now once we've written this down, we can replace these by their expressions from the equations of motion. And a system is Hamiltonian if the expression we have is something that is antisymmetric between F and H, the Hamiltonian. So H is the total energy. So that's the general form. So how does that look like for us? So first we will write down the Hamiltonian. So I've written it already. So this is kinetic energy. And when you do the least action principle, it turns out that the momentum, which is the variable conjugate to velocity, to position, is actually the absolute velocity. That is the velocity relative to Earth plus the rotation of the Earth. So to get kinetic energy, I must first subtract the rotation of the Earth, and then apply my metric tensor. And everything else is standard. And then when I look at how my energy evolves in time, here I find the mass flux here, density times contravalent component. Here I find something which looks like a Bernoulli function. So that's kinetic energy plus geopotential plus Gibbs function. That's a thermodynamic potential. And here I get temperature, right? Because when I derive internal energy with respect to entropy, I get temperature. So if I'm able to cast the equations of motion only in terms of these quantities, then I will find my Hamiltonian formulation. And that's actually well known. So I just written down all the steps here. Fortunately, it looks pretty ugly on the screen. So this is the standard way that we usually write the Euler equation. So dTU plus U grad U plus Coyolis plus gravity plus pressure gradient is 0. And then you can do this standard transformation. U grad U equals this, right? That's a very well known trick. You can also use the definition of the Gibbs potential to do this transformation here. And in the end, you end up with exactly what you wanted, which is that other terms can be identified as functional derivatives of the Hamiltonian. So now we get this other but still completed general form of the equations of motion. And now the conservation of energy only depends on the ability to do integration by parts. So when you do the total budget, the only operation you need to get zero at the end is an integration by parts. That's something that's relatively easy to get in a numerical method. So here I've shown how it works for Eulerian coordinates. In practice, what we use are non-Eulerian vertical coordinates. So I have extra slides if you want to know more about this. But basically, that's the idea. So we got a slightly more complicated system. But the basic idea is the same. If you can do integration by parts, then energy doesn't change in time. So now we know from what form to start in order to conserve energy. And we just need to discretize. Because everything I showed was still continuous. At some point, we need to go discrete. So we're working in this thin shell. And I will separate the horizontal coordinates, the two, the psi 1 and psi 2, which essentially are coordinates on the sphere. So I will actually write this now in a more intrinsic way, with n being some point on the unit sphere. And I have my other coordinates, which is an abstract vertical coordinate, going from 0 to 1. And I will discretize separately each of these. So the idea is that you always want in atmospheric or ocean models. As far as I know, everybody who has tried not to have the model mesh aligned with gravity has had terrible problems. So we didn't try to do it differently. So doing this guarantees that you're working with columns, which are perfectly aligned with gravity. So in the horizontal, we will use, so in the atmospheric nomenclature, it's called C grid. I know in other context, it's called MAC. So the idea is that you have, these are control volumes for mass. And we will have mass fluxes across cell edges. And we will describe velocity. We will have one degree of freedom of velocity per point here, per intersection point. So that's staggered mesh. In the vertical, we also have some staggering. So we have model layers. We will put some, the mass and thermodynamic quantities here. But we will, so in time, the attitude of these layers evolves. And we need to keep track of them. And so the attitude is defined at interfaces. So we have like one staggering in the vertical. So if you put this together, you have this nice columns. And then you need to decide how to represent your fields on this discrete mesh. And the choice that worked well for us is very simple in a way. So we will associate, we will attach quantities to the geometric object that they are, they naturally coupled to. So for instance, mass wants to be integrated over a cell volume so that you can do fine volume formulation. That's what we do here. That's an M. So we just integrate this over some cell. And that is the degree of freedom. That's all we know about the mass field is how much mass there is in each cell. To do our fine volume mass budget, we will look at mass fluxes across surfaces. So some other quantities are integrated over 2D surfaces. And the momentum, because it's a conjugate variable, it works well to integrate it along 1D lines. So what we do is we define our degree of freedom for velocity as the integral of momentum along these lines here. And the immediate benefit of this is that you can do very easily the curve. So you do the circulation around that triangle and you get the curl of momentum, which is involved in the equation of motion. So this idea, you can find it into the literature. It's actually very old, but the idea to cast it as a discrete sphere calculus is relatively recent. So the idea is that all these equations I've described, they are actually exterior derivatives. So they are derivatives that don't depend on the metric. And that's why in this context, they are exact. So when you do, if I give you exactly the mass here and the fluxes, then you can do an exact mass budget. That's because the geometry fits together. And another thing is because you compute the curl as a circulation like this, and the grad is just the difference between two points here, then the curl of grad is zero, which is nice. That means that your gradient term will not produce vertices while it shouldn't. So even if you have errors in the term, inside the gradient, the vertices will not be affected. So that's important because the large scale circulation is essentially a vertical circulation. And finally, you can show that you have this integration by parts property. So it's straightforward to derive an energy conserving 3D software. Okay, I have like five minutes left. Three, no, five, because I... Okay. So this is what we've been working on in the last few years. And now we are in a project where we want to apply this model to different applications. So I talked about paleoclimate at the beginning. So we work with paleoclimate colleagues to integrate the atmosphere solver into Earth system model. But the first application, which was more straightforward to obtain was with our planetary climate colleagues, who now work on giant planets. So now this shows a vorticity initially, and now it was showing the wind magnitude. And this is a preliminary attempt to model the circulation of Saturn. It would maybe remind you Jupiter, because it's less known that Saturn has this bounded structure too, but it's actually Jupiter. It's supposed to be Jupiter. So the model has some physics, so it can't really do the equatorial jet right, but it has some reasonable looking features. I just show it again. So this was a pretty big simulation. It's typically... We have more mesh points than what we expect to use in climate in the foreseeable future. I think you see the mesh at some point. I think this one, I think, was one-eighth of a degree, which would mean on Earth something like 12 or 15 kilometers. And we want more to be able to do like 15 kilometers or 25. But this was... So of course, this is driven by heating from the sun. So the colleagues put some... So they say it's simple. So that's a simple radiative transfer scheme. So on Earth, it's much more complicated than this. And this radiative transfer created the heating, which created the density differences, which created the pressure gradient, which put the whole atmosphere in motion. So it was a pretty big simulation. We were also using... oops, sorry, parallel input output in order to get the system also right, the computation that it was doing. So that's leveraging on work by Jan Mordeswa that has been going on for several years to actually offload the work of writing to disk to separate processes. And that's... Yeah, what you don't see... So that's a scaling curve. That's pretty horrible. But... Oh yeah, so you can see the number of cores here. So it's... So what you don't see is how big the problem was. So it's not strong scaling. So this is for... Blue is for, I think, one half of a degree. And I don't really see very well. So that's, say, a small problem. So half of a degree would be 50 kilometers on Earth. That's what we want to do for climate. That's quarter of a degree and that's eighth of a degree. So we can scale to, say, a few tens of thousands of processors on the current machines. Okay, so I will finish with a few words about really what our project here at SamRex is about. So at this point, we are relatively happy with this model. So we want to make it operational, do climate and so on. But we still think about the future. And yes, and so the point is that because the mesh is irregular, because every cell is different from the neighbor, you don't have the consolation of errors that you typically have on a regular mesh. So on a regular mesh, if you do a centered, finite difference in second order, but on an irregular mesh, it's first order. So typically, the model we have now is probably less accurate than the previous one in terms of order of convergence. So that's a bit... That's not a major concern, but we would prefer it possible to have a slightly more accurate model. So and probably find element methods can do the job. So of course, I think it's well known now that you can do high order find elements like spectral elements. So the question of feasibility is not really there. So the main question is whether you can still do this energy conservation exercise. So I will just try to illustrate that on a much simpler set of equations, which is just a nonlinear wave equation. So it's completely manufactured set, but it looks a little bit like the Euler equation I've shown. So here, this is just mass, a mass budget with some mass flux here. And this is some F of H. So it looks a little bit like the gravitational potential, for instance. So you can show that there is a conserved energy in this system where G is the potential for F. And this is a kind of kinetic energy here. So you can recast this nonlinear wave equation in something that looks more or less like a Hamiltonian formulation. And again, the conservation of energy is just the result of an integration by this. So now with this starting point, can we still conserve energy? We find elements. So first, of course, we find elements method, you start with picking spaces. So we have, if we, in the most general setting, we have four kinds of quantities, the time-evolving variables H and U, and the functional derivatives of which we take the gradient. So we have four, potentially four spaces. And if you look at the definition of the functional derivatives from energy, this tells you that these should be obtained by projection. I'm sorry, this really looks ugly on screen. So these terms here, first, you should obtain them as a projection. So that means that these two spaces, H and B, and U and F, they should be in duality, right? Because you need this, the mass matrix here to be well conditioned, to be invertible and well conditioned. And maybe the easiest is to take these spaces identical. So you have only two spaces to choose. And then, since you take gradients, you need at least one of these spaces to have square integrable gradients, to be H1. So let's say it's the flux space, right? So we can take the gradient of the flux here, which is actually divergence, right, but it's not. So now, once you have decided that this space here is H1, this equation will be tested precisely against this space. So since your test function has a gradient, you can do integration by parts and use a weak gradient here. So you don't need this space to have square integrable gradients. And by using a weak derivative, you actually enforce the discrete integration by parts property. Essentially, it's contained in your definition of the gradient. So it's actually pretty straightforward to have discrete integration by parts in a fine element context, too, to conserve energy. So the conservation of energy is not really a big issue. What is a bigger issue is to pick good spaces, especially to have stability and good condition. So that's what we... In principle, it should work and it really should give a very nice numerical method with all the qualitative properties that we are looking for. And the only downside is that typically, a FEM method will be more expensive to compute, right? In a fine difference method, it's very simple. If you look at the code, you just... For a gradient, you just take the difference between two points and you're done. A divergence is just some sum of our neighboring cells. So it's very... The model now, I'm sorry, this is falling, the fine difference model is very cheap, it's very simple, but it's low order. The question is how much more expensive... So we know our FEM method will be more expensive. The question is whether it's twice or 10 times or 100 times. And that can depend also on how you implement all these costly operations like assembling your matrices or computing right-hand sides. And so that's... So Chris, who's in the AHA project with Manel, who are there in the back. I saw them. Yes, I saw Chris. I don't see... Manel is there. So they are looking at a technique proposed by Kirby to rearrange the computation and use a dense matrix, matrix multiplication to make full use of externally optimized libraries. And I hope they will know and you will know in two weeks how much that can help in making this kind of FEM method sufficient. So that's my summary. So I hope I've given you some better knowledge of what we actually solve in a climate model, the kind of equations we actually solve. In our group, we have put a strong emphasis on conservation properties and some work to really go to the source of these conservation laws. And what we've tried to be is to be systematic. So everything I've shown, we can change, say, the equation of state or the geometry. And we know we can still have all the qualitative properties we want. So we're not stuck with a trick that would tie us to a specific set of equations. And that's thanks to the Hamiltonian formulation. So I haven't really showed much detail with the mimetic find difference method, but I have extra slides if you want. So it works very well with these mimetic find difference methods. And also probably we find element discretizations. But the end performance is an important issue. Thank you very much. So I'm rephrasing for the recording point of view. So given the size of the systems that you're considering, I would expect in finite elements that the assembling process, that is, of course, the stage when you are computing the matrices, should be negligible compared to the solution process, so the inversing of the system. Do we have any estimate of the respective times that it would require? We have only very preliminary experience, and the right person to ask is Chris about this. But I think assembly is still quite, well, these are two issues. But I think for high order, you do spend some, maybe not assembling the matrix, but at least computing the right-hand sides, the nonlinear terms. Maybe Chris, you can comment on that. I think that's very architecture and implementation. We think that we had kind of differences within a factor of, say, three to four. But I think it's going to be very implementation and architecture dependent. With respect to the number of degrees of freedom, right? Yeah. With respect to the number of elements, obviously, so it's directly proportional to the number of degrees of freedom. And then, of course, assembling is embarrassingly prone to parallel strategies, because you can distribute your assembling process over as many processes as you like. So it should be very easy to lower the cost of this phase, I think. I don't want to take too much. It should be very scalable. But we also think in terms of the overall resource that we are using. Other question? I have just one quick. Can you give us a little bit more insight about the lack of physics for the Jupiter case? Who was the... Ah, yeah. So, Jupiter... It's Saturn. It's actually Saturn. So on Earth, we have some heating from the surface, right? Because we have a solid surface and the radiation heats the surface and this causes convection. So on Saturn, you don't have a solid surface. It's a gaseous planet. But you still do have convection. And that's because there is some heat going from the interior. So actually, Saturn has a positive energy budget. It still radiates more energy than it receives from the Sun. So it's a belief. There are some arguments that in order to get the zonal... The jet at the equator in the correct direction, you need to have some input of energy from the inside of the planet. So here it's not there. And I don't think I've shown... So maybe I've shown that. The jet... No, I haven't. The jet is in the wrong direction at the equator. And I don't get the number of degrees of freedom that you have in the geopotential direction for your... So typically, for this one, eighth-degree simulation, I think we are about one million cells in the horizontal and a few tens in the vertical. When we discuss at the canum, we talk about parallelization in time. Is it still something you could imagine? Well, I think that's something very interesting. It will take a long time until it gets operational. But that's something that several people are looking at. Because probably we will reach the point where we have used all the available parallelism in space at some point in the next decade or so. I don't know exactly when. But so if it's feasible, and there are many obstacles, maybe not in the fluid solver, but in everything else, if it's feasible, then it would be a very interesting direction to use. So personally, I'm interested in thinking at least about it. Maybe the last one question for me. When we used to do energy conservation time discretization, it's most of the time implicit. Are you also... Yeah, so I didn't talk about time discretization at all. So in time, we use standard explicit rule-scruta time stepping, which means that even if the spatial discretization is energy conserving, you lose that once your discretization time. But you still get some stability properties. Like, you can show that linear stability is still there. And since in time, we are fully explicit. So we are limited by a current condition based on horizontal sand waves. So the time stepping is relatively small. So with a say third or fourth order runch-scruta, you have very small time discretization errors. Okay. Thank you very much. Other questions? Okay. So thanks again, Thomas. Thank you.
 The work I will present is the result of several years of work, including a former PhD student  Marine Torx, Sarvesh Dubey, who is now in the US, Yanmer Dessoif from LSE, another lab  from the Institut Pierre Simon la Place Federation, and Eva Gellos who will be here in two weeks.  And I will talk about climate modeling with an ion high performance.  I will explain how we designed an atmospheric flow solver based on mimetic finite differences.  And I will sketch ideas that we have to possibly go beyond these finite differences with similar  ideas but based on finite elements.  So first I will give some background on climate modeling and atmospheric modeling before  going into the details of the model that we have been developing the last few years and  finish with these finite element ideas.  So I think last week you had a presentation by my colleague Fabrice Voitus who is working  at Meteor France and doing weather forecasting.  So forecasting is really you model the atmospheric flow but it is an initial value problem.  You try, you do your best to know the flow tomorrow given an initial condition today.  So it's a short term in our time scales.  It's a very short term prediction.  And I think Fabrice has emphasized the very tight time constraints they have to do this  exercise.  So if you do the calculation typically they have one hour to give their forecast which  means the model must run about 100 times faster than reality to be useful.  So for climate we are not looking at the next few days.  We are interested in several decades or several centuries.  And it's not an initial value problem.  It's more a boundary value problem where we look for the statistical equilibrium of the  system given some external conditions like solar forcing or the composition of the atmosphere  things like this.  So we look for the statistics of the model but we do run an initial condition problem,  initial value problem for a long time.  But the initial value is not an important daytime in this exercise.  And because we are interested in these long time scales we need an even faster model.  So for modern climate if you have a patient, a PhD student who waits for one month until  the simulation gets out of the machine you need to run a thousand times faster than reality.  So modern climate is very interesting and important but you may also be interested in  paleoclimate.  So typically even so-called modern paleoclimate it's how climate changes over a few thousands  of years under say changing insulation conditions.  And even if you're very patient if you want to have 10,000 years of simulation you need  an even faster model.  So we will aim at high performance to reach these kind of time scales but we don't want  to do this at the expense of the quality of the simulation.  So I will explain what we focus on later on.  And also in order to reach these high throughput rates we can't really achieve the kind of  resolution that you will have in a weather forecasting model which resolves a few kilometers  now.  In a typical climate model today, sorry, it will be more in maybe 100 kilometers maybe  now it's a little bit below that depending on which group is running its model.  It used to be 500 kilometers 30, 40 years ago.  That means also that the fluid solver is only a small part of the model because very important  things happen at subgrid scale.  So what we really call a climate model now it's actually what better described as an  earth system model.  So it's a modeling platform which couples the atmosphere, the ocean, sea ice, possibly  biogeochemistry, plankton in the ocean, possibly usually you need for climate continental surfaces  that's subsurface hydrology and vegetation.  All of these different physics coupled together.  So what I will talk about is only the atmospheric component and in fact what I talk about is  just the fluid solver inside this atmospheric component.  And 90% of an atmospheric model is actually devoted to everything that's subgrid which  I won't talk at all about.  And this kind of model can be used for terrestrial climate but also for planetary climates like  Mars, Venus which are also very interesting scientific objects.  Okay, so now more to the specifics.  So I'll talk now about the flow solver.  So the atmosphere, it's just, you know, the loads of physics are valid for the atmosphere  as well.  So it's a compressible flow and it obeys Navier-Stokes equations for compressible flow.  I will talk only about the inviscid part of the equations for the talk.  One first difference compared to standard CFD is that we work in a rotating frame.  So we work in a frame attached to the Earth which is rotating.  And this implies that in addition to the usual pressure gradient and say gravity, so that  would be if the Earth was not rotating, you would have here the gradient of the gravitational  potential of the planet.  So because of rotation, we have this Coriolis term here.  So r is the velocity of the planet itself, the velocity field, a solid body rotation.  And the curl of r is twice the rotation rate of the planet.  So that's the Coriolis term, 2 omega cross velocity.  That's the fluid velocity.  And in addition to this Coriolis term, we have a centrifugal term which also derives  from a potential.  So you can put it together with the gravitational potential to form what's called the geopotential.  So even if the Earth were perfectly spherical, and even if V had a spherical symmetry, because  of that centrifugal term, the isosurfaces of geopotential would be slightly elliptical.  And in fact, V itself isn't perfectly spherical because the Earth is slightly flat.  Okay, so that's what we want to solve.  And we also have the mass budget here and some thermodynamic equations.  So s is entropy.  And for what I'd be interested in, this heating term will be zero.  So that's the second principle of thermodynamics.  Okay.  But we don't solve it, the atmosphere is not just a random standard flow.  It's a flow in a very specific regime.  And to understand that, it's good to have some scales in mind.  So there are some velocity scales given by the speed of sound and wind.  And some time scales.  So the buoyancy oscillation is the oscillation of a piece of air.  That imagine that you have some piece of air here.  And the air above is slightly warmer.  And if I bring this up, this cold air will be in a warmer environment and be heavier  and it will sink back.  And this creates oscillations at a certain frequency, which is this one.  And we have another time scale, which is given by the rotation rate of the planet.  So the velocity together with the compressibility of the gas measured by the speed of sound  gives rise to a very important scale, which is the high scale, which is essentially the  thickness of the atmosphere.  So if you go above by a number of these scales, essentially you're in outer space.  You can construct this other scale, which is given by the speed of sound on the Coriolis  force, which is about 1,000 kilometers on Earth.  And it's essentially the size of the cyclones and anti-cyclones.  So this is from Jupiter, actually.  But that would be the cyclones and anti-cyclones that you find typically at the mid-latitudes.  So on Earth, there are some small numbers.  Like the Mach number is quite small, the centrifugal acceleration compared to the gravity of the  Earth is small, the atmosphere is quite shallow, a few tens of kilometers compared to thousands  of kilometers.  So there's also a scale separation between the scale height and the Rosby radius.  So that means that typically a climate model will capture all these large scales, but in  current conditions it will not resolve these smaller scales here, which are widely separated.  So this means that there will be room for approximations in order to simplify a little  bit the equations we're solving.  And the most important one probably is the hydrostatic approximation.  So if I take the vertical momentum budget, typically this term here, the vertical acceleration,  is negligible at the scales at which a climate model is operating.  So you can do some standard estimates based on physics.  That's the physical way of doing estimates.  And you end up figuring out that this acceleration term is negligible compared to the other ones  if the horizontal scale of your flow is large enough compared to the scale height.  So for a climate model, this term is small and we typically will neglect vertical acceleration,  which has profound consequences on the structure of the equations we saw.  But for weather forecasting, they don't do that.  Also because the atmosphere is a very thin shell, we don't work in, say, Cartesian coordinating  in outer space.  We use coordinates that are fitted to the atmosphere.  And especially what has been used for a while is to use coordinates which follow the geopotential  surfaces.  So that means that if you stay on this surface over which phi is a constant, you don't feel  gravity because you're on an outer surface of the total potential, the geopotential.  And gravity by definition, the vertical by definition is the direction along the gradient  of phi.  So by using this kind of coordinates, you separate the gravity, which is the dominant  force, from the other forces.  That also means we need to choose some coordinate system in the horizontal, which cannot be  done everywhere on a sphere.  So either you use, say, lat-long coordinates, which have a singularity at the pole, or you  use more this kind of coordinates, but you need to connect together several patches,  like as is shown here.  That's the so-called cube sphere.  OK, so typically in the next slides, I will be working in curvilinear coordinates, but  it doesn't mean that everything becomes more complicated, especially for transport.  If you, so when you change from Cartesian to curvilinear coordinates, and if you work  with the right quantities, you don't even see that you're working in curvilinear coordinates.  So first, your velocity has to be, I'm not sure you see that very well, your velocity  is, so the motion of the fluid has to be described by the so-called contravariant velocity component,  which are just the Lagrangian derivative of a fluid parcel in your coordinate system.  And then you have to take into account the fact that the volume is not given by just  dx dy, the volume is given by dx dy dz, but because you have this curvilinear coordinate,  there is a Jacobian here.  But if you take, if you look at the so-called pseudo density, which is the density of the  fluid times the Jacobian, which is essentially the amount of mass per unit psi coordinates,  then your transport equations look exactly the same in curvilinear coordinates as in  Cartesian coordinates.  You don't see very well here.  So this would be lambda dot.  So lambda is longitude, lambda dot, phi is latitude, phi dot, and r is the distance from  the center of the planet.  So here I've just taken as an example, latitude, longitude, radius, so three coordinates.  So for transport, this looks exactly like Udxs plus, so the three components of Cartesian  velocities times the three components of the Cartesian gradient.  So if you do things right, you don't really see the coordinate system.  So where you do see the coordinate system is in the dynamics.  In the momentum equation, here you have to take into account this metric tensor.  The metric tensor tells you the distances between neighboring points, and only the dynamics  really need to know about that.  So you have nice Christopher symbols all over the place, and the rest doesn't change  too much.  That also means that once you've done this, once you've changed to these curvilinear coordinates,  you don't really have to be, well, you can decide how accurate you want to be on these  metric terms.  And typically what is done in atmospheric models is that because these geopotential surfaces  are very close to perfect spheres, we just pretend they are perfect spheres.  So in the end, the metric factor we have here are just the metric factors of the spherical  coordinates.  So you don't even see that the earth is not completely flat.  When you forget about this idea that you need to incorporate together gravitational and  centrifugal potentials.  So that's called the spherical geoid approximation.  You can also take into account the fact that the atmospheric is shallow to further simplify  the metric terms.  In the end, so the real equations of motion are up there, the compressible Euler equation  in 3D space.  And you can simplify the metric like this as I explained, spherical geoid.  You can further simplify by pretending basically that any air parcel is at a fixed distance  from the center of the earth.  And you can also do this hydrostatic approximation I talked about.  So in the end, climate models actually do solve these equations here, which look a lot  like standard Euler equations, but they are slightly different.  And an important consequence of the hydrostatic approximation is that because it introduces  a constraint in the system, because the balance, so the vertical momentum balance, which was  an equation used to predict vertical velocity, is now changed into purely balanced equations.  That means the system has fewer degrees of freedom when you do the hydrostatic approximation.  And this suppresses acoustic waves, which is good because we are not, we don't care  about acoustic waves.  And from a numerical point of view, they cause additional difficulties because they  are fast oscillations.  So it's easier to treat these systems without these fast oscillations rather than the full  system that Fabrice has been talking about.  So that's basically what all climate models do regarding the flow solver.  So what do we do to solve it numerically?  So what we did before is we had a formulation of these equations in lat-long coordinates,  which is nice because it has a Cartesian structure.  So you can do centered finite differences and have easily second order numerics with  some conservation properties.  And for a while, we've had finite volume transport.  So as I explained, the transport part doesn't really see anything special.  So you can use standard methods like finite volumes or your preferred transport method.  So we've had positive definite finite volume transport for a while, basically van Lier  schemes, and some more ancient finite difference schemes, but which had some interesting conservation  properties, which I will develop now.  The only problem with this is that we have this singularity at the poles, which is a  feature of the spherical topology.  And so from a computational point of view, you can't really distribute the work.  So we need to put filters at the pole for stability to remove the fastest oscillations  there to have a decent Courant stability condition.  So that means we need to apply some Fourier transforms back and forth there, which prevents  parallelization along the zonal direction.  So you can't really distribute the work in parallel as much as you would like with this  kind of model.  Essentially, you distribute along the latitude but not longitude.  So it creates a bottleneck in scalability for higher resolutions.  So the solution, in a way, is well known.  You want to go away from this singular mesh and have a more uniform mesh.  Then what you lose is the regularity.  So essentially, you have a mesh where every cell is different from its neighbors, and  it's harder to get second order or even first order accuracy there.  So what you get easily is scalability, because then it's just straightforward to distribute  the work.  What you need to work a little bit for is to preserve this consistency, so properties,  essentially halting discrete conservation laws at the discrete level.  And also, that's another thing that we wanted.  We didn't want to be locked in forever with this hydrostatic primitive equation set.  So we wanted to have an approach which would be adaptable to different sets of equations  of motion.  So especially, we will conserve energy.  So just what I call discrete conservation law is the idea that after discretization,  you have some quantity that is exactly preserved in the numerical model.  So for instance, it's pretty easy to conserve total mass.  If you use a finite-volume method, mass is going from one cell to the neighboring cell.  And so if you do your communications properly, everything that leaves one cell goes into  another cell.  And since we don't have real boundaries or we have no flux boundary conditions, total  mass is exactly preserved in the discrete model.  So that's pretty easy to do this for, say, total mass, total mass of water, total entropy,  because you can cast them as flux form.  d dt plus divergence of something is zero.  What's more difficult is energy, because energy is typically a complicated expression of the  variables that you're involving in time.  You have internal energy, which is very nonlinear, kinetic energy, which is quadratic, et cetera.  So what I will discuss is essentially how to conserve energy at the discrete level in  the model.  So first, maybe we want a good reason to make that effort, maybe, because if it's not a worthy  effort, we're not going to do it.  So there are physical reasons to do it.  If we do climate, climate modeling is essentially looking at how energy flows in the system.  The atmosphere receives energy from the sun, and then this creates motion, and this motion  will transport energy here and there from the equator to the pole, for instance.  And this will be sent back to space by radiation.  So we want to have a good energy budget in the model.  So that's a physical motivation to have an energy-conserving numeric.  But also from a purely numerical point of view, I think it's interesting to stress that  if you conserve energy, which is a convex function of the prognostic variables, then  you gain some stability guarantees.  So that's energy, right?  Kinetic internal potential.  And it's typically will prognose rho, u, and s, and this is all convex.  So you have a convex functional, and it's conserved.  So if you are close to a minimum, you should stay close to that minimum.  The thing is, this thing doesn't really have an interesting minimum, but if you have other  conserved quantities like total mass or total entropy, you can form what is called an absurdo  energy, which now has an interesting minimum.  So if you preserve all these three quantities, you can demonstrate that any isothermal state  of rest is stable.  So physically it's stable, but numerically it will be stable.  So it gives you some stability guarantee to your model, which means that conservation,  real conservation laws will create limits to what the model can do.  Even if it makes errors, these errors will remain within certain physically realizable  states.  So if now we're convinced we should conserve energy, we can reflect on where this conservation  of energy comes from.  And it comes from the fact that our equations of motion have a variational structure.  So they come from a least action principle.  And if you have equations of motion deriving from a least action principle with Lagrangian  being time independent, then you conserve energy.  So to see this, I rewrote these equations of motion in Cartesian vector form.  And you see that if I write Lagrangian, which is kinetic energy plus this Coriolis term,  plus internal and potential energy, then I can write the equations exactly in this form.  So dL dx dot is just the acceleration, the only term here, which depends on x dot is  this one.  So if I derive L with respect to x dot, I get x dot and so on.  So this internal energy will produce your pressure term, potential energy will produce  your gravity term and so on.  So you can show that all the equations that I've shown derive from an least action principle.  And you can rewrite this in curvilinear coordinates too, but it's just a rewriting of the same  thing.  So we have a good start.  We know that all our equations of motion can be written from a least action principle.  So we know why they have the conservation of energy.  So that's more a physical interesting thing.  Numerically, you will not discretize the least action principle.  Some people do that, but it's difficult.  Another approach is to use the Hamiltonian formulation.  So if you have in classical mechanics, if you have a system which derives from a least  action principle, you can derive a Hamiltonian formulation for it.  So there is a very elegant algebraic way to do it.  I will use a much lower level presentation.  So what does it mean that we have a Hamiltonian formulation?  So the Hamiltonian formulation says that instead of looking at individual evolution equations  for each variable of the system, we will derive a single equation which is valid for any quantity  depending on the current state of the system.  So if we do have those equations of motion, the kind of which I've already shown, and  we have some expression F, which depends on velocity, density, etc., typically some integral,  then we know how to compute this time derivative.  We just do the derivation under the integral and we can compute everything.  So what that means is in order to know how F evolves in time, we only need to know how  velocity, density, and entropy evolve in time, and how F depends on all these quantities.  That's just a chain rule.  There's nothing here but the chain rule.  So now once we've written this down, we can replace these by their expressions from the  equations of motion.  And a system is Hamiltonian if the expression we have is something that is antisymmetric  between F and H, the Hamiltonian.  So H is the total energy.  So that's the general form.  So how does that look like for us?  So first we will write down the Hamiltonian.  So I've written it already.  So this is kinetic energy.  And when you do the least action principle, it turns out that the momentum, which is the  variable conjugate to velocity, to position, is actually the absolute velocity.  That is the velocity relative to Earth plus the rotation of the Earth.  So to get kinetic energy, I must first subtract the rotation of the Earth, and then apply  my metric tensor.  And everything else is standard.  And then when I look at how my energy evolves in time, here I find the mass flux here, density  times contravalent component.  Here I find something which looks like a Bernoulli function.  So that's kinetic energy plus geopotential plus Gibbs function.  That's a thermodynamic potential.  And here I get temperature, right?  Because when I derive internal energy with respect to entropy, I get temperature.  So if I'm able to cast the equations of motion only in terms of these quantities, then I  will find my Hamiltonian formulation.  And that's actually well known.  So I just written down all the steps here.  Fortunately, it looks pretty ugly on the screen.  So this is the standard way that we usually write the Euler equation.  So dTU plus U grad U plus Coyolis plus gravity plus pressure gradient is 0.  And then you can do this standard transformation.  U grad U equals this, right?  That's a very well known trick.  You can also use the definition of the Gibbs potential to do this transformation here.  And in the end, you end up with exactly what you wanted, which is that other terms can  be identified as functional derivatives of the Hamiltonian.  So now we get this other but still completed general form of the equations of motion.  And now the conservation of energy only depends on the ability to do integration by parts.  So when you do the total budget, the only operation you need to get zero at the end  is an integration by parts.  That's something that's relatively easy to get in a numerical method.  So here I've shown how it works for Eulerian coordinates.  In practice, what we use are non-Eulerian vertical coordinates.  So I have extra slides if you want to know more about this.  But basically, that's the idea.  So we got a slightly more complicated system.  But the basic idea is the same.  If you can do integration by parts, then energy doesn't change in time.  So now we know from what form to start in order to conserve energy.  And we just need to discretize.  Because everything I showed was still continuous.  At some point, we need to go discrete.  So we're working in this thin shell.  And I will separate the horizontal coordinates, the two, the psi 1 and psi 2, which essentially  are coordinates on the sphere.  So I will actually write this now in a more intrinsic way, with n being some point on  the unit sphere.  And I have my other coordinates, which is an abstract vertical coordinate, going from  0 to 1.  And I will discretize separately each of these.  So the idea is that you always want in atmospheric or ocean models.  As far as I know, everybody who has tried not to have the model mesh aligned with gravity  has had terrible problems.  So we didn't try to do it differently.  So doing this guarantees that you're working with columns, which are perfectly aligned  with gravity.  So in the horizontal, we will use, so in the atmospheric nomenclature, it's called  C grid.  I know in other context, it's called MAC.  So the idea is that you have, these are control volumes for mass.  And we will have mass fluxes across cell edges.  And we will describe velocity.  We will have one degree of freedom of velocity per point here, per intersection point.  So that's staggered mesh.  In the vertical, we also have some staggering.  So we have model layers.  We will put some, the mass and thermodynamic quantities here.  But we will, so in time, the attitude of these layers evolves.  And we need to keep track of them.  And so the attitude is defined at interfaces.  So we have like one staggering in the vertical.  So if you put this together, you have this nice columns.  And then you need to decide how to represent your fields on this discrete mesh.  And the choice that worked well for us is very simple in a way.  So we will associate, we will attach quantities to the geometric object that they are, they  naturally coupled to.  So for instance, mass wants to be integrated over a cell volume so that you can do fine  volume formulation.  That's what we do here.  That's an M. So we just integrate this over some cell.  And that is the degree of freedom.  That's all we know about the mass field is how much mass there is in each cell.  To do our fine volume mass budget, we will look at mass fluxes across surfaces.  So some other quantities are integrated over 2D surfaces.  And the momentum, because it's a conjugate variable, it works well to integrate it along  1D lines.  So what we do is we define our degree of freedom for velocity as the integral of momentum  along these lines here.  And the immediate benefit of this is that you can do very easily the curve.  So you do the circulation around that triangle and you get the curl of momentum, which is  involved in the equation of motion.  So this idea, you can find it into the literature.  It's actually very old, but the idea to cast it as a discrete sphere calculus is relatively  recent.  So the idea is that all these equations I've described, they are actually exterior derivatives.  So they are derivatives that don't depend on the metric.  And that's why in this context, they are exact.  So when you do, if I give you exactly the mass here and the fluxes, then you can do  an exact mass budget.  That's because the geometry fits together.  And another thing is because you compute the curl as a circulation like this, and the  grad is just the difference between two points here, then the curl of grad is zero, which  is nice.  That means that your gradient term will not produce vertices while it shouldn't.  So even if you have errors in the term, inside the gradient, the vertices will not be affected.  So that's important because the large scale circulation is essentially a vertical circulation.  And finally, you can show that you have this integration by parts property.  So it's straightforward to derive an energy conserving 3D software.  Okay, I have like five minutes left.  Three, no, five, because I...  Okay.  So this is what we've been working on in the last few years.  And now we are in a project where we want to apply this model to different applications.  So I talked about paleoclimate at the beginning.  So we work with paleoclimate colleagues to integrate the atmosphere solver into Earth  system model.  But the first application, which was more straightforward to obtain was with our planetary  climate colleagues, who now work on giant planets.  So now this shows a vorticity initially, and now it was showing the wind magnitude.  And this is a preliminary attempt to model the circulation of Saturn.  It would maybe remind you Jupiter, because it's less known that Saturn has this bounded  structure too, but it's actually Jupiter.  It's supposed to be Jupiter.  So the model has some physics, so it can't really do the equatorial jet right, but it  has some reasonable looking features.  I just show it again.  So this was a pretty big simulation.  It's typically...  We have more mesh points than what we expect to use in climate in the foreseeable future.  I think you see the mesh at some point.  I think this one, I think, was one-eighth of a degree, which would mean on Earth something  like 12 or 15 kilometers.  And we want more to be able to do like 15 kilometers or 25.  But this was...  So of course, this is driven by heating from the sun.  So the colleagues put some...  So they say it's simple.  So that's a simple radiative transfer scheme.  So on Earth, it's much more complicated than this.  And this radiative transfer created the heating, which created the density differences, which  created the pressure gradient, which put the whole atmosphere in motion.  So it was a pretty big simulation.  We were also using... oops, sorry, parallel input output in order to get the system also  right, the computation that it was doing.  So that's leveraging on work by Jan Mordeswa that has been going on for several years to  actually offload the work of writing to disk to separate processes.  And that's...  Yeah, what you don't see...  So that's a scaling curve.  That's pretty horrible.  But...  Oh yeah, so you can see the number of cores here.  So it's...  So what you don't see is how big the problem was.  So it's not strong scaling.  So this is for...  Blue is for, I think, one half of a degree.  And I don't really see very well.  So that's, say, a small problem.  So half of a degree would be 50 kilometers on Earth.  That's what we want to do for climate.  That's quarter of a degree and that's eighth of a degree.  So we can scale to, say, a few tens of thousands of processors on the current machines.  Okay, so I will finish with a few words about really what our project here at SamRex is  about.  So at this point, we are relatively happy with this model.  So we want to make it operational, do climate and so on.  But we still think about the future.  And yes, and so the point is that because the mesh is irregular, because every cell  is different from the neighbor, you don't have the consolation of errors that you typically  have on a regular mesh.  So on a regular mesh, if you do a centered, finite difference in second order, but on  an irregular mesh, it's first order.  So typically, the model we have now is probably less accurate than the previous one in terms  of order of convergence.  So that's a bit...  That's not a major concern, but we would prefer it possible to have a slightly more accurate  model.  So and probably find element methods can do the job.  So of course, I think it's well known now that you can do high order find elements like  spectral elements.  So the question of feasibility is not really there.  So the main question is whether you can still do this energy conservation exercise.  So I will just try to illustrate that on a much simpler set of equations, which is just  a nonlinear wave equation.  So it's completely manufactured set, but it looks a little bit like the Euler equation  I've shown.  So here, this is just mass, a mass budget with some mass flux here.  And this is some F of H. So it looks a little bit like the gravitational potential, for  instance.  So you can show that there is a conserved energy in this system where G is the potential  for F. And this is a kind of kinetic energy here.  So you can recast this nonlinear wave equation in something that looks more or less like  a Hamiltonian formulation.  And again, the conservation of energy is just the result of an integration by this.  So now with this starting point, can we still conserve energy?  We find elements.  So first, of course, we find elements method, you start with picking spaces.  So we have, if we, in the most general setting, we have four kinds of quantities, the time-evolving  variables H and U, and the functional derivatives of which we take the gradient.  So we have four, potentially four spaces.  And if you look at the definition of the functional derivatives from energy, this tells you that  these should be obtained by projection.  I'm sorry, this really looks ugly on screen.  So these terms here, first, you should obtain them as a projection.  So that means that these two spaces, H and B, and U and F, they should be in duality,  right?  Because you need this, the mass matrix here to be well conditioned, to be invertible and  well conditioned.  And maybe the easiest is to take these spaces identical.  So you have only two spaces to choose.  And then, since you take gradients, you need at least one of these spaces to have square  integrable gradients, to be H1.  So let's say it's the flux space, right?  So we can take the gradient of the flux here, which is actually divergence, right, but it's  not.  So now, once you have decided that this space here is H1, this equation will be tested precisely  against this space.  So since your test function has a gradient, you can do integration by parts and use a  weak gradient here.  So you don't need this space to have square integrable gradients.  And by using a weak derivative, you actually enforce the discrete integration by parts  property.  Essentially, it's contained in your definition of the gradient.  So it's actually pretty straightforward to have discrete integration by parts in a fine  element context, too, to conserve energy.  So the conservation of energy is not really a big issue.  What is a bigger issue is to pick good spaces, especially to have stability and good condition.  So that's what we...  In principle, it should work and it really should give a very nice numerical method with  all the qualitative properties that we are looking for.  And the only downside is that typically, a FEM method will be more expensive to compute,  right?  In a fine difference method, it's very simple.  If you look at the code, you just...  For a gradient, you just take the difference between two points and you're done.  A divergence is just some sum of our neighboring cells.  So it's very...  The model now, I'm sorry, this is falling, the fine difference model is very cheap, it's  very simple, but it's low order.  The question is how much more expensive...  So we know our FEM method will be more expensive.  The question is whether it's twice or 10 times or 100 times.  And that can depend also on how you implement all these costly operations like assembling  your matrices or computing right-hand sides.  And so that's...  So Chris, who's in the AHA project with Manel, who are there in the back.  I saw them.  Yes, I saw Chris.  I don't see...  Manel is there.  So they are looking at a technique proposed by Kirby to rearrange the computation and  use a dense matrix, matrix multiplication to make full use of externally optimized libraries.  And I hope they will know and you will know in two weeks how much that can help in making  this kind of FEM method sufficient.  So that's my summary.  So I hope I've given you some better knowledge of what we actually solve in a climate model,  the kind of equations we actually solve.  In our group, we have put a strong emphasis on conservation properties and some work to  really go to the source of these conservation laws.  And what we've tried to be is to be systematic.  So everything I've shown, we can change, say, the equation of state or the geometry.  And we know we can still have all the qualitative properties we want.  So we're not stuck with a trick that would tie us to a specific set of equations.  And that's thanks to the Hamiltonian formulation.  So I haven't really showed much detail with the mimetic find difference method, but I  have extra slides if you want.  So it works very well with these mimetic find difference methods.  And also probably we find element discretizations.  But the end performance is an important issue.  Thank you very much.  So I'm rephrasing for the recording point of view.  So given the size of the systems that you're considering, I would expect in finite elements  that the assembling process, that is, of course, the stage when you are computing the matrices,  should be negligible compared to the solution process, so the inversing of the system.  Do we have any estimate of the respective times that it would require?  We have only very preliminary experience, and the right person to ask is Chris about  this.  But I think assembly is still quite, well, these are two issues.  But I think for high order, you do spend some, maybe not assembling the matrix, but at least  computing the right-hand sides, the nonlinear terms.  Maybe Chris, you can comment on that.  I think that's very architecture and implementation.  We think that we had kind of differences within a factor of, say, three to four.  But I think it's going to be very implementation and architecture dependent.  With respect to the number of degrees of freedom, right?  Yeah.  With respect to the number of elements, obviously, so it's directly proportional to the number  of degrees of freedom.  And then, of course, assembling is embarrassingly prone to parallel strategies, because you  can distribute your assembling process over as many processes as you like.  So it should be very easy to lower the cost of this phase, I think.  I don't want to take too much.  It should be very scalable.  But we also think in terms of the overall resource that we are using.  Other question?  I have just one quick.  Can you give us a little bit more insight about the lack of physics for the Jupiter case?  Who was the...  Ah, yeah.  So, Jupiter...  It's Saturn.  It's actually Saturn.  So on Earth, we have some heating from the surface, right?  Because we have a solid surface and the radiation heats the surface and this causes convection.  So on Saturn, you don't have a solid surface.  It's a gaseous planet.  But you still do have convection.  And that's because there is some heat going from the interior.  So actually, Saturn has a positive energy budget.  It still radiates more energy than it receives from the Sun.  So it's a belief.  There are some arguments that in order to get the zonal...  The jet at the equator in the correct direction, you need to have some input of energy from  the inside of the planet.  So here it's not there.  And I don't think I've shown...  So maybe I've shown that.  The jet...  No, I haven't.  The jet is in the wrong direction at the equator.  And I don't get the number of degrees of freedom that you have in the geopotential direction  for your...  So typically, for this one, eighth-degree simulation, I think we are about one million  cells in the horizontal and a few tens in the vertical.  When we discuss at the canum, we talk about parallelization in time.  Is it still something you could imagine?  Well, I think that's something very interesting.  It will take a long time until it gets operational.  But that's something that several people are looking at.  Because probably we will reach the point where we have used all the available parallelism  in space at some point in the next decade or so.  I don't know exactly when.  But so if it's feasible, and there are many obstacles, maybe not in the fluid solver, but  in everything else, if it's feasible, then it would be a very interesting direction to  use.  So personally, I'm interested in thinking at least about it.  Maybe the last one question for me.  When we used to do energy conservation time discretization, it's most of the time implicit.  Are you also...  Yeah, so I didn't talk about time discretization at all.  So in time, we use standard explicit rule-scruta time stepping, which means that even if the  spatial discretization is energy conserving, you lose that once your discretization time.  But you still get some stability properties.  Like, you can show that linear stability is still there.  And since in time, we are fully explicit.  So we are limited by a current condition based on horizontal sand waves.  So the time stepping is relatively small.  So with a say third or fourth order runch-scruta, you have very small time discretization errors.  Okay.  Thank you very much.  Other questions?  Okay.  So thanks again, Thomas.  Thank you.
2016-08-09_Dubos_0_47.png 2016-08-09_Dubos_48_95.png 2016-08-09_Dubos_96_155.png 2016-08-09_Dubos_156_203.png 2016-08-09_Dubos_204_275.png 2016-08-09_Dubos_276_395.png 2016-08-09_Dubos_396_515.png 2016-08-09_Dubos_516_1835.png 2016-08-09_Dubos_1836_2339.png 2016-08-09_Dubos_2340_3767.png 2016-08-09_Dubos_3768_6563.png 2016-08-09_Dubos_6564_7055.png 2016-08-09_Dubos_7056_7655.png 2016-08-09_Dubos_7656_9527.png 2016-08-09_Dubos_9528_13055.png 2016-08-09_Dubos_13056_16823.png 2016-08-09_Dubos_16824_18203.png 2016-08-09_Dubos_18204_20171.png 2016-08-09_Dubos_20172_21347.png 2016-08-09_Dubos_21348_23495.png 2016-08-09_Dubos_23496_23639.png 2016-08-09_Dubos_23640_24167.png 2016-08-09_Dubos_24168_24239.png 2016-08-09_Dubos_24240_24599.png 2016-08-09_Dubos_24600_24947.png 2016-08-09_Dubos_24948_25127.png 2016-08-09_Dubos_25128_25307.png 2016-08-09_Dubos_25308_25571.png 2016-08-09_Dubos_25572_27731.png 2016-08-09_Dubos_27732_28355.png 2016-08-09_Dubos_28356_28739.png 2016-08-09_Dubos_28740_32951.png 2016-08-09_Dubos_32952_33107.png 2016-08-09_Dubos_33108_34919.png 2016-08-09_Dubos_34920_37967.png 2016-08-09_Dubos_37968_41939.png 2016-08-09_Dubos_41940_44903.png 2016-08-09_Dubos_44904_46907.png 2016-08-09_Dubos_46908_49019.png 2016-08-09_Dubos_49020_49607.png 2016-08-09_Dubos_49608_50867.png 2016-08-09_Dubos_50868_52607.png 2016-08-09_Dubos_52608_53375.png 2016-08-09_Dubos_53376_55199.png 2016-08-09_Dubos_55200_58295.png 2016-08-09_Dubos_58296_58463.png 2016-08-09_Dubos_58464_58523.png 2016-08-09_Dubos_58524_58631.png 2016-08-09_Dubos_58632_58703.png 2016-08-09_Dubos_58704_59003.png 2016-08-09_Dubos_59004_59051.png 2016-08-09_Dubos_59052_60143.png 2016-08-09_Dubos_60144_60191.png 2016-08-09_Dubos_60192_60347.png 2016-08-09_Dubos_60348_60407.png 2016-08-09_Dubos_60408_60515.png 2016-08-09_Dubos_60516_60587.png 2016-08-09_Dubos_60588_60731.png 2016-08-09_Dubos_60732_60791.png 2016-08-09_Dubos_60792_60863.png 2016-08-09_Dubos_60864_60923.png 2016-08-09_Dubos_60924_60971.png 2016-08-09_Dubos_60972_61763.png 2016-08-09_Dubos_61764_61895.png 2016-08-09_Dubos_61896_62495.png 2016-08-09_Dubos_62496_63647.png 2016-08-09_Dubos_63648_64187.png 2016-08-09_Dubos_64188_64235.png 2016-08-09_Dubos_64236_71555.png 2016-08-09_Dubos_71556_74531.png 2016-08-09_Dubos_74532_81755.png 2016-08-09_Dubos_81756_81803.png 2016-08-09_Dubos_81804_81851.png 2016-08-09_Dubos_81852_85679.png 2016-08-09_Dubos_85680_85751.png 2016-08-09_Dubos_85752_85799.png 2016-08-09_Dubos_85800_85847.png 2016-08-09_Dubos_85848_85917.png
The work I will present is the result of several years of work, including a former PhD student Marine Torx, Sarvesh Dubey, who is now in the US, Yanmer Dessoif from LSE, another lab from the Institut Pierre Simon la Place Federation, and Eva Gellos who will be here in two weeks. And I will talk about climate modeling with an ion high performance. I will explain how we designed an atmospheric flow solver based on mimetic finite differences. And I will sketch ideas that we have to possibly go beyond these finite differences with similar ideas but based on finite elements. So first I will give some background on climate modeling and atmospheric modeling before going into the details of the model that we have been developing the last few years and finish with these finite element ideas. So I think last week you had a presentation by my colleague Fabrice Voitus who is working at Meteor France and doing weather forecasting. So forecasting is really you model the atmospheric flow but it is an initial value problem. You try, you do your best to know the flow tomorrow given an initial condition today. So it's a short term in our time scales. It's a very short term prediction. And I think Fabrice has emphasized the very tight time constraints they have to do this exercise. So if you do the calculation typically they have one hour to give their forecast which means the model must run about 100 times faster than reality to be useful. So for climate we are not looking at the next few days. We are interested in several decades or several centuries. And it's not an initial value problem. It's more a boundary value problem where we look for the statistical equilibrium of the system given some external conditions like solar forcing or the composition of the atmosphere things like this. So we look for the statistics of the model but we do run an initial condition problem, initial value problem for a long time. But the initial value is not an important daytime in this exercise. And because we are interested in these long time scales we need an even faster model. So for modern climate if you have a patient, a PhD student who waits for one month until the simulation gets out of the machine you need to run a thousand times faster than reality. So modern climate is very interesting and important but you may also be interested in paleoclimate. So typically even so-called modern paleoclimate it's how climate changes over a few thousands of years under say changing insulation conditions. And even if you're very patient if you want to have 10,000 years of simulation you need an even faster model. So we will aim at high performance to reach these kind of time scales but we don't want to do this at the expense of the quality of the simulation. So I will explain what we focus on later on. And also in order to reach these high throughput rates we can't really achieve the kind of resolution that you will have in a weather forecasting model which resolves a few kilometers now. In a typical climate model today, sorry, it will be more in maybe 100 kilometers maybe now it's a little bit below that depending on which group is running its model. It used to be 500 kilometers 30, 40 years ago. That means also that the fluid solver is only a small part of the model because very important things happen at subgrid scale. So what we really call a climate model now it's actually what better described as an earth system model. So it's a modeling platform which couples the atmosphere, the ocean, sea ice, possibly biogeochemistry, plankton in the ocean, possibly usually you need for climate continental surfaces that's subsurface hydrology and vegetation. All of these different physics coupled together. So what I will talk about is only the atmospheric component and in fact what I talk about is just the fluid solver inside this atmospheric component. And 90% of an atmospheric model is actually devoted to everything that's subgrid which I won't talk at all about. And this kind of model can be used for terrestrial climate but also for planetary climates like Mars, Venus which are also very interesting scientific objects. Okay, so now more to the specifics. So I'll talk now about the flow solver. So the atmosphere, it's just, you know, the loads of physics are valid for the atmosphere as well. So it's a compressible flow and it obeys Navier-Stokes equations for compressible flow. I will talk only about the inviscid part of the equations for the talk. One first difference compared to standard CFD is that we work in a rotating frame. So we work in a frame attached to the Earth which is rotating. And this implies that in addition to the usual pressure gradient and say gravity, so that would be if the Earth was not rotating, you would have here the gradient of the gravitational potential of the planet. So because of rotation, we have this Coriolis term here. So r is the velocity of the planet itself, the velocity field, a solid body rotation. And the curl of r is twice the rotation rate of the planet. So that's the Coriolis term, 2 omega cross velocity. That's the fluid velocity. And in addition to this Coriolis term, we have a centrifugal term which also derives from a potential. So you can put it together with the gravitational potential to form what's called the geopotential. So even if the Earth were perfectly spherical, and even if V had a spherical symmetry, because of that centrifugal term, the isosurfaces of geopotential would be slightly elliptical. And in fact, V itself isn't perfectly spherical because the Earth is slightly flat. Okay, so that's what we want to solve. And we also have the mass budget here and some thermodynamic equations. So s is entropy. And for what I'd be interested in, this heating term will be zero. So that's the second principle of thermodynamics. Okay. But we don't solve it, the atmosphere is not just a random standard flow. It's a flow in a very specific regime. And to understand that, it's good to have some scales in mind. So there are some velocity scales given by the speed of sound and wind. And some time scales. So the buoyancy oscillation is the oscillation of a piece of air. That imagine that you have some piece of air here. And the air above is slightly warmer. And if I bring this up, this cold air will be in a warmer environment and be heavier and it will sink back. And this creates oscillations at a certain frequency, which is this one. And we have another time scale, which is given by the rotation rate of the planet. So the velocity together with the compressibility of the gas measured by the speed of sound gives rise to a very important scale, which is the high scale, which is essentially the thickness of the atmosphere. So if you go above by a number of these scales, essentially you're in outer space. You can construct this other scale, which is given by the speed of sound on the Coriolis force, which is about 1,000 kilometers on Earth. And it's essentially the size of the cyclones and anti-cyclones. So this is from Jupiter, actually. But that would be the cyclones and anti-cyclones that you find typically at the mid-latitudes. So on Earth, there are some small numbers. Like the Mach number is quite small, the centrifugal acceleration compared to the gravity of the Earth is small, the atmosphere is quite shallow, a few tens of kilometers compared to thousands of kilometers. So there's also a scale separation between the scale height and the Rosby radius. So that means that typically a climate model will capture all these large scales, but in current conditions it will not resolve these smaller scales here, which are widely separated. So this means that there will be room for approximations in order to simplify a little bit the equations we're solving. And the most important one probably is the hydrostatic approximation. So if I take the vertical momentum budget, typically this term here, the vertical acceleration, is negligible at the scales at which a climate model is operating. So you can do some standard estimates based on physics. That's the physical way of doing estimates. And you end up figuring out that this acceleration term is negligible compared to the other ones if the horizontal scale of your flow is large enough compared to the scale height. So for a climate model, this term is small and we typically will neglect vertical acceleration, which has profound consequences on the structure of the equations we saw. But for weather forecasting, they don't do that. Also because the atmosphere is a very thin shell, we don't work in, say, Cartesian coordinating in outer space. We use coordinates that are fitted to the atmosphere. And especially what has been used for a while is to use coordinates which follow the geopotential surfaces. So that means that if you stay on this surface over which phi is a constant, you don't feel gravity because you're on an outer surface of the total potential, the geopotential. And gravity by definition, the vertical by definition is the direction along the gradient of phi. So by using this kind of coordinates, you separate the gravity, which is the dominant force, from the other forces. That also means we need to choose some coordinate system in the horizontal, which cannot be done everywhere on a sphere. So either you use, say, lat-long coordinates, which have a singularity at the pole, or you use more this kind of coordinates, but you need to connect together several patches, like as is shown here. That's the so-called cube sphere. OK, so typically in the next slides, I will be working in curvilinear coordinates, but it doesn't mean that everything becomes more complicated, especially for transport. If you, so when you change from Cartesian to curvilinear coordinates, and if you work with the right quantities, you don't even see that you're working in curvilinear coordinates. So first, your velocity has to be, I'm not sure you see that very well, your velocity is, so the motion of the fluid has to be described by the so-called contravariant velocity component, which are just the Lagrangian derivative of a fluid parcel in your coordinate system. And then you have to take into account the fact that the volume is not given by just dx dy, the volume is given by dx dy dz, but because you have this curvilinear coordinate, there is a Jacobian here. But if you take, if you look at the so-called pseudo density, which is the density of the fluid times the Jacobian, which is essentially the amount of mass per unit psi coordinates, then your transport equations look exactly the same in curvilinear coordinates as in Cartesian coordinates. You don't see very well here. So this would be lambda dot. So lambda is longitude, lambda dot, phi is latitude, phi dot, and r is the distance from the center of the planet. So here I've just taken as an example, latitude, longitude, radius, so three coordinates. So for transport, this looks exactly like Udxs plus, so the three components of Cartesian velocities times the three components of the Cartesian gradient. So if you do things right, you don't really see the coordinate system. So where you do see the coordinate system is in the dynamics. In the momentum equation, here you have to take into account this metric tensor. The metric tensor tells you the distances between neighboring points, and only the dynamics really need to know about that. So you have nice Christopher symbols all over the place, and the rest doesn't change too much. That also means that once you've done this, once you've changed to these curvilinear coordinates, you don't really have to be, well, you can decide how accurate you want to be on these metric terms. And typically what is done in atmospheric models is that because these geopotential surfaces are very close to perfect spheres, we just pretend they are perfect spheres. So in the end, the metric factor we have here are just the metric factors of the spherical coordinates. So you don't even see that the earth is not completely flat. When you forget about this idea that you need to incorporate together gravitational and centrifugal potentials. So that's called the spherical geoid approximation. You can also take into account the fact that the atmospheric is shallow to further simplify the metric terms. In the end, so the real equations of motion are up there, the compressible Euler equation in 3D space. And you can simplify the metric like this as I explained, spherical geoid. You can further simplify by pretending basically that any air parcel is at a fixed distance from the center of the earth. And you can also do this hydrostatic approximation I talked about. So in the end, climate models actually do solve these equations here, which look a lot like standard Euler equations, but they are slightly different. And an important consequence of the hydrostatic approximation is that because it introduces a constraint in the system, because the balance, so the vertical momentum balance, which was an equation used to predict vertical velocity, is now changed into purely balanced equations. That means the system has fewer degrees of freedom when you do the hydrostatic approximation. And this suppresses acoustic waves, which is good because we are not, we don't care about acoustic waves. And from a numerical point of view, they cause additional difficulties because they are fast oscillations. So it's easier to treat these systems without these fast oscillations rather than the full system that Fabrice has been talking about. So that's basically what all climate models do regarding the flow solver. So what do we do to solve it numerically? So what we did before is we had a formulation of these equations in lat-long coordinates, which is nice because it has a Cartesian structure. So you can do centered finite differences and have easily second order numerics with some conservation properties. And for a while, we've had finite volume transport. So as I explained, the transport part doesn't really see anything special. So you can use standard methods like finite volumes or your preferred transport method. So we've had positive definite finite volume transport for a while, basically van Lier schemes, and some more ancient finite difference schemes, but which had some interesting conservation properties, which I will develop now. The only problem with this is that we have this singularity at the poles, which is a feature of the spherical topology. And so from a computational point of view, you can't really distribute the work. So we need to put filters at the pole for stability to remove the fastest oscillations there to have a decent Courant stability condition. So that means we need to apply some Fourier transforms back and forth there, which prevents parallelization along the zonal direction. So you can't really distribute the work in parallel as much as you would like with this kind of model. Essentially, you distribute along the latitude but not longitude. So it creates a bottleneck in scalability for higher resolutions. So the solution, in a way, is well known. You want to go away from this singular mesh and have a more uniform mesh. Then what you lose is the regularity. So essentially, you have a mesh where every cell is different from its neighbors, and it's harder to get second order or even first order accuracy there. So what you get easily is scalability, because then it's just straightforward to distribute the work. What you need to work a little bit for is to preserve this consistency, so properties, essentially halting discrete conservation laws at the discrete level. And also, that's another thing that we wanted. We didn't want to be locked in forever with this hydrostatic primitive equation set. So we wanted to have an approach which would be adaptable to different sets of equations of motion. So especially, we will conserve energy. So just what I call discrete conservation law is the idea that after discretization, you have some quantity that is exactly preserved in the numerical model. So for instance, it's pretty easy to conserve total mass. If you use a finite-volume method, mass is going from one cell to the neighboring cell. And so if you do your communications properly, everything that leaves one cell goes into another cell. And since we don't have real boundaries or we have no flux boundary conditions, total mass is exactly preserved in the discrete model. So that's pretty easy to do this for, say, total mass, total mass of water, total entropy, because you can cast them as flux form. d dt plus divergence of something is zero. What's more difficult is energy, because energy is typically a complicated expression of the variables that you're involving in time. You have internal energy, which is very nonlinear, kinetic energy, which is quadratic, et cetera. So what I will discuss is essentially how to conserve energy at the discrete level in the model. So first, maybe we want a good reason to make that effort, maybe, because if it's not a worthy effort, we're not going to do it. So there are physical reasons to do it. If we do climate, climate modeling is essentially looking at how energy flows in the system. The atmosphere receives energy from the sun, and then this creates motion, and this motion will transport energy here and there from the equator to the pole, for instance. And this will be sent back to space by radiation. So we want to have a good energy budget in the model. So that's a physical motivation to have an energy-conserving numeric. But also from a purely numerical point of view, I think it's interesting to stress that if you conserve energy, which is a convex function of the prognostic variables, then you gain some stability guarantees. So that's energy, right? Kinetic internal potential. And it's typically will prognose rho, u, and s, and this is all convex. So you have a convex functional, and it's conserved. So if you are close to a minimum, you should stay close to that minimum. The thing is, this thing doesn't really have an interesting minimum, but if you have other conserved quantities like total mass or total entropy, you can form what is called an absurdo energy, which now has an interesting minimum. So if you preserve all these three quantities, you can demonstrate that any isothermal state of rest is stable. So physically it's stable, but numerically it will be stable. So it gives you some stability guarantee to your model, which means that conservation, real conservation laws will create limits to what the model can do. Even if it makes errors, these errors will remain within certain physically realizable states. So if now we're convinced we should conserve energy, we can reflect on where this conservation of energy comes from. And it comes from the fact that our equations of motion have a variational structure. So they come from a least action principle. And if you have equations of motion deriving from a least action principle with Lagrangian being time independent, then you conserve energy. So to see this, I rewrote these equations of motion in Cartesian vector form. And you see that if I write Lagrangian, which is kinetic energy plus this Coriolis term, plus internal and potential energy, then I can write the equations exactly in this form. So dL dx dot is just the acceleration, the only term here, which depends on x dot is this one. So if I derive L with respect to x dot, I get x dot and so on. So this internal energy will produce your pressure term, potential energy will produce your gravity term and so on. So you can show that all the equations that I've shown derive from an least action principle. And you can rewrite this in curvilinear coordinates too, but it's just a rewriting of the same thing. So we have a good start. We know that all our equations of motion can be written from a least action principle. So we know why they have the conservation of energy. So that's more a physical interesting thing. Numerically, you will not discretize the least action principle. Some people do that, but it's difficult. Another approach is to use the Hamiltonian formulation. So if you have in classical mechanics, if you have a system which derives from a least action principle, you can derive a Hamiltonian formulation for it. So there is a very elegant algebraic way to do it. I will use a much lower level presentation. So what does it mean that we have a Hamiltonian formulation? So the Hamiltonian formulation says that instead of looking at individual evolution equations for each variable of the system, we will derive a single equation which is valid for any quantity depending on the current state of the system. So if we do have those equations of motion, the kind of which I've already shown, and we have some expression F, which depends on velocity, density, etc., typically some integral, then we know how to compute this time derivative. We just do the derivation under the integral and we can compute everything. So what that means is in order to know how F evolves in time, we only need to know how velocity, density, and entropy evolve in time, and how F depends on all these quantities. That's just a chain rule. There's nothing here but the chain rule. So now once we've written this down, we can replace these by their expressions from the equations of motion. And a system is Hamiltonian if the expression we have is something that is antisymmetric between F and H, the Hamiltonian. So H is the total energy. So that's the general form. So how does that look like for us? So first we will write down the Hamiltonian. So I've written it already. So this is kinetic energy. And when you do the least action principle, it turns out that the momentum, which is the variable conjugate to velocity, to position, is actually the absolute velocity. That is the velocity relative to Earth plus the rotation of the Earth. So to get kinetic energy, I must first subtract the rotation of the Earth, and then apply my metric tensor. And everything else is standard. And then when I look at how my energy evolves in time, here I find the mass flux here, density times contravalent component. Here I find something which looks like a Bernoulli function. So that's kinetic energy plus geopotential plus Gibbs function. That's a thermodynamic potential. And here I get temperature, right? Because when I derive internal energy with respect to entropy, I get temperature. So if I'm able to cast the equations of motion only in terms of these quantities, then I will find my Hamiltonian formulation. And that's actually well known. So I just written down all the steps here. Fortunately, it looks pretty ugly on the screen. So this is the standard way that we usually write the Euler equation. So dTU plus U grad U plus Coyolis plus gravity plus pressure gradient is 0. And then you can do this standard transformation. U grad U equals this, right? That's a very well known trick. You can also use the definition of the Gibbs potential to do this transformation here. And in the end, you end up with exactly what you wanted, which is that other terms can be identified as functional derivatives of the Hamiltonian. So now we get this other but still completed general form of the equations of motion. And now the conservation of energy only depends on the ability to do integration by parts. So when you do the total budget, the only operation you need to get zero at the end is an integration by parts. That's something that's relatively easy to get in a numerical method. So here I've shown how it works for Eulerian coordinates. In practice, what we use are non-Eulerian vertical coordinates. So I have extra slides if you want to know more about this. But basically, that's the idea. So we got a slightly more complicated system. But the basic idea is the same. If you can do integration by parts, then energy doesn't change in time. So now we know from what form to start in order to conserve energy. And we just need to discretize. Because everything I showed was still continuous. At some point, we need to go discrete. So we're working in this thin shell. And I will separate the horizontal coordinates, the two, the psi 1 and psi 2, which essentially are coordinates on the sphere. So I will actually write this now in a more intrinsic way, with n being some point on the unit sphere. And I have my other coordinates, which is an abstract vertical coordinate, going from 0 to 1. And I will discretize separately each of these. So the idea is that you always want in atmospheric or ocean models. As far as I know, everybody who has tried not to have the model mesh aligned with gravity has had terrible problems. So we didn't try to do it differently. So doing this guarantees that you're working with columns, which are perfectly aligned with gravity. So in the horizontal, we will use, so in the atmospheric nomenclature, it's called C grid. I know in other context, it's called MAC. So the idea is that you have, these are control volumes for mass. And we will have mass fluxes across cell edges. And we will describe velocity. We will have one degree of freedom of velocity per point here, per intersection point. So that's staggered mesh. In the vertical, we also have some staggering. So we have model layers. We will put some, the mass and thermodynamic quantities here. But we will, so in time, the attitude of these layers evolves. And we need to keep track of them. And so the attitude is defined at interfaces. So we have like one staggering in the vertical. So if you put this together, you have this nice columns. And then you need to decide how to represent your fields on this discrete mesh. And the choice that worked well for us is very simple in a way. So we will associate, we will attach quantities to the geometric object that they are, they naturally coupled to. So for instance, mass wants to be integrated over a cell volume so that you can do fine volume formulation. That's what we do here. That's an M. So we just integrate this over some cell. And that is the degree of freedom. That's all we know about the mass field is how much mass there is in each cell. To do our fine volume mass budget, we will look at mass fluxes across surfaces. So some other quantities are integrated over 2D surfaces. And the momentum, because it's a conjugate variable, it works well to integrate it along 1D lines. So what we do is we define our degree of freedom for velocity as the integral of momentum along these lines here. And the immediate benefit of this is that you can do very easily the curve. So you do the circulation around that triangle and you get the curl of momentum, which is involved in the equation of motion. So this idea, you can find it into the literature. It's actually very old, but the idea to cast it as a discrete sphere calculus is relatively recent. So the idea is that all these equations I've described, they are actually exterior derivatives. So they are derivatives that don't depend on the metric. And that's why in this context, they are exact. So when you do, if I give you exactly the mass here and the fluxes, then you can do an exact mass budget. That's because the geometry fits together. And another thing is because you compute the curl as a circulation like this, and the grad is just the difference between two points here, then the curl of grad is zero, which is nice. That means that your gradient term will not produce vertices while it shouldn't. So even if you have errors in the term, inside the gradient, the vertices will not be affected. So that's important because the large scale circulation is essentially a vertical circulation. And finally, you can show that you have this integration by parts property. So it's straightforward to derive an energy conserving 3D software. Okay, I have like five minutes left. Three, no, five, because I... Okay. So this is what we've been working on in the last few years. And now we are in a project where we want to apply this model to different applications. So I talked about paleoclimate at the beginning. So we work with paleoclimate colleagues to integrate the atmosphere solver into Earth system model. But the first application, which was more straightforward to obtain was with our planetary climate colleagues, who now work on giant planets. So now this shows a vorticity initially, and now it was showing the wind magnitude. And this is a preliminary attempt to model the circulation of Saturn. It would maybe remind you Jupiter, because it's less known that Saturn has this bounded structure too, but it's actually Jupiter. It's supposed to be Jupiter. So the model has some physics, so it can't really do the equatorial jet right, but it has some reasonable looking features. I just show it again. So this was a pretty big simulation. It's typically... We have more mesh points than what we expect to use in climate in the foreseeable future. I think you see the mesh at some point. I think this one, I think, was one-eighth of a degree, which would mean on Earth something like 12 or 15 kilometers. And we want more to be able to do like 15 kilometers or 25. But this was... So of course, this is driven by heating from the sun. So the colleagues put some... So they say it's simple. So that's a simple radiative transfer scheme. So on Earth, it's much more complicated than this. And this radiative transfer created the heating, which created the density differences, which created the pressure gradient, which put the whole atmosphere in motion. So it was a pretty big simulation. We were also using... oops, sorry, parallel input output in order to get the system also right, the computation that it was doing. So that's leveraging on work by Jan Mordeswa that has been going on for several years to actually offload the work of writing to disk to separate processes. And that's... Yeah, what you don't see... So that's a scaling curve. That's pretty horrible. But... Oh yeah, so you can see the number of cores here. So it's... So what you don't see is how big the problem was. So it's not strong scaling. So this is for... Blue is for, I think, one half of a degree. And I don't really see very well. So that's, say, a small problem. So half of a degree would be 50 kilometers on Earth. That's what we want to do for climate. That's quarter of a degree and that's eighth of a degree. So we can scale to, say, a few tens of thousands of processors on the current machines. Okay, so I will finish with a few words about really what our project here at SamRex is about. So at this point, we are relatively happy with this model. So we want to make it operational, do climate and so on. But we still think about the future. And yes, and so the point is that because the mesh is irregular, because every cell is different from the neighbor, you don't have the consolation of errors that you typically have on a regular mesh. So on a regular mesh, if you do a centered, finite difference in second order, but on an irregular mesh, it's first order. So typically, the model we have now is probably less accurate than the previous one in terms of order of convergence. So that's a bit... That's not a major concern, but we would prefer it possible to have a slightly more accurate model. So and probably find element methods can do the job. So of course, I think it's well known now that you can do high order find elements like spectral elements. So the question of feasibility is not really there. So the main question is whether you can still do this energy conservation exercise. So I will just try to illustrate that on a much simpler set of equations, which is just a nonlinear wave equation. So it's completely manufactured set, but it looks a little bit like the Euler equation I've shown. So here, this is just mass, a mass budget with some mass flux here. And this is some F of H. So it looks a little bit like the gravitational potential, for instance. So you can show that there is a conserved energy in this system where G is the potential for F. And this is a kind of kinetic energy here. So you can recast this nonlinear wave equation in something that looks more or less like a Hamiltonian formulation. And again, the conservation of energy is just the result of an integration by this. So now with this starting point, can we still conserve energy? We find elements. So first, of course, we find elements method, you start with picking spaces. So we have, if we, in the most general setting, we have four kinds of quantities, the time-evolving variables H and U, and the functional derivatives of which we take the gradient. So we have four, potentially four spaces. And if you look at the definition of the functional derivatives from energy, this tells you that these should be obtained by projection. I'm sorry, this really looks ugly on screen. So these terms here, first, you should obtain them as a projection. So that means that these two spaces, H and B, and U and F, they should be in duality, right? Because you need this, the mass matrix here to be well conditioned, to be invertible and well conditioned. And maybe the easiest is to take these spaces identical. So you have only two spaces to choose. And then, since you take gradients, you need at least one of these spaces to have square integrable gradients, to be H1. So let's say it's the flux space, right? So we can take the gradient of the flux here, which is actually divergence, right, but it's not. So now, once you have decided that this space here is H1, this equation will be tested precisely against this space. So since your test function has a gradient, you can do integration by parts and use a weak gradient here. So you don't need this space to have square integrable gradients. And by using a weak derivative, you actually enforce the discrete integration by parts property. Essentially, it's contained in your definition of the gradient. So it's actually pretty straightforward to have discrete integration by parts in a fine element context, too, to conserve energy. So the conservation of energy is not really a big issue. What is a bigger issue is to pick good spaces, especially to have stability and good condition. So that's what we... In principle, it should work and it really should give a very nice numerical method with all the qualitative properties that we are looking for. And the only downside is that typically, a FEM method will be more expensive to compute, right? In a fine difference method, it's very simple. If you look at the code, you just... For a gradient, you just take the difference between two points and you're done. A divergence is just some sum of our neighboring cells. So it's very... The model now, I'm sorry, this is falling, the fine difference model is very cheap, it's very simple, but it's low order. The question is how much more expensive... So we know our FEM method will be more expensive. The question is whether it's twice or 10 times or 100 times. And that can depend also on how you implement all these costly operations like assembling your matrices or computing right-hand sides. And so that's... So Chris, who's in the AHA project with Manel, who are there in the back. I saw them. Yes, I saw Chris. I don't see... Manel is there. So they are looking at a technique proposed by Kirby to rearrange the computation and use a dense matrix, matrix multiplication to make full use of externally optimized libraries. And I hope they will know and you will know in two weeks how much that can help in making this kind of FEM method sufficient. So that's my summary. So I hope I've given you some better knowledge of what we actually solve in a climate model, the kind of equations we actually solve. In our group, we have put a strong emphasis on conservation properties and some work to really go to the source of these conservation laws. And what we've tried to be is to be systematic. So everything I've shown, we can change, say, the equation of state or the geometry. And we know we can still have all the qualitative properties we want. So we're not stuck with a trick that would tie us to a specific set of equations. And that's thanks to the Hamiltonian formulation. So I haven't really showed much detail with the mimetic find difference method, but I have extra slides if you want. So it works very well with these mimetic find difference methods. And also probably we find element discretizations. But the end performance is an important issue. Thank you very much. So I'm rephrasing for the recording point of view. So given the size of the systems that you're considering, I would expect in finite elements that the assembling process, that is, of course, the stage when you are computing the matrices, should be negligible compared to the solution process, so the inversing of the system. Do we have any estimate of the respective times that it would require? We have only very preliminary experience, and the right person to ask is Chris about this. But I think assembly is still quite, well, these are two issues. But I think for high order, you do spend some, maybe not assembling the matrix, but at least computing the right-hand sides, the nonlinear terms. Maybe Chris, you can comment on that. I think that's very architecture and implementation. We think that we had kind of differences within a factor of, say, three to four. But I think it's going to be very implementation and architecture dependent. With respect to the number of degrees of freedom, right? Yeah. With respect to the number of elements, obviously, so it's directly proportional to the number of degrees of freedom. And then, of course, assembling is embarrassingly prone to parallel strategies, because you can distribute your assembling process over as many processes as you like. So it should be very easy to lower the cost of this phase, I think. I don't want to take too much. It should be very scalable. But we also think in terms of the overall resource that we are using. Other question? I have just one quick. Can you give us a little bit more insight about the lack of physics for the Jupiter case? Who was the... Ah, yeah. So, Jupiter... It's Saturn. It's actually Saturn. So on Earth, we have some heating from the surface, right? Because we have a solid surface and the radiation heats the surface and this causes convection. So on Saturn, you don't have a solid surface. It's a gaseous planet. But you still do have convection. And that's because there is some heat going from the interior. So actually, Saturn has a positive energy budget. It still radiates more energy than it receives from the Sun. So it's a belief. There are some arguments that in order to get the zonal... The jet at the equator in the correct direction, you need to have some input of energy from the inside of the planet. So here it's not there. And I don't think I've shown... So maybe I've shown that. The jet... No, I haven't. The jet is in the wrong direction at the equator. And I don't get the number of degrees of freedom that you have in the geopotential direction for your... So typically, for this one, eighth-degree simulation, I think we are about one million cells in the horizontal and a few tens in the vertical. When we discuss at the canum, we talk about parallelization in time. Is it still something you could imagine? Well, I think that's something very interesting. It will take a long time until it gets operational. But that's something that several people are looking at. Because probably we will reach the point where we have used all the available parallelism in space at some point in the next decade or so. I don't know exactly when. But so if it's feasible, and there are many obstacles, maybe not in the fluid solver, but in everything else, if it's feasible, then it would be a very interesting direction to use. So personally, I'm interested in thinking at least about it. Maybe the last one question for me. When we used to do energy conservation time discretization, it's most of the time implicit. Are you also... Yeah, so I didn't talk about time discretization at all. So in time, we use standard explicit rule-scruta time stepping, which means that even if the spatial discretization is energy conserving, you lose that once your discretization time. But you still get some stability properties. Like, you can show that linear stability is still there. And since in time, we are fully explicit. So we are limited by a current condition based on horizontal sand waves. So the time stepping is relatively small. So with a say third or fourth order runch-scruta, you have very small time discretization errors. Okay. Thank you very much. Other questions? Okay. So thanks again, Thomas. Thank you.
Institut Pierre Simon la Place Federation, and performance. I will explain how we designed an atmospheric flow solver based on mimetic finite with similar ideas but based on finite elements. we are interested in these long time scales we to do this at the expense of the quality of the will talk only about the inviscid part of the and say gravity, so that would be if the Earth consequences on the structure of the equations we equations, but they are slightly different. And do centered finite differences and have easily properties. And for a while, we've had finite motion. So especially, we will conserve energy. easy to do this for, say, total mass, total mass to conserve energy at the discrete level in the that our equations of motion have a variational internal and potential energy, then I can write difficult. Another approach is to use the expressions from the equations of motion. And a I will find my Hamiltonian formulation. And conserve energy. And we just need to discretize. involved in the equation of motion. So this idea, current machines. Okay, so I will finish with a discrete integration by parts property. parts in a fine element context, too, to conserve energy. So the conservation of energy is not expensive... So we know our FEM method will be more expensive. The question is whether it's find difference methods. And also probably we matrices, should be negligible compared to the processes as you like. So it should be very easy