diff --git a/.buildinfo b/.buildinfo index d56e401..a19252e 100644 --- a/.buildinfo +++ b/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 2fc1b2936d785829fdbc462f16cc33bc -tags: fbb0d17656682115ca4d033fb2f83ba1 +config: a38950419e00abfe3c37662591290520 +tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/_sources/advection_diffusion.rst.txt b/_sources/advection_diffusion.rst.txt new file mode 100644 index 0000000..4b9b7ae --- /dev/null +++ b/_sources/advection_diffusion.rst.txt @@ -0,0 +1,215 @@ +The advection-diffusion-reaction equation +----------------------------------------- + +The advection-diffusion-reaction equation (also called the continuity equation in semiconductor physics) in flux form, is given by, + +.. math:: + u_t = (\mathcal{-F})_x + s(t,x,u) + +where :math:`\mathcal{F} = au - du_x`. This expression holds over the whole domain, therefore it must also hold over a subdomain :math:`\Omega_j=(x_{j-1/2}, x_{j+1/2})`. Integrating over the subdomain, or cell, gives, + +.. math:: + \int_{x_{j-1/2}}^{x_{j+1/2}} u_t~dx = \int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx + \int_{x_{j-1/2}}^{x_{j+1/2}} s(x,t,u)~dx + +If we now use :math:`w` and :math:`\bar{s}` to represent the cell averages of the solution variable and reaction term then we don't have to perform any actual integrations. The fundamental quantity we will use in the discretised equations is a cell average, + +.. math:: + w_j^{\prime} = -\frac{\mathcal{F}_{j+1/2}}{h_j} + \frac{\mathcal{F}_{j-1/2}}{h_{j}} + \bar{s}_j + +.. note:: + In the above we *can* perform the integration of the flux term. Integrating the divergence of the flux over the unit cell simply gives use the flux at the cell faces, :math:`\int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx = \left[ \mathcal{-F} \right]_{x_{j-1/2}}^{x_{j+1/2}}`. Because we want the cell average we divide by the cell width :math:`h_j`. In three dimension we would divide by the volume of the cell, hence the name, finite volume method. + +For the advection component of the flux we will use a linear interpolation (see :ref:`label-interpolation-aside`) to determine the contribution at the cell faces, and for the diffusion component with will use a simple cell average, + +.. math:: + \mathcal{F}_{j+\frac{1}{2}} = a_{j+\frac{1}{2}}\left( \frac{h_{j+1}}{2h_{+}} w_j + \frac{h_j}{2h_{+}} w_{j+1} \right) - d_{j+\frac{1}{2}} \frac{w_{j+1}-w_j}{h_{+}} + +.. math:: + \mathcal{F}_{j-\frac{1}{2}} = a_{j-\frac{1}{2}}\left( \frac{h_{j}}{2h_{-}} w_{j-1} + \frac{h_{j-1}}{2h_{-}} w_{j} \right) - d_{j-\frac{1}{2}} \frac{w_{j}-w_{j-1}}{h_{-}} + +The meaning of the "h" terms is shown in the figure below, + +.. figure:: img/cells_FVM.png + :scale: 100 % + :alt: Finite volume cells with important distances and positions labelled. + :align: center + + Finite volume cells with important distances and positions labelled. + +Substituting the definition of the fluxes into the above equation yields, + +.. math:: + w_j^{\prime} = \frac{w_{j-1}}{h_j} \left( \frac{a(x_{j-\frac{1}{2}}) h_j}{2h_{-}} + \frac{d(x_{j-\frac{1}{2}})}{h_{-}}\right) + \frac{w_j}{h_j}\left( \frac{a(x_{j-\frac{1}{2}})h_{j-1}}{2h_{-}} - \frac{a(x_{j+\frac{1}{2}})h_{j+1}}{2h_{+}} - \frac{d(x_{j-\frac{1}{2}})}{h_{-}} - \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right) + \frac{w_{j+1}}{h_j} \left( \frac{-a(x_{j+\frac{1}{2}})h_j}{2h_{+}} + \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right) + :label: eq1 + +This equation is the *semi-discrete* form because the equation has not been discretised in time. + +For convenience we will simply write the r.h.s. of :eq:`eq1` as, + +.. math:: + w_j^{\prime} = F(w) + +where it is understood that :math:`w = (w_{j-1}, w_j, w_{j+1})` stands for the discretisation stencil. + +Adaptive upwinding & exponential fitting +**************************************** + +Hundsdorfer (on pg. 266) notes that a adaptive upwinding scheme can be introduced by altering the diffusion coefficient, + +.. math:: + d(x_{j}) \rightarrow d(x_{j}) + \frac{1}{2}\kappa_{j} h_{j} a(x_{j}) + + +.. math:: + \kappa=\text{sgn}(a(x_{j})). + +Here the nature of the discretisation is defined by the sign of the :math:`a`, when :math:`a\neq0` the scheme results in an upwind discretisation. Only if :math:`a=0` (corresponding to the diffusion equation limit) does the scheme become based on a central difference (see the red line in the Figure below). A generalisation is possible by choosing a function for :math:`\kappa` which automatically adjusts the discretisation method depending on the local value of Peclet's number, :math:`\mu=ah/d`. Exponential fitting has shown to be a robust scheme for steady-state approaches. Exponential fitting can be introduced by, + +.. math:: + \kappa = \frac{e^{\mu}+1}{e^{\mu}-1} - \frac{2}{\mu} + +An approximation to the above, which does not require the evaluation of exponentials is commonly used, + +.. math:: + \kappa = \begin{cases} + \text{max}(0, 1-2/\mu) & \text{when}~ \mu>0 \\ + \text{min}(0, -1-2/\mu) & \text{when}~ \mu<0 + \end{cases} + +.. figure:: img/Exponential_fitting_and_approximation.png + :scale: 100 % + :alt: Exponential fitting and approximation + :align: center + + +This is a very elegant want to introduce adaptive upwinding or exponential fitting because it can be brought into the discretisation in an ad hoc manner. Note that for all of the adaptive expressions :math:`\kappa\rightarrow\pm1` as :math:`\mu\rightarrow\pm\infty`, and :math:`\kappa\rightarrow 0` as :math:`\mu\rightarrow 0`. This means the discretisation will automatically weight in the favour of the upwind scheme in locations where advection dominates. Conversely, where diffusion dominates the weighting factor shifts in favour of a central difference scheme. + +Explicit and implicit forms +*************************** + +In the discussion so far we have been ignoring one important deal: at what *time* is the r.h.s of the discretised equation evaluated? + +Moreover, if we choose the r.h.s to be evaluated at the *present* time step, :math:`t_n` this is known as an *explicit* method, + +.. math:: + w_j^{\prime} = F(w^n) + +Explicit methods are very simple. Starting with initial conditions at time :math:`t_n`, the above equation can be rearranged to find the solution variable :math:`w_j^{n+1}` at the future time step, :math:`t_{n+1}`. + +However the downside of using explicit methods is that they are often numerically unstable unless the time step is exceptionally (sometime unrealistically) small. + +Fortunately there is an second alternative, we can choice to write the r.h.s of the equation at the *future* time step, :math:`t_{n+1}`, this is known as an *implicit* method, + +.. math:: + w_j^{\prime} = F(w^{n+1}) + +Explicit methods are significantly more numerically robust, however they pose a problem, how does one write the equations because the solution variable at the future time step :math:`w^{n+1}` is unknown? The answer is that at each time step we must solve a linear system of equation to find :math:`w^{n+1}`. + +.. note:: + For the **linear** advection-diffusion-reaction equation implicit methods are simply to implement even though the computation cost is increases. One must simply write the equation in the linear form :math:`A\cdot x = d` and solve for :math:`x` which is the solution variable at the future time step. However if the equations are non-linear then implicit methods pose problem because the equation **cannot** be written in linear form. In these situations are there are a number of techniques that are used but they all most use a iterative procedure, such as a Newton-Raphson method, to solve the equations. + +The following section assumes that the equation is linear. + +The :math:`\theta`-method +************************* + +The :math:`\theta`-method is an approach which combines implicit and explicit forms into one method. It consists of writing the equation as the average of the implicit and explicit forms. If we let :math:`F_{w^{n}}` and :math:`F_{w^{n+1}}` stand for the r.h.s of the explicit and implicit forms of the equations then the :math:`\theta`-method gives, + +.. math:: + w_j^{\prime} = \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n}) + +Setting :math:`\theta=0` recovers a fully implicit scheme while :math:`\theta=1` gives a fully explicit discretisation. The value of :math:`\theta` is not limited to just 0 or 1. It is common to set :math:`\theta=1/2`, this is called the Crank-Nicolson method. It is particularly popular for diffusion problem because it preserves the stability of the implicit form but also increases the accuracy of the time integration from first to second order (because two points in time are being averaged). For advection diffusion problems the Crank-Nicolson method is also unconditionally stable. + +In the above equation, + +.. math:: + F(w_j^{n}) = r_a w_{j-1}^{n} + r_b w_{j}^{n} + r_c w_{j+1}^{n} \\ + F(w_j^{n+1}) = r_a w_{j-1}^{n+1} + r_b w_{j}^{n+1} + r_c w_{j+1}^{n+1} + +and the coefficients are, + +.. math:: + r_a & = \frac{1}{h_j} \left( \frac{a(x_{j-\frac{1}{2}}) h_j}{2h_{-}} + \frac{d(x_{j-\frac{1}{2}})}{h_{-}}\right) \\ + r_b & = \frac{1}{h_j}\left( \frac{a(x_{j-\frac{1}{2}})h_{j-1}}{2h_{-}} - \frac{a(x_{j+\frac{1}{2}})h_{j+1}}{2h_{+}} - \frac{d(x_{j-\frac{1}{2}})}{h_{-}} - \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right)\\ + r_c & = \frac{1}{h_j} \left( \frac{-a(x_{j+\frac{1}{2}})h_j}{2h_{+}} + \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right) + +We have written the coefficients without dependence on time to simplify the notation, but the coefficients should be calculated at the same time points as their solution variable. However, the coefficients must be linear, they should not depend on the values of the solution variable. + +Discretised equation in matrix form +*********************************** + +Dropping the spatial subscripts and writing in vector form the equations becomes, + +.. math:: + \frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n}) + +where, + +.. math:: + F(w^{n+1}) = M^{n+1} w^{n+1} + s^{n+1} \\ + F(w^{n}) = M^{n}w^{n} + s^{n} + +with the coefficient matrix, + +.. math:: + M = + \begin{align} + \begin{pmatrix} + r_b & r_c & & & 0 \\ + r_a & r_b & r_c & & \\ + & \ddots & \ddots & \ddots& \\ + & & r_a & r_b & r_c \\ + 0 & & & r_a & r_b + \end{pmatrix} + \end{align} + +The :math:`r`-terms have be previously defined as the coefficients that result from the discretisation method. Note that the subscripts for the :math:`r`-terms have been dropped to simplify the notation, but they are functions of space. For example, terms in the first row should be calculated with :math:`j=1`, in the second with :math:`j=2` etc. Solving linear equations is discussed late in section XX. + +.. _label-interpolation-aside: + +Aside: Linear interpolation between cell centre and face values +========================================================================= + +In general, linear interpolation between two points :math:`(x_0, x_1)` can be used to find the value of a function at :math:`f(x)`, + +.. math:: + f(x) = \frac{x - x_1}{x_0 - x_1}f(x_0) + \frac{x - x_0}{x_1 - x_0}f(x_1) + +In a cell centred grid we know the value of the variable :math:`w` at difference points, :math:`w_j` and :math:`w_{j+1}`. We can apply the linear interpolation formulae above to determine value at cell face :math:`w_{j+1/2}`. + +.. math:: + w_{j+1/2} = \frac{x_{j+1/2} - x_{j+1}}{x_{j} - x_{j+1}} w_j + \frac{x_{j+1/2} - x_j}{x_{j+1} - x_j} w_{j+1} + +This can be simplified firstly by using function to represent the distance between cell centres, + +.. math:: + h_{-} = x_j - x_{j-1} \quad h_{+} = x_{j+1} - x_{j} + +to give, + +.. math:: + w_{j+1/2} = \frac{x_{j+1} - x_{j+1/2}}{h_{+}} w_j + \frac{x_{j+1/2} - x_j}{h_{+}} w_{j+1} + +This expression still contains :math:`x_{j+1/2}` which we can simplify further by using an expression for the position of cell centres, + +.. math:: + x_j = \frac{1}{2} \left( x_{j-\frac{1}{2}} + x_{j+\frac{1}{2}} \right) \quad x_{j+1} = \frac{1}{2} \left( x_{j+\frac{1}{2}} + x_{j+\frac{3}{2}} \right) + + +Note, this expression is still valid of non-uniform grids, it simply says that cell centres are always equidistant from two faces. Rearranging the above expression and substituting in for :math:`x_{j}` and :math:`x_{j+1}` terms gives, + +.. math:: + w_{j+1/2} = \frac{\frac{1}{2} \left( x_{j+\frac{1}{2}} + x_{j+\frac{3}{2}} \right) - x_{j+1/2}}{h_{+}} w_j + \frac{x_{j+1/2} - \frac{1}{2} \left( x_{j-\frac{1}{2}} + x_{j+\frac{1}{2}} \right)}{h_{+}} w_{j+1} + + +Finally, by defining the distance between vertices as, :math:`h_j = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}`, we can simplify to the following expression, + +.. math:: + w_{j+1/2} = \frac{h_{j+1}}{2h_{+}} w_j + \frac{h_j}{2h_{+}} w_{j+1} + + +Similarly the :math:`w_{j-1/2}` can be found, + +.. math:: + w_{j-1/2} = \frac{h_{j}}{2h_{-}} w_{j-1} + \frac{h_{j-1}}{2h_{-}} w_{j} + diff --git a/_sources/advection_diffusion_boundary_conditions.rst.txt b/_sources/advection_diffusion_boundary_conditions.rst.txt new file mode 100644 index 0000000..32763bc --- /dev/null +++ b/_sources/advection_diffusion_boundary_conditions.rst.txt @@ -0,0 +1,174 @@ +Boundary conditions for the advection-diffusion-reaction equation +----------------------------------------------------------------- + +The advection-diffusion-reaction equation is a particularly good equation to explore apply boundary conditions because it is a more general version of other equations. For example, the diffusion equation, the transport equation and the Poisson equation can all be recovered from this basic form. Moreover, by developing a general scheme for boundary conditions of the advection-reaction-diffusion equation we automatically get a system for imposing boundary conditions on all equations of a similar form. + +Robin boundary conditions (known flux) +************************************** + +Robin specify a known *total flux* comprised of a diffusion and advection component. Note that in the diffusion equation limit (where :math:`a=0`) these boundary conditions reduce to Neumann boundary conditions. + +`Within the finite volume method Robin boundary conditions are naturally resolved `_. This means that there is no need for interpolation or ghost point substitution (although these approaches remain possible) to include the boundary conditions because the flux at the boundary appears naturally in the semi-discretised equation. For example consider the semi-discretised equation evaluated at the boundary cell :math:`\Omega_1`, + +.. figure:: img/boundary_cell_FVM.png + :scale: 100 % + :alt: Finite volume boundary cell at the left hand side. + :align: center + + Finite volume boundary cell at the left hand side. + +.. math:: + w_1^{\prime} = -\frac{\mathcal{F}_{3/2}}{h_1} + \frac{\mathcal{F}_{1/2}}{h_{1}} + \bar{s}_1 + +The Robin boundary condition specifies the flux at :math:`x_{1/2}`, + +.. math:: + \mathcal{F}_{1/2} = g_{R}(x_{1/2}) + +Therefore the boundary condition can be incorporated without invoking any information regarding the ghost cell, + +.. math:: + w_1^{\prime} = \frac{w_1}{h_1}\left( \frac{-a(x_{3/2})h_{2}}{2h_{+}} - \frac{d(x_{3/2})}{h_{+}} \right) + \frac{w_{2}}{h_1} \left( \frac{-a(x_{3/2}) h_1}{2h_{+}} + \frac{d(x_{3/2})}{h_{+}} \right) + \frac{g_{R}(x_L)}{h_1} + \bar{s}_1 + +Similarly applying the same procedure to the :math:`\Omega_J` cell at the right hand side boundary, + +.. figure:: img/boundary_cell_FVM_rhs.png + :scale: 100 % + :alt: Finite volume boundary cell at the right hand side. + :align: center + + Finite volume boundary cell at the right hand side. + +.. math:: + w_J^{\prime} = -\frac{\mathcal{F}_{J+1/2}}{h_J} + \frac{\mathcal{F}_{J-1/2}}{h_J} + \bar{s}_J + +The Robin boundary condition at the right hand side is, + +.. math:: + g_{R}(x_R) = \mathcal{F}_{J+1/2} + +Therefore the naturally resolved boundary condition on the right hand side becomes, + +.. math:: + w_J^{\prime} = \frac{w_{J-1}}{h_J}\left( \frac{a(x_{J-1/2})h_{J}}{2h_{-}} + \frac{d(x_{J-1/2})}{h_{-}} \right) + \frac{w_{J}}{h_J} \left( \frac{a(x_{J-1/2}) h_{J-1}}{2h_{-}} - \frac{d(x_{J-1/2})}{h_{-}} \right) - \frac{g_{R}(x_R)}{h_J} + \bar{s}_J + +Dirichlet boundary conditions (known value) +******************************************* + +.. warning:: + This is only valid provided the initial conditions also satisfy the boundary conditions + +Dirichlet (known value) boundary conditions are trivial to implement. All that is required is to fix the boundary term to be a constant value. If the initial conditions satisfy the boundary conditions (as they should) then all that is needed to to make sure that the time-derivative (the change with time) of the boundary points is always zero. For example, to impose Dirichlet boundary conditions + +.. math:: + w_1^{n+1} = g_D(x_1) \quad w_{J}^{n+1} = g_D(x_J) + +We simply need to specify, + +.. math:: + F_1(w_1^{n}) = 0 \quad F_1(w_1^{n+1}) = 0 + +This can be achieved by multiplying the r.h.s. by the modified identity matrix where the first and last elements are zero, + +.. math:: + w^{\prime} = \text{diag}(0,1,\cdots,1,0) \left( \theta F(w^{n+1}) + (1-\theta)F(w^{n}) \right) + +If the boundary conditions are time dependent, + +.. math:: + w_1^{n+1} = g_D(x_1, t_{n+1}) \quad w_{J}^{n+1} = g_D(x_J, t_{n+1}) + +The we should also add a matrix to the r.h.s. which is contains the changes to the boundary value between the time-steps, + +.. math:: + w^{\prime} = \text{diag}(0,1,\cdots,1,0) \left( \theta F(w^{n+1}) + (1-\theta)F(w^{n}) \right) + \text{diag}(\Delta g_D(x_1), 0, \cdots, 0, \Delta g_D(x_J) ) + +where, + +.. math:: + \Delta g_D(x_1) = g_D(x_1, t_{n+1}) - g_D(x_1, t_{n}) \\ + \Delta g_D(x_J) = g_D(x_J, t_{n+1}) - g_D(x_J, t_{n}) + +The following assumes time-invariant boundary conditions. + +General matrix form with a transient term +***************************************** + +The semi-discretised advection-diffusion-reaction equation can be written in the form below using the :math:`\theta`-method, + +.. math:: + w^{\prime} = \alpha\left( \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n}) \right) + \beta + +where :math:`F(w)` contains a matrix :math:`M`, a vector of the solution variable :math:`w` and a vector of the reaction term :math:`r`, + +.. math:: + F(w) = Mw + r + +We have introduced a new matrix, :math:`\alpha` and vector, :math:`\beta`. These are used to incorporate the boundary conditions, they have the form, + +.. math:: + \alpha & = \text{diag}\left( \alpha_1, 1, \cdots, 1, \alpha_J \right) \\ + \beta & = \left( \beta_1, 0, \cdots, 0, \beta_J \right) + +The modified coefficient matrix becomes, + +.. math:: + M = + \begin{align} + \begin{pmatrix} + b_1 & c_1 & & & 0 \\ + r_a & r_b & r_c & & \\ + && \ddots & \ddots & \ddots& \\ + && & r_a & r_b & r_c \\ + 0 && & & a_J & b_J + \end{pmatrix} + \end{align} + +Table showing coefficients which should be altered to apply either Robin or Dirichlet boundary conditions. + +.. tabularcolumns:: |m{5cm}|m{5cm}|m{5cm}| + ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| Symbol | Robin | Dirichlet | ++===================+==============================================================================================================+=====================+ +| :math:`b_1` | :math:`\frac{1}{h_1}\left( \frac{-a(x_{3/2})h_{2}}{2h_{+}} - \frac{d(x_{3/2})}{h_{+}} \right)` | :math:`1` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`c_1` | :math:`\frac{1}{h_1} \left( \frac{-a(x_{3/2}) h_1}{2h_{+}} + \frac{d(x_{3/2})}{h_{+}} \right)` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`a_J` | :math:`\frac{1}{h_J}\left( \frac{a(x_{J-1/2})h_{J}}{2h_{-}} + \frac{d(x_{J-1/2})}{h_{-}} \right)` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`b_J` | :math:`\frac{1}{h_J} \left( \frac{a(x_{J-1/2}) h_{J-1}}{2h_{-}} - \frac{d(x_{J-1/2})}{h_{-}} \right)` | :math:`1` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`\alpha_1` | :math:`1` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`\alpha_J` | :math:`1` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`\beta_1` | :math:`\frac{g_R(x_1)}{h_1}` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ +| :math:`\beta_J` | :math:`-\frac{g_R(x_J)}{h_J}` | :math:`0` | ++-------------------+--------------------------------------------------------------------------------------------------------------+---------------------+ + +General matrix form without a transient term +******************************************** + +The boundary conditions scheme discussed above is only valid for initial value problems; problems where an initial vector of the solution variable :math:`w^0` is specified along with boundary conditions. PDEs of the advection-diffusion-reaction form that **do not** contain a time derivative are an important class of equations they are called *boundary value problems*, + +.. math:: + 0 = \alpha\left( \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n}) \right) + \beta + +Boundary value problems do not require initial conditions for the solution variable. For this type of equations we need to use a different general scheme to implement boundary condition values. + +The element of the :math:`M` matrix remain the same, however for Dirichlet boundary conditions the elements must be modified in for the :math:`\alpha` and :math:`\beta` terms, + ++-------------------+----------------------------------------------------------------------------+---------------------------+ +| Symbol | Robin | Dirichlet | ++===================+============================================================================+===========================+ +| :math:`\alpha_1` | :math:`1` | :math:`1` | ++-------------------+----------------------------------------------------------------------------+---------------------------+ +| :math:`\alpha_J` | :math:`1` | :math:`1` | ++-------------------+----------------------------------------------------------------------------+---------------------------+ +| :math:`\beta_1` | :math:`\frac{g_R(x_1)}{h_1}` | :math:`-r(x_1)- g_D(x_1)` | ++-------------------+----------------------------------------------------------------------------+---------------------------+ +| :math:`\beta_J` | :math:`-\frac{g_R(x_J)}{h_J}` | :math:`-r(x_J)- g_D(x_J)` | ++-------------------+----------------------------------------------------------------------------+---------------------------+ + diff --git a/_sources/examples.rst.txt b/_sources/examples.rst.txt new file mode 100644 index 0000000..83edd85 --- /dev/null +++ b/_sources/examples.rst.txt @@ -0,0 +1,38 @@ +Example calculations +-------------------- + +.. note:: + For all simulations :math:`a=1`, :math:`d=0.001` on a domain where :math:`0\leq x \leq 1`, with mesh spacing :math:`h=0.02` and time step :math:`\tau=0.01`. The equations are solver with the boundary conditions :math:`u(0,t)=1` and :math:`u(1,t)=0` and initial conditions :math:`u(x,0)=sin(\pi x)^{100}`. With these boundary conditions a boundary layer will form near :math:`x=1`. + +Uniform grid +************ + +First we test the finite-volume method using a standard uniform grid. Note that the Peclet number for the above parameters is :math:`\mu=20` so the central discretisation scheme is not stable as illustrated by the oscillations in the solution. The upwind scheme does not have a stability criteria related to the Peclet number so the solutions for the *upwind* case are smooth. Finally, the *exponentially fitted* scheme has automatically weighted in favour of the upwind discretisation, the value of :math:`\kappa\approx 0.9`. + +.. raw:: html + +
+ +
+ +Random grid +*********** + +Although a random grid is of no practical use it is a good test of the code because bugs are more likely to show up when symmetry has been reduced. I the follow simulations :math:`\text{min}(\kappa)` =0.004 and :math:`\text{max}(\kappa)` =0.98 so the discretisation scheme is abruptly changing from cell to cell. + +.. raw:: html + +
+ +
+ +Nonuniform grid +*************** + +Nonuniform grids can be used to reduce to *improve* the solution as shown here. The following simulation contains the same number of cells as the previous simulations however the cell centres are clustered towards the right hand boundary. The increased density of cells allow the boundary layer to be resolved. The transient solution computed from the *central* scheme is still significantly affected by the high Peclet number but it is interesting to observe that the steady-state solution of all three methods are very similar. Furthermore, :math:`\kappa\approx0.01-0.1` in the region of the boundary layer which implies the local value of Peclet number as been reduced enough so allow the exponentially fitted scheme to weight in favour of the higher accuracy central discretisation. + +.. raw:: html + +
+ +
diff --git a/_sources/implementation_advection_diffusion.rst.txt b/_sources/implementation_advection_diffusion.rst.txt new file mode 100644 index 0000000..c30b0fc --- /dev/null +++ b/_sources/implementation_advection_diffusion.rst.txt @@ -0,0 +1,199 @@ +Solving the advection-diffusion-reaction equation in Python +----------------------------------------------------------- + +Here we discuss how to implement a solver for the advection-diffusion equation in Python. The notes will consider how to design a solver which minimises code complexity and maximise readability. + +Overview +******** + +At a high-level usage of the code looks like the following, + +.. code-block:: python + :emphasize-lines: 3,6,7,11 + + # Define a mesh + faces = np.linspace(-0.5, 1, 100) + mesh = Mesh(faces) + + # Define coefficients + a = CellVariable(0.01, mesh=mesh) # Advection velocity + d = CellVariable(1e-3, mesh=mesh) # Diffusion coefficient + + # Make a 'model' and apply boundary conditions + k = 1 # Time step + model = Model(faces, a, d, k) + model.set_boundary_conditions(left_value=1., right_value=0.) + + # Ask for the discretised system... + M = model.coefficient_matrix() + alpha = model.alpha_matrix() + beta = model.beta_vector() + + # Solve... + +There are a number of classes (highlighted above) which abstract away details of dealing with finite-volume equations: + + 1. :code:`Mesh` + 2. :code:`CellVariable` + 3. :code:`Model` + +The the :code:`Mesh` and :code:`CellVariable` classes have been inspired by the approach of Fipy_. + +.. _Fipy: http://www.ctcms.nist.gov/fipy/ + +In the following we will highlight the main features of each class and point out how they are useful. The classes do not attempt at doing "*too much*", they simply aid in the creation of the linear system. + +The :code:`Mesh` class +********************** + +:code:`Mesh` objects are initialised with a list of faces locations, which can be non-uniformly distributed if desired. A mesh is completely defined by locations of cell **faces**. Some useful methods, + +.. code:: python + + def h(self, i): + ... + +:code:`h()` returns the cell width for cell with index :code:`i`. This is function is vectorisable by passing an array of the required indices but it *does not* accept "fancy indexing". The reason being, that it is very easy to make mistakes with subscript indexing, the goal here is to make the user be explicit when requesting elements. Note the :code:`self.cell_widths` instance variable returns the numpy array of cell widths is a second way of accessing this data. + +.. code:: python + + def hm(self, i): + ... + +:code:`hm()` returns the distance between cell centroids for the cell at index :code:`i` and :code:`i-1`, that is in the *backwards* or **minus** direction. + +.. code:: python + + def hp(self, i): + ... + +Similarly :code:`hp()` returns the distance between cell centroids for the cell at index :code:`i` and :code:`i+1`, that is in the *forwards* or **plus** direction. + +In addition to the above method, the class contains the instance variables, :code:`self.cells` which returns an array of cell centroid locations, :code:`self.J` which contains the number of cells, and also a copy of the face locations (an array) via :code:`self.faces`. + +The :code:`CellVariable` class +****************************** + +The goal of the :code:`CellVariable` class is to provide a elegant way of automatically interpolating between the cell value and the face value. The class holes values which correspond to the **cell average**. Internally, this class is a subclass of :code:`numpy.ndarray` so it is a fully functioning numpy array. It has a new constructor and additional method which return interpolated values at the cell surfaces. + +A :code:`CellVariable` is initialised with a value for the cell average (this can be a constant or an array-like quantity) and the :code:`Mesh` on which the cell variable is defined. My coupling the cell variable with the mesh the class can perform interpolation between the cell and face values using the methods, + +.. code:: python + + def p(self, i): + ... + + def m(self, i) + ... + +Again :code:`self.p(i)` stands for the *plus* direction and :code:`self.m(i)` stands for the *minus* direction, as such they return values at the right and left face of the cell. The mesh can be returned via the instance variable :code:`cell_variable.mesh`. + + +The :code:`Model` class +*********************** + +The model class is where the creating of the matrices occurs and where boundary conditions can be applied to the problem. For these reasons the class is fairly complicated. + +There are method which return different element of the final matrix. The interior elements are fairly homogenous, the only real difference is where there are spatially varying coefficient of cell widths. For this reason the the method :code:`_interior_matrix_elements()` returns elements which correspond to the lower, central and upper diagonals for a specific index. For example, to calculate the interior matrix elements for mesh point at value :code:`index` one would do the following, + +.. code:: python + + # Return the the interior matrix elements (the r-terms) + # for a particular spatial index + model = Model(...) + index = ... + # The lower, central and, upper diagonal terms of the stencil + ra, rb, rc = model._interior_matrix_elements(index) + +The function names here correspond to the matrix element in the previous section. + +Note that the function is prefixed with an underscore this is because are private, 'users' should have no need to call these method. It is called internally when constructing the finite-volume matrices. However, as this is a overview of how to implement this is an readable and useful way we include this detail. + +The following methods play a similar role, + +.. code:: python + + def _robin_boundary_condition_matrix_elements_left(self): + ... + + def _robin_boundary_condition_matrix_elements_right(self): + ... + + def _dirichlet_boundary_condition_matrix_elements_left(self): + ... + + def _dirichlet_boundary_condition_matrix_elements_right(self): + ... + +They return a list of index-value pairs :code:`([(1,1), a11], [(i,2), b12] ...)`. The functions return the value of element which need to change (with respect to the interior values) in order include boundary conditions. The index-value pair facilitates automatic insertion of the values into the correct matrix element. As we will see later, rather than hard coding the position of the various element if the index and value are specified it makes the destination of the element unambiguous. It also allows the value of the matrix element to be defined at the same point in the code as the location. This is beneficial for providing context and should reduce bugs and complexity. + +The elements for the :math:`\beta` boundary condition vector (which is added to the linear system) are generated from the functions below, + +.. code:: python + + def _robin_boundary_condition_vector_elements_left(self): + ... + + def _robin_boundary_condition_vector_elements_right(self): + ... + + def _dirichlet_boundary_condition_vector_elements_left(self): + ... + + def _dirichlet_boundary_condition_vector_elements_right(self): + ... + +Again, these method should return *index-values* pairs. + +The :code:`Model` class also include some convenience function for checking the value of the Peclet number and the CFL conditions which can be called via, + +.. code:: python + + def peclet_number(self): + return self.a * self.mesh.cell_widths / self.d + + def CFL_condition(self): + return self.a * self.k / self.mesh.cell_widths + + +The method which are intended for the user to actually call when constructing the linear system are, + + +.. code:: python + + def coefficient_matrix(self): + ... + + def alpha_matrix(self): + ... + + def beta_vector(self): + ... + +The linear system for time-stepping can be constructed easily, + +.. code:: python + + # In pseudo-code + model = AdvectionDiffusionModel(...) + M = model.coefficient_matrix() + alpha = model.alpha_matrix() + beta = model.beta_vector() + I = sparse.identity(model.mesh.J) + + # Apply boundary conditions + u_init = ... # + ... + + tau = 0.01 # time step + theta = 0.5 # Implicit/explicit parameter + + u = u_init + for i in range(...): + # time step the linear system, A.x = d + A = I - tau * theta * alpha * M + d = (I + tau * (1-theta) * alpha * M) * u + u = spsolve(A,d) + + +Finally, when initialising a :code:`Model` object the keyword argument :code:`, discretisation` is important. Is can be set to one of the following :code:`'upwind', 'central', 'exponential'`. The :code:`upwind` option uses the classic *first order upwind* discretisation, :code:`central` uses *second-order central* and setting to :code:`exponential` uses an adaptive scheme which will use weight between the central and upwind scheme depending on the local value of the Peclet number. This is the classic 'exponential fitting' or 'Scharfetter-Gummel' discretisation. **N.B.** Scharfetter-Gummel also refers to a method of solving the advection-diffusion equation is a non-coupled manner, this is not the case here where it only refers to the the discretisation method. diff --git a/_sources/index.rst.txt b/_sources/index.rst.txt new file mode 100644 index 0000000..335f28b --- /dev/null +++ b/_sources/index.rst.txt @@ -0,0 +1,31 @@ +.. FVM Docs documentation master file, created by + sphinx-quickstart on Fri Jun 21 16:17:12 2013. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Notes on implementing the finite-volume method for physical simulations +======================================================================= + +Contents: + +.. toctree:: + :maxdepth: 2 + + overview.rst + advection_diffusion.rst + advection_diffusion_boundary_conditions.rst + solving_linear_equations.rst + solving_nonlinear_equations.rst + poisson.rst + implementation_advection_diffusion.rst + examples.rst + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + diff --git a/_sources/overview.rst.txt b/_sources/overview.rst.txt new file mode 100644 index 0000000..640bec8 --- /dev/null +++ b/_sources/overview.rst.txt @@ -0,0 +1,29 @@ +Overview +-------- + +These are notes that I am making while learning the finite-volume method (FVM). The methods discussed are based on the the book by W. Hundsdorfer and J. G. Verwer, Numerical solutions of time-dependent advection-diffusion reaction equations. + +See the the accompanying project on Github for implementation of the solvers in Python branch, http://danieljfarrell.github.come/FVM/. + + +Finite volumes vs. finite differences +************************************* + +When first starting to solve (partial) differential equation problems a student will first encounter the finite difference method (FDM). The FDM relies on approximating the a differential operator as a series of distinct points, + +.. figure:: img/grid_FDM.png + :scale: 100 % + :alt: Finite difference grid + :align: center + + +In contrast, the FVM uses a fundamental unit of averages over small spatial elements, or cells, + +.. figure:: img/cells_FVM.png + :scale: 100 % + :alt: Finite volume grid. + :align: center + +Although this distinction may seem trivial at first or even over complicated it is a particularly power technique when applied to conservation equations because the conservative properties of the governing equation are maintained even at the discretised level without introducing any approximations. This property is achieved by symbolic integrating the equations over a finite volume, then using solve the problem in terms of the integrated variables. + +This sounds much more complicated that it is, let's look at an example. diff --git a/_sources/poisson.rst.txt b/_sources/poisson.rst.txt new file mode 100644 index 0000000..cd6ce72 --- /dev/null +++ b/_sources/poisson.rst.txt @@ -0,0 +1,101 @@ +The Poisson equation +-------------------- + +The Poisson equation forms the basis of electrostatics and is of the form, + +.. math:: + 0 = (\epsilon(x) \phi_x)_x + \rho(x,t) + +where :math:`\phi` is the electric potential, :math:`\epsilon(x)` is the materials dielectric constant, :math:`\rho(x, t)` is a charge distribution (possibly varying with time), and the :math:`x` subscripts indicate a spatial partial derivative. Note that this equation is reaction-diffusion equation with the diffusion flux :math:`\mathcal{F} = -\epsilon(x) \phi_x`. As with the advection-diffusion equation, this expression holds over the whole domain, therefore it must also hold over a subdomain :math:`\Omega_j=(x_{j-1/2}, x_{j+1/2})`. Integrating over the subdomain, + +.. math:: + 0 = \int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx + \int_{x_{j-1/2}}^{x_{j+1/2}} \rho(x,t)~dx + +The cell averaged values are therefore, + +.. math:: + 0 = -\frac{\mathcal{F}_{j+1/2}}{h_j} + \frac{\mathcal{F}_{j-1/2}}{h_{j}} + \bar{\rho}_j + +To determine the flux at the cell faces we will take the mean of the fluxes in cells adjacent to the interfaces, + +.. math:: + \mathcal{F}_{j+\frac{1}{2}} = - \epsilon_{j+\frac{1}{2}} \frac{\phi_{j+1}-\phi_j}{h_{+}} + +.. math:: + \mathcal{F}_{j-\frac{1}{2}} = - \epsilon_{j-\frac{1}{2}} \frac{\phi_{j}-\phi_{j-1}}{h_{-}} + +Substituting the definition of the fluxes into the integrated equation and factoring for the electrical potential terms yields, + +.. math:: + 0 = \frac{\phi_{j-1}}{h_j} \left( \frac{\epsilon(x_{j-\frac{1}{2}})}{h_{-}} \right) - \frac{\phi_{j}}{h_j} \left( \frac{\epsilon(x_{j-\frac{1}{2}})}{h_{-}} + \frac{\epsilon(x_{j+\frac{1}{2}})}{h_{+}} \right) + \frac{\phi_{j+1}}{h_j} \left( \frac{\epsilon(x_{j+\frac{1}{2}})}{h_{+}} \right) + \bar{\rho}_j + +Naturally resolved Neumann boundary conditions +********************************************** + +Neumann boundary conditions applied to the finite volume form of the Poisson equation are naturally resolved; that is to say that ghost cell or extrapolation approaches are not needed. Consider the integral form of the Poisson equation at the last cell of the domain on the right hand side, + +.. math:: + 0 = \frac{1}{h_J}\left( \epsilon(x_{R}) \phi_x\bigg|_{x=x_R} - \epsilon(x_{J-1/2}) \phi_x\bigg|_{J-1/2} \right) + \bar{\rho}_J + +The Neumann boundary condition defines the derivative of :math:`\phi` at the the right hand boundary, + +.. math:: + \phi_x\bigg|_{x=x_R} = \sigma_R + +This boundary condition appears naturally as a term in the integral equation. Therefore the boundary condition can be satisfied by substitution into the integral form of the equation at the right hand boundary, + +.. math:: + \frac{1}{h_J}\left( \epsilon(x_{R}) \sigma_R - \epsilon(x_{J-1/2}) \phi_x\bigg|_{J-1/2} \right) + \bar{\rho}_J + +Dirichlet boundary conditions +***************************** + +Fixed value, or Dirichlet, boundary conditions are not naturally resolved and can only be applied to the discretised form of the equation. We need to apply a Dirichlet condition to the left hand side, such that full equation for the cell :math:`j=1` reads, + +.. math:: + \phi(x_L) = \omega_L + +Clearly, this is trivial as all that is required is fix the first row of the matrix equation to be the boundary condition. + +Poisson equation in matrix form +******************************* + +Defining the (vector) coefficients, + +.. math:: + \ell_a(j) & = \frac{1}{h_j}\frac{\epsilon(x_{j-1/2})}{h_{-}} \\ + \ell_b(j) & = -\frac{1}{h_j}\left( \frac{\epsilon_{j-1/2}}{h_{-}} + \frac{\epsilon_{j+1/2}}{h_{+}} \right) \\ + \ell_c(j) & = \frac{1}{h_j}\frac{\epsilon(x_{j+1/2})}{h+} \\ + +With Dirichlet boundary condition on the left hand side and a Neumann boundary condition on the right hand side the Poisson equation can be written in matrix form, + +.. math:: + \begin{align} + \begin{pmatrix} + 1 & 0 & & & 0 \\ + \ell_a(2) & \ell_b(2) & \ell_c(2) & & \\ + & \ddots & \ddots & \ddots & \\ + & & \ell_a(J-1) & \ell_b(J-1) & \ell_c(J-1) \\ + 0 & & & \ell_a(J) & -\ell_a(J) + \end{pmatrix} + \begin{pmatrix} + \phi_1 \\ + \phi_2 \\ + \vdots \\ + \phi_{J-1} \\ + \phi_{J} \\ + \end{pmatrix} + + + \begin{pmatrix} + -\omega_L \\ + \bar{\rho}_2 \\ + \vdots \\ + \bar{\rho}_{J-1} \\ + 1/h_j\epsilon(x_J)\sigma_R + \bar{\rho}_{J} \\ + \end{pmatrix} + = 0 + \end{align} + + + + diff --git a/_sources/solving_linear_equations.rst.txt b/_sources/solving_linear_equations.rst.txt new file mode 100644 index 0000000..a5d7160 --- /dev/null +++ b/_sources/solving_linear_equations.rst.txt @@ -0,0 +1,32 @@ +Solving linear equations +------------------------ + +Provided the equation is linear, meaning that the coefficients, nor the reaction term depend on the solution variable a time-stepping approach can be used to solve the equation. + +Time stepping the advection-diffusion equation +********************************************** + +Here we will rearrange the :math:`\theta`-method form of the advection-diffusion equation, + +.. math:: + \frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n}) + +into the form of a linear system :math:`A\cdot x = d`. Firstly lets move all terms involving future time points the l.h.s, + +.. math:: + w^{n+1} - \theta \tau F(w^{n+1}) = w^{n} + (1-\theta) F(w^{n+1}) + +Replacing the :math:`F` terms with the full matrix expressions and factoring yields, + +.. math:: + \underbrace{(I - \theta\tau M^{n+1})}_{A}\underbrace{w^{n+1}}_{x} = \underbrace{(I + (1-\theta)\tau M^{n})w^{n}}_{d} + +Time-stepping involves solving this equation iteratively. First the initial conditions :math:`w^0` is used to calculate the solution variable at the next point in time :math:`w^1`, then values of the solution variable are updated such that :math:`w^1` is used to calculate :math:`w^2`, and so on and so forth. This procedure is shown in the Figure below, + +.. figure:: img/time_stepping_linear_advection_diffusion.png + :scale: 50 % + :alt: Time stepping the linear advection-diffusion equation. + :align: center + +.. note:: + Parameters used for this simulation. :math:`a=1`, :math:`d=1\times 10^{-3}`, :math:`h=5 \times 10^{-3}`, :math:`\tau=5 \times 10^{-4}` with exponentially fitted discretisation. The initial and boundary conditions where, :math:`u(x,t_{0}) =sin(\pi x)^{100}` with :math:`u(x_1)=1` and :math:`u_x(x_J)=0`. The equation was time stepped for 400 iterations, the plots so the solution at times :math:`t=0,0.05,0.1,0.15,0.2`. \ No newline at end of file diff --git a/_sources/solving_nonlinear_equations.rst.txt b/_sources/solving_nonlinear_equations.rst.txt new file mode 100644 index 0000000..9cac420 --- /dev/null +++ b/_sources/solving_nonlinear_equations.rst.txt @@ -0,0 +1,108 @@ +Solving nonlinear equations +--------------------------- + +.. note:: + Nonlinear equations are defined as those having coefficients which are functions of the solution variable. + +In the previous section we discussed how to solve the linear advection-diffusion-reaction equation with method of time stepping. We used the Crank-Nicolson scheme: the implicit terms provides stability (unrestricted time step) and by also using the explicit term the accuracy of the time integration is improved to second order. However, this approach is only valid for linear equations. + +Nonlinear equations (by definition) cannot be written in linear form, as such the time-stepping approach cannot be used. It remains possible to use a fully explicit method, however in practice this is usually not done because the time step is too restrictive or the explicit form is *unconditionally unstable*. + +Here we discuss three approach which can be used to solve nonlinear equations: + + 1. IMEX (implicit/explicit) time-stepping + 2. Newton iteration + 3. Method of lines + +Fisher's equation +***************** + +We will use Fishers' nonlinear reaction-diffusion equation as an example because it has an analytical solution, + +.. math:: + u_t = d u_{xx} + \beta u(1 - u) + +We can directly apply the discretisation scheme already developed, but we must set :math:`a=0` and note the nonlinear reaction term :math:`\beta u(1 - u)`. With semi-open boundary condition :math:`u(x_0,t)=1` Fishers' equation has the analytical form, + +.. math:: + u(x,t) = \frac{1}{\left(1 + e^{c^{1/2}\left(x - Ut\right)} \right)^2} \\ + +where :math:`U = 5\left(\frac{1}{c}\right)^{1/2}`, and, :math:`c=\frac{6}{d\beta}` + + +IMEX time-stepping +****************** + +The IMEX scheme treats the linear parts of the equation using implicit (or semi-implicit) scheme and the non-linear parts full explicitly. Therefore for Fishers' equation this becomes, + +.. math:: + w^{\prime} = \theta M^{n+1}u^{n+1} + (1-\theta) M^{n}u^{n} + \beta u^n(1-u^n) + +This can be rearranged into the form :math:`A\cdot x = d` and solved by standard time stepping techniques. First we discretised the time derivative :math:`u^{\prime} = \frac{u^{n+1} - u^{n}}{\tau}`, then we rearrange the equation so that future time points are on the l.h.s. and known time points are on the right, + +.. math:: + u^{n+1} - \theta \tau M^{n+1}u^{n+1} = u^{n} + (1-\theta) \tau M^{n}u^{n} + \tau\beta u^n(1-u^n) + +Using the identity matrix, :math:`I` gives the desired result, + +.. math:: + (I - \theta \tau M^{n+1})u^{n+1} = (I + (1-\theta) \tau M^{n} ) u^{n} + \tau \beta u^n(1-u^n) + +This can be readily solved using any linear solver, the results of this approach are shown in the Figure below. Note the poor agreement with the analytical solution for large values of times. IMEX method is known to introduce a global error that grows exponentially with iteration number. The error in the above Figure has been accentuated by deliberately choosing a fairly large time step. A small time step can be used suppress errors. + +.. figure:: img/IMEX_solution_fvm.png + :scale: 50 % + :alt: IMEX solution of Fishers' equation. + :align: center + + Solution of Fishers' equation with an IMEX scheme. + +Newton iteration +**************** + +We have a semi-discretised system of the form :math:`u^{\prime}(t) = F(u(t))` with the application of the :math:`\theta`-method this gives, + +.. math:: + w^{n+1} - w^{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) = 0 + +with unknown vector :math:`w^{n+1}`. This equation can be solved using Newton iteration. + +.. math:: + \nu^{k+1} = \nu^{k} - (I - \theta\tau A^{n})^{-1} \left( \nu^{k} - u_{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) \right) + +where :math:`k` is the iteration index (:math:`k\geq0`) and :math:`A^{n}` is the Jacobian matrix of :math:`F(w^n)`. We use the symbol :math:`\nu^{k}` for iteration variables such that they are distinguished from solution of the equation at a real time point :math:`u^n`. The iteration needs a starting value, it is perfectly value to choose :math:`\nu^0 = w^n` or for a better estimate we can precompute one iteration :math:`\nu^0 = w^n + \tau F(w^{n})`. Here we have described the so-called *modified* Newton iteration because the Jacobian is not updated during the iteration, this is known to work well for stiff PDEs. + +Applying the modified Newton iteration to Fisher's equation yields the results shown in the Figure below. The results are much better than for IMEX scheme as the error (particularly for long time points) is much lower. However, the accumulation of global error can still be observed. This is a known feature of so called *single step methods*. In ODE terminology a single step method is one that uses only the most recent past data to evaluate the equation at the future time step. Note even if the equations are implicit (or semi-implicit as with the :math:`\theta`-method) the future state is predicted only from one historic data point. In contrast a *multistep* method would used many previous data points to achieve highly accurate time integration. + + +.. figure:: img/Newton_solution_fvm.png + :scale: 50 % + :alt: IMEX solution of Fishers' equation. + :align: center + + Solution of Fishers' equation with a modified Newton iteration. + +Method of lines +*************** + +.. note:: + For an excellent introduction to the Method of Lines see `the article at scholarpedia.org `_ + +The method of lines is not really a method, it is a way of writing PDEs in such a way that they can be solved by black box *ODE* solvers. The time integration of ODEs is a mature topic. Robust and well tested solvers are readily available on the internet and integrated into many open-source projects for easy use. The method of lines is a way of letting PDE equation take advantage of the mature nature of ODE solvers. + +We have a semi-discretised system of the form, :math:`u^{\prime}(t) = F(u(t))` to use the method of lines we leave the PDE is the semi-discristed form such that is has a continuous time derivative on the l.h.s., i.e. :math:`u^{\prime}(t) = \frac{du(t)}{dt}`. In the semi-discretised form what we really have is no longer a PDE but a system of coupled ODEs. Moreover, a PDE equation is a differential equation which contains differential terms with respect to more that one variable. In this semi-discretised form we have replaced the spatial differential operator with a matrix so there is only one derivative. The equation has become a system of ODEs with number of equation proportional to the number of grid points in our discretisation. + +.. note:: + The discretised boundary conditions are expected to be included in :math:`F(u(t))` such that the problem is well posed. + +There are a wide variety of ODE solvers we are interested in the *implicit* type as we know many PDE problems are unstable in explicit form. The result in the Figure below have been calculated using an implicit multistep Adams-Moulton method (the algorithm used is the highly popular *vode* solver from `netlib.org `_, although we used actually used a Python wrapper from `scipy `_ to perform the calculations shown). This can be thought of as a generalised of the :math:`\theta`-method to include a multiple number of future time steps which are solved simultaneously with a Newtom iteration. Unlike single step method, multistep methods do not show a global error accumulation, so the time integration is a *much* better approximation to analytical solution. Fisher's equation is not particularly stiff so this approach worked well. However for stiff equations a technique called *Backward differentiation formulas* (BDF) is used. The vode solver has the ability to dynamically adapt between stiff and non stiff mode (i.e. Adams-Moulton to BDF) which is one of the reasons for it's popularity. + +The method of lines is a powerful technique for solving semi-discretised PDE problems as much of the details of numerically solving the equation can be off loaded to mature and robust solvers. I do not yet have enough experience with the technique to point out short coming or failings. But it seems that provided the spatial discretisation is stable the details of the time integration can be almost ignored! + +.. figure:: img/MOL_solution_fvm.png + :scale: 50 % + :alt: MOL solution of Fishers' equation. + :align: center + + Solution of Fishers' equation by the method of lines. + diff --git a/_static/basic.css b/_static/basic.css index 43e8baf..f316efc 100644 --- a/_static/basic.css +++ b/_static/basic.css @@ -4,7 +4,7 @@ * * Sphinx stylesheet -- basic theme. * - * :copyright: Copyright 2007-2011 by the Sphinx team, see AUTHORS. + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. * :license: BSD, see LICENSE for details. * */ @@ -15,6 +15,12 @@ div.clearer { clear: both; } +div.section::after { + display: block; + content: ''; + clear: left; +} + /* -- relbar ---------------------------------------------------------------- */ div.related { @@ -52,6 +58,8 @@ div.sphinxsidebar { width: 230px; margin-left: -100%; font-size: 90%; + word-wrap: break-word; + overflow-wrap : break-word; } div.sphinxsidebar ul { @@ -79,16 +87,29 @@ div.sphinxsidebar input { font-size: 1em; } +div.sphinxsidebar #searchbox form.search { + overflow: hidden; +} + div.sphinxsidebar #searchbox input[type="text"] { - width: 170px; + float: left; + width: 80%; + padding: 0.25em; + box-sizing: border-box; } div.sphinxsidebar #searchbox input[type="submit"] { - width: 30px; + float: left; + width: 20%; + border-left: none; + padding: 0.25em; + box-sizing: border-box; } + img { border: 0; + max-width: 100%; } /* -- search page ----------------------------------------------------------- */ @@ -109,7 +130,7 @@ ul.search li a { font-weight: bold; } -ul.search li div.context { +ul.search li p.context { color: #888; margin: 2px 0 0 30px; text-align: left; @@ -123,6 +144,8 @@ ul.keywordmatches li.goodmatch a { table.contentstable { width: 90%; + margin-left: auto; + margin-right: auto; } table.contentstable p.biglink { @@ -150,9 +173,14 @@ table.indextable td { vertical-align: top; } -table.indextable dl, table.indextable dd { +table.indextable ul { margin-top: 0; margin-bottom: 0; + list-style-type: none; +} + +table.indextable > tbody > tr > td > ul { + padding-left: 0em; } table.indextable tr.pcap { @@ -184,19 +212,45 @@ div.genindex-jumpbox { padding: 0.4em; } +/* -- domain module index --------------------------------------------------- */ + +table.modindextable td { + padding: 2px; + border-collapse: collapse; +} + /* -- general body styles --------------------------------------------------- */ +div.body { + min-width: 360px; + max-width: 800px; +} + +div.body p, div.body dd, div.body li, div.body blockquote { + -moz-hyphens: auto; + -ms-hyphens: auto; + -webkit-hyphens: auto; + hyphens: auto; +} + a.headerlink { visibility: hidden; } +a:visited { + color: #551A8B; +} + h1:hover > a.headerlink, h2:hover > a.headerlink, h3:hover > a.headerlink, h4:hover > a.headerlink, h5:hover > a.headerlink, h6:hover > a.headerlink, -dt:hover > a.headerlink { +dt:hover > a.headerlink, +caption:hover > a.headerlink, +p.caption:hover > a.headerlink, +div.code-block-caption:hover > a.headerlink { visibility: visible; } @@ -208,10 +262,6 @@ div.body td { text-align: left; } -.field-list ul { - padding-left: 1em; -} - .first { margin-top: 0 !important; } @@ -221,19 +271,25 @@ p.rubric { font-weight: bold; } -img.align-left, .figure.align-left, object.align-left { +img.align-left, figure.align-left, .figure.align-left, object.align-left { clear: left; float: left; margin-right: 1em; } -img.align-right, .figure.align-right, object.align-right { +img.align-right, figure.align-right, .figure.align-right, object.align-right { clear: right; float: right; margin-left: 1em; } -img.align-center, .figure.align-center, object.align-center { +img.align-center, figure.align-center, .figure.align-center, object.align-center { + display: block; + margin-left: auto; + margin-right: auto; +} + +img.align-default, figure.align-default, .figure.align-default { display: block; margin-left: auto; margin-right: auto; @@ -247,30 +303,45 @@ img.align-center, .figure.align-center, object.align-center { text-align: center; } +.align-default { + text-align: center; +} + .align-right { text-align: right; } /* -- sidebars -------------------------------------------------------------- */ -div.sidebar { +div.sidebar, +aside.sidebar { margin: 0 0 0.5em 1em; border: 1px solid #ddb; - padding: 7px 7px 0 7px; + padding: 7px; background-color: #ffe; width: 40%; float: right; + clear: right; + overflow-x: auto; } p.sidebar-title { font-weight: bold; } +nav.contents, +aside.topic, +div.admonition, div.topic, blockquote { + clear: left; +} + /* -- topics ---------------------------------------------------------------- */ +nav.contents, +aside.topic, div.topic { border: 1px solid #ccc; - padding: 7px 7px 0 7px; + padding: 7px; margin: 10px 0 10px 0; } @@ -292,10 +363,6 @@ div.admonition dt { font-weight: bold; } -div.admonition dl { - margin-bottom: 0; -} - p.admonition-title { margin: 0px 10px 5px 0px; font-weight: bold; @@ -306,13 +373,55 @@ div.body p.centered { margin-top: 25px; } +/* -- content of sidebars/topics/admonitions -------------------------------- */ + +div.sidebar > :last-child, +aside.sidebar > :last-child, +nav.contents > :last-child, +aside.topic > :last-child, +div.topic > :last-child, +div.admonition > :last-child { + margin-bottom: 0; +} + +div.sidebar::after, +aside.sidebar::after, +nav.contents::after, +aside.topic::after, +div.topic::after, +div.admonition::after, +blockquote::after { + display: block; + content: ''; + clear: both; +} + /* -- tables ---------------------------------------------------------------- */ table.docutils { + margin-top: 10px; + margin-bottom: 10px; border: 0; border-collapse: collapse; } +table.align-center { + margin-left: auto; + margin-right: auto; +} + +table.align-default { + margin-left: auto; + margin-right: auto; +} + +table caption span.caption-number { + font-style: italic; +} + +table caption span.caption-text { +} + table.docutils td, table.docutils th { padding: 1px 8px 1px 5px; border-top: 0; @@ -321,14 +430,6 @@ table.docutils td, table.docutils th { border-bottom: 1px solid #aaa; } -table.field-list td, table.field-list th { - border: 0 !important; -} - -table.footnote td, table.footnote th { - border: 0 !important; -} - th { text-align: left; padding-right: 5px; @@ -343,6 +444,126 @@ table.citation td { border-bottom: none; } +th > :first-child, +td > :first-child { + margin-top: 0px; +} + +th > :last-child, +td > :last-child { + margin-bottom: 0px; +} + +/* -- figures --------------------------------------------------------------- */ + +div.figure, figure { + margin: 0.5em; + padding: 0.5em; +} + +div.figure p.caption, figcaption { + padding: 0.3em; +} + +div.figure p.caption span.caption-number, +figcaption span.caption-number { + font-style: italic; +} + +div.figure p.caption span.caption-text, +figcaption span.caption-text { +} + +/* -- field list styles ----------------------------------------------------- */ + +table.field-list td, table.field-list th { + border: 0 !important; +} + +.field-list ul { + margin: 0; + padding-left: 1em; +} + +.field-list p { + margin: 0; +} + +.field-name { + -moz-hyphens: manual; + -ms-hyphens: manual; + -webkit-hyphens: manual; + hyphens: manual; +} + +/* -- hlist styles ---------------------------------------------------------- */ + +table.hlist { + margin: 1em 0; +} + +table.hlist td { + vertical-align: top; +} + +/* -- object description styles --------------------------------------------- */ + +.sig { + font-family: 'Consolas', 'Menlo', 'DejaVu Sans Mono', 'Bitstream Vera Sans Mono', monospace; +} + +.sig-name, code.descname { + background-color: transparent; + font-weight: bold; +} + +.sig-name { + font-size: 1.1em; +} + +code.descname { + font-size: 1.2em; +} + +.sig-prename, code.descclassname { + background-color: transparent; +} + +.optional { + font-size: 1.3em; +} + +.sig-paren { + font-size: larger; +} + +.sig-param.n { + font-style: italic; +} + +/* C++ specific styling */ + +.sig-inline.c-texpr, +.sig-inline.cpp-texpr { + font-family: unset; +} + +.sig.c .k, .sig.c .kt, +.sig.cpp .k, .sig.cpp .kt { + color: #0033B3; +} + +.sig.c .m, +.sig.cpp .m { + color: #1750EB; +} + +.sig.c .s, .sig.c .sc, +.sig.cpp .s, .sig.cpp .sc { + color: #067D17; +} + + /* -- other body styles ----------------------------------------------------- */ ol.arabic { @@ -365,11 +586,81 @@ ol.upperroman { list-style: upper-roman; } +:not(li) > ol > li:first-child > :first-child, +:not(li) > ul > li:first-child > :first-child { + margin-top: 0px; +} + +:not(li) > ol > li:last-child > :last-child, +:not(li) > ul > li:last-child > :last-child { + margin-bottom: 0px; +} + +ol.simple ol p, +ol.simple ul p, +ul.simple ol p, +ul.simple ul p { + margin-top: 0; +} + +ol.simple > li:not(:first-child) > p, +ul.simple > li:not(:first-child) > p { + margin-top: 0; +} + +ol.simple p, +ul.simple p { + margin-bottom: 0; +} + +aside.footnote > span, +div.citation > span { + float: left; +} +aside.footnote > span:last-of-type, +div.citation > span:last-of-type { + padding-right: 0.5em; +} +aside.footnote > p { + margin-left: 2em; +} +div.citation > p { + margin-left: 4em; +} +aside.footnote > p:last-of-type, +div.citation > p:last-of-type { + margin-bottom: 0em; +} +aside.footnote > p:last-of-type:after, +div.citation > p:last-of-type:after { + content: ""; + clear: both; +} + +dl.field-list { + display: grid; + grid-template-columns: fit-content(30%) auto; +} + +dl.field-list > dt { + font-weight: bold; + word-break: break-word; + padding-left: 0.5em; + padding-right: 5px; +} + +dl.field-list > dd { + padding-left: 0.5em; + margin-top: 0em; + margin-left: 0em; + margin-bottom: 0em; +} + dl { margin-bottom: 15px; } -dd p { +dd > :first-child { margin-top: 0px; } @@ -383,30 +674,32 @@ dd { margin-left: 30px; } -dt:target, .highlighted { - background-color: #fbe54e; +.sig dd { + margin-top: 0px; + margin-bottom: 0px; } -dl.glossary dt { - font-weight: bold; - font-size: 1.1em; +.sig dl { + margin-top: 0px; + margin-bottom: 0px; } -.field-list ul { - margin: 0; - padding-left: 1em; +dl > dd:last-child, +dl > dd:last-child > :last-child { + margin-bottom: 0; } -.field-list p { - margin: 0; +dt:target, span.highlighted { + background-color: #fbe54e; } -.refcount { - color: #060; +rect.highlighted { + fill: #fbe54e; } -.optional { - font-size: 1.3em; +dl.glossary dt { + font-weight: bold; + font-size: 1.1em; } .versionmodified { @@ -447,11 +740,26 @@ dl.glossary dt { font-style: oblique; } +.classifier:before { + font-style: normal; + margin: 0 0.5em; + content: ":"; + display: inline-block; +} + abbr, acronym { border-bottom: dotted 1px; cursor: help; } +.translated { + background-color: rgba(207, 255, 207, 0.2) +} + +.untranslated { + background-color: rgba(255, 207, 207, 0.2) +} + /* -- code displays --------------------------------------------------------- */ pre { @@ -459,37 +767,105 @@ pre { overflow-y: hidden; /* fixes display issues on Chrome browsers */ } +pre, div[class*="highlight-"] { + clear: both; +} + +span.pre { + -moz-hyphens: none; + -ms-hyphens: none; + -webkit-hyphens: none; + hyphens: none; + white-space: nowrap; +} + +div[class*="highlight-"] { + margin: 1em 0; +} + td.linenos pre { - padding: 5px 0px; border: 0; background-color: transparent; color: #aaa; } table.highlighttable { - margin-left: 0.5em; + display: block; +} + +table.highlighttable tbody { + display: block; +} + +table.highlighttable tr { + display: flex; } table.highlighttable td { - padding: 0 0.5em 0 0.5em; + margin: 0; + padding: 0; } -tt.descname { - background-color: transparent; - font-weight: bold; - font-size: 1.2em; +table.highlighttable td.linenos { + padding-right: 0.5em; +} + +table.highlighttable td.code { + flex: 1; + overflow: hidden; +} + +.highlight .hll { + display: block; +} + +div.highlight pre, +table.highlighttable pre { + margin: 0; +} + +div.code-block-caption + div { + margin-top: 0; +} + +div.code-block-caption { + margin-top: 1em; + padding: 2px 5px; + font-size: small; } -tt.descclassname { +div.code-block-caption code { background-color: transparent; } -tt.xref, a tt { +table.highlighttable td.linenos, +span.linenos, +div.highlight span.gp { /* gp: Generic.Prompt */ + user-select: none; + -webkit-user-select: text; /* Safari fallback only */ + -webkit-user-select: none; /* Chrome/Safari */ + -moz-user-select: none; /* Firefox */ + -ms-user-select: none; /* IE10+ */ +} + +div.code-block-caption span.caption-number { + padding: 0.1em 0.3em; + font-style: italic; +} + +div.code-block-caption span.caption-text { +} + +div.literal-block-wrapper { + margin: 1em 0; +} + +code.xref, a code { background-color: transparent; font-weight: bold; } -h1 tt, h2 tt, h3 tt, h4 tt, h5 tt, h6 tt { +h1 code, h2 code, h3 code, h4 code, h5 code, h6 code { background-color: transparent; } @@ -521,6 +897,15 @@ span.eqno { float: right; } +span.eqno a.headerlink { + position: absolute; + z-index: 1; +} + +div.math:hover a.headerlink { + visibility: visible; +} + /* -- printout stylesheet --------------------------------------------------- */ @media print { diff --git a/_static/classic.css b/_static/classic.css new file mode 100644 index 0000000..5530147 --- /dev/null +++ b/_static/classic.css @@ -0,0 +1,269 @@ +/* + * classic.css_t + * ~~~~~~~~~~~~~ + * + * Sphinx stylesheet -- classic theme. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * + */ + +@import url("basic.css"); + +/* -- page layout ----------------------------------------------------------- */ + +html { + /* CSS hack for macOS's scrollbar (see #1125) */ + background-color: #FFFFFF; +} + +body { + font-family: sans-serif; + font-size: 100%; + background-color: #11303d; + color: #000; + margin: 0; + padding: 0; +} + +div.document { + display: flex; + background-color: #1c4e63; +} + +div.documentwrapper { + float: left; + width: 100%; +} + +div.bodywrapper { + margin: 0 0 0 230px; +} + +div.body { + background-color: #ffffff; + color: #000000; + padding: 0 20px 30px 20px; +} + +div.footer { + color: #ffffff; + width: 100%; + padding: 9px 0 9px 0; + text-align: center; + font-size: 75%; +} + +div.footer a { + color: #ffffff; + text-decoration: underline; +} + +div.related { + background-color: #133f52; + line-height: 30px; + color: #ffffff; +} + +div.related a { + color: #ffffff; +} + +div.sphinxsidebar { +} + +div.sphinxsidebar h3 { + font-family: 'Trebuchet MS', sans-serif; + color: #ffffff; + font-size: 1.4em; + font-weight: normal; + margin: 0; + padding: 0; +} + +div.sphinxsidebar h3 a { + color: #ffffff; +} + +div.sphinxsidebar h4 { + font-family: 'Trebuchet MS', sans-serif; + color: #ffffff; + font-size: 1.3em; + font-weight: normal; + margin: 5px 0 0 0; + padding: 0; +} + +div.sphinxsidebar p { + color: #ffffff; +} + +div.sphinxsidebar p.topless { + margin: 5px 10px 10px 10px; +} + +div.sphinxsidebar ul { + margin: 10px; + padding: 0; + color: #ffffff; +} + +div.sphinxsidebar a { + color: #98dbcc; +} + +div.sphinxsidebar input { + border: 1px solid #98dbcc; + font-family: sans-serif; + font-size: 1em; +} + + + +/* -- hyperlink styles ------------------------------------------------------ */ + +a { + color: #355f7c; + text-decoration: none; +} + +a:visited { + color: #551a8b; + text-decoration: none; +} + +a:hover { + text-decoration: underline; +} + + + +/* -- body styles ----------------------------------------------------------- */ + +div.body h1, +div.body h2, +div.body h3, +div.body h4, +div.body h5, +div.body h6 { + font-family: 'Trebuchet MS', sans-serif; + background-color: #f2f2f2; + font-weight: normal; + color: #20435c; + border-bottom: 1px solid #ccc; + margin: 20px -20px 10px -20px; + padding: 3px 0 3px 10px; +} + +div.body h1 { margin-top: 0; font-size: 200%; } +div.body h2 { font-size: 160%; } +div.body h3 { font-size: 140%; } +div.body h4 { font-size: 120%; } +div.body h5 { font-size: 110%; } +div.body h6 { font-size: 100%; } + +a.headerlink { + color: #c60f0f; + font-size: 0.8em; + padding: 0 4px 0 4px; + text-decoration: none; +} + +a.headerlink:hover { + background-color: #c60f0f; + color: white; +} + +div.body p, div.body dd, div.body li, div.body blockquote { + text-align: justify; + line-height: 130%; +} + +div.admonition p.admonition-title + p { + display: inline; +} + +div.admonition p { + margin-bottom: 5px; +} + +div.admonition pre { + margin-bottom: 5px; +} + +div.admonition ul, div.admonition ol { + margin-bottom: 5px; +} + +div.note { + background-color: #eee; + border: 1px solid #ccc; +} + +div.seealso { + background-color: #ffc; + border: 1px solid #ff6; +} + +nav.contents, +aside.topic, +div.topic { + background-color: #eee; +} + +div.warning { + background-color: #ffe4e4; + border: 1px solid #f66; +} + +p.admonition-title { + display: inline; +} + +p.admonition-title:after { + content: ":"; +} + +pre { + padding: 5px; + background-color: unset; + color: unset; + line-height: 120%; + border: 1px solid #ac9; + border-left: none; + border-right: none; +} + +code { + background-color: #ecf0f3; + padding: 0 1px 0 1px; + font-size: 0.95em; +} + +th, dl.field-list > dt { + background-color: #ede; +} + +.warning code { + background: #efc2c2; +} + +.note code { + background: #d6d6d6; +} + +.viewcode-back { + font-family: sans-serif; +} + +div.viewcode-block:target { + background-color: #f4debf; + border-top: 1px solid #ac9; + border-bottom: 1px solid #ac9; +} + +div.code-block-caption { + color: #efefef; + background-color: #1c4e63; +} \ No newline at end of file diff --git a/_static/default.css b/_static/default.css new file mode 100644 index 0000000..81b9363 --- /dev/null +++ b/_static/default.css @@ -0,0 +1 @@ +@import url("classic.css"); diff --git a/_static/doctools.js b/_static/doctools.js index d4619fd..4d67807 100644 --- a/_static/doctools.js +++ b/_static/doctools.js @@ -2,246 +2,155 @@ * doctools.js * ~~~~~~~~~~~ * - * Sphinx JavaScript utilities for all documentation. + * Base JavaScript utilities for all Sphinx HTML documentation. * - * :copyright: Copyright 2007-2011 by the Sphinx team, see AUTHORS. + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. * :license: BSD, see LICENSE for details. * */ - -/** - * select a different prefix for underscore - */ -$u = _.noConflict(); - -/** - * make the code below compatible with browsers without - * an installed firebug like debugger -if (!window.console || !console.firebug) { - var names = ["log", "debug", "info", "warn", "error", "assert", "dir", - "dirxml", "group", "groupEnd", "time", "timeEnd", "count", "trace", - "profile", "profileEnd"]; - window.console = {}; - for (var i = 0; i < names.length; ++i) - window.console[names[i]] = function() {}; -} - */ - -/** - * small helper function to urldecode strings - */ -jQuery.urldecode = function(x) { - return decodeURIComponent(x).replace(/\+/g, ' '); -} - -/** - * small helper function to urlencode strings - */ -jQuery.urlencode = encodeURIComponent; - -/** - * This function returns the parsed url parameters of the - * current request. Multiple values per key are supported, - * it will always return arrays of strings for the value parts. - */ -jQuery.getQueryParameters = function(s) { - if (typeof s == 'undefined') - s = document.location.search; - var parts = s.substr(s.indexOf('?') + 1).split('&'); - var result = {}; - for (var i = 0; i < parts.length; i++) { - var tmp = parts[i].split('=', 2); - var key = jQuery.urldecode(tmp[0]); - var value = jQuery.urldecode(tmp[1]); - if (key in result) - result[key].push(value); - else - result[key] = [value]; - } - return result; -}; - -/** - * small function to check if an array contains - * a given item. - */ -jQuery.contains = function(arr, item) { - for (var i = 0; i < arr.length; i++) { - if (arr[i] == item) - return true; - } - return false; -}; - -/** - * highlight a given string on a jquery object by wrapping it in - * span elements with the given class name. - */ -jQuery.fn.highlightText = function(text, className) { - function highlight(node) { - if (node.nodeType == 3) { - var val = node.nodeValue; - var pos = val.toLowerCase().indexOf(text); - if (pos >= 0 && !jQuery(node.parentNode).hasClass(className)) { - var span = document.createElement("span"); - span.className = className; - span.appendChild(document.createTextNode(val.substr(pos, text.length))); - node.parentNode.insertBefore(span, node.parentNode.insertBefore( - document.createTextNode(val.substr(pos + text.length)), - node.nextSibling)); - node.nodeValue = val.substr(0, pos); - } - } - else if (!jQuery(node).is("button, select, textarea")) { - jQuery.each(node.childNodes, function() { - highlight(this); - }); - } +"use strict"; + +const BLACKLISTED_KEY_CONTROL_ELEMENTS = new Set([ + "TEXTAREA", + "INPUT", + "SELECT", + "BUTTON", +]); + +const _ready = (callback) => { + if (document.readyState !== "loading") { + callback(); + } else { + document.addEventListener("DOMContentLoaded", callback); } - return this.each(function() { - highlight(this); - }); }; /** * Small JavaScript module for the documentation. */ -var Documentation = { - - init : function() { - this.fixFirefoxAnchorBug(); - this.highlightSearchWords(); - this.initIndexTable(); +const Documentation = { + init: () => { + Documentation.initDomainIndexTable(); + Documentation.initOnKeyListeners(); }, /** * i18n support */ - TRANSLATIONS : {}, - PLURAL_EXPR : function(n) { return n == 1 ? 0 : 1; }, - LOCALE : 'unknown', + TRANSLATIONS: {}, + PLURAL_EXPR: (n) => (n === 1 ? 0 : 1), + LOCALE: "unknown", // gettext and ngettext don't access this so that the functions // can safely bound to a different name (_ = Documentation.gettext) - gettext : function(string) { - var translated = Documentation.TRANSLATIONS[string]; - if (typeof translated == 'undefined') - return string; - return (typeof translated == 'string') ? translated : translated[0]; - }, - - ngettext : function(singular, plural, n) { - var translated = Documentation.TRANSLATIONS[singular]; - if (typeof translated == 'undefined') - return (n == 1) ? singular : plural; - return translated[Documentation.PLURALEXPR(n)]; - }, - - addTranslations : function(catalog) { - for (var key in catalog.messages) - this.TRANSLATIONS[key] = catalog.messages[key]; - this.PLURAL_EXPR = new Function('n', 'return +(' + catalog.plural_expr + ')'); - this.LOCALE = catalog.locale; - }, - - /** - * add context elements like header anchor links - */ - addContextElements : function() { - $('div[id] > :header:first').each(function() { - $('\u00B6'). - attr('href', '#' + this.id). - attr('title', _('Permalink to this headline')). - appendTo(this); - }); - $('dt[id]').each(function() { - $('\u00B6'). - attr('href', '#' + this.id). - attr('title', _('Permalink to this definition')). - appendTo(this); - }); + gettext: (string) => { + const translated = Documentation.TRANSLATIONS[string]; + switch (typeof translated) { + case "undefined": + return string; // no translation + case "string": + return translated; // translation exists + default: + return translated[0]; // (singular, plural) translation tuple exists + } }, - /** - * workaround a firefox stupidity - */ - fixFirefoxAnchorBug : function() { - if (document.location.hash && $.browser.mozilla) - window.setTimeout(function() { - document.location.href += ''; - }, 10); + ngettext: (singular, plural, n) => { + const translated = Documentation.TRANSLATIONS[singular]; + if (typeof translated !== "undefined") + return translated[Documentation.PLURAL_EXPR(n)]; + return n === 1 ? singular : plural; }, - /** - * highlight the search words provided in the url in the text - */ - highlightSearchWords : function() { - var params = $.getQueryParameters(); - var terms = (params.highlight) ? params.highlight[0].split(/\s+/) : []; - if (terms.length) { - var body = $('div.body'); - window.setTimeout(function() { - $.each(terms, function() { - body.highlightText(this.toLowerCase(), 'highlighted'); - }); - }, 10); - $('

' + _('Hide Search Matches') + '

') - .appendTo($('#searchbox')); - } + addTranslations: (catalog) => { + Object.assign(Documentation.TRANSLATIONS, catalog.messages); + Documentation.PLURAL_EXPR = new Function( + "n", + `return (${catalog.plural_expr})` + ); + Documentation.LOCALE = catalog.locale; }, /** - * init the domain index toggle buttons + * helper function to focus on search bar */ - initIndexTable : function() { - var togglers = $('img.toggler').click(function() { - var src = $(this).attr('src'); - var idnum = $(this).attr('id').substr(7); - $('tr.cg-' + idnum).toggle(); - if (src.substr(-9) == 'minus.png') - $(this).attr('src', src.substr(0, src.length-9) + 'plus.png'); - else - $(this).attr('src', src.substr(0, src.length-8) + 'minus.png'); - }).css('display', ''); - if (DOCUMENTATION_OPTIONS.COLLAPSE_INDEX) { - togglers.click(); - } + focusSearchBar: () => { + document.querySelectorAll("input[name=q]")[0]?.focus(); }, /** - * helper function to hide the search marks again + * Initialise the domain index toggle buttons */ - hideSearchWords : function() { - $('#searchbox .highlight-link').fadeOut(300); - $('span.highlighted').removeClass('highlighted'); + initDomainIndexTable: () => { + const toggler = (el) => { + const idNumber = el.id.substr(7); + const toggledRows = document.querySelectorAll(`tr.cg-${idNumber}`); + if (el.src.substr(-9) === "minus.png") { + el.src = `${el.src.substr(0, el.src.length - 9)}plus.png`; + toggledRows.forEach((el) => (el.style.display = "none")); + } else { + el.src = `${el.src.substr(0, el.src.length - 8)}minus.png`; + toggledRows.forEach((el) => (el.style.display = "")); + } + }; + + const togglerElements = document.querySelectorAll("img.toggler"); + togglerElements.forEach((el) => + el.addEventListener("click", (event) => toggler(event.currentTarget)) + ); + togglerElements.forEach((el) => (el.style.display = "")); + if (DOCUMENTATION_OPTIONS.COLLAPSE_INDEX) togglerElements.forEach(toggler); }, - /** - * make the url absolute - */ - makeURL : function(relativeURL) { - return DOCUMENTATION_OPTIONS.URL_ROOT + '/' + relativeURL; - }, + initOnKeyListeners: () => { + // only install a listener if it is really needed + if ( + !DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS && + !DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS + ) + return; + + document.addEventListener("keydown", (event) => { + // bail for input elements + if (BLACKLISTED_KEY_CONTROL_ELEMENTS.has(document.activeElement.tagName)) return; + // bail with special keys + if (event.altKey || event.ctrlKey || event.metaKey) return; + + if (!event.shiftKey) { + switch (event.key) { + case "ArrowLeft": + if (!DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) break; + + const prevLink = document.querySelector('link[rel="prev"]'); + if (prevLink && prevLink.href) { + window.location.href = prevLink.href; + event.preventDefault(); + } + break; + case "ArrowRight": + if (!DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) break; + + const nextLink = document.querySelector('link[rel="next"]'); + if (nextLink && nextLink.href) { + window.location.href = nextLink.href; + event.preventDefault(); + } + break; + } + } - /** - * get the current relative url - */ - getCurrentURL : function() { - var path = document.location.pathname; - var parts = path.split(/\//); - $.each(DOCUMENTATION_OPTIONS.URL_ROOT.split(/\//), function() { - if (this == '..') - parts.pop(); + // some keyboard layouts may need Shift to get / + switch (event.key) { + case "/": + if (!DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS) break; + Documentation.focusSearchBar(); + event.preventDefault(); + } }); - var url = parts.join('/'); - return path.substring(url.lastIndexOf('/') + 1, path.length - 1); - } + }, }; // quick alias for translations -_ = Documentation.gettext; +const _ = Documentation.gettext; -$(document).ready(function() { - Documentation.init(); -}); +_ready(Documentation.init); diff --git a/_static/documentation_options.js b/_static/documentation_options.js new file mode 100644 index 0000000..e21c068 --- /dev/null +++ b/_static/documentation_options.js @@ -0,0 +1,13 @@ +const DOCUMENTATION_OPTIONS = { + VERSION: '0.1', + LANGUAGE: 'en', + COLLAPSE_INDEX: false, + BUILDER: 'html', + FILE_SUFFIX: '.html', + LINK_SUFFIX: '.html', + HAS_SOURCE: true, + SOURCELINK_SUFFIX: '.txt', + NAVIGATION_WITH_KEYS: false, + SHOW_SEARCH_SUMMARY: true, + ENABLE_SEARCH_SHORTCUTS: true, +}; \ No newline at end of file diff --git a/_static/file.png b/_static/file.png index d18082e..a858a41 100644 Binary files a/_static/file.png and b/_static/file.png differ diff --git a/_static/language_data.js b/_static/language_data.js new file mode 100644 index 0000000..367b8ed --- /dev/null +++ b/_static/language_data.js @@ -0,0 +1,199 @@ +/* + * language_data.js + * ~~~~~~~~~~~~~~~~ + * + * This script contains the language-specific data used by searchtools.js, + * namely the list of stopwords, stemmer, scorer and splitter. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * + */ + +var stopwords = ["a", "and", "are", "as", "at", "be", "but", "by", "for", "if", "in", "into", "is", "it", "near", "no", "not", "of", "on", "or", "such", "that", "the", "their", "then", "there", "these", "they", "this", "to", "was", "will", "with"]; + + +/* Non-minified version is copied as a separate JS file, if available */ + +/** + * Porter Stemmer + */ +var Stemmer = function() { + + var step2list = { + ational: 'ate', + tional: 'tion', + enci: 'ence', + anci: 'ance', + izer: 'ize', + bli: 'ble', + alli: 'al', + entli: 'ent', + eli: 'e', + ousli: 'ous', + ization: 'ize', + ation: 'ate', + ator: 'ate', + alism: 'al', + iveness: 'ive', + fulness: 'ful', + ousness: 'ous', + aliti: 'al', + iviti: 'ive', + biliti: 'ble', + logi: 'log' + }; + + var step3list = { + icate: 'ic', + ative: '', + alize: 'al', + iciti: 'ic', + ical: 'ic', + ful: '', + ness: '' + }; + + var c = "[^aeiou]"; // consonant + var v = "[aeiouy]"; // vowel + var C = c + "[^aeiouy]*"; // consonant sequence + var V = v + "[aeiou]*"; // vowel sequence + + var mgr0 = "^(" + C + ")?" + V + C; // [C]VC... is m>0 + var meq1 = "^(" + C + ")?" + V + C + "(" + V + ")?$"; // [C]VC[V] is m=1 + var mgr1 = "^(" + C + ")?" + V + C + V + C; // [C]VCVC... is m>1 + var s_v = "^(" + C + ")?" + v; // vowel in stem + + this.stemWord = function (w) { + var stem; + var suffix; + var firstch; + var origword = w; + + if (w.length < 3) + return w; + + var re; + var re2; + var re3; + var re4; + + firstch = w.substr(0,1); + if (firstch == "y") + w = firstch.toUpperCase() + w.substr(1); + + // Step 1a + re = /^(.+?)(ss|i)es$/; + re2 = /^(.+?)([^s])s$/; + + if (re.test(w)) + w = w.replace(re,"$1$2"); + else if (re2.test(w)) + w = w.replace(re2,"$1$2"); + + // Step 1b + re = /^(.+?)eed$/; + re2 = /^(.+?)(ed|ing)$/; + if (re.test(w)) { + var fp = re.exec(w); + re = new RegExp(mgr0); + if (re.test(fp[1])) { + re = /.$/; + w = w.replace(re,""); + } + } + else if (re2.test(w)) { + var fp = re2.exec(w); + stem = fp[1]; + re2 = new RegExp(s_v); + if (re2.test(stem)) { + w = stem; + re2 = /(at|bl|iz)$/; + re3 = new RegExp("([^aeiouylsz])\\1$"); + re4 = new RegExp("^" + C + v + "[^aeiouwxy]$"); + if (re2.test(w)) + w = w + "e"; + else if (re3.test(w)) { + re = /.$/; + w = w.replace(re,""); + } + else if (re4.test(w)) + w = w + "e"; + } + } + + // Step 1c + re = /^(.+?)y$/; + if (re.test(w)) { + var fp = re.exec(w); + stem = fp[1]; + re = new RegExp(s_v); + if (re.test(stem)) + w = stem + "i"; + } + + // Step 2 + re = /^(.+?)(ational|tional|enci|anci|izer|bli|alli|entli|eli|ousli|ization|ation|ator|alism|iveness|fulness|ousness|aliti|iviti|biliti|logi)$/; + if (re.test(w)) { + var fp = re.exec(w); + stem = fp[1]; + suffix = fp[2]; + re = new RegExp(mgr0); + if (re.test(stem)) + w = stem + step2list[suffix]; + } + + // Step 3 + re = /^(.+?)(icate|ative|alize|iciti|ical|ful|ness)$/; + if (re.test(w)) { + var fp = re.exec(w); + stem = fp[1]; + suffix = fp[2]; + re = new RegExp(mgr0); + if (re.test(stem)) + w = stem + step3list[suffix]; + } + + // Step 4 + re = /^(.+?)(al|ance|ence|er|ic|able|ible|ant|ement|ment|ent|ou|ism|ate|iti|ous|ive|ize)$/; + re2 = /^(.+?)(s|t)(ion)$/; + if (re.test(w)) { + var fp = re.exec(w); + stem = fp[1]; + re = new RegExp(mgr1); + if (re.test(stem)) + w = stem; + } + else if (re2.test(w)) { + var fp = re2.exec(w); + stem = fp[1] + fp[2]; + re2 = new RegExp(mgr1); + if (re2.test(stem)) + w = stem; + } + + // Step 5 + re = /^(.+?)e$/; + if (re.test(w)) { + var fp = re.exec(w); + stem = fp[1]; + re = new RegExp(mgr1); + re2 = new RegExp(meq1); + re3 = new RegExp("^" + C + v + "[^aeiouwxy]$"); + if (re.test(stem) || (re2.test(stem) && !(re3.test(stem)))) + w = stem; + } + re = /ll$/; + re2 = new RegExp(mgr1); + if (re.test(w) && re2.test(w)) { + re = /.$/; + w = w.replace(re,""); + } + + // and turn initial Y back to y + if (firstch == "y") + w = firstch.toLowerCase() + w.substr(1); + return w; + } +} + diff --git a/_static/minus.png b/_static/minus.png index da1c562..d96755f 100644 Binary files a/_static/minus.png and b/_static/minus.png differ diff --git a/_static/plus.png b/_static/plus.png index b3cb374..7107cec 100644 Binary files a/_static/plus.png and b/_static/plus.png differ diff --git a/_static/pygments.css b/_static/pygments.css index d79caa1..0d49244 100644 --- a/_static/pygments.css +++ b/_static/pygments.css @@ -1,15 +1,23 @@ +pre { line-height: 125%; } +td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } +span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } +td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } +span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } .highlight .hll { background-color: #ffffcc } -.highlight { background: #eeffcc; } +.highlight { background: #eeffcc; } .highlight .c { color: #408090; font-style: italic } /* Comment */ .highlight .err { border: 1px solid #FF0000 } /* Error */ .highlight .k { color: #007020; font-weight: bold } /* Keyword */ .highlight .o { color: #666666 } /* Operator */ +.highlight .ch { color: #408090; font-style: italic } /* Comment.Hashbang */ .highlight .cm { color: #408090; font-style: italic } /* Comment.Multiline */ .highlight .cp { color: #007020 } /* Comment.Preproc */ +.highlight .cpf { color: #408090; font-style: italic } /* Comment.PreprocFile */ .highlight .c1 { color: #408090; font-style: italic } /* Comment.Single */ .highlight .cs { color: #408090; background-color: #fff0f0 } /* Comment.Special */ .highlight .gd { color: #A00000 } /* Generic.Deleted */ .highlight .ge { font-style: italic } /* Generic.Emph */ +.highlight .ges { font-weight: bold; font-style: italic } /* Generic.EmphStrong */ .highlight .gr { color: #FF0000 } /* Generic.Error */ .highlight .gh { color: #000080; font-weight: bold } /* Generic.Heading */ .highlight .gi { color: #00A000 } /* Generic.Inserted */ @@ -40,12 +48,15 @@ .highlight .nv { color: #bb60d5 } /* Name.Variable */ .highlight .ow { color: #007020; font-weight: bold } /* Operator.Word */ .highlight .w { color: #bbbbbb } /* Text.Whitespace */ +.highlight .mb { color: #208050 } /* Literal.Number.Bin */ .highlight .mf { color: #208050 } /* Literal.Number.Float */ .highlight .mh { color: #208050 } /* Literal.Number.Hex */ .highlight .mi { color: #208050 } /* Literal.Number.Integer */ .highlight .mo { color: #208050 } /* Literal.Number.Oct */ +.highlight .sa { color: #4070a0 } /* Literal.String.Affix */ .highlight .sb { color: #4070a0 } /* Literal.String.Backtick */ .highlight .sc { color: #4070a0 } /* Literal.String.Char */ +.highlight .dl { color: #4070a0 } /* Literal.String.Delimiter */ .highlight .sd { color: #4070a0; font-style: italic } /* Literal.String.Doc */ .highlight .s2 { color: #4070a0 } /* Literal.String.Double */ .highlight .se { color: #4070a0; font-weight: bold } /* Literal.String.Escape */ @@ -56,7 +67,9 @@ .highlight .s1 { color: #4070a0 } /* Literal.String.Single */ .highlight .ss { color: #517918 } /* Literal.String.Symbol */ .highlight .bp { color: #007020 } /* Name.Builtin.Pseudo */ +.highlight .fm { color: #06287e } /* Name.Function.Magic */ .highlight .vc { color: #bb60d5 } /* Name.Variable.Class */ .highlight .vg { color: #bb60d5 } /* Name.Variable.Global */ .highlight .vi { color: #bb60d5 } /* Name.Variable.Instance */ +.highlight .vm { color: #bb60d5 } /* Name.Variable.Magic */ .highlight .il { color: #208050 } /* Literal.Number.Integer.Long */ \ No newline at end of file diff --git a/_static/searchtools.js b/_static/searchtools.js index 663be4c..92da3f8 100644 --- a/_static/searchtools.js +++ b/_static/searchtools.js @@ -1,560 +1,619 @@ /* - * searchtools.js_t + * searchtools.js * ~~~~~~~~~~~~~~~~ * - * Sphinx JavaScript utilties for the full-text search. + * Sphinx JavaScript utilities for the full-text search. * - * :copyright: Copyright 2007-2011 by the Sphinx team, see AUTHORS. + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. * :license: BSD, see LICENSE for details. * */ +"use strict"; /** - * helper function to return a node containing the - * search summary for a given text. keywords is a list - * of stemmed words, hlwords is the list of normal, unstemmed - * words. the first one is used to find the occurance, the - * latter for highlighting it. + * Simple result scoring code. */ - -jQuery.makeSearchSummary = function(text, keywords, hlwords) { - var textLower = text.toLowerCase(); - var start = 0; - $.each(keywords, function() { - var i = textLower.indexOf(this.toLowerCase()); - if (i > -1) - start = i; - }); - start = Math.max(start - 120, 0); - var excerpt = ((start > 0) ? '...' : '') + - $.trim(text.substr(start, 240)) + - ((start + 240 - text.length) ? '...' : ''); - var rv = $('
').text(excerpt); - $.each(hlwords, function() { - rv = rv.highlightText(this, 'highlighted'); - }); - return rv; +if (typeof Scorer === "undefined") { + var Scorer = { + // Implement the following function to further tweak the score for each result + // The function takes a result array [docname, title, anchor, descr, score, filename] + // and returns the new score. + /* + score: result => { + const [docname, title, anchor, descr, score, filename] = result + return score + }, + */ + + // query matches the full name of an object + objNameMatch: 11, + // or matches in the last dotted part of the object name + objPartialMatch: 6, + // Additive scores depending on the priority of the object + objPrio: { + 0: 15, // used to be importantResults + 1: 5, // used to be objectResults + 2: -5, // used to be unimportantResults + }, + // Used when the priority is not in the mapping. + objPrioDefault: 0, + + // query found in title + title: 15, + partialTitle: 7, + // query found in terms + term: 5, + partialTerm: 2, + }; } +const _removeChildren = (element) => { + while (element && element.lastChild) element.removeChild(element.lastChild); +}; /** - * Porter Stemmer + * See https://developer.mozilla.org/en-US/docs/Web/JavaScript/Guide/Regular_Expressions#escaping */ -var Stemmer = function() { - - var step2list = { - ational: 'ate', - tional: 'tion', - enci: 'ence', - anci: 'ance', - izer: 'ize', - bli: 'ble', - alli: 'al', - entli: 'ent', - eli: 'e', - ousli: 'ous', - ization: 'ize', - ation: 'ate', - ator: 'ate', - alism: 'al', - iveness: 'ive', - fulness: 'ful', - ousness: 'ous', - aliti: 'al', - iviti: 'ive', - biliti: 'ble', - logi: 'log' - }; - - var step3list = { - icate: 'ic', - ative: '', - alize: 'al', - iciti: 'ic', - ical: 'ic', - ful: '', - ness: '' - }; - - var c = "[^aeiou]"; // consonant - var v = "[aeiouy]"; // vowel - var C = c + "[^aeiouy]*"; // consonant sequence - var V = v + "[aeiou]*"; // vowel sequence - - var mgr0 = "^(" + C + ")?" + V + C; // [C]VC... is m>0 - var meq1 = "^(" + C + ")?" + V + C + "(" + V + ")?$"; // [C]VC[V] is m=1 - var mgr1 = "^(" + C + ")?" + V + C + V + C; // [C]VCVC... is m>1 - var s_v = "^(" + C + ")?" + v; // vowel in stem - - this.stemWord = function (w) { - var stem; - var suffix; - var firstch; - var origword = w; - - if (w.length < 3) - return w; - - var re; - var re2; - var re3; - var re4; - - firstch = w.substr(0,1); - if (firstch == "y") - w = firstch.toUpperCase() + w.substr(1); - - // Step 1a - re = /^(.+?)(ss|i)es$/; - re2 = /^(.+?)([^s])s$/; - - if (re.test(w)) - w = w.replace(re,"$1$2"); - else if (re2.test(w)) - w = w.replace(re2,"$1$2"); - - // Step 1b - re = /^(.+?)eed$/; - re2 = /^(.+?)(ed|ing)$/; - if (re.test(w)) { - var fp = re.exec(w); - re = new RegExp(mgr0); - if (re.test(fp[1])) { - re = /.$/; - w = w.replace(re,""); - } - } - else if (re2.test(w)) { - var fp = re2.exec(w); - stem = fp[1]; - re2 = new RegExp(s_v); - if (re2.test(stem)) { - w = stem; - re2 = /(at|bl|iz)$/; - re3 = new RegExp("([^aeiouylsz])\\1$"); - re4 = new RegExp("^" + C + v + "[^aeiouwxy]$"); - if (re2.test(w)) - w = w + "e"; - else if (re3.test(w)) { - re = /.$/; - w = w.replace(re,""); - } - else if (re4.test(w)) - w = w + "e"; - } - } - - // Step 1c - re = /^(.+?)y$/; - if (re.test(w)) { - var fp = re.exec(w); - stem = fp[1]; - re = new RegExp(s_v); - if (re.test(stem)) - w = stem + "i"; - } - - // Step 2 - re = /^(.+?)(ational|tional|enci|anci|izer|bli|alli|entli|eli|ousli|ization|ation|ator|alism|iveness|fulness|ousness|aliti|iviti|biliti|logi)$/; - if (re.test(w)) { - var fp = re.exec(w); - stem = fp[1]; - suffix = fp[2]; - re = new RegExp(mgr0); - if (re.test(stem)) - w = stem + step2list[suffix]; - } +const _escapeRegExp = (string) => + string.replace(/[.*+\-?^${}()|[\]\\]/g, "\\$&"); // $& means the whole matched string + +const _displayItem = (item, searchTerms, highlightTerms) => { + const docBuilder = DOCUMENTATION_OPTIONS.BUILDER; + const docFileSuffix = DOCUMENTATION_OPTIONS.FILE_SUFFIX; + const docLinkSuffix = DOCUMENTATION_OPTIONS.LINK_SUFFIX; + const showSearchSummary = DOCUMENTATION_OPTIONS.SHOW_SEARCH_SUMMARY; + const contentRoot = document.documentElement.dataset.content_root; + + const [docName, title, anchor, descr, score, _filename] = item; + + let listItem = document.createElement("li"); + let requestUrl; + let linkUrl; + if (docBuilder === "dirhtml") { + // dirhtml builder + let dirname = docName + "/"; + if (dirname.match(/\/index\/$/)) + dirname = dirname.substring(0, dirname.length - 6); + else if (dirname === "index/") dirname = ""; + requestUrl = contentRoot + dirname; + linkUrl = requestUrl; + } else { + // normal html builders + requestUrl = contentRoot + docName + docFileSuffix; + linkUrl = docName + docLinkSuffix; + } + let linkEl = listItem.appendChild(document.createElement("a")); + linkEl.href = linkUrl + anchor; + linkEl.dataset.score = score; + linkEl.innerHTML = title; + if (descr) { + listItem.appendChild(document.createElement("span")).innerHTML = + " (" + descr + ")"; + // highlight search terms in the description + if (SPHINX_HIGHLIGHT_ENABLED) // set in sphinx_highlight.js + highlightTerms.forEach((term) => _highlightText(listItem, term, "highlighted")); + } + else if (showSearchSummary) + fetch(requestUrl) + .then((responseData) => responseData.text()) + .then((data) => { + if (data) + listItem.appendChild( + Search.makeSearchSummary(data, searchTerms, anchor) + ); + // highlight search terms in the summary + if (SPHINX_HIGHLIGHT_ENABLED) // set in sphinx_highlight.js + highlightTerms.forEach((term) => _highlightText(listItem, term, "highlighted")); + }); + Search.output.appendChild(listItem); +}; +const _finishSearch = (resultCount) => { + Search.stopPulse(); + Search.title.innerText = _("Search Results"); + if (!resultCount) + Search.status.innerText = Documentation.gettext( + "Your search did not match any documents. Please make sure that all words are spelled correctly and that you've selected enough categories." + ); + else + Search.status.innerText = _( + "Search finished, found ${resultCount} page(s) matching the search query." + ).replace('${resultCount}', resultCount); +}; +const _displayNextItem = ( + results, + resultCount, + searchTerms, + highlightTerms, +) => { + // results left, load the summary and display it + // this is intended to be dynamic (don't sub resultsCount) + if (results.length) { + _displayItem(results.pop(), searchTerms, highlightTerms); + setTimeout( + () => _displayNextItem(results, resultCount, searchTerms, highlightTerms), + 5 + ); + } + // search finished, update title and status message + else _finishSearch(resultCount); +}; +// Helper function used by query() to order search results. +// Each input is an array of [docname, title, anchor, descr, score, filename]. +// Order the results by score (in opposite order of appearance, since the +// `_displayNextItem` function uses pop() to retrieve items) and then alphabetically. +const _orderResultsByScoreThenName = (a, b) => { + const leftScore = a[4]; + const rightScore = b[4]; + if (leftScore === rightScore) { + // same score: sort alphabetically + const leftTitle = a[1].toLowerCase(); + const rightTitle = b[1].toLowerCase(); + if (leftTitle === rightTitle) return 0; + return leftTitle > rightTitle ? -1 : 1; // inverted is intentional + } + return leftScore > rightScore ? 1 : -1; +}; - // Step 3 - re = /^(.+?)(icate|ative|alize|iciti|ical|ful|ness)$/; - if (re.test(w)) { - var fp = re.exec(w); - stem = fp[1]; - suffix = fp[2]; - re = new RegExp(mgr0); - if (re.test(stem)) - w = stem + step3list[suffix]; - } +/** + * Default splitQuery function. Can be overridden in ``sphinx.search`` with a + * custom function per language. + * + * The regular expression works by splitting the string on consecutive characters + * that are not Unicode letters, numbers, underscores, or emoji characters. + * This is the same as ``\W+`` in Python, preserving the surrogate pair area. + */ +if (typeof splitQuery === "undefined") { + var splitQuery = (query) => query + .split(/[^\p{Letter}\p{Number}_\p{Emoji_Presentation}]+/gu) + .filter(term => term) // remove remaining empty strings +} - // Step 4 - re = /^(.+?)(al|ance|ence|er|ic|able|ible|ant|ement|ment|ent|ou|ism|ate|iti|ous|ive|ize)$/; - re2 = /^(.+?)(s|t)(ion)$/; - if (re.test(w)) { - var fp = re.exec(w); - stem = fp[1]; - re = new RegExp(mgr1); - if (re.test(stem)) - w = stem; - } - else if (re2.test(w)) { - var fp = re2.exec(w); - stem = fp[1] + fp[2]; - re2 = new RegExp(mgr1); - if (re2.test(stem)) - w = stem; +/** + * Search Module + */ +const Search = { + _index: null, + _queued_query: null, + _pulse_status: -1, + + htmlToText: (htmlString, anchor) => { + const htmlElement = new DOMParser().parseFromString(htmlString, 'text/html'); + for (const removalQuery of [".headerlinks", "script", "style"]) { + htmlElement.querySelectorAll(removalQuery).forEach((el) => { el.remove() }); } + if (anchor) { + const anchorContent = htmlElement.querySelector(`[role="main"] ${anchor}`); + if (anchorContent) return anchorContent.textContent; - // Step 5 - re = /^(.+?)e$/; - if (re.test(w)) { - var fp = re.exec(w); - stem = fp[1]; - re = new RegExp(mgr1); - re2 = new RegExp(meq1); - re3 = new RegExp("^" + C + v + "[^aeiouwxy]$"); - if (re.test(stem) || (re2.test(stem) && !(re3.test(stem)))) - w = stem; + console.warn( + `Anchored content block not found. Sphinx search tries to obtain it via DOM query '[role=main] ${anchor}'. Check your theme or template.` + ); } - re = /ll$/; - re2 = new RegExp(mgr1); - if (re.test(w) && re2.test(w)) { - re = /.$/; - w = w.replace(re,""); - } - - // and turn initial Y back to y - if (firstch == "y") - w = firstch.toLowerCase() + w.substr(1); - return w; - } -} + // if anchor not specified or not found, fall back to main content + const docContent = htmlElement.querySelector('[role="main"]'); + if (docContent) return docContent.textContent; -/** - * Search Module - */ -var Search = { - - _index : null, - _queued_query : null, - _pulse_status : -1, - - init : function() { - var params = $.getQueryParameters(); - if (params.q) { - var query = params.q[0]; - $('input[name="q"]')[0].value = query; - this.performSearch(query); - } + console.warn( + "Content block not found. Sphinx search tries to obtain it via DOM query '[role=main]'. Check your theme or template." + ); + return ""; }, - loadIndex : function(url) { - $.ajax({type: "GET", url: url, data: null, success: null, - dataType: "script", cache: true}); + init: () => { + const query = new URLSearchParams(window.location.search).get("q"); + document + .querySelectorAll('input[name="q"]') + .forEach((el) => (el.value = query)); + if (query) Search.performSearch(query); }, - setIndex : function(index) { - var q; - this._index = index; - if ((q = this._queued_query) !== null) { - this._queued_query = null; - Search.query(q); + loadIndex: (url) => + (document.body.appendChild(document.createElement("script")).src = url), + + setIndex: (index) => { + Search._index = index; + if (Search._queued_query !== null) { + const query = Search._queued_query; + Search._queued_query = null; + Search.query(query); } }, - hasIndex : function() { - return this._index !== null; - }, + hasIndex: () => Search._index !== null, - deferQuery : function(query) { - this._queued_query = query; - }, + deferQuery: (query) => (Search._queued_query = query), - stopPulse : function() { - this._pulse_status = 0; - }, + stopPulse: () => (Search._pulse_status = -1), - startPulse : function() { - if (this._pulse_status >= 0) - return; - function pulse() { + startPulse: () => { + if (Search._pulse_status >= 0) return; + + const pulse = () => { Search._pulse_status = (Search._pulse_status + 1) % 4; - var dotString = ''; - for (var i = 0; i < Search._pulse_status; i++) - dotString += '.'; - Search.dots.text(dotString); - if (Search._pulse_status > -1) - window.setTimeout(pulse, 500); + Search.dots.innerText = ".".repeat(Search._pulse_status); + if (Search._pulse_status >= 0) window.setTimeout(pulse, 500); }; pulse(); }, /** - * perform a search for something + * perform a search for something (or wait until index is loaded) */ - performSearch : function(query) { + performSearch: (query) => { // create the required interface elements - this.out = $('#search-results'); - this.title = $('

' + _('Searching') + '

').appendTo(this.out); - this.dots = $('').appendTo(this.title); - this.status = $('

').appendTo(this.out); - this.output = $('
    ').appendTo(this.out); - - $('#search-progress').text(_('Preparing search...')); - this.startPulse(); + const searchText = document.createElement("h2"); + searchText.textContent = _("Searching"); + const searchSummary = document.createElement("p"); + searchSummary.classList.add("search-summary"); + searchSummary.innerText = ""; + const searchList = document.createElement("ul"); + searchList.classList.add("search"); + + const out = document.getElementById("search-results"); + Search.title = out.appendChild(searchText); + Search.dots = Search.title.appendChild(document.createElement("span")); + Search.status = out.appendChild(searchSummary); + Search.output = out.appendChild(searchList); + + const searchProgress = document.getElementById("search-progress"); + // Some themes don't use the search progress node + if (searchProgress) { + searchProgress.innerText = _("Preparing search..."); + } + Search.startPulse(); // index already loaded, the browser was quick! - if (this.hasIndex()) - this.query(query); - else - this.deferQuery(query); + if (Search.hasIndex()) Search.query(query); + else Search.deferQuery(query); }, - query : function(query) { - var stopwords = ["and","then","into","it","as","are","in","if","for","no","there","their","was","is","be","to","that","but","they","not","such","with","by","a","on","these","of","will","this","near","the","or","at"]; - - // Stem the searchterms and add them to the correct list - var stemmer = new Stemmer(); - var searchterms = []; - var excluded = []; - var hlterms = []; - var tmp = query.split(/\s+/); - var objectterms = []; - for (var i = 0; i < tmp.length; i++) { - if (tmp[i] != "") { - objectterms.push(tmp[i].toLowerCase()); - } + _parseQuery: (query) => { + // stem the search terms and add them to the correct list + const stemmer = new Stemmer(); + const searchTerms = new Set(); + const excludedTerms = new Set(); + const highlightTerms = new Set(); + const objectTerms = new Set(splitQuery(query.toLowerCase().trim())); + splitQuery(query.trim()).forEach((queryTerm) => { + const queryTermLower = queryTerm.toLowerCase(); + + // maybe skip this "word" + // stopwords array is from language_data.js + if ( + stopwords.indexOf(queryTermLower) !== -1 || + queryTerm.match(/^\d+$/) + ) + return; - if ($u.indexOf(stopwords, tmp[i]) != -1 || tmp[i].match(/^\d+$/) || - tmp[i] == "") { - // skip this "word" - continue; - } // stem the word - var word = stemmer.stemWord(tmp[i]).toLowerCase(); + let word = stemmer.stemWord(queryTermLower); // select the correct list - if (word[0] == '-') { - var toAppend = excluded; - word = word.substr(1); - } + if (word[0] === "-") excludedTerms.add(word.substr(1)); else { - var toAppend = searchterms; - hlterms.push(tmp[i].toLowerCase()); + searchTerms.add(word); + highlightTerms.add(queryTermLower); } - // only add if not already in the list - if (!$.contains(toAppend, word)) - toAppend.push(word); - }; - var highlightstring = '?highlight=' + $.urlencode(hlterms.join(" ")); + }); - // console.debug('SEARCH: searching for:'); - // console.info('required: ', searchterms); - // console.info('excluded: ', excluded); + if (SPHINX_HIGHLIGHT_ENABLED) { // set in sphinx_highlight.js + localStorage.setItem("sphinx_highlight_terms", [...highlightTerms].join(" ")) + } - // prepare search - var filenames = this._index.filenames; - var titles = this._index.titles; - var terms = this._index.terms; - var fileMap = {}; - var files = null; - // different result priorities - var importantResults = []; - var objectResults = []; - var regularResults = []; - var unimportantResults = []; - $('#search-progress').empty(); + // console.debug("SEARCH: searching for:"); + // console.info("required: ", [...searchTerms]); + // console.info("excluded: ", [...excludedTerms]); - // lookup as object - for (var i = 0; i < objectterms.length; i++) { - var others = [].concat(objectterms.slice(0,i), - objectterms.slice(i+1, objectterms.length)) - var results = this.performObjectSearch(objectterms[i], others); - // Assume first word is most likely to be the object, - // other words more likely to be in description. - // Therefore put matches for earlier words first. - // (Results are eventually used in reverse order). - objectResults = results[0].concat(objectResults); - importantResults = results[1].concat(importantResults); - unimportantResults = results[2].concat(unimportantResults); - } + return [query, searchTerms, excludedTerms, highlightTerms, objectTerms]; + }, - // perform the search on the required terms - for (var i = 0; i < searchterms.length; i++) { - var word = searchterms[i]; - // no match but word was a required one - if ((files = terms[word]) == null) - break; - if (files.length == undefined) { - files = [files]; - } - // create the mapping - for (var j = 0; j < files.length; j++) { - var file = files[j]; - if (file in fileMap) - fileMap[file].push(word); - else - fileMap[file] = [word]; + /** + * execute search (requires search index to be loaded) + */ + _performSearch: (query, searchTerms, excludedTerms, highlightTerms, objectTerms) => { + const filenames = Search._index.filenames; + const docNames = Search._index.docnames; + const titles = Search._index.titles; + const allTitles = Search._index.alltitles; + const indexEntries = Search._index.indexentries; + + // Collect multiple result groups to be sorted separately and then ordered. + // Each is an array of [docname, title, anchor, descr, score, filename]. + const normalResults = []; + const nonMainIndexResults = []; + + _removeChildren(document.getElementById("search-progress")); + + const queryLower = query.toLowerCase().trim(); + for (const [title, foundTitles] of Object.entries(allTitles)) { + if (title.toLowerCase().trim().includes(queryLower) && (queryLower.length >= title.length/2)) { + for (const [file, id] of foundTitles) { + let score = Math.round(100 * queryLower.length / title.length) + normalResults.push([ + docNames[file], + titles[file] !== title ? `${titles[file]} > ${title}` : title, + id !== null ? "#" + id : "", + null, + score, + filenames[file], + ]); + } } } - // now check if the files don't contain excluded terms - for (var file in fileMap) { - var valid = true; - - // check if all requirements are matched - if (fileMap[file].length != searchterms.length) - continue; - - // ensure that none of the excluded terms is in the - // search result. - for (var i = 0; i < excluded.length; i++) { - if (terms[excluded[i]] == file || - $.contains(terms[excluded[i]] || [], file)) { - valid = false; - break; + // search for explicit entries in index directives + for (const [entry, foundEntries] of Object.entries(indexEntries)) { + if (entry.includes(queryLower) && (queryLower.length >= entry.length/2)) { + for (const [file, id, isMain] of foundEntries) { + const score = Math.round(100 * queryLower.length / entry.length); + const result = [ + docNames[file], + titles[file], + id ? "#" + id : "", + null, + score, + filenames[file], + ]; + if (isMain) { + normalResults.push(result); + } else { + nonMainIndexResults.push(result); + } } } + } + + // lookup as object + objectTerms.forEach((term) => + normalResults.push(...Search.performObjectSearch(term, objectTerms)) + ); + + // lookup as search terms in fulltext + normalResults.push(...Search.performTermsSearch(searchTerms, excludedTerms)); - // if we have still a valid result we can add it - // to the result list - if (valid) - regularResults.push([filenames[file], titles[file], '', null]); + // let the scorer override scores with a custom scoring function + if (Scorer.score) { + normalResults.forEach((item) => (item[4] = Scorer.score(item))); + nonMainIndexResults.forEach((item) => (item[4] = Scorer.score(item))); } - // delete unused variables in order to not waste - // memory until list is retrieved completely - delete filenames, titles, terms; + // Sort each group of results by score and then alphabetically by name. + normalResults.sort(_orderResultsByScoreThenName); + nonMainIndexResults.sort(_orderResultsByScoreThenName); + + // Combine the result groups in (reverse) order. + // Non-main index entries are typically arbitrary cross-references, + // so display them after other results. + let results = [...nonMainIndexResults, ...normalResults]; + + // remove duplicate search results + // note the reversing of results, so that in the case of duplicates, the highest-scoring entry is kept + let seen = new Set(); + results = results.reverse().reduce((acc, result) => { + let resultStr = result.slice(0, 4).concat([result[5]]).map(v => String(v)).join(','); + if (!seen.has(resultStr)) { + acc.push(result); + seen.add(resultStr); + } + return acc; + }, []); - // now sort the regular results descending by title - regularResults.sort(function(a, b) { - var left = a[1].toLowerCase(); - var right = b[1].toLowerCase(); - return (left > right) ? -1 : ((left < right) ? 1 : 0); - }); + return results.reverse(); + }, + + query: (query) => { + const [searchQuery, searchTerms, excludedTerms, highlightTerms, objectTerms] = Search._parseQuery(query); + const results = Search._performSearch(searchQuery, searchTerms, excludedTerms, highlightTerms, objectTerms); - // combine all results - var results = unimportantResults.concat(regularResults) - .concat(objectResults).concat(importantResults); + // for debugging + //Search.lastresults = results.slice(); // a copy + // console.info("search results:", Search.lastresults); // print the results - var resultCount = results.length; - function displayNextItem() { - // results left, load the summary and display it - if (results.length) { - var item = results.pop(); - var listItem = $('
  • '); - if (DOCUMENTATION_OPTIONS.FILE_SUFFIX == '') { - // dirhtml builder - var dirname = item[0] + '/'; - if (dirname.match(/\/index\/$/)) { - dirname = dirname.substring(0, dirname.length-6); - } else if (dirname == 'index/') { - dirname = ''; - } - listItem.append($('').attr('href', - DOCUMENTATION_OPTIONS.URL_ROOT + dirname + - highlightstring + item[2]).html(item[1])); - } else { - // normal html builders - listItem.append($('').attr('href', - item[0] + DOCUMENTATION_OPTIONS.FILE_SUFFIX + - highlightstring + item[2]).html(item[1])); - } - if (item[3]) { - listItem.append($(' (' + item[3] + ')')); - Search.output.append(listItem); - listItem.slideDown(5, function() { - displayNextItem(); + _displayNextItem(results, results.length, searchTerms, highlightTerms); + }, + + /** + * search for object names + */ + performObjectSearch: (object, objectTerms) => { + const filenames = Search._index.filenames; + const docNames = Search._index.docnames; + const objects = Search._index.objects; + const objNames = Search._index.objnames; + const titles = Search._index.titles; + + const results = []; + + const objectSearchCallback = (prefix, match) => { + const name = match[4] + const fullname = (prefix ? prefix + "." : "") + name; + const fullnameLower = fullname.toLowerCase(); + if (fullnameLower.indexOf(object) < 0) return; + + let score = 0; + const parts = fullnameLower.split("."); + + // check for different match types: exact matches of full name or + // "last name" (i.e. last dotted part) + if (fullnameLower === object || parts.slice(-1)[0] === object) + score += Scorer.objNameMatch; + else if (parts.slice(-1)[0].indexOf(object) > -1) + score += Scorer.objPartialMatch; // matches in last name + + const objName = objNames[match[1]][2]; + const title = titles[match[0]]; + + // If more than one term searched for, we require other words to be + // found in the name/title/description + const otherTerms = new Set(objectTerms); + otherTerms.delete(object); + if (otherTerms.size > 0) { + const haystack = `${prefix} ${name} ${objName} ${title}`.toLowerCase(); + if ( + [...otherTerms].some((otherTerm) => haystack.indexOf(otherTerm) < 0) + ) + return; + } + + let anchor = match[3]; + if (anchor === "") anchor = fullname; + else if (anchor === "-") anchor = objNames[match[1]][1] + "-" + fullname; + + const descr = objName + _(", in ") + title; + + // add custom score for some objects according to scorer + if (Scorer.objPrio.hasOwnProperty(match[2])) + score += Scorer.objPrio[match[2]]; + else score += Scorer.objPrioDefault; + + results.push([ + docNames[match[0]], + fullname, + "#" + anchor, + descr, + score, + filenames[match[0]], + ]); + }; + Object.keys(objects).forEach((prefix) => + objects[prefix].forEach((array) => + objectSearchCallback(prefix, array) + ) + ); + return results; + }, + + /** + * search for full-text terms in the index + */ + performTermsSearch: (searchTerms, excludedTerms) => { + // prepare search + const terms = Search._index.terms; + const titleTerms = Search._index.titleterms; + const filenames = Search._index.filenames; + const docNames = Search._index.docnames; + const titles = Search._index.titles; + + const scoreMap = new Map(); + const fileMap = new Map(); + + // perform the search on the required terms + searchTerms.forEach((word) => { + const files = []; + const arr = [ + { files: terms[word], score: Scorer.term }, + { files: titleTerms[word], score: Scorer.title }, + ]; + // add support for partial matches + if (word.length > 2) { + const escapedWord = _escapeRegExp(word); + if (!terms.hasOwnProperty(word)) { + Object.keys(terms).forEach((term) => { + if (term.match(escapedWord)) + arr.push({ files: terms[term], score: Scorer.partialTerm }); }); - } else if (DOCUMENTATION_OPTIONS.HAS_SOURCE) { - $.get(DOCUMENTATION_OPTIONS.URL_ROOT + '_sources/' + - item[0] + '.txt', function(data) { - if (data != '') { - listItem.append($.makeSearchSummary(data, searchterms, hlterms)); - Search.output.append(listItem); - } - listItem.slideDown(5, function() { - displayNextItem(); - }); - }, "text"); - } else { - // no source available, just display title - Search.output.append(listItem); - listItem.slideDown(5, function() { - displayNextItem(); + } + if (!titleTerms.hasOwnProperty(word)) { + Object.keys(titleTerms).forEach((term) => { + if (term.match(escapedWord)) + arr.push({ files: titleTerms[term], score: Scorer.partialTitle }); }); } } - // search finished, update title and status message - else { - Search.stopPulse(); - Search.title.text(_('Search Results')); - if (!resultCount) - Search.status.text(_('Your search did not match any documents. Please make sure that all words are spelled correctly and that you\'ve selected enough categories.')); - else - Search.status.text(_('Search finished, found %s page(s) matching the search query.').replace('%s', resultCount)); - Search.status.fadeIn(500); - } + + // no match but word was a required one + if (arr.every((record) => record.files === undefined)) return; + + // found search word in contents + arr.forEach((record) => { + if (record.files === undefined) return; + + let recordFiles = record.files; + if (recordFiles.length === undefined) recordFiles = [recordFiles]; + files.push(...recordFiles); + + // set score for the word in each file + recordFiles.forEach((file) => { + if (!scoreMap.has(file)) scoreMap.set(file, {}); + scoreMap.get(file)[word] = record.score; + }); + }); + + // create the mapping + files.forEach((file) => { + if (!fileMap.has(file)) fileMap.set(file, [word]); + else if (fileMap.get(file).indexOf(word) === -1) fileMap.get(file).push(word); + }); + }); + + // now check if the files don't contain excluded terms + const results = []; + for (const [file, wordList] of fileMap) { + // check if all requirements are matched + + // as search terms with length < 3 are discarded + const filteredTermCount = [...searchTerms].filter( + (term) => term.length > 2 + ).length; + if ( + wordList.length !== searchTerms.size && + wordList.length !== filteredTermCount + ) + continue; + + // ensure that none of the excluded terms is in the search result + if ( + [...excludedTerms].some( + (term) => + terms[term] === file || + titleTerms[term] === file || + (terms[term] || []).includes(file) || + (titleTerms[term] || []).includes(file) + ) + ) + break; + + // select one (max) score for the file. + const score = Math.max(...wordList.map((w) => scoreMap.get(file)[w])); + // add result to the result list + results.push([ + docNames[file], + titles[file], + "", + null, + score, + filenames[file], + ]); } - displayNextItem(); + return results; }, - performObjectSearch : function(object, otherterms) { - var filenames = this._index.filenames; - var objects = this._index.objects; - var objnames = this._index.objnames; - var titles = this._index.titles; - - var importantResults = []; - var objectResults = []; - var unimportantResults = []; - - for (var prefix in objects) { - for (var name in objects[prefix]) { - var fullname = (prefix ? prefix + '.' : '') + name; - if (fullname.toLowerCase().indexOf(object) > -1) { - var match = objects[prefix][name]; - var objname = objnames[match[1]][2]; - var title = titles[match[0]]; - // If more than one term searched for, we require other words to be - // found in the name/title/description - if (otherterms.length > 0) { - var haystack = (prefix + ' ' + name + ' ' + - objname + ' ' + title).toLowerCase(); - var allfound = true; - for (var i = 0; i < otherterms.length; i++) { - if (haystack.indexOf(otherterms[i]) == -1) { - allfound = false; - break; - } - } - if (!allfound) { - continue; - } - } - var descr = objname + _(', in ') + title; - anchor = match[3]; - if (anchor == '') - anchor = fullname; - else if (anchor == '-') - anchor = objnames[match[1]][1] + '-' + fullname; - result = [filenames[match[0]], fullname, '#'+anchor, descr]; - switch (match[2]) { - case 1: objectResults.push(result); break; - case 0: importantResults.push(result); break; - case 2: unimportantResults.push(result); break; - } - } - } - } + /** + * helper function to return a node containing the + * search summary for a given text. keywords is a list + * of stemmed words. + */ + makeSearchSummary: (htmlText, keywords, anchor) => { + const text = Search.htmlToText(htmlText, anchor); + if (text === "") return null; - // sort results descending - objectResults.sort(function(a, b) { - return (a[1] > b[1]) ? -1 : ((a[1] < b[1]) ? 1 : 0); - }); + const textLower = text.toLowerCase(); + const actualStartPosition = [...keywords] + .map((k) => textLower.indexOf(k.toLowerCase())) + .filter((i) => i > -1) + .slice(-1)[0]; + const startWithContext = Math.max(actualStartPosition - 120, 0); - importantResults.sort(function(a, b) { - return (a[1] > b[1]) ? -1 : ((a[1] < b[1]) ? 1 : 0); - }); + const top = startWithContext === 0 ? "" : "..."; + const tail = startWithContext + 240 < text.length ? "..." : ""; - unimportantResults.sort(function(a, b) { - return (a[1] > b[1]) ? -1 : ((a[1] < b[1]) ? 1 : 0); - }); + let summary = document.createElement("p"); + summary.classList.add("context"); + summary.textContent = top + text.substr(startWithContext, 240).trim() + tail; - return [importantResults, objectResults, unimportantResults] - } -} + return summary; + }, +}; -$(document).ready(function() { - Search.init(); -}); \ No newline at end of file +_ready(Search.init); diff --git a/_static/sidebar.js b/_static/sidebar.js new file mode 100644 index 0000000..f28c206 --- /dev/null +++ b/_static/sidebar.js @@ -0,0 +1,70 @@ +/* + * sidebar.js + * ~~~~~~~~~~ + * + * This script makes the Sphinx sidebar collapsible. + * + * .sphinxsidebar contains .sphinxsidebarwrapper. This script adds + * in .sphixsidebar, after .sphinxsidebarwrapper, the #sidebarbutton + * used to collapse and expand the sidebar. + * + * When the sidebar is collapsed the .sphinxsidebarwrapper is hidden + * and the width of the sidebar and the margin-left of the document + * are decreased. When the sidebar is expanded the opposite happens. + * This script saves a per-browser/per-session cookie used to + * remember the position of the sidebar among the pages. + * Once the browser is closed the cookie is deleted and the position + * reset to the default (expanded). + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * + */ + +const initialiseSidebar = () => { + + + + + // global elements used by the functions. + const bodyWrapper = document.getElementsByClassName("bodywrapper")[0] + const sidebar = document.getElementsByClassName("sphinxsidebar")[0] + const sidebarWrapper = document.getElementsByClassName('sphinxsidebarwrapper')[0] + const sidebarButton = document.getElementById("sidebarbutton") + const sidebarArrow = sidebarButton.querySelector('span') + + // for some reason, the document has no sidebar; do not run into errors + if (typeof sidebar === "undefined") return; + + const flipArrow = element => element.innerText = (element.innerText === "»") ? "«" : "»" + + const collapse_sidebar = () => { + bodyWrapper.style.marginLeft = ".8em"; + sidebar.style.width = ".8em" + sidebarWrapper.style.display = "none" + flipArrow(sidebarArrow) + sidebarButton.title = _('Expand sidebar') + window.localStorage.setItem("sidebar", "collapsed") + } + + const expand_sidebar = () => { + bodyWrapper.style.marginLeft = "" + sidebar.style.removeProperty("width") + sidebarWrapper.style.display = "" + flipArrow(sidebarArrow) + sidebarButton.title = _('Collapse sidebar') + window.localStorage.setItem("sidebar", "expanded") + } + + sidebarButton.addEventListener("click", () => { + (sidebarWrapper.style.display === "none") ? expand_sidebar() : collapse_sidebar() + }) + + if (!window.localStorage.getItem("sidebar")) return + const value = window.localStorage.getItem("sidebar") + if (value === "collapsed") collapse_sidebar(); + else if (value === "expanded") expand_sidebar(); +} + +if (document.readyState !== "loading") initialiseSidebar() +else document.addEventListener("DOMContentLoaded", initialiseSidebar) \ No newline at end of file diff --git a/_static/sphinx_highlight.js b/_static/sphinx_highlight.js new file mode 100644 index 0000000..8a96c69 --- /dev/null +++ b/_static/sphinx_highlight.js @@ -0,0 +1,154 @@ +/* Highlighting utilities for Sphinx HTML documentation. */ +"use strict"; + +const SPHINX_HIGHLIGHT_ENABLED = true + +/** + * highlight a given string on a node by wrapping it in + * span elements with the given class name. + */ +const _highlight = (node, addItems, text, className) => { + if (node.nodeType === Node.TEXT_NODE) { + const val = node.nodeValue; + const parent = node.parentNode; + const pos = val.toLowerCase().indexOf(text); + if ( + pos >= 0 && + !parent.classList.contains(className) && + !parent.classList.contains("nohighlight") + ) { + let span; + + const closestNode = parent.closest("body, svg, foreignObject"); + const isInSVG = closestNode && closestNode.matches("svg"); + if (isInSVG) { + span = document.createElementNS("http://www.w3.org/2000/svg", "tspan"); + } else { + span = document.createElement("span"); + span.classList.add(className); + } + + span.appendChild(document.createTextNode(val.substr(pos, text.length))); + const rest = document.createTextNode(val.substr(pos + text.length)); + parent.insertBefore( + span, + parent.insertBefore( + rest, + node.nextSibling + ) + ); + node.nodeValue = val.substr(0, pos); + /* There may be more occurrences of search term in this node. So call this + * function recursively on the remaining fragment. + */ + _highlight(rest, addItems, text, className); + + if (isInSVG) { + const rect = document.createElementNS( + "http://www.w3.org/2000/svg", + "rect" + ); + const bbox = parent.getBBox(); + rect.x.baseVal.value = bbox.x; + rect.y.baseVal.value = bbox.y; + rect.width.baseVal.value = bbox.width; + rect.height.baseVal.value = bbox.height; + rect.setAttribute("class", className); + addItems.push({ parent: parent, target: rect }); + } + } + } else if (node.matches && !node.matches("button, select, textarea")) { + node.childNodes.forEach((el) => _highlight(el, addItems, text, className)); + } +}; +const _highlightText = (thisNode, text, className) => { + let addItems = []; + _highlight(thisNode, addItems, text, className); + addItems.forEach((obj) => + obj.parent.insertAdjacentElement("beforebegin", obj.target) + ); +}; + +/** + * Small JavaScript module for the documentation. + */ +const SphinxHighlight = { + + /** + * highlight the search words provided in localstorage in the text + */ + highlightSearchWords: () => { + if (!SPHINX_HIGHLIGHT_ENABLED) return; // bail if no highlight + + // get and clear terms from localstorage + const url = new URL(window.location); + const highlight = + localStorage.getItem("sphinx_highlight_terms") + || url.searchParams.get("highlight") + || ""; + localStorage.removeItem("sphinx_highlight_terms") + url.searchParams.delete("highlight"); + window.history.replaceState({}, "", url); + + // get individual terms from highlight string + const terms = highlight.toLowerCase().split(/\s+/).filter(x => x); + if (terms.length === 0) return; // nothing to do + + // There should never be more than one element matching "div.body" + const divBody = document.querySelectorAll("div.body"); + const body = divBody.length ? divBody[0] : document.querySelector("body"); + window.setTimeout(() => { + terms.forEach((term) => _highlightText(body, term, "highlighted")); + }, 10); + + const searchBox = document.getElementById("searchbox"); + if (searchBox === null) return; + searchBox.appendChild( + document + .createRange() + .createContextualFragment( + '

    ' + + '' + + _("Hide Search Matches") + + "

    " + ) + ); + }, + + /** + * helper function to hide the search marks again + */ + hideSearchWords: () => { + document + .querySelectorAll("#searchbox .highlight-link") + .forEach((el) => el.remove()); + document + .querySelectorAll("span.highlighted") + .forEach((el) => el.classList.remove("highlighted")); + localStorage.removeItem("sphinx_highlight_terms") + }, + + initEscapeListener: () => { + // only install a listener if it is really needed + if (!DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS) return; + + document.addEventListener("keydown", (event) => { + // bail for input elements + if (BLACKLISTED_KEY_CONTROL_ELEMENTS.has(document.activeElement.tagName)) return; + // bail with special keys + if (event.shiftKey || event.altKey || event.ctrlKey || event.metaKey) return; + if (DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS && (event.key === "Escape")) { + SphinxHighlight.hideSearchWords(); + event.preventDefault(); + } + }); + }, +}; + +_ready(() => { + /* Do not call highlightSearchWords() when we are on the search page. + * It will highlight words from the *previous* search query. + */ + if (typeof Search === "undefined") SphinxHighlight.highlightSearchWords(); + SphinxHighlight.initEscapeListener(); +}); diff --git a/advection_diffusion.html b/advection_diffusion.html index 5b9e363..743360f 100644 --- a/advection_diffusion.html +++ b/advection_diffusion.html @@ -1,306 +1,150 @@ - - - - - + - - The advection-diffusion-reaction equation — FVM Docs 0.1 documentation - - - - + + The advection-diffusion-reaction equation — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Overview

    -

    Next topic

    -

    Boundary conditions for the advection-diffusion-reaction equation

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Overview

    -

    Next topic

    -

    Boundary conditions for the advection-diffusion-reaction equation

    -

    This Page

    - - -
    -
    - + + -
    -
    -
    -
    -
    - -
    -

    The advection-diffusion-reaction equation

    +
    +
    +
    +
    + +
    +

    The advection-diffusion-reaction equation

    The advection-diffusion-reaction equation (also called the continuity equation in semiconductor physics) in flux form, is given by,

    -
    +
    \[u_t = (\mathcal{-F})_x + s(t,x,u)\]
    -

    where \(\mathcal{F} = au - du_x\). This expression holds over the whole domain, therefore it must also hold over a subdomain \(\Omega_j=(x_{j-1/2}, x_{j+1/2})\). Integrating over the subdomain, or cell, gives,

    -
    +

    where \(\mathcal{F} = au - du_x\). This expression holds over the whole domain, therefore it must also hold over a subdomain \(\Omega_j=(x_{j-1/2}, x_{j+1/2})\). Integrating over the subdomain, or cell, gives,

    +
    \[\int_{x_{j-1/2}}^{x_{j+1/2}} u_t~dx = \int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx + \int_{x_{j-1/2}}^{x_{j+1/2}} s(x,t,u)~dx\]
    -

    If we now use \(w\) and \(\bar{s}\) to represent the cell averages of the solution variable and reaction term then we don’t have to perform any actual integrations. The fundamental quantity we will use in the discretised equations is a cell average,

    -
    +

    If we now use \(w\) and \(\bar{s}\) to represent the cell averages of the solution variable and reaction term then we don’t have to perform any actual integrations. The fundamental quantity we will use in the discretised equations is a cell average,

    +
    \[w_j^{\prime} = -\frac{\mathcal{F}_{j+1/2}}{h_j} + \frac{\mathcal{F}_{j-1/2}}{h_{j}} + \bar{s}_j\]
    -

    Note

    -

    In the above we can perform the integration of the flux term. Integrating the divergence of the flux over the unit cell simply gives use the flux at the cell faces, \(\int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx = \left[ \mathcal{-F} \right]_{x_{j-1/2}}^{x_{j+1/2}}\). Because we want the cell average we divide by the cell width \(h_j\). In three dimension we would divide by the volume of the cell, hence the name, finite volume method.

    +

    Note

    +

    In the above we can perform the integration of the flux term. Integrating the divergence of the flux over the unit cell simply gives use the flux at the cell faces, \(\int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx = \left[ \mathcal{-F} \right]_{x_{j-1/2}}^{x_{j+1/2}}\). Because we want the cell average we divide by the cell width \(h_j\). In three dimension we would divide by the volume of the cell, hence the name, finite volume method.

    -

    For the advection component of the flux we will use a linear interpolation (see Aside: Linear interpolation between cell centre and face values) to determine the contribution at the cell faces, and for the diffusion component with will use a simple cell average,

    -
    +

    For the advection component of the flux we will use a linear interpolation (see Aside: Linear interpolation between cell centre and face values) to determine the contribution at the cell faces, and for the diffusion component with will use a simple cell average,

    +
    \[\mathcal{F}_{j+\frac{1}{2}} = a_{j+\frac{1}{2}}\left( \frac{h_{j+1}}{2h_{+}} w_j + \frac{h_j}{2h_{+}} w_{j+1} \right) - d_{j+\frac{1}{2}} \frac{w_{j+1}-w_j}{h_{+}}\]
    -
    +
    \[\mathcal{F}_{j-\frac{1}{2}} = a_{j-\frac{1}{2}}\left( \frac{h_{j}}{2h_{-}} w_{j-1} + \frac{h_{j-1}}{2h_{-}} w_{j} \right) - d_{j-\frac{1}{2}} \frac{w_{j}-w_{j-1}}{h_{-}}\]
    -

    The meaning of the “h” terms is shown in the figure below,

    -
    -Finite volume cells with important distances and positions labelled. -

    Finite volume cells with important distances and positions labelled.

    -
    +

    The meaning of the “h” terms is shown in the figure below,

    +
    +Finite volume cells with important distances and positions labelled. + +
    +

    Finite volume cells with important distances and positions labelled.

    +
    +

    Substituting the definition of the fluxes into the above equation yields,

    -
    -(1)\[w_j^{\prime} = \frac{w_{j-1}}{h_j} \left( \frac{a(x_{j-\frac{1}{2}}) h_j}{2h_{-}} + \frac{d(x_{j-\frac{1}{2}})}{h_{-}}\right) + \frac{w_j}{h_j}\left( \frac{a(x_{j-\frac{1}{2}})h_{j-1}}{2h_{-}} - \frac{a(x_{j+\frac{1}{2}})h_{j+1}}{2h_{+}} - \frac{d(x_{j-\frac{1}{2}})}{h_{-}} - \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right) + \frac{w_{j+1}}{h_j} \left( \frac{-a(x_{j+\frac{1}{2}})h_j}{2h_{+}} + \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right)\]
    +
    +(1)\[w_j^{\prime} = \frac{w_{j-1}}{h_j} \left( \frac{a(x_{j-\frac{1}{2}}) h_j}{2h_{-}} + \frac{d(x_{j-\frac{1}{2}})}{h_{-}}\right) + \frac{w_j}{h_j}\left( \frac{a(x_{j-\frac{1}{2}})h_{j-1}}{2h_{-}} - \frac{a(x_{j+\frac{1}{2}})h_{j+1}}{2h_{+}} - \frac{d(x_{j-\frac{1}{2}})}{h_{-}} - \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right) + \frac{w_{j+1}}{h_j} \left( \frac{-a(x_{j+\frac{1}{2}})h_j}{2h_{+}} + \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right)\]

    This equation is the semi-discrete form because the equation has not been discretised in time.

    -

    For convenience we will simply write the r.h.s. of (1) as,

    -
    +

    For convenience we will simply write the r.h.s. of (1) as,

    +
    \[w_j^{\prime} = F(w)\]
    -

    where it is understood that \(w = (w_{j-1}, w_j, w_{j+1})\) stands for the discretisation stencil.

    -
    -

    Adaptive upwinding & exponential fitting

    +

    where it is understood that \(w = (w_{j-1}, w_j, w_{j+1})\) stands for the discretisation stencil.

    +
    +

    Adaptive upwinding & exponential fitting

    Hundsdorfer (on pg. 266) notes that a adaptive upwinding scheme can be introduced by altering the diffusion coefficient,

    -
    +
    \[d(x_{j}) \rightarrow d(x_{j}) + \frac{1}{2}\kappa_{j} h_{j} a(x_{j})\]
    -
    +
    \[\kappa=\text{sgn}(a(x_{j})).\]
    -

    Here the nature of the discretisation is defined by the sign of the \(a\), when \(a\neq0\) the scheme results in an upwind discretisation. Only if \(a=0\) (corresponding to the diffusion equation limit) does the scheme become based on a central difference (see the red line in the Figure below). A generalisation is possible by choosing a function for \(\kappa\) which automatically adjusts the discretisation method depending on the local value of Peclet’s number, \(\mu=ah/d\). Exponential fitting has shown to be a robust scheme for steady-state approaches. Exponential fitting can be introduced by,

    -
    +

    Here the nature of the discretisation is defined by the sign of the \(a\), when \(a\neq0\) the scheme results in an upwind discretisation. Only if \(a=0\) (corresponding to the diffusion equation limit) does the scheme become based on a central difference (see the red line in the Figure below). A generalisation is possible by choosing a function for \(\kappa\) which automatically adjusts the discretisation method depending on the local value of Peclet’s number, \(\mu=ah/d\). Exponential fitting has shown to be a robust scheme for steady-state approaches. Exponential fitting can be introduced by,

    +
    \[\kappa = \frac{e^{\mu}+1}{e^{\mu}-1} - \frac{2}{\mu}\]

    An approximation to the above, which does not require the evaluation of exponentials is commonly used,

    -
    +
    \[\begin{split}\kappa = \begin{cases} \text{max}(0, 1-2/\mu) & \text{when}~ \mu>0 \\ \text{min}(0, -1-2/\mu) & \text{when}~ \mu<0 \end{cases}\end{split}\]
    -
    -Exponential fitting and approximation -
    -

    This is a very elegant want to introduce adaptive upwinding or exponential fitting because it can be brought into the discretisation in an ad hoc manner. Note that for all of the adaptive expressions \(\kappa\rightarrow\pm1\) as \(\mu\rightarrow\pm\infty\), and \(\kappa\rightarrow 0\) as \(\mu\rightarrow 0\). This means the discretisation will automatically weight in the favour of the upwind scheme in locations where advection dominates. Conversely, where diffusion dominates the weighting factor shifts in favour of a central difference scheme.

    -
    -
    -

    Explicit and implicit forms

    +
    +Exponential fitting and approximation + +
    +

    This is a very elegant want to introduce adaptive upwinding or exponential fitting because it can be brought into the discretisation in an ad hoc manner. Note that for all of the adaptive expressions \(\kappa\rightarrow\pm1\) as \(\mu\rightarrow\pm\infty\), and \(\kappa\rightarrow 0\) as \(\mu\rightarrow 0\). This means the discretisation will automatically weight in the favour of the upwind scheme in locations where advection dominates. Conversely, where diffusion dominates the weighting factor shifts in favour of a central difference scheme.

    +
    +
    +

    Explicit and implicit forms

    In the discussion so far we have been ignoring one important deal: at what time is the r.h.s of the discretised equation evaluated?

    -

    Moreover, if we choose the r.h.s to be evaluated at the present time step, \(t_n\) this is known as an explicit method,

    -
    +

    Moreover, if we choose the r.h.s to be evaluated at the present time step, \(t_n\) this is known as an explicit method,

    +
    \[w_j^{\prime} = F(w^n)\]
    -

    Explicit methods are very simple. Starting with initial conditions at time \(t_n\), the above equation can be rearranged to find the solution variable \(w_j^{n+1}\) at the future time step, \(t_{n+1}\).

    +

    Explicit methods are very simple. Starting with initial conditions at time \(t_n\), the above equation can be rearranged to find the solution variable \(w_j^{n+1}\) at the future time step, \(t_{n+1}\).

    However the downside of using explicit methods is that they are often numerically unstable unless the time step is exceptionally (sometime unrealistically) small.

    -

    Fortunately there is an second alternative, we can choice to write the r.h.s of the equation at the future time step, \(t_{n+1}\), this is known as an implicit method,

    -
    +

    Fortunately there is an second alternative, we can choice to write the r.h.s of the equation at the future time step, \(t_{n+1}\), this is known as an implicit method,

    +
    \[w_j^{\prime} = F(w^{n+1})\]
    -

    Explicit methods are significantly more numerically robust, however they pose a problem, how does one write the equations because the solution variable at the future time step \(w^{n+1}\) is unknown? The answer is that at each time step we must solve a linear system of equation to find \(w^{n+1}\).

    +

    Explicit methods are significantly more numerically robust, however they pose a problem, how does one write the equations because the solution variable at the future time step \(w^{n+1}\) is unknown? The answer is that at each time step we must solve a linear system of equation to find \(w^{n+1}\).

    -

    Note

    -

    For the linear advection-diffusion-reaction equation implicit methods are simply to implement even though the computation cost is increases. One must simply write the equation in the linear form \(A\cdot x = d\) and solve for \(x\) which is the solution variable at the future time step. However if the equations are non-linear then implicit methods pose problem because the equation cannot be written in linear form. In these situations are there are a number of techniques that are used but they all most use a iterative procedure, such as a Newton-Raphson method, to solve the equations.

    +

    Note

    +

    For the linear advection-diffusion-reaction equation implicit methods are simply to implement even though the computation cost is increases. One must simply write the equation in the linear form \(A\cdot x = d\) and solve for \(x\) which is the solution variable at the future time step. However if the equations are non-linear then implicit methods pose problem because the equation cannot be written in linear form. In these situations are there are a number of techniques that are used but they all most use a iterative procedure, such as a Newton-Raphson method, to solve the equations.

    The following section assumes that the equation is linear.

    -
    -
    -

    The \(\theta\)-method

    -

    The \(\theta\)-method is an approach which combines implicit and explicit forms into one method. It consists of writing the equation as the average of the implicit and explicit forms. If we let \(F_{w^{n}}\) and \(F_{w^{n+1}}\) stand for the r.h.s of the explicit and implicit forms of the equations then the \(\theta\)-method gives,

    -
    +
    +
    +

    The \(\theta\)-method

    +

    The \(\theta\)-method is an approach which combines implicit and explicit forms into one method. It consists of writing the equation as the average of the implicit and explicit forms. If we let \(F_{w^{n}}\) and \(F_{w^{n+1}}\) stand for the r.h.s of the explicit and implicit forms of the equations then the \(\theta\)-method gives,

    +
    \[w_j^{\prime} = \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n})\]
    -

    Setting \(\theta=0\) recovers a fully implicit scheme while \(\theta=1\) gives a fully explicit discretisation. The value of \(\theta\) is not limited to just 0 or 1. It is common to set \(\theta=1/2\), this is called the Crank-Nicolson method. It is particularly popular for diffusion problem because it preserves the stability of the implicit form but also increases the accuracy of the time integration from first to second order (because two points in time are being averaged). For advection diffusion problems the Crank-Nicolson method is also unconditionally stable.

    +

    Setting \(\theta=0\) recovers a fully implicit scheme while \(\theta=1\) gives a fully explicit discretisation. The value of \(\theta\) is not limited to just 0 or 1. It is common to set \(\theta=1/2\), this is called the Crank-Nicolson method. It is particularly popular for diffusion problem because it preserves the stability of the implicit form but also increases the accuracy of the time integration from first to second order (because two points in time are being averaged). For advection diffusion problems the Crank-Nicolson method is also unconditionally stable.

    In the above equation,

    -
    +
    \[\begin{split}F(w_j^{n}) = r_a w_{j-1}^{n} + r_b w_{j}^{n} + r_c w_{j+1}^{n} \\ F(w_j^{n+1}) = r_a w_{j-1}^{n+1} + r_b w_{j}^{n+1} + r_c w_{j+1}^{n+1}\end{split}\]

    and the coefficients are,

    -
    +
    \[\begin{split}r_a & = \frac{1}{h_j} \left( \frac{a(x_{j-\frac{1}{2}}) h_j}{2h_{-}} + \frac{d(x_{j-\frac{1}{2}})}{h_{-}}\right) \\ r_b & = \frac{1}{h_j}\left( \frac{a(x_{j-\frac{1}{2}})h_{j-1}}{2h_{-}} - \frac{a(x_{j+\frac{1}{2}})h_{j+1}}{2h_{+}} - \frac{d(x_{j-\frac{1}{2}})}{h_{-}} - \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right)\\ r_c & = \frac{1}{h_j} \left( \frac{-a(x_{j+\frac{1}{2}})h_j}{2h_{+}} + \frac{d(x_{j+\frac{1}{2}})}{h_{+}} \right)\end{split}\]

    We have written the coefficients without dependence on time to simplify the notation, but the coefficients should be calculated at the same time points as their solution variable. However, the coefficients must be linear, they should not depend on the values of the solution variable.

    -
    -
    -

    Discretised equation in matrix form

    +
    +
    +

    Discretised equation in matrix form

    Dropping the spatial subscripts and writing in vector form the equations becomes,

    -
    -\[\begin{split}\frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n})\end{split}\]
    +
    +\[\frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n})\]

    where,

    -
    +
    \[\begin{split}F(w^{n+1}) = M^{n+1} w^{n+1} + s^{n+1} \\ F(w^{n}) = M^{n}w^{n} + s^{n}\end{split}\]

    with the coefficient matrix,

    -
    +
    \[\begin{split}M = \begin{align} \begin{pmatrix} @@ -311,72 +155,110 @@

    Discretised equation in matrix form\(r\)-terms have be previously defined as the coefficients that result from the discretisation method. Note that the subscripts for the \(r\)-terms have been dropped to simplify the notation, but they are functions of space. For example, terms in the first row should be calculated with \(j=1\), in the second with \(j=2\) etc. Solving linear equations is discussed late in section XX.

    -
    -

    Aside: Linear interpolation between cell centre and face values

    -

    In general, linear interpolation between two points \((x_0, x_1)\) can be used to find the value of a function at \(f(x)\),

    -
    +

    The \(r\)-terms have be previously defined as the coefficients that result from the discretisation method. Note that the subscripts for the \(r\)-terms have been dropped to simplify the notation, but they are functions of space. For example, terms in the first row should be calculated with \(j=1\), in the second with \(j=2\) etc. Solving linear equations is discussed late in section XX.

    +
    +

    Aside: Linear interpolation between cell centre and face values

    +

    In general, linear interpolation between two points \((x_0, x_1)\) can be used to find the value of a function at \(f(x)\),

    +
    \[f(x) = \frac{x - x_1}{x_0 - x_1}f(x_0) + \frac{x - x_0}{x_1 - x_0}f(x_1)\]
    -

    In a cell centred grid we know the value of the variable \(w\) at difference points, \(w_j\) and \(w_{j+1}\). We can apply the linear interpolation formulae above to determine value at cell face \(w_{j+1/2}\).

    -
    +

    In a cell centred grid we know the value of the variable \(w\) at difference points, \(w_j\) and \(w_{j+1}\). We can apply the linear interpolation formulae above to determine value at cell face \(w_{j+1/2}\).

    +
    \[w_{j+1/2} = \frac{x_{j+1/2} - x_{j+1}}{x_{j} - x_{j+1}} w_j + \frac{x_{j+1/2} - x_j}{x_{j+1} - x_j} w_{j+1}\]

    This can be simplified firstly by using function to represent the distance between cell centres,

    -
    +
    \[h_{-} = x_j - x_{j-1} \quad h_{+} = x_{j+1} - x_{j}\]

    to give,

    -
    +
    \[w_{j+1/2} = \frac{x_{j+1} - x_{j+1/2}}{h_{+}} w_j + \frac{x_{j+1/2} - x_j}{h_{+}} w_{j+1}\]
    -

    This expression still contains \(x_{j+1/2}\) which we can simplify further by using an expression for the position of cell centres,

    -
    +

    This expression still contains \(x_{j+1/2}\) which we can simplify further by using an expression for the position of cell centres,

    +
    \[x_j = \frac{1}{2} \left( x_{j-\frac{1}{2}} + x_{j+\frac{1}{2}} \right) \quad x_{j+1} = \frac{1}{2} \left( x_{j+\frac{1}{2}} + x_{j+\frac{3}{2}} \right)\]
    -

    Note, this expression is still valid of non-uniform grids, it simply says that cell centres are always equidistant from two faces. Rearranging the above expression and substituting in for \(x_{j}\) and \(x_{j+1}\) terms gives,

    -
    +

    Note, this expression is still valid of non-uniform grids, it simply says that cell centres are always equidistant from two faces. Rearranging the above expression and substituting in for \(x_{j}\) and \(x_{j+1}\) terms gives,

    +
    \[w_{j+1/2} = \frac{\frac{1}{2} \left( x_{j+\frac{1}{2}} + x_{j+\frac{3}{2}} \right) - x_{j+1/2}}{h_{+}} w_j + \frac{x_{j+1/2} - \frac{1}{2} \left( x_{j-\frac{1}{2}} + x_{j+\frac{1}{2}} \right)}{h_{+}} w_{j+1}\]
    -

    Finally, by defining the distance between vertices as, \(h_j = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}\), we can simplify to the following expression,

    -
    +

    Finally, by defining the distance between vertices as, \(h_j = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}\), we can simplify to the following expression,

    +
    \[w_{j+1/2} = \frac{h_{j+1}}{2h_{+}} w_j + \frac{h_j}{2h_{+}} w_{j+1}\]
    -

    Similarly the \(w_{j-1/2}\) can be found,

    -
    +

    Similarly the \(w_{j-1/2}\) can be found,

    +
    \[w_{j-1/2} = \frac{h_{j}}{2h_{-}} w_{j-1} + \frac{h_{j-1}}{2h_{-}} w_{j}\]
    -
    -
    -
    +
    +

    +
    -
    -
    -
    +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/advection_diffusion_boundary_conditions.html b/advection_diffusion_boundary_conditions.html index cbeb00e..8c2e704 100644 --- a/advection_diffusion_boundary_conditions.html +++ b/advection_diffusion_boundary_conditions.html @@ -1,274 +1,126 @@ - - - - - + - - Boundary conditions for the advection-diffusion-reaction equation — FVM Docs 0.1 documentation - - - - + + Boundary conditions for the advection-diffusion-reaction equation — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - + + - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    The advection-diffusion-reaction equation

    -

    Next topic

    -

    Solving linear equations

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    The advection-diffusion-reaction equation

    -

    Next topic

    -

    Solving linear equations

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    Boundary conditions for the advection-diffusion-reaction equation

    +
    +
    +
    +
    + +
    +

    Boundary conditions for the advection-diffusion-reaction equation

    The advection-diffusion-reaction equation is a particularly good equation to explore apply boundary conditions because it is a more general version of other equations. For example, the diffusion equation, the transport equation and the Poisson equation can all be recovered from this basic form. Moreover, by developing a general scheme for boundary conditions of the advection-reaction-diffusion equation we automatically get a system for imposing boundary conditions on all equations of a similar form.

    -
    -

    Robin boundary conditions (known flux)

    -

    Robin specify a known total flux comprised of a diffusion and advection component. Note that in the diffusion equation limit (where \(a=0\)) these boundary conditions reduce to Neumann boundary conditions.

    -

    Within the finite volume method Robin boundary conditions are naturally resolved. This means that there is no need for interpolation or ghost point substitution (although these approaches remain possible) to include the boundary conditions because the flux at the boundary appears naturally in the semi-discretised equation. For example consider the semi-discretised equation evaluated at the boundary cell \(\Omega_1\),

    -
    -Finite volume boundary cell at the left hand side. -

    Finite volume boundary cell at the left hand side.

    -
    -
    +
    +

    Robin boundary conditions (known flux)

    +

    Robin specify a known total flux comprised of a diffusion and advection component. Note that in the diffusion equation limit (where \(a=0\)) these boundary conditions reduce to Neumann boundary conditions.

    +

    Within the finite volume method Robin boundary conditions are naturally resolved. This means that there is no need for interpolation or ghost point substitution (although these approaches remain possible) to include the boundary conditions because the flux at the boundary appears naturally in the semi-discretised equation. For example consider the semi-discretised equation evaluated at the boundary cell \(\Omega_1\),

    +
    +Finite volume boundary cell at the left hand side. + +
    +

    Finite volume boundary cell at the left hand side.

    +
    +
    +
    \[w_1^{\prime} = -\frac{\mathcal{F}_{3/2}}{h_1} + \frac{\mathcal{F}_{1/2}}{h_{1}} + \bar{s}_1\]
    -

    The Robin boundary condition specifies the flux at \(x_{1/2}\),

    -
    +

    The Robin boundary condition specifies the flux at \(x_{1/2}\),

    +
    \[\mathcal{F}_{1/2} = g_{R}(x_{1/2})\]

    Therefore the boundary condition can be incorporated without invoking any information regarding the ghost cell,

    -
    +
    \[w_1^{\prime} = \frac{w_1}{h_1}\left( \frac{-a(x_{3/2})h_{2}}{2h_{+}} - \frac{d(x_{3/2})}{h_{+}} \right) + \frac{w_{2}}{h_1} \left( \frac{-a(x_{3/2}) h_1}{2h_{+}} + \frac{d(x_{3/2})}{h_{+}} \right) + \frac{g_{R}(x_L)}{h_1} + \bar{s}_1\]
    -

    Similarly applying the same procedure to the \(\Omega_J\) cell at the right hand side boundary,

    -
    -Finite volume boundary cell at the right hand side. -

    Finite volume boundary cell at the right hand side.

    -
    -
    +

    Similarly applying the same procedure to the \(\Omega_J\) cell at the right hand side boundary,

    +
    +Finite volume boundary cell at the right hand side. + +
    +

    Finite volume boundary cell at the right hand side.

    +
    +
    +
    \[w_J^{\prime} = -\frac{\mathcal{F}_{J+1/2}}{h_J} + \frac{\mathcal{F}_{J-1/2}}{h_J} + \bar{s}_J\]

    The Robin boundary condition at the right hand side is,

    -
    +
    \[g_{R}(x_R) = \mathcal{F}_{J+1/2}\]

    Therefore the naturally resolved boundary condition on the right hand side becomes,

    -
    +
    \[w_J^{\prime} = \frac{w_{J-1}}{h_J}\left( \frac{a(x_{J-1/2})h_{J}}{2h_{-}} + \frac{d(x_{J-1/2})}{h_{-}} \right) + \frac{w_{J}}{h_J} \left( \frac{a(x_{J-1/2}) h_{J-1}}{2h_{-}} - \frac{d(x_{J-1/2})}{h_{-}} \right) - \frac{g_{R}(x_R)}{h_J} + \bar{s}_J\]
    -
    -
    -

    Dirichlet boundary conditions (known value)

    +
    +
    +

    Dirichlet boundary conditions (known value)

    -

    Warning

    -

    This is only valid provided the initial conditions also satisfy the boundary conditions

    +

    Warning

    +

    This is only valid provided the initial conditions also satisfy the boundary conditions

    Dirichlet (known value) boundary conditions are trivial to implement. All that is required is to fix the boundary term to be a constant value. If the initial conditions satisfy the boundary conditions (as they should) then all that is needed to to make sure that the time-derivative (the change with time) of the boundary points is always zero. For example, to impose Dirichlet boundary conditions

    -
    +
    \[w_1^{n+1} = g_D(x_1) \quad w_{J}^{n+1} = g_D(x_J)\]

    We simply need to specify,

    -
    +
    \[F_1(w_1^{n}) = 0 \quad F_1(w_1^{n+1}) = 0\]

    This can be achieved by multiplying the r.h.s. by the modified identity matrix where the first and last elements are zero,

    -
    +
    \[w^{\prime} = \text{diag}(0,1,\cdots,1,0) \left( \theta F(w^{n+1}) + (1-\theta)F(w^{n}) \right)\]

    If the boundary conditions are time dependent,

    -
    +
    \[w_1^{n+1} = g_D(x_1, t_{n+1}) \quad w_{J}^{n+1} = g_D(x_J, t_{n+1})\]

    The we should also add a matrix to the r.h.s. which is contains the changes to the boundary value between the time-steps,

    -
    +
    \[w^{\prime} = \text{diag}(0,1,\cdots,1,0) \left( \theta F(w^{n+1}) + (1-\theta)F(w^{n}) \right) + \text{diag}(\Delta g_D(x_1), 0, \cdots, 0, \Delta g_D(x_J) )\]

    where,

    -
    +
    \[\begin{split}\Delta g_D(x_1) = g_D(x_1, t_{n+1}) - g_D(x_1, t_{n}) \\ \Delta g_D(x_J) = g_D(x_J, t_{n+1}) - g_D(x_J, t_{n})\end{split}\]

    The following assumes time-invariant boundary conditions.

    -
    -
    -

    General matrix form with a transient term

    -

    The semi-discretised advection-diffusion-reaction equation can be written in the form below using the \(\theta\)-method,

    -
    +
    +
    +

    General matrix form with a transient term

    +

    The semi-discretised advection-diffusion-reaction equation can be written in the form below using the \(\theta\)-method,

    +
    \[w^{\prime} = \alpha\left( \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n}) \right) + \beta\]
    -

    where \(F(w)\) contains a matrix \(M\), a vector of the solution variable \(w\) and a vector of the reaction term \(r\),

    -
    +

    where \(F(w)\) contains a matrix \(M\), a vector of the solution variable \(w\) and a vector of the reaction term \(r\),

    +
    \[F(w) = Mw + r\]
    -

    We have introduced a new matrix, \(\alpha\) and vector, \(\beta\). These are used to incorporate the boundary conditions, they have the form,

    -
    +

    We have introduced a new matrix, \(\alpha\) and vector, \(\beta\). These are used to incorporate the boundary conditions, they have the form,

    +
    \[\begin{split}\alpha & = \text{diag}\left( \alpha_1, 1, \cdots, 1, \alpha_J \right) \\ \beta & = \left( \beta_1, 0, \cdots, 0, \beta_J \right)\end{split}\]

    The modified coefficient matrix becomes,

    -
    +
    \[\begin{split}M = \begin{align} \begin{pmatrix} @@ -280,130 +132,155 @@

    General matrix form with a transient term - - - - - - -Symbol -Robin -Dirichlet + + + + + - - - - + + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + +

    Symbol

    Robin

    Dirichlet
    \(b_1\)\(\frac{1}{h_1}\left( \frac{-a(x_{3/2})h_{2}}{2h_{+}} - \frac{d(x_{3/2})}{h_{+}} \right)\)\(1\)

    \(b_1\)

    \(\frac{1}{h_1}\left( \frac{-a(x_{3/2})h_{2}}{2h_{+}} - \frac{d(x_{3/2})}{h_{+}} \right)\)

    \(1\)
    \(c_1\)\(\frac{1}{h_1} \left( \frac{-a(x_{3/2}) h_1}{2h_{+}} + \frac{d(x_{3/2})}{h_{+}} \right)\)\(0\)

    \(c_1\)

    \(\frac{1}{h_1} \left( \frac{-a(x_{3/2}) h_1}{2h_{+}} + \frac{d(x_{3/2})}{h_{+}} \right)\)

    \(0\)
    \(a_J\)\(\frac{1}{h_J}\left( \frac{a(x_{J-1/2})h_{J}}{2h_{-}} + \frac{d(x_{J-1/2})}{h_{-}} \right)\)\(0\)

    \(a_J\)

    \(\frac{1}{h_J}\left( \frac{a(x_{J-1/2})h_{J}}{2h_{-}} + \frac{d(x_{J-1/2})}{h_{-}} \right)\)

    \(0\)
    \(b_J\)\(\frac{1}{h_J} \left( \frac{a(x_{J-1/2}) h_{J-1}}{2h_{-}} - \frac{d(x_{J-1/2})}{h_{-}} \right)\)\(1\)

    \(b_J\)

    \(\frac{1}{h_J} \left( \frac{a(x_{J-1/2}) h_{J-1}}{2h_{-}} - \frac{d(x_{J-1/2})}{h_{-}} \right)\)

    \(1\)
    \(\alpha_1\)\(1\)\(0\)

    \(\alpha_1\)

    \(1\)

    \(0\)
    \(\alpha_J\)\(1\)\(0\)

    \(\alpha_J\)

    \(1\)

    \(0\)
    \(\beta_1\)\(\frac{g_R(x_1)}{h_1}\)\(0\)

    \(\beta_1\)

    \(\frac{g_R(x_1)}{h_1}\)

    \(0\)
    \(\beta_J\)\(-\frac{g_R(x_J)}{h_J}\)\(0\)

    \(\beta_J\)

    \(-\frac{g_R(x_J)}{h_J}\)

    \(0\)
    -

    -
    -

    General matrix form without a transient term

    -

    The boundary conditions scheme discussed above is only valid for initial value problems; problems where an initial vector of the solution variable \(w^0\) is specified along with boundary conditions. PDEs of the advection-diffusion-reaction form that do not contain a time derivative are an important class of equations they are called boundary value problems,

    -
    +
    +
    +

    General matrix form without a transient term

    +

    The boundary conditions scheme discussed above is only valid for initial value problems; problems where an initial vector of the solution variable \(w^0\) is specified along with boundary conditions. PDEs of the advection-diffusion-reaction form that do not contain a time derivative are an important class of equations they are called boundary value problems,

    +
    \[0 = \alpha\left( \theta F(w_j^{n+1}) + (1-\theta)F(w_j^{n}) \right) + \beta\]

    Boundary value problems do not require initial conditions for the solution variable. For this type of equations we need to use a different general scheme to implement boundary condition values.

    -

    The element of the \(M\) matrix remain the same, however for Dirichlet boundary conditions the elements must be modified in for the \(\alpha\) and \(\beta\) terms,

    - ----- - - - - +

    The element of the \(M\) matrix remain the same, however for Dirichlet boundary conditions the elements must be modified in for the \(\alpha\) and \(\beta\) terms,

    +
    SymbolRobinDirichlet
    + + + + - - - - + + + + - - - + + + - - - + + + - - - + + +

    Symbol

    Robin

    Dirichlet
    \(\alpha_1\)\(1\)\(1\)

    \(\alpha_1\)

    \(1\)

    \(1\)
    \(\alpha_J\)\(1\)\(1\)

    \(\alpha_J\)

    \(1\)

    \(1\)
    \(\beta_1\)\(\frac{g_R(x_1)}{h_1}\)\(-r(x_1)- g_D(x_1)\)

    \(\beta_1\)

    \(\frac{g_R(x_1)}{h_1}\)

    \(-r(x_1)- g_D(x_1)\)
    \(\beta_J\)\(-\frac{g_R(x_J)}{h_J}\)\(-r(x_J)- g_D(x_J)\)

    \(\beta_J\)

    \(-\frac{g_R(x_J)}{h_J}\)

    \(-r(x_J)- g_D(x_J)\)
    -
    -
    +
    + -
    -
    -
    +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/doctrees/advection_diffusion.doctree b/doctrees/advection_diffusion.doctree index 353d1f6..5ab7050 100644 Binary files a/doctrees/advection_diffusion.doctree and b/doctrees/advection_diffusion.doctree differ diff --git a/doctrees/advection_diffusion_boundary_conditions.doctree b/doctrees/advection_diffusion_boundary_conditions.doctree index b4a08dc..1662543 100644 Binary files a/doctrees/advection_diffusion_boundary_conditions.doctree and b/doctrees/advection_diffusion_boundary_conditions.doctree differ diff --git a/doctrees/environment.pickle b/doctrees/environment.pickle index 32bd07b..36367bf 100644 Binary files a/doctrees/environment.pickle and b/doctrees/environment.pickle differ diff --git a/doctrees/examples.doctree b/doctrees/examples.doctree index 7287fb4..ad54c3c 100644 Binary files a/doctrees/examples.doctree and b/doctrees/examples.doctree differ diff --git a/doctrees/implementation_advection_diffusion.doctree b/doctrees/implementation_advection_diffusion.doctree index 893daa2..370f73c 100644 Binary files a/doctrees/implementation_advection_diffusion.doctree and b/doctrees/implementation_advection_diffusion.doctree differ diff --git a/doctrees/index.doctree b/doctrees/index.doctree index ca71818..0d1b738 100644 Binary files a/doctrees/index.doctree and b/doctrees/index.doctree differ diff --git a/doctrees/overview.doctree b/doctrees/overview.doctree index 76ca403..e30225e 100644 Binary files a/doctrees/overview.doctree and b/doctrees/overview.doctree differ diff --git a/doctrees/poisson.doctree b/doctrees/poisson.doctree index 12d8d02..c66b0c7 100644 Binary files a/doctrees/poisson.doctree and b/doctrees/poisson.doctree differ diff --git a/doctrees/solving_linear_equations.doctree b/doctrees/solving_linear_equations.doctree index d5c7f54..6bc2586 100644 Binary files a/doctrees/solving_linear_equations.doctree and b/doctrees/solving_linear_equations.doctree differ diff --git a/doctrees/solving_nonlinear_equations.doctree b/doctrees/solving_nonlinear_equations.doctree index 22a6be5..a145641 100644 Binary files a/doctrees/solving_nonlinear_equations.doctree and b/doctrees/solving_nonlinear_equations.doctree differ diff --git a/examples.html b/examples.html index ed9da87..c54d062 100644 --- a/examples.html +++ b/examples.html @@ -1,248 +1,129 @@ - - - - - + - - Example calculations — FVM Docs 0.1 documentation - - - - + + Example calculations — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving the advection-diffusion-reaction equation in Python

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving the advection-diffusion-reaction equation in Python

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    Example calculations

    + + + +
    +
    +
    +
    + +
    +

    Example calculations

    -

    Note

    -

    For all simulations \(a=1\), \(d=0.001\) on a domain where \(0\leq x \leq 1\), with mesh spacing \(h=0.02\) and time step \(\tau=0.01\). The equations are solver with the boundary conditions \(u(0,t)=1\) and \(u(1,t)=0\) and initial conditions \(u(x,0)=sin(\pi x)^{100}\). With these boundary conditions a boundary layer will form near \(x=1\).

    +

    Note

    +

    For all simulations \(a=1\), \(d=0.001\) on a domain where \(0\leq x \leq 1\), with mesh spacing \(h=0.02\) and time step \(\tau=0.01\). The equations are solver with the boundary conditions \(u(0,t)=1\) and \(u(1,t)=0\) and initial conditions \(u(x,0)=sin(\pi x)^{100}\). With these boundary conditions a boundary layer will form near \(x=1\).

    -
    -

    Uniform grid

    -

    First we test the finite-volume method using a standard uniform grid. Note that the Peclet number for the above parameters is \(\mu=20\) so the central discretisation scheme is not stable as illustrated by the oscillations in the solution. The upwind scheme does not have a stability criteria related to the Peclet number so the solutions for the upwind case are smooth. Finally, the exponentially fitted scheme has automatically weighted in favour of the upwind discretisation, the value of \(\kappa\approx 0.9\).

    +
    +

    Uniform grid

    +

    First we test the finite-volume method using a standard uniform grid. Note that the Peclet number for the above parameters is \(\mu=20\) so the central discretisation scheme is not stable as illustrated by the oscillations in the solution. The upwind scheme does not have a stability criteria related to the Peclet number so the solutions for the upwind case are smooth. Finally, the exponentially fitted scheme has automatically weighted in favour of the upwind discretisation, the value of \(\kappa\approx 0.9\).

    -
    -
    -

    Random grid

    -

    Although a random grid is of no practical use it is a good test of the code because bugs are more likely to show up when symmetry has been reduced. I the follow simulations \(\text{min}(\kappa)\) =0.004 and \(\text{max}(\kappa)\) =0.98 so the discretisation scheme is abruptly changing from cell to cell.

    +
    +
    +

    Random grid

    +

    Although a random grid is of no practical use it is a good test of the code because bugs are more likely to show up when symmetry has been reduced. I the follow simulations \(\text{min}(\kappa)\) =0.004 and \(\text{max}(\kappa)\) =0.98 so the discretisation scheme is abruptly changing from cell to cell.

    -
    -
    -

    Nonuniform grid

    -

    Nonuniform grids can be used to reduce to improve the solution as shown here. The following simulation contains the same number of cells as the previous simulations however the cell centres are clustered towards the right hand boundary. The increased density of cells allow the boundary layer to be resolved. The transient solution computed from the central scheme is still significantly affected by the high Peclet number but it is interesting to observe that the steady-state solution of all three methods are very similar. Furthermore, \(\kappa\approx0.01-0.1\) in the region of the boundary layer which implies the local value of Peclet number as been reduced enough so allow the exponentially fitted scheme to weight in favour of the higher accuracy central discretisation.

    +
    +
    +

    Nonuniform grid

    +

    Nonuniform grids can be used to reduce to improve the solution as shown here. The following simulation contains the same number of cells as the previous simulations however the cell centres are clustered towards the right hand boundary. The increased density of cells allow the boundary layer to be resolved. The transient solution computed from the central scheme is still significantly affected by the high Peclet number but it is interesting to observe that the steady-state solution of all three methods are very similar. Furthermore, \(\kappa\approx0.01-0.1\) in the region of the boundary layer which implies the local value of Peclet number as been reduced enough so allow the exponentially fitted scheme to weight in favour of the higher accuracy central discretisation.

    -
    -
    +
    + -
    -
    -
    +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/genindex.html b/genindex.html index 0f855c7..15f3790 100644 --- a/genindex.html +++ b/genindex.html @@ -1,156 +1,36 @@ - - - - - - - + + + Index — FVM Docs 0.1 documentation + + - Index — FVM Docs 0.1 documentation - - - - - - - - - - - - - - - - - - - - - - - - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    + + + - - -
    - - - -
    - -
    - Open Table Of Contents -
    - - - - -
    -
    -
    - - - -
    - -
    -
    - - - - -
    -
    - - -
    -
    -
    -
    -
    - + + + + + +
    +
    +
    +
    +

    Index

    @@ -159,38 +39,39 @@

    Index

    -
    -
    -
    +
    - - -
    - - -
    -
    -
    - -
    -
    -
    - - -
    - © Copyright 2013, Daniel J Farrell. - Created using Sphinx 1.1.3. -
    - - +
    + +
    +
    + + - - \ No newline at end of file diff --git a/implementation_advection_diffusion.html b/implementation_advection_diffusion.html index 4f474dd..242196d 100644 --- a/implementation_advection_diffusion.html +++ b/implementation_advection_diffusion.html @@ -1,283 +1,133 @@ - - - - - + - - Solving the advection-diffusion-reaction equation in Python — FVM Docs 0.1 documentation - - - - + + Solving the advection-diffusion-reaction equation in Python — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    The Poisson equation

    -

    Next topic

    -

    Example calculations

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    The Poisson equation

    -

    Next topic

    -

    Example calculations

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    Solving the advection-diffusion-reaction equation in Python

    + + + +
    +
    +
    +
    + +
    +

    Solving the advection-diffusion-reaction equation in Python

    Here we discuss how to implement a solver for the advection-diffusion equation in Python. The notes will consider how to design a solver which minimises code complexity and maximise readability.

    -
    -

    Overview

    +
    +

    Overview

    At a high-level usage of the code looks like the following,

    -
    # Define a mesh
+
    # Define a mesh
 faces = np.linspace(-0.5, 1, 100)
 mesh =  Mesh(faces)
 
-# Define coefficients
-a = CellVariable(0.01, mesh=mesh) # Advection velocity
-d = CellVariable(1e-3, mesh=mesh) # Diffusion coefficient
+# Define coefficients
+a = CellVariable(0.01, mesh=mesh) # Advection velocity
+d = CellVariable(1e-3, mesh=mesh) # Diffusion coefficient
 
-# Make a 'model' and apply boundary conditions
-k = 1 # Time step
+# Make a 'model' and apply boundary conditions
+k = 1 # Time step
 model = Model(faces, a, d, k)
 model.set_boundary_conditions(left_value=1., right_value=0.)
 
-# Ask for the discretised system...
+# Ask for the discretised system...
 M = model.coefficient_matrix()
 alpha = model.alpha_matrix()
 beta = model.beta_vector()
 
-# Solve...
+# Solve...
 
    
 
    
 
    There are a number of classes (highlighted above) which abstract away details of dealing with finite-volume equations:
    
 
    
 
      
-
    1. Mesh
      2. 
-
    2. CellVariable
      3. 
-
    3. Model
      4. 
+
    4. Mesh
      5. 
+
    5. CellVariable
      6. 
+
    6. Model
      7. 
 
    
 
    
-
    The the Mesh and CellVariable classes have been inspired by the approach of Fipy.
    
-
    In the following we will highlight the main features of each class and point out how they are useful. The classes do not attempt at doing “too much”, they simply aid in the creation of the linear system.
    
-
    -
    -

    The Mesh class

    -

    Mesh objects are initialised with a list of faces locations, which can be non-uniformly distributed if desired. A mesh is completely defined by locations of cell faces. Some useful methods,

    -
    def h(self, i):
+
    The the Mesh and CellVariable classes have been inspired by the approach of Fipy.
    
+
    In the following we will highlight the main features of each class and point out how they are useful. The classes do not attempt at doing “too much”, they simply aid in the creation of the linear system.
    
+
    +
    +

    The Mesh class

    +

    Mesh objects are initialised with a list of faces locations, which can be non-uniformly distributed if desired. A mesh is completely defined by locations of cell faces. Some useful methods,

    +
    def h(self, i):
         ...
    -

    h() returns the cell width for cell with index i. This is function is vectorisable by passing an array of the required indices but it does not accept “fancy indexing”. The reason being, that it is very easy to make mistakes with subscript indexing, the goal here is to make the user be explicit when requesting elements. Note the self.cell_widths instance variable returns the numpy array of cell widths is a second way of accessing this data.

    -
    def hm(self, i):
-...
    -
    -

    hm() returns the distance between cell centroids for the cell at index i and i-1, that is in the backwards or minus direction.

    -
    def hp(self, i):
-...
    +

    h() returns the cell width for cell with index i. This is function is vectorisable by passing an array of the required indices but it does not accept “fancy indexing”. The reason being, that it is very easy to make mistakes with subscript indexing, the goal here is to make the user be explicit when requesting elements. Note the self.cell_widths instance variable returns the numpy array of cell widths is a second way of accessing this data.

    +
    def hm(self, i):
+...
+
    -

    Similarly hp() returns the distance between cell centroids for the cell at index i and i+1, that is in the forwards or plus direction.

    -

    In addition to the above method, the class contains the instance variables, self.cells which returns an array of cell centroid locations, self.J which contains the number of cells, and also a copy of the face locations (an array) via self.faces.

    +

    hm() returns the distance between cell centroids for the cell at index i and i-1, that is in the backwards or minus direction.

    +
    def hp(self, i):
+...
+
    -
    -

    The CellVariable class

    -

    The goal of the CellVariable class is to provide a elegant way of automatically interpolating between the cell value and the face value. The class holes values which correspond to the cell average. Internally, this class is a subclass of numpy.ndarray so it is a fully functioning numpy array. It has a new constructor and additional method which return interpolated values at the cell surfaces.

    -

    A CellVariable is initialised with a value for the cell average (this can be a constant or an array-like quantity) and the Mesh on which the cell variable is defined. My coupling the cell variable with the mesh the class can perform interpolation between the cell and face values using the methods,

    -
    def p(self, i):
-...
+
    Similarly hp() returns the distance between cell centroids for the cell at index i and i+1, that is in the forwards or plus direction.
    
+
    In addition to the above method, the class contains the instance variables, self.cells which returns an array of cell centroid locations, self.J which contains the number of cells, and also a copy of the face locations (an array) via self.faces.
    
+
    +
    +

    The CellVariable class

    +

    The goal of the CellVariable class is to provide a elegant way of automatically interpolating between the cell value and the face value. The class holes values which correspond to the cell average. Internally, this class is a subclass of numpy.ndarray so it is a fully functioning numpy array. It has a new constructor and additional method which return interpolated values at the cell surfaces.

    +

    A CellVariable is initialised with a value for the cell average (this can be a constant or an array-like quantity) and the Mesh on which the cell variable is defined. My coupling the cell variable with the mesh the class can perform interpolation between the cell and face values using the methods,

    +
    def p(self, i):
+...
 
-def m(self, i)
-...
    -
    -

    Again self.p(i) stands for the plus direction and self.m(i) stands for the minus direction, as such they return values at the right and left face of the cell. The mesh can be returned via the instance variable cell_variable.mesh.

    +def m(self, i) +... +
    -
    -

    The Model class

    +

    Again self.p(i) stands for the plus direction and self.m(i) stands for the minus direction, as such they return values at the right and left face of the cell. The mesh can be returned via the instance variable cell_variable.mesh.

    +
    +
    +

    The Model class

    The model class is where the creating of the matrices occurs and where boundary conditions can be applied to the problem. For these reasons the class is fairly complicated.

    -

    There are method which return different element of the final matrix. The interior elements are fairly homogenous, the only real difference is where there are spatially varying coefficient of cell widths. For this reason the the method _interior_matrix_elements() returns elements which correspond to the lower, central and upper diagonals for a specific index. For example, to calculate the interior matrix elements for mesh point at value index one would do the following,

    -
    # Return the the interior matrix elements (the r-terms)
-# for a particular spatial index
-    model = Model(...)
-index = ...
-# The lower, central and, upper diagonal terms of the stencil
-    ra, rb, rc = model._interior_matrix_elements(index)
    +

    There are method which return different element of the final matrix. The interior elements are fairly homogenous, the only real difference is where there are spatially varying coefficient of cell widths. For this reason the the method _interior_matrix_elements() returns elements which correspond to the lower, central and upper diagonals for a specific index. For example, to calculate the interior matrix elements for mesh point at value index one would do the following,

    +
    # Return the the interior matrix elements (the r-terms)
+# for a particular spatial index
+    model = Model(...)
+index = ...
+# The lower, central and, upper diagonal terms of the stencil
+    ra, rb, rc = model._interior_matrix_elements(index)
+

    The function names here correspond to the matrix element in the previous section.

    -

    Note that the function is prefixed with an underscore this is because are private, ‘users’ should have no need to call these method. It is called internally when constructing the finite-volume matrices. However, as this is a overview of how to implement this is an readable and useful way we include this detail.

    +

    Note that the function is prefixed with an underscore this is because are private, ‘users’ should have no need to call these method. It is called internally when constructing the finite-volume matrices. However, as this is a overview of how to implement this is an readable and useful way we include this detail.

    The following methods play a similar role,

    -
    def _robin_boundary_condition_matrix_elements_left(self):
+
    def _robin_boundary_condition_matrix_elements_left(self):
         ...
 
 def _robin_boundary_condition_matrix_elements_right(self):
@@ -290,9 +140,9 @@ 
    The Model cl
         ...
 
    
 
    
-
    They return a list of index-value pairs ([(1,1), a11], [(i,2), b12] ...). The functions return the value of element which need to change (with respect to the interior values) in order include boundary conditions. The index-value pair facilitates automatic insertion of the values into the correct matrix element. As we will see later, rather than hard coding the position of the various element if the index and value are specified it makes the destination of the element unambiguous. It also allows the value of the matrix element to be defined at the same point in the code as the location. This is beneficial for providing context and should reduce bugs and complexity.
    
-
    The elements for the \(\beta\) boundary condition vector (which is added to the linear system) are generated from the functions below,
    
-
    def _robin_boundary_condition_vector_elements_left(self):
+
    They return a list of index-value pairs ([(1,1), a11], [(i,2), b12] ...). The functions return the value of element which need to change (with respect to the interior values) in order include boundary conditions. The index-value pair facilitates automatic insertion of the values into the correct matrix element. As we will see later, rather than hard coding the position of the various element if the index and value are specified it makes the destination of the element unambiguous. It also allows the value of the matrix element to be defined at the same point in the code as the location. This is beneficial for providing context and should reduce bugs and complexity.
    
+
    The elements for the \(\beta\) boundary condition vector (which is added to the linear system) are generated from the functions below,
    
+
    def _robin_boundary_condition_vector_elements_left(self):
     ...
 
 def _robin_boundary_condition_vector_elements_right(self):
@@ -306,8 +156,8 @@ 
    The Model cl
 
    
 
    
 
    Again, these method should return index-values pairs.
    
-
    The Model class also include some convenience function for checking the value of the Peclet number and the CFL conditions which can be called via,
    
-
    def peclet_number(self):
+
    The Model class also include some convenience function for checking the value of the Peclet number and the CFL conditions which can be called via,
    
+
    def peclet_number(self):
         return self.a * self.mesh.cell_widths / self.d
 
 def CFL_condition(self):
@@ -315,7 +165,7 @@ 
    The Model cl
 
    
 
    
 
    The method which are intended for the user to actually call when constructing the linear system are,
    
-
    def coefficient_matrix(self):
+
    def coefficient_matrix(self):
     ...
 
 def alpha_matrix(self):
@@ -326,67 +176,102 @@ 
    The Model cl
 
    
 
    
 
    The linear system for time-stepping can be constructed easily,
    
-
    # In pseudo-code
+
    # In pseudo-code
 model = AdvectionDiffusionModel(...)
 M = model.coefficient_matrix()
 alpha = model.alpha_matrix()
 beta = model.beta_vector()
 I = sparse.identity(model.mesh.J)
 
-# Apply boundary conditions
-u_init = ... #
+# Apply boundary conditions
+u_init = ... #
 ...
 
-tau = 0.01   # time step
-theta = 0.5  # Implicit/explicit parameter
+tau = 0.01   # time step
+theta = 0.5  # Implicit/explicit parameter
 
 u = u_init
 for i in range(...):
-    # time step the linear system, A.x = d
+    # time step the linear system, A.x = d
     A = I - tau * theta * alpha * M
     d = (I + tau * (1-theta) * alpha * M) * u
     u = spsolve(A,d)
 
    
 
    
-
    Finally, when initialising a Model object the keyword argument , discretisation is important. Is can be set to one of the following 'upwind', 'central', 'exponential'. The upwind option uses the classic first order upwind discretisation, central uses second-order central and setting to exponential uses an adaptive scheme which will use weight between the central and upwind scheme depending on the local value of the Peclet number. This is the classic ‘exponential fitting’ or ‘Scharfetter-Gummel’ discretisation. N.B. Scharfetter-Gummel also refers to a method of solving the advection-diffusion equation is a non-coupled manner, this is not the case here where it only refers to the the discretisation method.
    
-
    
-
    
+
    Finally, when initialising a Model object the keyword argument , discretisation is important. Is can be set to one of the following 'upwind', 'central', 'exponential'. The upwind option uses the classic first order upwind discretisation, central uses second-order central and setting to exponential uses an adaptive scheme which will use weight between the central and upwind scheme depending on the local value of the Peclet number. This is the classic ‘exponential fitting’ or ‘Scharfetter-Gummel’ discretisation. N.B. Scharfetter-Gummel also refers to a method of solving the advection-diffusion equation is a non-coupled manner, this is not the case here where it only refers to the the discretisation method.
    
+
    + -
    -
    -
    +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/index.html b/index.html index 64de78e..d1f412a 100644 --- a/index.html +++ b/index.html @@ -1,186 +1,44 @@ - - - - - + - - Notes on implementing the finite-volume method for physical simulations — FVM Docs 0.1 documentation - - - - + + Notes on implementing the finite-volume method for physical simulations — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Next topic

    -

    Overview

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Next topic

    -

    Overview

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    Notes on implementing the finite-volume method for physical simulations

    + + + +
    +
    +
    +
    + +
    +

    Notes on implementing the finite-volume method for physical simulations

    Contents:

    -
    - + -
    -
    -
    +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/objects.inv b/objects.inv index ac64426..b721504 100644 Binary files a/objects.inv and b/objects.inv differ diff --git a/overview.html b/overview.html index ab0f278..6a87749 100644 --- a/overview.html +++ b/overview.html @@ -1,248 +1,133 @@ - - - - - + - - Overview — FVM Docs 0.1 documentation - - - - + + Overview — FVM Docs 0.1 documentation + + - - - - - - - - + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - + + + +
    +
    +
    +
    + +
    +

    Overview

    +

    These are notes that I am making while learning the finite-volume method (FVM). The methods discussed are based on the the book by W. Hundsdorfer and J. G. Verwer, Numerical solutions of time-dependent advection-diffusion reaction equations.

    +

    See the the accompanying project on Github for implementation of the solvers in Python branch, http://danieljfarrell.github.come/FVM/.

    +
    +

    Finite volumes vs. finite differences

    +

    When first starting to solve (partial) differential equation problems a student will first encounter the finite difference method (FDM). The FDM relies on approximating the a differential operator as a series of distinct points,

    +
    +Finite difference grid + +
    +

    In contrast, the FVM uses a fundamental unit of averages over small spatial elements, or cells,

    +
    +Finite volume grid. + +
    +

    Although this distinction may seem trivial at first or even over complicated it is a particularly power technique when applied to conservation equations because the conservative properties of the governing equation are maintained even at the discretised level without introducing any approximations. This property is achieved by symbolic integrating the equations over a finite volume, then using solve the problem in terms of the integrated variables.

    +

    This sounds much more complicated that it is, let’s look at an example.

    +
    +
    + + +
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    -
    + + - - \ No newline at end of file diff --git a/poisson.html b/poisson.html index 4ed5d10..8c9df21 100644 --- a/poisson.html +++ b/poisson.html @@ -1,245 +1,93 @@ - - - - - + - - The Poisson equation — FVM Docs 0.1 documentation - - - - + + The Poisson equation — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving nonlinear equations

    -

    Next topic

    -

    Solving the advection-diffusion-reaction equation in Python

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving nonlinear equations

    -

    Next topic

    -

    Solving the advection-diffusion-reaction equation in Python

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    The Poisson equation

    + + + +
    +
    +
    +
    + +
    +

    The Poisson equation

    The Poisson equation forms the basis of electrostatics and is of the form,

    -
    +
    \[0 = (\epsilon(x) \phi_x)_x + \rho(x,t)\]
    -

    where \(\phi\) is the electric potential, \(\epsilon(x)\) is the materials dielectric constant, \(\rho(x, t)\) is a charge distribution (possibly varying with time), and the \(x\) subscripts indicate a spatial partial derivative. Note that this equation is reaction-diffusion equation with the diffusion flux \(\mathcal{F} = -\epsilon(x) \phi_x\). As with the advection-diffusion equation, this expression holds over the whole domain, therefore it must also hold over a subdomain \(\Omega_j=(x_{j-1/2}, x_{j+1/2})\). Integrating over the subdomain,

    -
    +

    where \(\phi\) is the electric potential, \(\epsilon(x)\) is the materials dielectric constant, \(\rho(x, t)\) is a charge distribution (possibly varying with time), and the \(x\) subscripts indicate a spatial partial derivative. Note that this equation is reaction-diffusion equation with the diffusion flux \(\mathcal{F} = -\epsilon(x) \phi_x\). As with the advection-diffusion equation, this expression holds over the whole domain, therefore it must also hold over a subdomain \(\Omega_j=(x_{j-1/2}, x_{j+1/2})\). Integrating over the subdomain,

    +
    \[0 = \int_{x_{j-1/2}}^{x_{j+1/2}} (\mathcal{-F})_x~dx + \int_{x_{j-1/2}}^{x_{j+1/2}} \rho(x,t)~dx\]

    The cell averaged values are therefore,

    -
    +
    \[0 = -\frac{\mathcal{F}_{j+1/2}}{h_j} + \frac{\mathcal{F}_{j-1/2}}{h_{j}} + \bar{\rho}_j\]

    To determine the flux at the cell faces we will take the mean of the fluxes in cells adjacent to the interfaces,

    -
    +
    \[\mathcal{F}_{j+\frac{1}{2}} = - \epsilon_{j+\frac{1}{2}} \frac{\phi_{j+1}-\phi_j}{h_{+}}\]
    -
    +
    \[\mathcal{F}_{j-\frac{1}{2}} = - \epsilon_{j-\frac{1}{2}} \frac{\phi_{j}-\phi_{j-1}}{h_{-}}\]

    Substituting the definition of the fluxes into the integrated equation and factoring for the electrical potential terms yields,

    -
    +
    \[0 = \frac{\phi_{j-1}}{h_j} \left( \frac{\epsilon(x_{j-\frac{1}{2}})}{h_{-}} \right) - \frac{\phi_{j}}{h_j} \left( \frac{\epsilon(x_{j-\frac{1}{2}})}{h_{-}} + \frac{\epsilon(x_{j+\frac{1}{2}})}{h_{+}} \right) + \frac{\phi_{j+1}}{h_j} \left( \frac{\epsilon(x_{j+\frac{1}{2}})}{h_{+}} \right) + \bar{\rho}_j\]
    -
    -

    Naturally resolved Neumann boundary conditions

    +
    +

    Naturally resolved Neumann boundary conditions

    Neumann boundary conditions applied to the finite volume form of the Poisson equation are naturally resolved; that is to say that ghost cell or extrapolation approaches are not needed. Consider the integral form of the Poisson equation at the last cell of the domain on the right hand side,

    -
    +
    \[0 = \frac{1}{h_J}\left( \epsilon(x_{R}) \phi_x\bigg|_{x=x_R} - \epsilon(x_{J-1/2}) \phi_x\bigg|_{J-1/2} \right) + \bar{\rho}_J\]
    -

    The Neumann boundary condition defines the derivative of \(\phi\) at the the right hand boundary,

    -
    +

    The Neumann boundary condition defines the derivative of \(\phi\) at the the right hand boundary,

    +
    \[\phi_x\bigg|_{x=x_R} = \sigma_R\]

    This boundary condition appears naturally as a term in the integral equation. Therefore the boundary condition can be satisfied by substitution into the integral form of the equation at the right hand boundary,

    -
    +
    \[\frac{1}{h_J}\left( \epsilon(x_{R}) \sigma_R - \epsilon(x_{J-1/2}) \phi_x\bigg|_{J-1/2} \right) + \bar{\rho}_J\]
    -
    -
    -

    Dirichlet boundary conditions

    -

    Fixed value, or Dirichlet, boundary conditions are not naturally resolved and can only be applied to the discretised form of the equation. We need to apply a Dirichlet condition to the left hand side, such that full equation for the cell \(j=1\) reads,

    -
    +
    +
    +

    Dirichlet boundary conditions

    +

    Fixed value, or Dirichlet, boundary conditions are not naturally resolved and can only be applied to the discretised form of the equation. We need to apply a Dirichlet condition to the left hand side, such that full equation for the cell \(j=1\) reads,

    +
    \[\phi(x_L) = \omega_L\]

    Clearly, this is trivial as all that is required is fix the first row of the matrix equation to be the boundary condition.

    -
    -
    -

    Poisson equation in matrix form

    +
    +
    +

    Poisson equation in matrix form

    Defining the (vector) coefficients,

    -
    +
    \[\begin{split}\ell_a(j) & = \frac{1}{h_j}\frac{\epsilon(x_{j-1/2})}{h_{-}} \\ \ell_b(j) & = -\frac{1}{h_j}\left( \frac{\epsilon_{j-1/2}}{h_{-}} + \frac{\epsilon_{j+1/2}}{h_{+}} \right) \\ \ell_c(j) & = \frac{1}{h_j}\frac{\epsilon(x_{j+1/2})}{h+} \\\end{split}\]

    With Dirichlet boundary condition on the left hand side and a Neumann boundary condition on the right hand side the Poisson equation can be written in matrix form,

    -
    +
    \[\begin{split}\begin{align} \begin{pmatrix} 1 & 0 & & & 0 \\ @@ -265,44 +113,78 @@

    Poisson equation in matrix form

    - - -
    +
    + + + - - \ No newline at end of file diff --git a/search.html b/search.html index 23f3397..439f46c 100644 --- a/search.html +++ b/search.html @@ -1,184 +1,94 @@ - - - - - + + + Search — FVM Docs 0.1 documentation + + - Search — FVM Docs 0.1 documentation - - - - - - - - - - - - - + + + - - - - - - - - - - - - - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -
    -
    -
    - - - -
    - -
    -
    -
    -
    - - -
    -
    -
    -
    -
    - + + + + + + + + + + + +
    +
    +
    +
    +

    Search

    -
    - + + + +

    - From here you can search these documents. Enter your search - words into the box below and click "search". Note that the search - function will automatically search for all of the words. Pages - containing fewer words won't appear in the result list. + Searching for multiple words only shows matches that contain + all words.

    -
    -
    - - -
    + + + + +
    -
    -
    +
    + -
    -
    -
    +
    - - -
    - - -
    -
    -
    - -
    -
    -
    - - -
    - © Copyright 2013, Daniel J Farrell. - Created using Sphinx 1.1.3. -
    - - +
    + +
    +
    + + - - \ No newline at end of file diff --git a/searchindex.js b/searchindex.js index bc9895b..2a3d7ce 100644 --- a/searchindex.js +++ b/searchindex.js @@ -1 +1 @@ -Search.setIndex({objects:{},terms:{all:[8,5,3,6,7],code:[1,6],partial:[2,7],impli:6,illustr:6,global:4,"_robin_boundary_condition_matrix_elements_left":1,prefix:1,extrapol:7,alpha_matrix:1,abil:4,follow:[1,5,8,6],alt:[],row:[5,7],semiconductor:5,privat:1,accur:4,depend:[1,2,5,8,3],readabl:1,articl:4,present:5,larg:4,introduc:[4,2,5,8],"case":[1,6],sourc:4,far:5,netlib:4,ident:[1,4,8],volum:[0,1,2,5,6,7,8],govern:2,veri:[1,5,6],affect:6,solver:[1,2,4,6],list:1,iter:[4,0,5,3],adjust:5,danieljfarrel:2,red:5,small:[4,2,5],div:6,initialis:1,dimens:5,upper:1,cfl:1,fortun:5,natur:[4,0,5,8,7],cell_width:1,direct:1,sign:5,past:4,second:[1,4,5],cost:5,design:1,pass:1,minimis:1,even:[4,2,5],index:[1,0,4],what:[4,5],nonlinear:[4,0],appear:[8,7],compar:[],section:[1,4,5],asid:5,favour:[5,6],uniform:[0,5,6],access:1,version:8,"new":[1,8],method:[0,1,2,3,4,5,6,8],contrast:[4,2],thin:[],full:[4,3,7],deriv:[4,8,7],gener:[1,0,5,8],coeffici:[1,3,4,5,7,8],here:[1,4,5,3,6],satisfi:[8,7],let:[4,2,5,3],along:8,vertic:5,modifi:[4,8],implicit:[1,0,5,4],valu:[0,1,3,4,5,6,7,8],box:4,search:0,convers:5,howev:[1,4,5,8,6],shift:5,observ:[4,6],poisson:[0,8,7],implement:[1,0,5,8,2],commonli:5,via:1,set_boundary_condit:1,approxim:[4,2,5],raphson:5,readili:4,brought:5,total:8,unit:[2,5],highli:4,plot:3,from:[1,4,5,8,6],describ:4,would:[1,4,5],beta_vector:1,distinct:2,two:5,next:3,few:[],"_robin_boundary_condition_vector_elements_right":1,call:[1,4,5,8],centroid:1,criteria:6,type:[4,8],start:[4,2,5],more:[4,2,5,8,6],desir:[1,4],src:6,relat:6,finit:[0,1,2,5,6,7,8],warn:8,"transient":[0,8,6],accept:1,particular:1,known:[4,0,5,8],actual:[1,4,5],hole:1,hold:[5,7],must:[4,5,8,7],"_dirichlet_boundary_condition_matrix_elements_left":1,topic:4,player:6,generalis:[4,5],work:4,other:8,remain:[4,8],can:[1,3,4,5,6,7,8],learn:2,def:1,tau:1,explor:8,give:[4,5],predict:4,calcul:[0,1,3,4,5,6],indic:[1,0,7],high:[1,6],want:5,occur:1,alwai:[5,8],differenti:[4,2],multipl:4,hoc:5,divid:5,rather:1,write:[4,5],how:[1,4,5],answer:5,simpl:5,updat:[4,3],surfac:1,mesh:[1,0,6],mai:2,underscor:1,data:[1,4],grow:4,physic:[0,5],goal:1,stabil:[4,5,6],"short":4,attempt:1,practic:[4,6],explicit:[1,0,5,4],correspond:[1,5],element:[1,2,8],inform:8,maintain:2,combin:5,incorpor:8,order:[1,4,5],over:[2,5,7],move:3,becaus:[1,2,4,5,6,8],moulton:4,still:[4,5,6],vari:[1,7],dynam:4,paramet:[1,3,6],style:6,scholarpedia:4,img:[],fix:[8,7],better:4,main:1,alter:[5,8],non:[1,4,5],good:[8,6],"return":1,thei:[1,4,5,8],python:[1,0,4,2],initi:[8,5,3,6],number:[1,4,5,6],facilit:1,now:5,discuss:[1,2,5,4,8],nor:3,introduct:4,choic:5,term:[0,1,2,3,4,5,7,8],subdomain:[5,7],name:[1,5],stiff:4,linspac:1,easili:1,achiev:[4,2,8],upwind:[1,0,5,6],mode:4,compris:8,each:[1,5],electrostat:7,complet:1,side:[8,7],mean:[8,5,3,7],domain:[5,6,7],weight:[1,5,6],replac:[4,3],hard:1,continu:[4,5],procedur:[8,5,3],realli:4,significantli:[5,6],our:4,out:[1,4],variabl:[1,2,3,4,5,8],shown:[4,5,3,6],matrix:[0,1,3,4,5,7,8],space:[5,6],newton_solution_fvm:[],semi:[4,5,8],publish:[],content:0,vector:[1,4,5,8,7],adapt:[1,0,5,4],internet:4,formula:[4,5],forth:3,correct:1,matric:1,common:5,linear:[1,0,5,4,3],written:[4,5,8,7],situat:5,given:5,standard:[4,6],reason:[1,4],base:[2,5],ask:1,org:4,"_dirichlet_boundary_condition_matrix_elements_right":1,basi:7,dirichlet:[0,8,7],veloc:1,exponenti:[0,1,3,4,5,6],webkitallowfullscreen:6,recov:[5,8],caption:[],enforc:[],"5th":[],first:[1,2,3,4,5,6,7,8],oper:[4,2],rang:1,perfectli:4,directli:4,arrai:1,scene:[],sometim:5,restrict:4,unlik:4,alreadi:4,done:4,wrapper:4,symmetri:6,stabl:[4,5,6],open:4,fanci:1,differ:[1,0,5,8,2],stencil:[1,5],advect:[0,1,2,3,4,5,7,8],unknown:[4,5],top:6,unrestrict:4,system:[1,4,5,8,3],construct:1,downsid:5,too:[1,4],statement:[],similarli:[1,5,8],scheme:[1,4,5,8,6],"final":[1,5,6],left_valu:1,copi:1,specifi:[1,8],part:4,accompani:2,somewhat:[],conserv:2,hundsdorf:[2,5],domin:5,than:[1,4],png:[],wide:4,conveni:[1,5],keyword:1,provid:[1,4,8,3],eleg:[1,5],zero:8,exampl:[0,1,2,4,5,6,8],project:[4,2],posit:[1,5],video:6,toward:6,seri:2,sai:[5,7],comput:[5,6],ani:[4,2,5,8],manner:[1,5],have:[1,4,5,8,6],further:5,need:[1,4,8,7],seem:[4,2],minu:1,optimis:[],option:1,automat:[1,5,8,6],techniqu:[4,2,5],diverg:5,moreov:[4,5,8],accuraci:[4,5,6],note:[0,1,2,3,4,5,6,7,8],also:[1,4,5,8,7],discret:5,take:[4,7],which:[1,4,5,8,6],jacobian:4,interior:1,singl:4,coefficient_matrix:1,simplifi:5,sure:8,unless:5,allow:[1,6],though:5,multipli:8,object:1,most:[4,5],beta:1,bdf:4,dielectr:7,pair:1,alpha:1,mozallowfullscreen:6,"class":[1,0,8],charg:7,homogen:1,don:5,later:1,request:1,doe:[1,5,6],determin:[5,7],method_of_lin:[],left:[1,8,7],steadi:[5,6],fact:[],show:[4,8,6],random:[0,6],particularli:[4,2,5,8],vectoris:1,trivial:[2,8,7],find:5,involv:3,cell:[1,2,5,6,7,8],onli:[1,4,5,8,7],explicitli:4,locat:[1,5],solut:[2,3,4,5,6,8],state:[4,5,6],should:[1,5,8],black:4,factor:[5,3,7],analyt:4,local:[1,5,6],contribut:5,variou:1,get:8,express:[5,3,7],cannot:[4,5],theta:1,increas:[5,6],subscript:[1,5,7],requir:[1,5,8,7],yield:[4,5,3,7],integr:[4,2,5,7],contain:[1,4,5,8,6],where:[1,3,4,5,6,7,8],set:[1,4,5],smooth:6,see:[1,2,5,4],result:[4,5],fail:4,boundari:[0,1,3,4,6,7,8],label:5,scharfett:1,behind:[],between:[1,4,5,8],"import":[1,5,8],awai:1,experi:4,approach:[1,3,4,5,7,8],altern:5,spars:1,numer:[4,2,5],averag:[1,2,5,7],abruptli:6,fipi:1,uncondition:[4,5],solv:[0,1,2,3,4,5],come:[4,2],addit:1,reaction:[0,1,2,3,4,5,7,8],last:[8,7],fit:[1,0,5,3,6],region:6,etc:5,instanc:1,multistep:4,context:1,improv:[4,6],pde:[4,8],whole:[5,7],load:4,simpli:[1,5,8],pose:[4,5],figur:[4,5,3],overview:[1,0,2],inspir:1,unrealist:5,height:6,distinguish:4,"_robin_boundary_condition_matrix_elements_right":1,mistak:1,assum:[5,8],math:[],late:5,becom:[4,5,8],evalu:[4,5,8],coupl:[1,4],numpi:1,three:[4,5,6],been:[1,4,5,6],github:2,compon:[5,8],accumul:4,much:[1,2,4],treat:4,interest:[4,6],electr:7,unambigu:1,dure:4,suppress:4,argument:1,deliber:4,those:4,discretis:[0,1,2,3,4,5,6,7,8],sound:2,exception:5,therefor:[4,5,8,7],look:[1,2],quantiti:[1,5],vimeo:6,align:[],properti:2,histor:4,"_robin_boundary_condition_vector_elements_left":1,defin:[1,4,5,7],"while":[2,5],abov:[1,4,5,8,6],error:4,aid:1,margin:6,robin:[0,8],layer:6,advantag:4,almost:4,centr:[5,6],tabl:[0,8],henc:5,destin:1,blah:[],scipi:4,maximis:1,develop:[4,8],perform:[1,4,5],make:[1,2,8],same:[1,5,8,6],preserv:5,complex:1,document:[],higher:6,ifram:6,http:[2,6],hand:[8,6,7],fairli:[1,4],user:1,student:2,mani:4,robust:[4,5],chang:[1,8,6],recent:4,lower:[1,4],off:4,center:[],choos:[4,5],firstli:[5,3],com:6,well:4,thought:4,without:[0,5,8,2],thi:[1,2,3,4,5,7,8],interpol:[1,5,8],model:[1,0],self:1,usual:4,distanc:[1,5],just:5,newtom:4,when:[1,2,5,6],newton:[4,0,5],unfair:[],rest:[],simultan:4,yet:4,previous:5,fvm:2,cluster:6,easi:[1,4],point:[1,2,3,4,5,8],cellvari:[1,0],framebord:6,add:8,valid:[4,5,8],densiti:6,role:1,modul:0,build:[],real:[1,4],applic:4,insert:1,read:7,imex:[4,0],spatial:[1,2,5,4,7],grid:[4,0,5,6],know:[4,5],verwer:2,cfl_condit:1,like:[1,6],specif:1,resolv:[0,8,6,7],flux:[0,5,8,7],unstabl:[4,5],either:8,"_dirichlet_boundary_condition_vector_elements_left":1,popular:[4,5],page:0,advectiondiffusionmodel:1,encount:2,www:[],drop:5,often:5,deal:[1,5],creation:1,some:1,understood:5,intern:1,respect:[1,4],a11:1,transport:8,distribut:[1,7],basic:8,peclet_numb:1,definit:[4,5,7],subclass:1,substitut:[5,8,7],estim:4,leav:4,condit:[0,1,3,4,5,6,7,8],gummel:1,complic:[1,2],refer:1,plu:1,previou:[1,4,6],nonuniform:[0,6],power:[4,2],"10px":6,usag:1,agreement:4,step:[0,1,3,4,5,6,8],although:[4,2,8,6],found:5,appli:[1,2,4,5,7,8],impos:8,fulli:[1,4,5],comparison:[],central:[1,5,6],face:[1,5,7],materi:7,spsolv:1,plai:1,simul:[0,3,6],allowfullscreen:6,stand:[1,5],constructor:1,produc:[],invok:8,within:8,bound:[],terminolog:4,diagon:1,right:[1,4,8,6,7],reli:2,precomput:4,neumann:[0,8,7],wai:[1,4],vode:4,uniformli:1,"long":4,avail:4,width:[1,5,6],adjac:7,peclet:[1,5,6],interfac:7,includ:[1,4,8],forward:1,"function":[1,4,5],ghost:[8,7],form:[0,3,4,5,6,7,8],enough:[4,6],regard:8,oscil:6,"_dirichlet_boundary_condition_vector_elements_right":1,cell_vari:1,diffus:[0,1,2,3,4,5,7,8],crank:[4,5],line:[4,0,5],level:[1,2],"true":[],bug:[1,6],notat:5,furthermor:6,consist:5,possibl:[4,5,8,7],clearli:7,fdm:2,below:[1,4,5,8,3],limit:[5,8],fundament:[2,5],highlight:1,problem:[1,2,5,4,8],similar:[1,8,6],expect:4,featur:[1,4],constant:[1,8,7],creat:1,classic:1,"abstract":1,accentu:4,"_interior_matrix_el":1,repres:5,proport:4,rearrang:[4,5,3],check:1,again:1,readi:[],equat:[0,1,2,3,4,5,6,7,8],excel:4,detail:[1,4],right_valu:1,book:2,b12:1,futur:[4,5,3],branch:2,varieti:4,test:[4,6],nice:[],poor:4,matur:4,intend:1,benefici:1,symbol:[4,2,8],ndarrai:1,scale:[],consid:[1,8,7],fisher:[4,0],equidist:5,reduc:[1,8,6],invari:8,longer:4,algorithm:4,nicolson:[4,5],pseudo:1,u_init:1,ignor:[4,5],adam:4,potenti:7,time:[0,1,2,3,4,5,6,7,8],discrist:4,backward:[1,4]},objtypes:{},titles:["Notes on implementing the finite-volume method for physical simulations","Solving the advection-diffusion-reaction equation in Python","Overview","Solving linear equations","Solving nonlinear equations","The advection-diffusion-reaction equation","Example calculations","The Poisson equation","Boundary conditions for the advection-diffusion-reaction equation"],objnames:{},filenames:["index","implementation_advection_diffusion","overview","solving_linear_equations","solving_nonlinear_equations","advection_diffusion","examples","poisson","advection_diffusion_boundary_conditions"]}) \ No newline at end of file +Search.setIndex({"alltitles": {"Adaptive upwinding & exponential fitting": [[0, "adaptive-upwinding-exponential-fitting"]], "Aside: Linear interpolation between cell centre and face values": [[0, "aside-linear-interpolation-between-cell-centre-and-face-values"]], "Boundary conditions for the advection-diffusion-reaction equation": [[1, "boundary-conditions-for-the-advection-diffusion-reaction-equation"]], "Dirichlet boundary conditions": [[6, "dirichlet-boundary-conditions"]], "Dirichlet boundary conditions (known value)": [[1, "dirichlet-boundary-conditions-known-value"]], "Discretised equation in matrix form": [[0, "discretised-equation-in-matrix-form"]], "Example calculations": [[2, "example-calculations"]], "Explicit and implicit forms": [[0, "explicit-and-implicit-forms"]], "Finite volumes vs. finite differences": [[5, "finite-volumes-vs-finite-differences"]], "Fisher\u2019s equation": [[8, "fisher-s-equation"]], "General matrix form with a transient term": [[1, "general-matrix-form-with-a-transient-term"]], "General matrix form without a transient term": [[1, "general-matrix-form-without-a-transient-term"]], "IMEX time-stepping": [[8, "imex-time-stepping"]], "Indices and tables": [[4, "indices-and-tables"]], "Method of lines": [[8, "method-of-lines"]], "Naturally resolved Neumann boundary conditions": [[6, "naturally-resolved-neumann-boundary-conditions"]], "Newton iteration": [[8, "newton-iteration"]], "Nonuniform grid": [[2, "nonuniform-grid"]], "Notes on implementing the finite-volume method for physical simulations": [[4, "notes-on-implementing-the-finite-volume-method-for-physical-simulations"]], "Overview": [[3, "overview"], [5, "overview"]], "Poisson equation in matrix form": [[6, "poisson-equation-in-matrix-form"]], "Random grid": [[2, "random-grid"]], "Robin boundary conditions (known flux)": [[1, "robin-boundary-conditions-known-flux"]], "Solving linear equations": [[7, "solving-linear-equations"]], "Solving nonlinear equations": [[8, "solving-nonlinear-equations"]], "Solving the advection-diffusion-reaction equation in Python": [[3, "solving-the-advection-diffusion-reaction-equation-in-python"]], "The CellVariable class": [[3, "the-cellvariable-class"]], "The Mesh class": [[3, "the-mesh-class"]], "The Model class": [[3, "the-model-class"]], "The Poisson equation": [[6, "the-poisson-equation"]], "The \\theta-method": [[0, "the-theta-method"]], "The advection-diffusion-reaction equation": [[0, "the-advection-diffusion-reaction-equation"]], "Time stepping the advection-diffusion equation": [[7, "time-stepping-the-advection-diffusion-equation"]], "Uniform grid": [[2, "uniform-grid"]]}, "docnames": ["advection_diffusion", "advection_diffusion_boundary_conditions", "examples", "implementation_advection_diffusion", "index", "overview", "poisson", "solving_linear_equations", "solving_nonlinear_equations"], "envversion": {"sphinx": 61, "sphinx.domains.c": 3, "sphinx.domains.changeset": 1, "sphinx.domains.citation": 1, "sphinx.domains.cpp": 9, "sphinx.domains.index": 1, "sphinx.domains.javascript": 3, "sphinx.domains.math": 2, "sphinx.domains.python": 4, "sphinx.domains.rst": 2, "sphinx.domains.std": 2}, "filenames": ["advection_diffusion.rst", "advection_diffusion_boundary_conditions.rst", "examples.rst", "implementation_advection_diffusion.rst", "index.rst", "overview.rst", "poisson.rst", "solving_linear_equations.rst", "solving_nonlinear_equations.rst"], "indexentries": {}, "objects": {}, "objnames": {}, "objtypes": {}, "terms": {"": [0, 1, 4, 5, 7], "0": [0, 1, 2, 3, 6, 7, 8], "001": 2, "004": 2, "01": [2, 3], "02": 2, "05": 7, "1": [0, 1, 2, 3, 6, 7, 8], "10": 7, "100": [2, 3, 7], "15": 7, "1e": 3, "2": [0, 1, 3, 6, 7, 8], "20": 2, "266": 0, "2h_": [0, 1], "3": [0, 1, 3, 7], "4": 7, "400": 7, "5": [3, 7, 8], "6": 8, "9": 2, "98": 2, "A": [0, 3, 7, 8], "As": [3, 6], "At": 3, "But": 8, "For": [0, 1, 2, 3, 8], "If": [0, 1], "In": [0, 3, 5, 8], "It": [0, 3, 8], "One": 0, "The": [1, 2, 4, 5, 7, 8], "There": [3, 8], "These": [1, 5], "To": 6, "With": [2, 6, 8], "_": [0, 1, 6, 7], "_1": 1, "_2": 6, "_dirichlet_boundary_condition_matrix_elements_left": 3, "_dirichlet_boundary_condition_matrix_elements_right": 3, "_dirichlet_boundary_condition_vector_elements_left": 3, "_dirichlet_boundary_condition_vector_elements_right": 3, "_interior_matrix_el": 3, "_j": [0, 1, 6], "_robin_boundary_condition_matrix_elements_left": 3, "_robin_boundary_condition_matrix_elements_right": 3, "_robin_boundary_condition_vector_elements_left": 3, "_robin_boundary_condition_vector_elements_right": 3, "_x": [0, 6], "a11": 3, "a_": 0, "a_j": 1, "abil": 8, "abov": [0, 1, 2, 3, 8], "abruptli": 2, "abstract": 3, "accentu": 8, "accept": 3, "access": 3, "accompani": 5, "accumul": 8, "accur": 8, "accuraci": [0, 2, 8], "achiev": [1, 5, 8], "actual": [0, 3, 8], "ad": [0, 3], "adam": 8, "adapt": [3, 4, 8], "add": 1, "addit": 3, "adjac": 6, "adjust": 0, "advantag": 8, "advect": [4, 5, 6, 8], "advectiondiffusionmodel": 3, "affect": 2, "again": 3, "agreement": 8, "ah": 0, "aid": 3, "algorithm": 8, "align": [0, 1, 6], "all": [0, 1, 2, 6, 7], "allow": [2, 3], "almost": 8, "along": 1, "alpha": [1, 3], "alpha_1": 1, "alpha_j": 1, "alpha_matrix": 3, "alreadi": 8, "also": [0, 1, 3, 6, 8], "alter": [0, 1], "altern": 0, "although": [1, 2, 5, 8], "alwai": [0, 1], "am": 5, "an": [0, 1, 3, 5, 8], "analyt": 8, "ani": [0, 1, 5, 8], "answer": 0, "appear": [1, 6], "appli": [0, 1, 3, 5, 6, 8], "applic": 8, "approach": [0, 1, 3, 6, 7, 8], "approx": 2, "approx0": 2, "approxim": [0, 5, 8], "ar": [0, 1, 2, 3, 5, 6, 7, 8], "argument": 3, "arrai": 3, "articl": 8, "ask": 3, "assum": [0, 1], "attempt": 3, "au": 0, "automat": [0, 1, 2, 3], "avail": 8, "averag": [0, 3, 5, 6], "awai": 3, "b": 3, "b12": 3, "b_1": 1, "b_j": 1, "backward": [3, 8], "bar": [0, 1, 6], "base": [0, 5], "basi": 6, "basic": 1, "bdf": 8, "becaus": [0, 1, 2, 3, 5, 8], "becom": [0, 1, 8], "been": [0, 2, 3, 8], "begin": [0, 1, 6], "being": [0, 3], "below": [0, 1, 3, 7, 8], "benefici": 3, "beta": [1, 3, 8], "beta_1": 1, "beta_j": 1, "beta_vector": 3, "better": 8, "between": [1, 3, 8], "bigg": 6, "black": 8, "book": 5, "boundari": [2, 3, 4, 7, 8], "box": 8, "branch": 5, "brought": 0, "bug": [2, 3], "c": 8, "c_1": 1, "calcul": [0, 3, 4, 7, 8], "call": [0, 1, 3, 8], "can": [0, 1, 2, 3, 6, 7, 8], "cannot": [0, 8], "case": [0, 2, 3], "cdot": [0, 1, 7, 8], "cell": [1, 2, 3, 5, 6], "cell_vari": 3, "cell_width": 3, "cellvari": 4, "centr": 2, "central": [0, 2, 3], "centroid": 3, "cfl": 3, "cfl_condit": 3, "chang": [1, 2, 3], "charg": 6, "check": 3, "choic": 0, "choos": [0, 8], "class": [1, 4], "classic": 3, "clearli": 6, "cluster": 2, "code": [2, 3], "coeffici": [0, 1, 3, 6, 7, 8], "coefficient_matrix": 3, "combin": 0, "come": [5, 8], "common": 0, "commonli": 0, "complet": 3, "complex": 3, "complic": [3, 5], "compon": [0, 1], "compris": 1, "comput": [0, 2], "condit": [0, 2, 3, 4, 7, 8], "conserv": 5, "consid": [1, 3, 6], "consist": 0, "constant": [1, 3, 6], "construct": 3, "constructor": 3, "contain": [0, 1, 2, 3, 8], "content": 4, "context": 3, "continu": [0, 8], "contrast": [5, 8], "contribut": 0, "conveni": [0, 3], "convers": 0, "copi": 3, "correct": 3, "correspond": [0, 3], "cost": 0, "coupl": [3, 8], "crank": [0, 8], "creat": 3, "creation": 3, "criteria": 2, "d": [0, 1, 2, 3, 7, 8], "d_": 0, "danieljfarrel": 5, "data": [3, 8], "ddot": [0, 1, 6], "deal": [0, 3], "def": 3, "defin": [0, 3, 6, 8], "definit": [0, 6, 8], "deliber": 8, "delta": 1, "densiti": 2, "depend": [0, 1, 3, 5, 7], "deriv": [1, 6, 8], "describ": 8, "design": 3, "desir": [3, 8], "destin": 3, "detail": [3, 8], "determin": [0, 6], "develop": [1, 8], "diag": 1, "diagon": 3, "dielectr": 6, "differ": [0, 1, 3, 4], "differenti": [5, 8], "diffus": [4, 5, 6, 8], "dimens": 0, "direct": 3, "directli": 8, "dirichlet": 4, "discret": 0, "discretis": [1, 2, 3, 4, 5, 6, 7, 8], "discrist": 8, "discuss": [0, 1, 3, 5, 8], "distanc": [0, 3], "distinct": 5, "distinguish": 8, "distribut": [3, 6], "diverg": 0, "divid": 0, "do": [1, 3, 8], "doe": [0, 2, 3], "domain": [0, 2, 6], "domin": 0, "don": 0, "done": 8, "downsid": 0, "drop": 0, "dt": 8, "du": 8, "du_x": 0, "dure": 8, "dx": [0, 6], "dynam": 8, "e": [0, 8], "each": [0, 3], "easi": [3, 8], "easili": 3, "either": 1, "electr": 6, "electrostat": 6, "eleg": [0, 3], "element": [1, 3, 5], "ell_a": 6, "ell_b": 6, "ell_c": 6, "encount": 5, "end": [0, 1, 6], "enough": [2, 8], "epsilon": 6, "epsilon_": 6, "equat": [2, 4, 5], "equidist": 0, "error": 8, "estim": 8, "etc": 0, "evalu": [0, 1, 8], "even": [0, 5, 8], "exampl": [0, 1, 3, 4, 5, 8], "excel": 8, "exception": 0, "expect": 8, "experi": 8, "explicit": [3, 4, 8], "explicitli": 8, "explor": 1, "exponenti": [2, 3, 4, 7, 8], "express": [0, 6, 7], "extrapol": 6, "f": [0, 1, 6, 7, 8], "f_": 0, "f_1": 1, "face": [3, 6], "facilit": 3, "factor": [0, 6, 7], "fail": 8, "fairli": [3, 8], "fanci": 3, "far": 0, "favour": [0, 2], "fdm": 5, "featur": [3, 8], "figur": [0, 7, 8], "final": [0, 2, 3], "find": 0, "finit": [0, 1, 2, 3, 6], "fipi": 3, "first": [0, 1, 2, 3, 5, 6, 7, 8], "firstli": [0, 7], "fisher": 4, "fit": [2, 3, 4, 7], "fix": [1, 6], "flux": [0, 4, 6], "follow": [0, 1, 2, 3], "form": [2, 4, 7, 8], "formula": [0, 8], "forth": 7, "fortun": 0, "forward": 3, "found": 0, "frac": [0, 1, 6, 7, 8], "from": [0, 1, 2, 3, 8], "full": [6, 7, 8], "fulli": [0, 3, 8], "function": [0, 3, 8], "fundament": [0, 5], "further": 0, "furthermor": 2, "futur": [0, 7, 8], "fvm": 5, "g": 5, "g_": 1, "g_d": 1, "g_r": 1, "gener": [0, 3, 4], "generalis": [0, 8], "geq0": 8, "get": 1, "ghost": [1, 6], "github": 5, "give": [0, 8], "given": 0, "global": 8, "goal": 3, "good": [1, 2], "govern": 5, "grid": [0, 4, 8], "grow": 8, "gummel": 3, "h": [0, 1, 2, 3, 6, 7, 8], "h_": [0, 1, 6], "h_1": 1, "h_j": [0, 1, 6], "ha": [0, 2, 3, 8], "hand": [1, 2, 6], "hard": 3, "have": [0, 1, 2, 3, 8], "henc": 0, "here": [0, 2, 3, 7, 8], "high": [2, 3], "higher": 2, "highli": 8, "highlight": 3, "histor": 8, "hm": 3, "hoc": 0, "hold": [0, 6], "hole": 3, "homogen": 3, "how": [0, 3, 8], "howev": [0, 1, 2, 3, 8], "hp": 3, "http": 5, "hundsdorf": [0, 5], "i": [0, 1, 2, 3, 5, 6, 7, 8], "ident": [1, 3, 8], "ignor": [0, 8], "illustr": 2, "imex": 4, "implement": [0, 1, 3, 5], "impli": 2, "implicit": [3, 4, 8], "import": [0, 1, 3], "impos": 1, "improv": [2, 8], "includ": [1, 3, 8], "incorpor": 1, "increas": [0, 2], "index": [3, 4, 8], "indic": [3, 6], "inform": 1, "infti": 0, "initi": [0, 1, 2, 7], "initialis": 3, "insert": 3, "inspir": 3, "instanc": 3, "int_": [0, 6], "integr": [0, 5, 6, 8], "intend": 3, "interest": [2, 8], "interfac": 6, "interior": 3, "intern": 3, "internet": 8, "interpol": [1, 3], "introduc": [0, 1, 5, 8], "introduct": 8, "invari": 1, "invok": 1, "involv": 7, "iter": [0, 4, 7], "j": [0, 1, 3, 5, 6], "jacobian": 8, "just": 0, "k": [3, 8], "kappa": [0, 2], "kappa_": 0, "keyword": 3, "know": [0, 8], "known": [0, 4, 8], "l": [7, 8], "label": 0, "larg": 8, "last": [1, 6], "late": 0, "later": 3, "layer": 2, "learn": 5, "leav": 8, "left": [0, 1, 3, 6, 8], "left_valu": 3, "leq": 2, "let": [0, 5, 7, 8], "level": [3, 5], "like": [2, 3], "limit": [0, 1], "line": [0, 4], "linear": [3, 4, 8], "linspac": 3, "list": 3, "load": 8, "local": [0, 2, 3], "locat": [0, 3], "long": 8, "longer": 8, "look": [3, 5], "lower": [3, 8], "m": [0, 1, 3, 7, 8], "mai": 5, "main": 3, "maintain": 5, "make": [1, 3, 5], "mani": 8, "manner": [0, 3], "materi": 6, "mathcal": [0, 1, 6], "matric": 3, "matrix": [3, 4, 7, 8], "matur": 8, "max": [0, 2], "maximis": 3, "mean": [0, 1, 6, 7], "mesh": [2, 4], "method": [1, 2, 3, 5, 7], "min": [0, 2], "minimis": 3, "minu": 3, "mistak": 3, "mode": 8, "model": 4, "modifi": [1, 8], "modul": 4, "more": [0, 1, 2, 5, 8], "moreov": [0, 1, 8], "most": [0, 8], "moulton": 8, "move": 7, "mu": [0, 2], "much": [3, 5, 8], "multipl": 8, "multipli": 1, "multistep": 8, "must": [0, 1, 6, 8], "mw": 1, "my": 3, "n": [0, 1, 3, 7, 8], "name": [0, 3], "natur": [0, 1, 4, 8], "ndarrai": 3, "need": [1, 3, 6, 8], "neq0": 0, "netlib": 8, "neumann": [1, 4], "new": [1, 3], "newtom": 8, "newton": [0, 4], "next": 7, "nicolson": [0, 8], "non": [0, 3, 8], "nonlinear": 4, "nonuniform": 4, "nor": 7, "notat": 0, "note": [0, 1, 2, 3, 5, 6, 8], "now": 0, "np": 3, "nu": 8, "number": [0, 2, 3, 8], "numer": [0, 5, 8], "numpi": 3, "object": 3, "observ": [2, 8], "occur": 3, "od": 8, "off": 8, "often": 0, "omega_1": 1, "omega_j": [0, 1, 6], "omega_l": 6, "one": [0, 3, 8], "onli": [0, 1, 3, 6, 8], "open": 8, "oper": [5, 8], "option": 3, "order": [0, 3, 8], "org": 8, "oscil": 2, "other": 1, "our": 8, "out": [3, 8], "over": [0, 5, 6], "overview": 4, "p": 3, "page": 4, "pair": 3, "paramet": [2, 3, 7], "part": 8, "partial": [5, 6], "particular": 3, "particularli": [0, 1, 5, 8], "pass": 3, "past": 8, "pde": [1, 8], "peclet": [0, 2, 3], "peclet_numb": 3, "perfectli": 8, "perform": [0, 3, 8], "pg": 0, "phi": 6, "phi_": 6, "phi_1": 6, "phi_2": 6, "phi_j": 6, "phi_x": 6, "physic": 0, "pi": [2, 7], "plai": 3, "plot": 7, "plu": 3, "pm": 0, "pm1": 0, "pmatrix": [0, 1, 6], "point": [0, 1, 3, 5, 7, 8], "poisson": [1, 4], "poor": 8, "popular": [0, 8], "pose": [0, 8], "posit": [0, 3], "possibl": [0, 1, 8], "possibli": 6, "potenti": 6, "power": [5, 8], "practic": [2, 8], "precomput": 8, "predict": 8, "prefix": 3, "present": 0, "preserv": 0, "previou": [2, 3, 8], "previous": 0, "prime": [0, 1, 8], "privat": 3, "problem": [0, 1, 3, 5, 8], "procedur": [0, 1, 7], "project": [5, 8], "properti": 5, "proport": 8, "provid": [1, 3, 7, 8], "pseudo": 3, "python": [4, 5, 8], "quad": [0, 1], "quantiti": [0, 3], "r": [0, 1, 3, 6], "r_a": [0, 1], "r_b": [0, 1], "r_c": [0, 1], "ra": 3, "random": 4, "rang": 3, "raphson": 0, "rather": 3, "rb": 3, "rc": 3, "reaction": [4, 5, 6, 7, 8], "read": 6, "readabl": 3, "readili": 8, "real": [3, 8], "realli": 8, "rearrang": [0, 7, 8], "reason": [3, 8], "recent": 8, "recov": [0, 1], "red": 0, "reduc": [1, 2, 3], "refer": 3, "regard": 1, "region": 2, "relat": 2, "reli": 5, "remain": [1, 8], "replac": [7, 8], "repres": 0, "request": 3, "requir": [0, 1, 3, 6], "resolv": [1, 2, 4], "respect": [3, 8], "restrict": 8, "result": [0, 8], "return": 3, "rho": 6, "right": [0, 1, 2, 3, 6, 8], "right_valu": 3, "rightarrow": 0, "robin": 4, "robust": [0, 8], "role": 3, "row": [0, 6], "sai": [0, 6], "same": [0, 1, 2, 3], "satisfi": [1, 6], "scharfett": 3, "scheme": [0, 1, 2, 3, 8], "scholarpedia": 8, "scipi": 8, "search": 4, "second": [0, 3, 8], "section": [0, 3, 8], "see": [0, 3, 5, 8], "seem": [5, 8], "self": 3, "semi": [0, 1, 8], "semiconductor": 0, "seri": 5, "set": [0, 3, 8], "set_boundary_condit": 3, "sgn": 0, "shift": 0, "short": 8, "should": [0, 1, 3], "show": [1, 2, 8], "shown": [0, 2, 7, 8], "side": [1, 6], "sigma_r": 6, "sign": 0, "significantli": [0, 2], "similar": [1, 2, 3], "similarli": [0, 1, 3], "simpl": 0, "simpli": [0, 1, 3], "simplifi": 0, "simul": [2, 7], "simultan": 8, "sin": [2, 7], "singl": 8, "situat": 0, "small": [0, 5, 8], "smooth": 2, "so": [0, 2, 3, 7, 8], "solut": [0, 1, 2, 5, 7, 8], "solv": [0, 4, 5], "solver": [2, 3, 5, 8], "some": 3, "sometim": 0, "sound": 5, "sourc": 8, "space": [0, 2], "spars": 3, "spatial": [0, 3, 5, 6, 8], "specif": 3, "specifi": [1, 3], "spsolv": 3, "stabil": [0, 2, 8], "stabl": [0, 2, 8], "stand": [0, 3], "standard": [2, 8], "start": [0, 5, 8], "state": [0, 2, 8], "steadi": [0, 2], "stencil": [0, 3], "step": [0, 1, 2, 3, 4], "stiff": 8, "still": [0, 2, 8], "student": 5, "subclass": 3, "subdomain": [0, 6], "subscript": [0, 3, 6], "substitut": [0, 1, 6], "suppress": 8, "sure": 1, "surfac": 3, "symbol": [1, 5, 8], "symmetri": 2, "system": [0, 1, 3, 7, 8], "t": [0, 2, 6, 7, 8], "t_": [0, 1, 7], "t_n": 0, "tabl": 1, "take": [6, 8], "tau": [0, 2, 3, 7, 8], "techniqu": [0, 5, 8], "term": [0, 3, 4, 5, 6, 7, 8], "terminologi": 8, "test": [2, 8], "text": [0, 1, 2], "than": [3, 8], "thei": [0, 1, 3, 8], "therefor": [0, 1, 6, 8], "theta": [1, 3, 4, 7, 8], "thi": [0, 1, 3, 5, 6, 7, 8], "those": 8, "though": 0, "thought": 8, "three": [0, 2, 8], "time": [0, 1, 2, 3, 4, 5, 6], "too": [3, 8], "topic": 8, "total": 1, "toward": 2, "transient": [2, 4], "transport": 1, "treat": 8, "trivial": [1, 5, 6], "two": 0, "type": [1, 8], "u": [0, 2, 3, 7, 8], "u_": 8, "u_init": 3, "u_t": [0, 8], "u_x": 7, "unambigu": 3, "uncondition": [0, 8], "underbrac": 7, "underscor": 3, "understood": 0, "uniform": [0, 4], "uniformli": 3, "unit": [0, 5], "unknown": [0, 8], "unless": 0, "unlik": 8, "unrealist": 0, "unrestrict": 8, "unstabl": [0, 8], "up": 2, "updat": [7, 8], "upper": 3, "upwind": [2, 3, 4], "us": [0, 1, 2, 3, 5, 7, 8], "usag": 3, "user": 3, "usual": 8, "ut": 8, "v": 4, "valid": [0, 1, 8], "valu": [2, 3, 4, 6, 7, 8], "vari": [3, 6], "variabl": [0, 1, 3, 5, 7, 8], "varieti": 8, "variou": 3, "vdot": 6, "vector": [0, 1, 3, 6, 8], "vectoris": 3, "veloc": 3, "veri": [0, 2, 3], "version": 1, "vertic": 0, "verwer": 5, "via": 3, "vode": 8, "volum": [0, 1, 2, 3, 6], "w": [0, 1, 5, 7, 8], "w_": [0, 1], "w_1": 1, "w_j": [0, 1], "wa": 7, "wai": [3, 8], "want": 0, "we": [0, 1, 2, 3, 6, 7, 8], "weight": [0, 2, 3], "well": 8, "what": [0, 8], "when": [0, 2, 3, 5], "where": [0, 1, 2, 3, 6, 7, 8], "which": [0, 1, 2, 3, 8], "while": [0, 5], "whole": [0, 6], "wide": 8, "width": [0, 3], "within": 1, "without": [0, 4, 5], "work": 8, "would": [0, 3, 8], "wrapper": 8, "write": [0, 8], "written": [0, 1, 6, 8], "x": [0, 2, 3, 6, 7, 8], "x_": [0, 1, 6], "x_0": [0, 8], "x_1": [0, 1, 7], "x_j": [0, 1, 6, 7], "x_l": [1, 6], "x_r": [1, 6], "xx": [0, 8], "yet": 8, "yield": [0, 6, 7, 8], "zero": 1}, "titles": ["The advection-diffusion-reaction equation", "Boundary conditions for the advection-diffusion-reaction equation", "Example calculations", "Solving the advection-diffusion-reaction equation in Python", "Notes on implementing the finite-volume method for physical simulations", "Overview", "The Poisson equation", "Solving linear equations", "Solving nonlinear equations"], "titleterms": {"": 8, "The": [0, 3, 6], "adapt": 0, "advect": [0, 1, 3, 7], "asid": 0, "between": 0, "boundari": [1, 6], "calcul": 2, "cell": 0, "cellvari": 3, "centr": 0, "class": 3, "condit": [1, 6], "differ": 5, "diffus": [0, 1, 3, 7], "dirichlet": [1, 6], "discretis": 0, "equat": [0, 1, 3, 6, 7, 8], "exampl": 2, "explicit": 0, "exponenti": 0, "face": 0, "finit": [4, 5], "fisher": 8, "fit": 0, "flux": 1, "form": [0, 1, 6], "gener": 1, "grid": 2, "imex": 8, "implement": 4, "implicit": 0, "indic": 4, "interpol": 0, "iter": 8, "known": 1, "line": 8, "linear": [0, 7], "matrix": [0, 1, 6], "mesh": 3, "method": [0, 4, 8], "model": 3, "natur": 6, "neumann": 6, "newton": 8, "nonlinear": 8, "nonuniform": 2, "note": 4, "overview": [3, 5], "physic": 4, "poisson": 6, "python": 3, "random": 2, "reaction": [0, 1, 3], "resolv": 6, "robin": 1, "simul": 4, "solv": [3, 7, 8], "step": [7, 8], "tabl": 4, "term": 1, "theta": 0, "time": [7, 8], "transient": 1, "uniform": 2, "upwind": 0, "v": 5, "valu": [0, 1], "volum": [4, 5], "without": 1}}) \ No newline at end of file diff --git a/solving_linear_equations.html b/solving_linear_equations.html index ed3b856..e7fed8d 100644 --- a/solving_linear_equations.html +++ b/solving_linear_equations.html @@ -1,254 +1,139 @@ - - - - - + - - Solving linear equations — FVM Docs 0.1 documentation - - - - + + Solving linear equations — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - + + + +
    +
    +
    +
    + +
    +

    Solving linear equations

    +

    Provided the equation is linear, meaning that the coefficients, nor the reaction term depend on the solution variable a time-stepping approach can be used to solve the equation.

    +
    +

    Time stepping the advection-diffusion equation

    +

    Here we will rearrange the \(\theta\)-method form of the advection-diffusion equation,

    +
    +\[\frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n})\]
    +

    into the form of a linear system \(A\cdot x = d\). Firstly lets move all terms involving future time points the l.h.s,

    +
    +\[w^{n+1} - \theta \tau F(w^{n+1}) = w^{n} + (1-\theta) F(w^{n+1})\]
    +

    Replacing the \(F\) terms with the full matrix expressions and factoring yields,

    +
    +\[\underbrace{(I - \theta\tau M^{n+1})}_{A}\underbrace{w^{n+1}}_{x} = \underbrace{(I + (1-\theta)\tau M^{n})w^{n}}_{d}\]
    +

    Time-stepping involves solving this equation iteratively. First the initial conditions \(w^0\) is used to calculate the solution variable at the next point in time \(w^1\), then values of the solution variable are updated such that \(w^1\) is used to calculate \(w^2\), and so on and so forth. This procedure is shown in the Figure below,

    +
    +Time stepping the linear advection-diffusion equation. + +
    +
    +

    Note

    +

    Parameters used for this simulation. \(a=1\), \(d=1\times 10^{-3}\), \(h=5 \times 10^{-3}\), \(\tau=5 \times 10^{-4}\) with exponentially fitted discretisation. The initial and boundary conditions where, \(u(x,t_{0}) =sin(\pi x)^{100}\) with \(u(x_1)=1\) and \(u_x(x_J)=0\). The equation was time stepped for 400 iterations, the plots so the solution at times \(t=0,0.05,0.1,0.15,0.2\).

    +
    +
    +
    + + +
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    -
      + - - - -
      - -
      -
      -

      Table Of Contents

      - - -

      Previous topic

      -

      Boundary conditions for the advection-diffusion-reaction equation

      -

      Next topic

      -

      Solving nonlinear equations

      -

      This Page

      - - + -
      - - -
      -
      -
      -
      -
      - -
      -

      Solving linear equations

      -

      Provided the equation is linear, meaning that the coefficients, nor the reaction term depend on the solution variable a time-stepping approach can be used to solve the equation.

      -
      -

      Time stepping the advection-diffusion equation

      -

      Here we will rearrange the \(\theta\)-method form of the advection-diffusion equation,

      -
      -\[\begin{split}\frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n})\end{split}\]
      -

      into the form of a linear system \(A\cdot x = d\). Firstly lets move all terms involving future time points the l.h.s,

      -
      -\[w^{n+1} - \theta \tau F(w^{n+1}) = w^{n} + (1-\theta) F(w^{n+1})\]
      -

      Replacing the \(F\) terms with the full matrix expressions and factoring yields,

      -
      -\[\underbrace{(I - \theta\tau M^{n+1})}_{A}\underbrace{w^{n+1}}_{x} = \underbrace{(I + (1-\theta)\tau M^{n})w^{n}}_{d}\]
      -

      Time-stepping involves solving this equation iteratively. First the initial conditions \(w^0\) is used to calculate the solution variable at the next point in time \(w^1\), then values of the solution variable are updated such that \(w^1\) is used to calculate \(w^2\), and so on and so forth. This procedure is shown in the Figure below,

      -
      -Time stepping the linear advection-diffusion equation. -
      -
      -

      Note

      -

      Parameters used for this simulation. \(a=1\), \(d=1\times 10^{-3}\), \(h=5 \times 10^{-3}\), \(\tau=5 \times 10^{-4}\) with exponentially fitted discretisation. The initial and boundary conditions where, \(u(x,t_{0}) =sin(\pi x)^{100}\) with \(u(x_1)=1\) and \(u_x(x_J)=0\). The equation was time stepped for 400 iterations, the plots so the solution at times \(t=0,0.05,0.1,0.15,0.2\).

      -
      -
      -
      - - -
      -
      -
      -
      -
      - - -
      - - -
      -
      -
      - + -
      -
      - - -
      - © Copyright 2013, Daniel J Farrell. - Created using Sphinx 1.1.3. -
      - - +
      +

      This Page

      + +
      + + +
      +
      +
      +
    + + - - \ No newline at end of file diff --git a/solving_nonlinear_equations.html b/solving_nonlinear_equations.html index d5ce27d..9632396 100644 --- a/solving_nonlinear_equations.html +++ b/solving_nonlinear_equations.html @@ -1,314 +1,204 @@ - - - - - + - - Solving nonlinear equations — FVM Docs 0.1 documentation - - - - + + Solving nonlinear equations — FVM Docs 0.1 documentation + + - - - - - - - - + + + + - - - - - - - - + + - - -
    -
    -
    - - - - - - FVM Docs 0.1 documentation -
    - -
    -
    -
    -
    - - - -
    - - - -
    - -
    - Open Table Of Contents -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving linear equations

    -

    Next topic

    -

    The Poisson equation

    -

    This Page

    - - -
    -
    -
    - - - -
    - -
    -
    -

    Table Of Contents

    - - -

    Previous topic

    -

    Solving linear equations

    -

    Next topic

    -

    The Poisson equation

    -

    This Page

    - - -
    -
    - - -
    -
    -
    -
    -
    - -
    -

    Solving nonlinear equations

    + + + +
    +
    +
    +
    + +
    +

    Solving nonlinear equations

    -

    Note

    -

    Nonlinear equations are defined as those having coefficients which are functions of the solution variable.

    +

    Note

    +

    Nonlinear equations are defined as those having coefficients which are functions of the solution variable.

    In the previous section we discussed how to solve the linear advection-diffusion-reaction equation with method of time stepping. We used the Crank-Nicolson scheme: the implicit terms provides stability (unrestricted time step) and by also using the explicit term the accuracy of the time integration is improved to second order. However, this approach is only valid for linear equations.

    Nonlinear equations (by definition) cannot be written in linear form, as such the time-stepping approach cannot be used. It remains possible to use a fully explicit method, however in practice this is usually not done because the time step is too restrictive or the explicit form is unconditionally unstable.

    Here we discuss three approach which can be used to solve nonlinear equations:

      -
    1. IMEX (implicit/explicit) time-stepping
      2. -
    2. Newton iteration
      3. -
    3. Method of lines
      4. +

    4. IMEX (implicit/explicit) time-stepping

      5. +

    5. Newton iteration

      6. +

    6. Method of lines

    -
    -

    Fisher’s equation

    -

    We will use Fishers’ nonlinear reaction-diffusion equation as an example because it has an analytical solution,

    -
    +
    +

    Fisher’s equation

    +

    We will use Fishers’ nonlinear reaction-diffusion equation as an example because it has an analytical solution,

    +
    \[u_t = d u_{xx} + \beta u(1 - u)\]
    -

    We can directly apply the discretisation scheme already developed, but we must set \(a=0\) and note the nonlinear reaction term \(\beta u(1 - u)\). With semi-open boundary condition \(u(x_0,t)=1\) Fishers’ equation has the analytical form,

    -
    +

    We can directly apply the discretisation scheme already developed, but we must set \(a=0\) and note the nonlinear reaction term \(\beta u(1 - u)\). With semi-open boundary condition \(u(x_0,t)=1\) Fishers’ equation has the analytical form,

    +
    \[\begin{split}u(x,t) = \frac{1}{\left(1 + e^{c^{1/2}\left(x - Ut\right)} \right)^2} \\\end{split}\]
    -

    where \(U = 5\left(\frac{1}{c}\right)^{1/2}\), and, \(c=\frac{6}{d\beta}\)

    -
    -
    -

    IMEX time-stepping

    -

    The IMEX scheme treats the linear parts of the equation using implicit (or semi-implicit) scheme and the non-linear parts full explicitly. Therefore for Fishers’ equation this becomes,

    -
    +

    where \(U = 5\left(\frac{1}{c}\right)^{1/2}\), and, \(c=\frac{6}{d\beta}\)

    +
    +
    +

    IMEX time-stepping

    +

    The IMEX scheme treats the linear parts of the equation using implicit (or semi-implicit) scheme and the non-linear parts full explicitly. Therefore for Fishers’ equation this becomes,

    +
    \[w^{\prime} = \theta M^{n+1}u^{n+1} + (1-\theta) M^{n}u^{n} + \beta u^n(1-u^n)\]
    -

    This can be rearranged into the form \(A\cdot x = d\) and solved by standard time stepping techniques. First we discretised the time derivative \(u^{\prime} = \frac{u^{n+1} - u^{n}}{\tau}\), then we rearrange the equation so that future time points are on the l.h.s. and known time points are on the right,

    -
    +

    This can be rearranged into the form \(A\cdot x = d\) and solved by standard time stepping techniques. First we discretised the time derivative \(u^{\prime} = \frac{u^{n+1} - u^{n}}{\tau}\), then we rearrange the equation so that future time points are on the l.h.s. and known time points are on the right,

    +
    \[u^{n+1} - \theta \tau M^{n+1}u^{n+1} = u^{n} + (1-\theta) \tau M^{n}u^{n} + \tau\beta u^n(1-u^n)\]
    -

    Using the identity matrix, \(I\) gives the desired result,

    -
    +

    Using the identity matrix, \(I\) gives the desired result,

    +
    \[(I - \theta \tau M^{n+1})u^{n+1} = (I + (1-\theta) \tau M^{n} ) u^{n} + \tau \beta u^n(1-u^n)\]

    This can be readily solved using any linear solver, the results of this approach are shown in the Figure below. Note the poor agreement with the analytical solution for large values of times. IMEX method is known to introduce a global error that grows exponentially with iteration number. The error in the above Figure has been accentuated by deliberately choosing a fairly large time step. A small time step can be used suppress errors.

    -
    -IMEX solution of Fishers' equation. -

    Solution of Fishers’ equation with an IMEX scheme.

    -
    -
    -
    -

    Newton iteration

    -

    We have a semi-discretised system of the form \(u^{\prime}(t) = F(u(t))\) with the application of the \(\theta\)-method this gives,

    -
    +
    +IMEX solution of Fishers' equation. + +
    +

    Solution of Fishers’ equation with an IMEX scheme.

    +
    +
    +
    +
    +

    Newton iteration

    +

    We have a semi-discretised system of the form \(u^{\prime}(t) = F(u(t))\) with the application of the \(\theta\)-method this gives,

    +
    \[w^{n+1} - w^{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) = 0\]
    -

    with unknown vector \(w^{n+1}\). This equation can be solved using Newton iteration.

    -
    +

    with unknown vector \(w^{n+1}\). This equation can be solved using Newton iteration.

    +
    \[\nu^{k+1} = \nu^{k} - (I - \theta\tau A^{n})^{-1} \left( \nu^{k} - u_{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) \right)\]
    -

    where \(k\) is the iteration index (\(k\geq0\)) and \(A^{n}\) is the Jacobian matrix of \(F(w^n)\). We use the symbol \(\nu^{k}\) for iteration variables such that they are distinguished from solution of the equation at a real time point \(u^n\). The iteration needs a starting value, it is perfectly value to choose \(\nu^0 = w^n\) or for a better estimate we can precompute one iteration \(\nu^0 = w^n + \tau F(w^{n})\). Here we have described the so-called modified Newton iteration because the Jacobian is not updated during the iteration, this is known to work well for stiff PDEs.

    -

    Applying the modified Newton iteration to Fisher’s equation yields the results shown in the Figure below. The results are much better than for IMEX scheme as the error (particularly for long time points) is much lower. However, the accumulation of global error can still be observed. This is a known feature of so called single step methods. In ODE terminology a single step method is one that uses only the most recent past data to evaluate the equation at the future time step. Note even if the equations are implicit (or semi-implicit as with the \(\theta\)-method) the future state is predicted only from one historic data point. In contrast a multistep method would used many previous data points to achieve highly accurate time integration.

    -
    -IMEX solution of Fishers' equation. -

    Solution of Fishers’ equation with a modified Newton iteration.

    -
    -
    -
    -

    Method of lines

    +

    where \(k\) is the iteration index (\(k\geq0\)) and \(A^{n}\) is the Jacobian matrix of \(F(w^n)\). We use the symbol \(\nu^{k}\) for iteration variables such that they are distinguished from solution of the equation at a real time point \(u^n\). The iteration needs a starting value, it is perfectly value to choose \(\nu^0 = w^n\) or for a better estimate we can precompute one iteration \(\nu^0 = w^n + \tau F(w^{n})\). Here we have described the so-called modified Newton iteration because the Jacobian is not updated during the iteration, this is known to work well for stiff PDEs.

    +

    Applying the modified Newton iteration to Fisher’s equation yields the results shown in the Figure below. The results are much better than for IMEX scheme as the error (particularly for long time points) is much lower. However, the accumulation of global error can still be observed. This is a known feature of so called single step methods. In ODE terminology a single step method is one that uses only the most recent past data to evaluate the equation at the future time step. Note even if the equations are implicit (or semi-implicit as with the \(\theta\)-method) the future state is predicted only from one historic data point. In contrast a multistep method would used many previous data points to achieve highly accurate time integration.

    +
    +IMEX solution of Fishers' equation. + +
    +

    Solution of Fishers’ equation with a modified Newton iteration.

    +
    +
    +
    +
    +

    Method of lines

    -

    Note

    -

    For an excellent introduction to the Method of Lines see the article at scholarpedia.org

    +

    Note

    +

    For an excellent introduction to the Method of Lines see the article at scholarpedia.org

    The method of lines is not really a method, it is a way of writing PDEs in such a way that they can be solved by black box ODE solvers. The time integration of ODEs is a mature topic. Robust and well tested solvers are readily available on the internet and integrated into many open-source projects for easy use. The method of lines is a way of letting PDE equation take advantage of the mature nature of ODE solvers.

    -

    We have a semi-discretised system of the form, \(u^{\prime}(t) = F(u(t))\) to use the method of lines we leave the PDE is the semi-discristed form such that is has a continuous time derivative on the l.h.s., i.e. \(u^{\prime}(t) = \frac{du(t)}{dt}\). In the semi-discretised form what we really have is no longer a PDE but a system of coupled ODEs. Moreover, a PDE equation is a differential equation which contains differential terms with respect to more that one variable. In this semi-discretised form we have replaced the spatial differential operator with a matrix so there is only one derivative. The equation has become a system of ODEs with number of equation proportional to the number of grid points in our discretisation.

    +

    We have a semi-discretised system of the form, \(u^{\prime}(t) = F(u(t))\) to use the method of lines we leave the PDE is the semi-discristed form such that is has a continuous time derivative on the l.h.s., i.e. \(u^{\prime}(t) = \frac{du(t)}{dt}\). In the semi-discretised form what we really have is no longer a PDE but a system of coupled ODEs. Moreover, a PDE equation is a differential equation which contains differential terms with respect to more that one variable. In this semi-discretised form we have replaced the spatial differential operator with a matrix so there is only one derivative. The equation has become a system of ODEs with number of equation proportional to the number of grid points in our discretisation.

    -

    Note

    -

    The discretised boundary conditions are expected to be included in \(F(u(t))\) such that the problem is well posed.

    +

    Note

    +

    The discretised boundary conditions are expected to be included in \(F(u(t))\) such that the problem is well posed.

    -

    There are a wide variety of ODE solvers we are interested in the implicit type as we know many PDE problems are unstable in explicit form. The result in the Figure below have been calculated using an implicit multistep Adams-Moulton method (the algorithm used is the highly popular vode solver from netlib.org, although we used actually used a Python wrapper from scipy to perform the calculations shown). This can be thought of as a generalised of the \(\theta\)-method to include a multiple number of future time steps which are solved simultaneously with a Newtom iteration. Unlike single step method, multistep methods do not show a global error accumulation, so the time integration is a much better approximation to analytical solution. Fisher’s equation is not particularly stiff so this approach worked well. However for stiff equations a technique called Backward differentiation formulas (BDF) is used. The vode solver has the ability to dynamically adapt between stiff and non stiff mode (i.e. Adams-Moulton to BDF) which is one of the reasons for it’s popularity.

    +

    There are a wide variety of ODE solvers we are interested in the implicit type as we know many PDE problems are unstable in explicit form. The result in the Figure below have been calculated using an implicit multistep Adams-Moulton method (the algorithm used is the highly popular vode solver from netlib.org, although we used actually used a Python wrapper from scipy to perform the calculations shown). This can be thought of as a generalised of the \(\theta\)-method to include a multiple number of future time steps which are solved simultaneously with a Newtom iteration. Unlike single step method, multistep methods do not show a global error accumulation, so the time integration is a much better approximation to analytical solution. Fisher’s equation is not particularly stiff so this approach worked well. However for stiff equations a technique called Backward differentiation formulas (BDF) is used. The vode solver has the ability to dynamically adapt between stiff and non stiff mode (i.e. Adams-Moulton to BDF) which is one of the reasons for it’s popularity.

    The method of lines is a powerful technique for solving semi-discretised PDE problems as much of the details of numerically solving the equation can be off loaded to mature and robust solvers. I do not yet have enough experience with the technique to point out short coming or failings. But it seems that provided the spatial discretisation is stable the details of the time integration can be almost ignored!

    -
    -MOL solution of Fishers' equation. -

    Solution of Fishers’ equation by the method of lines.

    -
    -
    -
    - - -
    -
    -
    +
    +MOL solution of Fishers' equation. + +
    +

    Solution of Fishers’ equation by the method of lines.

    +
    +
    + + + + +
    - - -
    +
    + + + - - \ No newline at end of file diff --git a/source/advection_diffusion.rst b/source/advection_diffusion.rst index 4b9b7ae..a7e741e 100644 --- a/source/advection_diffusion.rst +++ b/source/advection_diffusion.rst @@ -141,7 +141,7 @@ Discretised equation in matrix form Dropping the spatial subscripts and writing in vector form the equations becomes, .. math:: - \frac{w^{n+1} - w^{n}}{\tau} & = \theta F(w^{n+1}) + (1-\theta)F(w^{n}) + \frac{w^{n+1} - w^{n}}{\tau} = \theta F(w^{n+1}) + (1-\theta)F(w^{n}) where, diff --git a/source/conf.py b/source/conf.py index 7d848c7..fbe55f2 100644 --- a/source/conf.py +++ b/source/conf.py @@ -25,8 +25,7 @@ # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.mathjax', - 'sphinxjp.themecore',] +extensions = ['sphinx.ext.mathjax'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -42,7 +41,7 @@ # General information about the project. project = u'FVM Docs' -copyright = u'2013, Daniel J Farrell' +copyright = u'2013-2024, Daniel J Farrell' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -95,8 +94,8 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. -#html_theme = 'default' -html_theme = 'basicstrap' +html_theme = 'default' +#html_theme = 'basicstrap' #html_theme = 'dotted' #html_theme = 'trstyle' #html_theme = 'solarized'