From 4f7c95f2877a67d4af527957709f3b0746c95bc8 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Sun, 7 Jan 2018 02:37:07 +0100 Subject: [PATCH] Extended step-26 documentation with regard to positivity preservation. --- examples/step-26/doc/results.dox | 87 ++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/examples/step-26/doc/results.dox b/examples/step-26/doc/results.dox index 6c45c86fe7..da1c4dccc8 100644 --- a/examples/step-26/doc/results.dox +++ b/examples/step-26/doc/results.dox @@ -120,3 +120,90 @@ current context we at least know what the right hand side is), and we should allow for more than one refinement level. Of course, all of this can be done using deal.II, it just requires a bit of algorithmic thinking in how to make this work! + + +

Positivity preservation

+ +To increase the accuracy and resolution of your simulation in time, one +typically decreases the time step size $k_n$. If you start playing around +with the time step in this particular example, you will notice that the +solution becomes partly negative, if $k_n$ is below a certain threshold. +This is not what we would expect to happen (in nature). + +To get an idea of this behavior mathematically, let us consider a general, +fully discrete problem: +@f{align*} + A u^{n} = B u^{n-1}. +@f} +The general form of the $i$th equation then reads: +@f{align*} + a_{ii} u^{n}_i &= b_{ii} u^{n-1}_i + + \sum\limits_{j \in S_i} \left( b_{ij} u^{n-1}_j - a_{ij} u^{n}_j \right), +@f} +where $S_i$ is the set of nearest neighbors to node $i$. If all coefficients +fulfill the following conditions: +@f{align*} + a_{ii} &> 0, & b_{ii} &\leq 0, & a_{ij} &\leq 0, & b_{ij} &\geq 0, + & + \forall j &\in S_i, +@f} +all solutions $u^{n}$ keep their sign from the previous ones $u^{n-1}$, and +consequently from the initial values $u^0$. See e.g. +Kuzmin, Hämäläinen +for more information on positivity preservation. + +Depending on the PDE to solve and the time integration scheme used, one is +able to deduce conditions for the time step $k_n$. For the heat equation with +the Crank-Nicolson scheme, +Schatz et. al. have +translated it to the following ones: +@f{align*} + (1 - \theta) k a_{ii} &\leq m_{ii},~ \forall i, + & + \theta k \left| a_{ij} \right| &\geq m_{ij},~ j \neq i, +@f} +where $M = m_{ij}$ denotes the mass matrix and $A = a_{ij}$ the stiffness +matrix with $a_{ij} \leq 0$ for $j \neq i$, respectively. With +$a_{ij} \leq 0$, we can formulate bounds for the global time step $k$ as +follows: +@f{align*} + k_{\text{max}} &= \frac{ 1 }{ 1 - \theta } + \min\left( \frac{ m_{ii} }{ a_{ii} } \right),~ \forall i, + & + k_{\text{min}} &= \frac{ 1 }{ \theta } + \max\left( \frac{ m_{ij} }{ \left|a_{ij}\right| } \right),~ j \neq i. +@f} + +Here, it is worth mentioning that the time step is constrained by both a lower +and upper bound in case of a Crank-Nicolson scheme. These bounds should be +considered along with the CFL condition to ensure significance of the performed +simulations. + +To ensure positivity preservation in this particular tutorial, we can use +the SparseMatrix objects for both mass and stiffness that are created via +the MatrixCreator functions. Iterating through each entry via SparseMatrixIterators +lets us check for diagonal and off-diagonal entries to set a proper time step +dynamically. For quadratic matrices, the diagonal element is stored as the +first member of a row (see SparseMatrix documentation). A exemplary code +snippet on how to grab the entries of interest from the mass_matrix +is shown below. + +@code +Assert (mass_matrix.m() == mass_matrix.n(), ExcNotQuadratic()); +const unsigned int& num_rows = mass_matrix.m(); +double mass_matrix_min_diag = 1e100, + mass_matrix_max_offdiag = 0.; + +SparseMatrixIterators::Iterator row_it (&mass_matrix, 0); + +for(unsigned int m = 0; mvalue(), mass_matrix_min_diag); + ++row_it; + + // check the off-diagonal elements + for(; row_it != mass_matrix.end(m); ++row_it) + mass_matrix_max_offdiag = std::max(row_it->value(), mass_matrix_max_offdiag); +} +@endcode -- 2.39.5