From: bangerth Date: Thu, 21 Aug 2008 04:34:38 +0000 (+0000) Subject: Comment on stability issues. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8e016365ab4ab3a7d91209a3e86f6d7b0786460b;p=dealii-svn.git Comment on stability issues. git-svn-id: https://svn.dealii.org/trunk@16619 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index a50ff82ba9..33ce876a6f 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -235,7 +235,9 @@ Courant-Friedrichs-Lewy (CFL) condition just like it was in @ref step_21 "step-21". (We wouldn't have had that stability condition if we treated the advection term implicitly since the BDF-2 scheme is A-stable, at the price that we needed to build a new temperature matrix at each time -step.) In particular this CFL condition means that the time step +step.) We will discuss the exact choice of time step in the results section, but for the moment of importance is that +this CFL condition means that the time step size k may change from time step to time step, and that we have to modify the above formula slightly. If $k_n,k_{n-1}$ are the time steps sizes of the current and previous time step, then we use the @@ -434,16 +436,17 @@ follows: @f} Here, $\beta$ is a stabilization constant (a dimensional analysis -reveals that it is unitless and therefore independent of scaling) and +reveals that it is unitless and therefore independent of scaling; we will +discuss its choice in the results section) and $c(\mathbf{u},T)$ is a normalization constant that must have units $\frac{m^{\alpha-1}K^\alpha}{s}$. We will choose it as $c(\mathbf{u},T) = - \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T) + c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T) \ |\mathrm{diam}(\Omega)|^{\alpha-2}$, where $\mathrm{var}(T)=\max_\Omega T - \min_\Omega T$ is the range of present temperature values (remember that buoyancy is driven by temperature -variations, not the absolute temperature). -To understand why this method works consider this: If on a particular +variations, not the absolute temperature) and $c_R$ is a dimensionless +constant. To understand why this method works consider this: If on a particular cell $K$ the temperature field is smooth, then we expect the residual to be small there (in fact to be on the order of ${\cal O}(h_K)$) and the stabilization term that injects artificial diffusion will there be diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox index cf98c5b102..5210a8fdb7 100644 --- a/deal.II/examples/step-31/doc/results.dox +++ b/deal.II/examples/step-31/doc/results.dox @@ -4,8 +4,146 @@

Numerical experiments to determine optimal parameters

-q1 vs q2 for temperature -q1/q1 stokes +The program as is has three parameters that we don't have much of a +theoretical handle on how to choose in an optimal way. These are: + +In all of these cases, we will have to expect that the correct choice of each +value depends on that of the others, and most likely also on the space +dimension and polynomial degree of the finite element used for the +temperature. Below we'll discuss a few numerical experiments to choose +constants. + + +
Choosing ck and β
+ +These two constants are definitely linked in some way. The reason is easy to +see: In the case of a pure advection problem, +$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \gamma$, any +explicit scheme has to satisfy a CFL condition of the form +$k\le \min_K \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. On the other hand, +for a pure diffusion problem, +$\frac{\partial T}{\partial t} + \nu \Delta T = \gamma$, +explicit schemes need to satisfy a condition +$k\le \min_K \frac{c_k^d h_K^2}{\nu}$. So given the form of $\nu$ above, an +advection diffusion problem like the one we have to solve here will result in +a condition of the form +$ +k\le \min_K \min \left\{ + \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}, + \frac{c_k^d h_K^2}{\beta \|mathbf{u}\|_{L^\infty(K)} h_K}\right\} + = + \min_K \left( \min \left\{ + c_k^a, + \frac{c_k^d}{\beta}\right\} + \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}} \right) +$. +It follows that we have to face the fact that we might want to choose $\beta$ +larger to improve the stability of the numerical scheme (by increasing the +amount of artificial diffusion), but we have to pay a price in the form of +smaller, and consequently more time steps. In practice, one would therefore +like to choose $\beta$ as small as possible to keep the transport problem +sufficiently stabilized while at the same time trying to choose the time step +as large as possible to reduce the overall amount of work. + +The find the right balance, the only way is to do a few computational +experiments. Here's what we did: We modified the program slightly to allow +less mesh refinement (so we don't always have to wait that long) and to choose +$ + \nu(T)|_K + = + \beta + \|\mathbf{u}\|_{L^\infty(K)} h_K +$ to eliminate the effect of of the constant $c_R$. We then run the program +for different values $c_k,\beta$ and observe maximal and minimal temperatures +in the domain. What we expect to see is this: If we choose the time step too +big (i.e. choose a $c_k$ bigger than theoretically allowed) then we will get +exponential growth of the temperature. If we choose $\beta$ too small, then +the transport stabilization becomes insufficient and the solution will show +significant oscillations but not exponential growth. Here is what we get for +$\beta=0.01, \beta=0.1$, and $\beta=0.5$, different choices of $c_k$, and +bilinear elements (temperature_degree=1) in 2d: + + + + + + + + + + + +
+ @image html "step-33.timestep.q1.beta=0.01.png" "" width=4cm + + @image html "step-33.timestep.q1.beta=0.03.png" "" width=4cm +
+ @image html "step-33.timestep.q1.beta=0.1.png" "" width=4cm + + @image html "step-33.timestep.q1.beta=0.5.png" "" width=4cm +
+ +The way to interpret these graphs goes like this: for $\beta=0.01$ and +$c_k=\frac 12,\frac 14$, we see exponential growth or at least large +variations, but if we choose +$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ +or smaller, then the scheme is +stable though a bit wobbly. For more artificial diffusion, we can choose +$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ +or smaller for $\beta=0.1$, and again need +$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ +for $\beta=0.5$ (this time because much diffusion requires a small time +step). + +So how to choose? If we were simply interested in a large time step, then we +would go with $\beta=...$ and +$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. +On the other hand, we're also interested in accuracy and here it may be of +interest to actually investigate what these curves show. To this end note that +we start with a zero temperature and that our sources are positive — so +we would intuitively expect that the temperature can never drop below +zero. But it does, a consequence of Gibb's phenomenon when using continuous +elements to approximate a discontinuous solution. We can therefore see that +choosing $\beta$ too small is bad: too little artificial diffusion leads to +over- and undershoots that aren't diffused away. On the other hand, for large +$\beta$, the minimum temperature drops below zero at the beginning but then +quickly diffuses back to zero. + +On the other hand, let's also look at the maximum temperature. Watching the +movie of the solution, we see that initially the fluid is at rest. The source +keeps heating the same volume of fluid whose temperature increases linearly at +the beginning until its buoyancy is able to move it upwards. The hottest part +of the fluid is therefore transported away from the solution and fluid taking +its place is heated for only a short time before being moved out of the source +region, therefore remaining cooler than the initial bubble. If $\kappa=0$ +(in the program it is nonzero but very small) then the hottest part of the +fluid should be advected along with the flow with its temperature +constant. That's what we can see in the graphs with the smallest $\beta$: Once +the maximum temperature is reached, it hardly changes any more. On the other +hand, the larger the artificial diffusion, the more the hot spot is +diffused. Note that for this criterion, the time step size does not play a +significant role.

Possible extensions