Time integration methods for finite element discretizations in weather forecasting
This is a modal window.
Das Video konnte nicht geladen werden, da entweder ein Server- oder Netzwerkfehler auftrat oder das Format nicht unterstützt wird.
Formale Metadaten
| Titel |
| |
| Serientitel | ||
| Anzahl der Teile | 17 | |
| Autor | ||
| Lizenz | CC-Namensnennung 3.0 Deutschland: Sie dürfen das Werk bzw. den Inhalt zu jedem legalen Zweck nutzen, verändern und in unveränderter oder veränderter Form vervielfältigen, verbreiten und öffentlich zugänglich machen, sofern Sie den Namen des Autors/Rechteinhabers in der von ihm festgelegten Weise nennen. | |
| Identifikatoren | 10.5446/40496 (DOI) | |
| Herausgeber | ||
| Erscheinungsjahr | ||
| Sprache |
Inhaltliche Metadaten
| Fachgebiet | ||
| Genre | ||
| Abstract |
|
Leibniz MMS Days 20197 / 17
00:00
Finite-Elemente-MethodeSpezifisches VolumenMinimumGüte der AnpassungEntscheidungstheorieElement <Gruppentheorie>MomentenproblemDiagrammVorlesung/Konferenz
00:20
Element <Gruppentheorie>Desintegration <Mathematik>UnendlichkeitMachsches PrinzipPrognoseverfahrenFinite-Elemente-MethodeExponentialgleichungt-TestMultiplikationsoperatorInelastischer StoßIntegral
01:01
PrognoseverfahrenElement <Gruppentheorie>Numerische MathematikMathematikVariationsprinzipStandardabweichungDualitätstheorieGlattheit <Mathematik>GitterverfeinerungIkosaederKugelUnendlichkeitVolumenGeschwindigkeitDruckverlaufVertikaleFrequenzFunktion <Mathematik>PolynomKomplex <Algebra>Stetige FunktionTensorProdukt <Mathematik>Numerische MathematikFinite-Elemente-MethodeStichprobenumfangMomentenproblemEntscheidungstheorieMeterMereologieSchlussregelElement <Gruppentheorie>ExponentialgleichungAffiner RaumGeschwindigkeitCoxeter-GruppeDruckverlaufFunktionalFlächeninhaltMultiplikationsoperatorRichtungt-TestPhysikalisches SystemMinkowski-MetrikZusammenhängender GraphStetige FunktionGüte der AnpassungBasis <Mathematik>Stützpunkt <Mathematik>AggregatzustandZentrische StreckungRuhmasseDifferenteProjektive EbeneKategorie <Mathematik>Nichtlineares GleichungssystemMathematikVorhersagbarkeitStabilitätstheorie <Logik>BestimmtheitsmaßBewegungsgleichungModulformOrtsoperatorKartesische KoordinatenFraktalgeometrieStandardabweichungDifferenzenrechnungDivergenz <Vektoranalysis>GravitationModelltheorieIntegralTensorproduktPolynomZweiKompakter RaumMinimalgradBoussinesq-ApproximationLebesgue-IntegralComputeranimation
08:04
Funktion <Mathematik>KonstanteLineare GleichungARCH-ProzessGeschwindigkeitMinkowski-MetrikDruckverlaufUnabhängige MengeUnendlichkeitDifferentialNichtlineares GleichungssystemExponentialgleichungElement <Gruppentheorie>FunktionalMereologieDruckverlaufIterationMultiplikationsoperatorGeschwindigkeitEinbettung <Mathematik>Zusammenhängender GraphAusreißer <Statistik>ModelltheoriePhysikalisches SystemGerichteter GraphRechter WinkelMinkowski-MetrikLebesgue-IntegralPunktTensorproduktExponentialgleichungDifferenteMathematikerinStetige FunktionStabilitätstheorie <Logik>Algebraische StrukturDistributionenraumZentrische StreckungGesetz <Physik>KonstanteMomentenproblemRichtungNichtlineares GleichungssystemGeradeRuhmasseDifferentialgleichungLineare GleichungHeegaard-ZerlegungTermFinite-Elemente-MethodeKategorie <Mathematik>Reelle ZahlSignifikanztestVektorraumFinitismusSpezifisches VolumenQuadratische FunktionStichprobenfehlerDivergenz <Vektoranalysis>Gewöhnliche Differentialgleichung
14:57
Matrizenrechnungp-BlockDiagonale <Geometrie>InverseSchwach besetzte MatrixZeitbereichKartesisches ProduktGravitationWellenlehreMachsches PrinzipGreen-FunktionTermSchwach besetzte MatrixMultiplikationsoperatorTermNichtlineares GleichungssystemEigenwertproblemMatrizenrechnungPhysikalisches SystemFinite-Elemente-MethodeStabilitätstheorie <Logik>RuhmasseMetrisches SystemExponentialgleichungSortierte LogikZweiPunktInverseAnalysisNichtlinearer OperatorForcingAlgebraische StrukturSupremum <Mathematik>ZeitbereichDifferenteGeschwindigkeitMinkowski-MetrikMathematikElement <Gruppentheorie>Lokales MinimumAbstimmung <Frequenz>Berechenbare FunktionEntscheidungstheorieIndexberechnungGruppenoperationRechter WinkelDiagonale <Geometrie>MengenlehreRichtungHypermatrixMatrixnormStörungstheorieHeegaard-ZerlegungInvertierbare MatrixParametersystemDruckverlaufMereologie
21:50
Vorlesung/Konferenz
Transkript: Englisch(automatisch erzeugt)
00:01
My first attempts to go from finite volumes to finite element discretization because it's, at the moment, an actual topic. And in my last days at our research institute, I want to learn a little bit about that to maybe do
00:21
something also later. And I'd ask this together with a bachelor's student, Catherine Lubashevsky. So it's a bachelor's student, so it does only first steps. And I want to explain a little bit what we have implemented and what are our experience.
00:41
So I want to give a short introduction and what I learned about finite elements, what we've done at the end that we only apply it to a linear problem. And then my special interests are time integration methods. And I want to explain what are the problems with this time
01:00
integration methods if we are confronted with finite element discretizations, some numerical example, and that's it. So at the moment, I have here listed some ideas what happens at the moment. And I think the mathematics for numerical weather prediction
01:23
reaches, in some sense, a really new frontiers. And there, you can look in different directions in this research. One is that you are going from standard latitude, longitude
01:43
grids, to unstructured grids. OK, unstructured grids we already have, maybe also already 20 years ago or something like that, or 30 years ago. But I think in the last 10 years together with that, we are changing from fully composable or non-hydric systems that these unstructured grids also
02:03
are getting some new interest. And we go to these grids in different directions. So we have different type of grids. And at the same time, we want grid refinement.
02:20
And we want to have also on this grid refined things. And for different grids, we want to have some type of unified methods. And we want to have stable methods and methods which preserve as much as possible from the continuous model.
02:40
And we also go with these grids to more sophisticated spatial subsidization methods. So we are looking for finite differences. And on this unstructured grid, we are going to something which people in other areas call mimetic or compacted finite element pairs.
03:03
Or we are looking for another type of presentation of the equations of motions like in the Hamilton or Nambu formulation and going then from this continuous formulation to a discrete one. And from that, we get something about the discretization.
03:22
And that we don't only do in space. We do only also the same thing in time. And there are now at the moment also appear a lot of things about new time discretization methods which follow from this Hamiltonian formulation.
03:42
And this year, for instance, there is a special workshop about this so-called variational integrators at the end of the year. And you see that at all the time, I think, in numerical weather prediction, there are now people are sitting which do deep, deep mathematics
04:03
in some sense. So what now I've done is no deep mathematics. So what we've done with the bachelor student is that we want to solve these linearized equations, which is a Boussinesq type equation where it's two dimensional.
04:24
You have two velocity components. We have a pressure component. And we have a buoyancy component. And we have in the system three time scales. So it's an advection time scale. We have this gravity time scale.
04:43
And we have this acoustic time scale, all the time scales which we also have in the full compositional model. And there was a paper, I think, last of this year in quarterly of mathematical society, where the people from the Gung Ho project
05:04
presented things about finite elements for this application. And we follow in the bachelor work to re-implement that and then to see how our time integration methods can be used for that. So here already I have a split equation in two parts.
05:27
So we have this low time scale. Think about maybe 20 meters per second. And you think about this fast time scale sound speed, 260 meters per second. And this one is somewhere in between.
05:43
And I split already the equation here in this linearized version in these two parts, which we also have later on. I call this the slow part and the fast part. So finite elements. So what we've done is we want to solve this simple equation
06:01
on a Cartesian grid. So we start with a student with simple one-dimensional finite elements and then go to tensor product elements and implement that as 41. And I couldn't also go into detail also not with the bachelor's student, because you have that knowledge
06:23
about all that, as I usually also not have in full case. We started with two types of finite elements. So on the one side, we use continuous finite elements, which are cell-wise or polynomial function
06:44
degrees of r plus 1. And we denote that by phi C r. And we have discontinuous elements, which are cell-wise polynomial functions, which also are of degree r plus 1. And we denote them by dG r.
07:01
And then we take spaces. We have to choose taken spaces for the two velocity components. We have to take a space for the buoyancy. And we have to take one for the pressure. And that what people first propose a little bit
07:20
is that you have to take three of these finite element spaces. You have to take one which is continuous and where you have here a quadratic elements. Then you have to take a space which is continuous.
07:40
And in the one-dimensional case, you have linear elements. And then that you have to take for the buoyancy, that for the velocity and for the pressure, you have to take this continuous space. And these three spaces then has this property that you can take the gradient of that. You come in here and take the divergence on one
08:02
to the equation again, you will end up somewhere here. But now we have to go to 2D. And in 2D, we have to choose the buoyancy, the velocity, and the pressure. And now for the buoyancy, we take this space here
08:21
as a tensor product. But for the velocity components, we have a tensor product which is linear for you. In this direction, we take a continuous space. And on the other side, we take a discontinuous space which one order less.
08:41
For v, it's vice versa. And for p again, we take a tensor product of two discontinuous spaces. And then these spaces, if you take it again, you have this type of embedding property. So you can, starting with the buoyancy, you compute this type of a curl.
09:03
You end up in this vector finite element space. And then you make divergence, and you end up in this space. So now, this is more know people. So what we have, we make only experiments
09:21
with low-order ansatz functions. So we have piecewise constant ansatz function, piecewise linear ansatz functions. And we have to work also with, oh, this is in German here, piecewise quadratic functions. You have to implement that all. So we have learned a little bit how
09:40
to do that in some efficient way. But it's also not easy if you go to finite elements. And here, you see how these finite elements look like. At the end, if you want to go to finite elements, you have to go from this to a weak formulation.
10:02
And we have to choose our trial functions, either the velocity components, and the pressure, and the b from this ansatz function, which I have showed. And again, we take test functions, which are from the same ansatz functions.
10:23
We multiply with these test functions, make here for the pressure term some partial integrations. And after that, we are looking for trial functions in the finite element space, which are fulfilled for all test functions, this equation.
10:46
And if you then compute all. So I want not to explain in detail how we compute the advection term, especially for the pressure, or here for it's discontinuous, or here for v. So we have to do something special.
11:04
At the end, we end up with a linear system. So we have a system of ordinary differential equations. And I have written again here the advection part. So if you cut this part, then you
11:22
see all the equations are each equation you can solve by itself. And we have here the pressure term, we have the buoyancy term, here's the acoustic term, and here we have again a buoyancy term. And we have to solve that equation in time.
11:43
So now comes in my favorite methods, which I use in the finite volume context, which are split explicit methods. So for the simple models, I want to solve this with fully explicit models, but with different time steps because of these different time scales.
12:03
So that's it. So we are going back to the system here. So I have a slow part and a fast part. And we have developed sometimes a combined special type
12:24
of Runge-Kutta methods, which we use for the slow part. So think that this fast part is 0. Then what we do is in each stage, we solve again a differential equation.
12:41
In our case, all equations are linear. But here, this part is frozen, so it's fixed. And only this part, which is the fast part, is a real right-hand side function, which we can compute in each point. And we can solve this linear system with this fixed source term.
13:02
And then at the end, we get the new stage value. We can compute again the slow part can form again this type of a right-hand side, solve this linear equation. And we have here maybe three stages with a fixed time step.
13:22
And then we take smaller time steps to solve this system with a special method. And if you are looking now for this fast part here, then we will see I have a little bit simplified.
13:42
If you take only momentum and pressure, that this fast part has this symplectic structure. So you see u depends on p, and p depends on u. And then what mathematicians call a symplectic Euler, I think, is a small error.
14:02
So we solve the first equation with a forward Euler. And then here should be m plus 1. Then we solve the second equation with a backward Euler. This is a stable method. And there is also something which
14:21
is a stabilized symplectic Euler, which is also called in the better literature something like of divergence dumping. So we make here again a forward Euler step and replace here this forward Euler step also
14:42
with some leg term that we go a little bit ahead. And this stabilize the whole procedure. That is also not fully understand at the moment. So now we are going to apply the split explicit methods.
15:00
And if we want to apply them, we have to bring the mass matrix, which is in front. What you have seen, if you work with finite elements, you have to bring this to the right hand side. So for each computation of both of these terms, you have to solve in some sense the linear systems. You can make before an LU decomposition,
15:21
but if you go to 3D, it's a minor problem. Now what people do in the finite element literature, they do some mass lumping. They replace this matrix by a diagonal or a block diagonal matrix, which I call here ML.
15:41
So that you have now further options. You replace both of them with this lump matrix or the other one. You only see because that I compute on a large time step. Here I take the exact mass matrix that I do on a small time step.
16:02
So I take the lump matrix where I have to solve a simple linear system. So it's diagonal or block diagonal. No problem. And there's a third idea. So what I do now, I subtract something. I take here an intermediate. I subtract that here from the full formulation,
16:25
and I add it again. So now this is something which I call my slow term. I have to compute it. And this is my fast term, but I have to compute it often. And we try that first.
16:40
And we see, OK, this was of them. OK, they will work with the split explicit methods, but this idea doesn't work. So it was unstable. We also not then I make some eigenvalue analysis and see that at least if you see that this operator has some eigenvalues which
17:00
are sitting not on the left hand side, and this makes the problems. And now the idea was not to take here this simple mass lump matrix, but to go to a matrix where we say we are looking for something which is an inverse matrix of the mass matrix
17:20
but has a sparse structure, something like an approximate inverse. And that we work here, that we add this term here. And again, so you don't change the equation. You only change what is the slow part and what is the fast part. And for this inverse matrix, we tested two ideas.
17:42
The first one is we are taking an M inverse matrix, which is the minimizer of this Frobenius norm. And M inverse has a prescribed sparsity pattern. So you couldn't fulfill this, so the minimum is not zero.
18:05
Or the other idea is that we compute the inverse and then take only a sparse pattern of this inverse matrix. And this sparsity pattern is usually the same as the sparsity pattern of M,
18:23
something like that. And with this idea, we make some experiments. So we take a Cartesian domain. So 150 kilometer to the east and to the west,
18:42
and 10 kilometers in the height, 3,000 seconds. And we start with an atmosphere at west, where the velocity and this pressure perturbation is zero. And the buoyancy has this type of a special structure. And this buoyancy perturbation then
19:02
will give us some movement then. So here are the parameters. And we take a spatial discretization. So it's in the x direction. And also in the z direction, it's one kilometer. So in the vertical direction, it's very coarse.
19:23
And we take a slow time step of 20 seconds. And we have to take something like 19 small time steps per large time step to fulfill the stability criteria for this. And this is the buoyancy at the starting point.
19:45
And then after 3,000 seconds, if you also include some horizontal movement here, you get this type of picture. And we are now looking somewhere here in the middle of this buoyancy in the case
20:01
where we have no movement to see what happens with this different mass lumping procedures. So I have plotted here now three type of lumping procedures. So I have first these two are the same.
20:21
So blue is we make no lumping, what we call the exact solution in some sense. Then we make a full lumping. So we have for both terms, we have the lumping matrix. And then we make green, where we have only this lumping on the first term. And there, I said all three things will work.
20:40
And then now, a little bit, and now we are going. Maybe this is the last one. Now what we now do is that we go to this sparse idea that we subtract a fast term.
21:03
And in the case that we take the inverse, which is our sparse, then you see first if we do it without this stabilization, we have here some problems. And then if you make also this small stabilization in this fast, what I call this stabilized split explicit,
21:21
then you see that the exact solution and that what we get here, you see now no differences. So this was this experiment. And that was something like the outcome of this bachelor work and what I want to present here. Now I also have already a finite element
21:42
code in 2D for the fully composable Euler equations. But no tests with that. So that's what I want to say. Thank you.