8. Conservation & mixed method#
\(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h}\)
A single boundary value problem can admit different weak formulations, possibly set in different spaces. For the Poisson equation, the primal weak form can be approximated by the Lagrange finite element method, while the dual weak form or the mixed weak form can be approximated by the pair of Raviart-Thomas and DG finite element spaces. This notebook is on the latter.
From the practical standpoint, one of the primary reasons for using this latter mixed method for the Poisson equation is conservation, a discrete structure-preservation property, which we define in this notebook. We will use the Raviart-Thomas elements we discussed earlier to implement and study a conservative mixed method in this notebook.
8.1. Two typical examples#
Many boundary value problems arise by combining physical conservation laws with empirical constitutive laws. Methods that preserve the conservation law discretely in some sense are referred to as conservative methods.
Here are two examples of how conservation laws and constitutive laws combine to produce boundary value problems.
8.1.1. Example A: Porous media flow#
For slow viscous flow through a permeable porous medium, the Darcy law is an empirical statement connecting the flux of the fluid \(q\) to the pressure (\(p\)) gradient by
where \(K\) is the permeability tensor of the rock or the porous formation. Conservation of mass then implies
where \(f\) represents injection sources, wells, or sinks.
Historically, Example A was studied extensively by the fossil fuel industry. However, it also finds ecological uses in groundwater remediation. And, perhaps even more fittingly for us here in Portland, in modeling the percolation of hot water through coffee grounds!
8.1.2. Example B: Diffusion of heat#
Consider the stationary heat dissipation example from a prior notebook. Recall that Fourier’s law of heat conduction is a constitutive law that states (in the absence of advective and radiative effects) that the flux of heat energy \(q\) (sometimes also called heat flux density) through a material is related to the gradient of the temperature \(T\) by
where \(\kappa\) is the measured thermal conductivity of the material. Conservation of energy implies that
where \(f\) represents the heat source.
Note that in both examples, elimination of flux results in an equation of the form
(with \(u\) equal to either \(T\) or \(p\), and \(\alpha\) equal to either \(\kappa\) or \(K\)), which we have treated previously using Lagrange elements. In this notebook, rather than eliminating \(q\), our focus is on approximating \(q\).
8.2. Conservation in the finite element context#
Suppose we use a mesh \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} \oh\) of a \(d\)-dimensional domain \(\om\) to compute a discrete approximation \(q_h\) to the exact flux \(q\). Let \(P_p(\oh) = \prod_{K \in \oh} P_p(K)\) denote the DG space, so piecewise polynomial vector fields of degree at most \(p\) are in \(P_p(\oh)^d.\)
Definition: We say that a vector field \(q_h\) in \(P_p(\oh)^d\) is a conservative flux approximation of an exact flux \(q\) if the following two conditions hold:
Condition (1): on all interior mesh facets, ⟦ \(q_h\cdot n\) ⟧ \( = 0\), and
Condition (2): on every mesh element \(K \in \oh\),
Condition (1) requires the normal flux to be single-valued. There we have used the notation ⟦ \(q_h\cdot n\) ⟧ for the jump of the normal component of the flux \(q_h\).
Condition (2) is an element-wise conservation condition. Note that by the divergence theorem the right hand side of the equality in Condition (2) also equals
with \(f = \dive q.\)
Also note that both Conditions (1) and (2) refer only to the values of \(q_h\) on the facets. (Indeed, there are methods that produce flux approximations only along the mesh facets, and the same definition can be used to decide if such fluxes are conservative or not. The mixed method however will produce fluxes on the whole domain.)
To understand why Conditions (1) and (2) are interesting properties for an approximation to have, consider one of the above-mentioned applications, say the diffusion of heat, and recall how the conservation equation arises. In physics, conservation of energy is often written in integral form over a control volume \(S\) as
i.e., the flux of heat leaving the subdomain \(S\) through its boundary (equalling the left hand side above) must equal the amount of heat produced by sources \(f\) within \(S\) (equal to the right hand side above). Because this conservation statement in integral form holds for any control volume \(S\), the divergence theorem implies the equivalent differential equation \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} \dive q = f\), which is one of the conservation equations stated previously.
Now, consider how we may obtain a discrete version of the integral form of the conservation law. Instead of arbitrary control volumes, we restrict ourselves to discrete control volumes, which are the unions of the elements used in the computation. Namely, consider any subset of mesh elements \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} O_h \subseteq \oh\) and put \(S_h = \cup\{ K: K \in O_h\}\). If the discrete flux out of the subdomain \(S_h\) satisfies
then it makes sense to declare that \(q_h\) satisfies a discrete conservation law, since equation (DC) is exactly the discrete analog of the exact integral form of the conservation law, namely the equation \( \int_{\partial S} q \cdot n = \int_S f \) discussed previously. The only difference is that while the exact form holds on any \(S\), equation (DC) only holds for the discrete control volumes \(S_h\).
One can immediately verify that the above equation (DC) is a direct consequence of Conditions (1) and (2) in our definition. Note how we need both Conditions (1) and (2) to accomplish the interelement cancellations of influx and outflux within \(S_h\).
One of the modern themes in numerical solution of partial differential equations, called structure preservation, studies questions like this: how can we construct methods that (not only converge optimally, but also) yield solutions that preserve certain critical features or structures of the exact solution? In the context of our model problem, methods that produce a conservative flux \(q_h\), as defined above, are structure-preserving, the preserved structure being a conservation law.
Another way of thinking about conservation is in terms of superconvergence of certain functionals. (The phenomena where the error, or some functional of the error, goes to zero at an unexpectedly high speed is called superconvergence. We have seen a superconvergence example in a prior notebook in a completely different context.) Condition (2) in the definition of conservation, \( \int_{\partial K} q_h \cdot n = \int_{\partial K} q\cdot n, \) can equivalently be written using the functional
as
For good methods, we expect the error \(q_h - q\) to go to zero as the mesh size \(h \to 0\), so we expect \(G_K(q) - G_K(q_h) \to 0\). But it is exceptional to get zero for a functional of the error. For a conservative method, that exceptional property holds and \(G_K(q) - G_K(q_h) = 0\).
8.3. A simple model problem#
Let us solve a specific example built from Example B. We want to simulate a steady heat flux \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} q: \om \to \R^2\) on a rectangular bar-shaped domain \(\om\) of length 6 units and width 2 units, divided into equal left and right halves. More specifics of the problem appear after we draw the geometry.
import ngsolve as ng
from ngsolve.webgui import Draw
from netgen.occ import X, Y, MoveTo, Rectangle, Glue, OCCGeometry
rgtbar = Rectangle(3, 2).Face()
lftbar = MoveTo(-3, 0).Rectangle(3, 2).Face()
lftbar.faces.name = 'lftbar'
rgtbar.faces.name = 'rgtbar'
lftbar.edges.Min(X).name = 'lft'
lftbar.edges.Min(Y).name = 'bot'
lftbar.edges.Max(Y).name = 'top'
lftbar.edges.Max(X).name = 'mid'
rgtbar.edges.Max(X).name = 'rgt'
rgtbar.edges.Min(Y).name = 'bot'
rgtbar.edges.Max(Y).name = 'top'
bar = Glue([lftbar, rgtbar])
Draw(bar);
geo = OCCGeometry(bar, dim=2)
mesh = ng.Mesh(geo.GenerateMesh(maxh=1/4))
Here are the further specifications of the thermal problem:
Materials whose isotropic thermal conductivity (\(\kappa\)) values equal 1 and 10 occupy the left and right halves, respectively.
The left boundary edge of \(\Omega\) is kept at temperature \(T=1\) while the right edge is kept at \(T=1/10\).
The top and bottom boundary edges are perfectly insulated, i.e, the outward flux \(q \cdot n\) vanishes there.
The bar is heated by a source which is modeled by $\( f(x, y) = 5 \exp(-10 [(x/5)^2 + (y-1)^2]). \)$
Here is a plot of \(f\):
from ngsolve import x, y, exp
f = 5 * exp(-10*( (x/5)**2 + (y-1)**2))
Draw(f, mesh);
Clearly, this heat source appears centered in the domain. However, the material is not homogeneous as seen from the thermal conductivity plot next:
kappa = mesh.MaterialCF({'lftbar': 1, 'rgtbar': 10})
Draw(kappa, mesh);
Finally, the other piece of given information on boundary conditions tells us that heat may enter the bar through the left end. Here is a visualization of the boundary data extension:
Tbdr = mesh.BoundaryCF({'lft': 1, 'rgt': 0.1})
Tboundary = ng.GridFunction(ng.H1(mesh)); Tboundary.Set(Tbdr, definedon=mesh.Boundaries('lft|rgt'));
Draw(Tboundary);
Given all the above bits of information, the problem is to compute the temperature \(T\) and the heat flux \(q\). We will solve it momentarily.
Questions for discussion:
Do you have any physical guesses on which side will get heated up more in steady state?
What features in the problem generate heat and which ones dissipate heat?
Will the centered source create a steady hot spot in the center?
8.4. The Raviart-Thomas mixed method#
Consider the boundary value problem of finding \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} u: \om\to \R\) and \(q: \om \to \R^d\) (\(d\ge 2\)), given \(f: \om \to \R\) and (pointwise) invertible \(\kappa : \om \to \R^{d \times d}\), satisfying
with Dirichlet boundary conditions \(u =0 \) on \(\partial \om\). From the lectures, we know that the mixed weak formulation of this problem reads as follows.
Find \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} u \in L^2(\om)\) and \(q \in H(\dive, \om)\) satisfying
\[\begin{split} \begin{aligned} \int_\om \kappa^{-1} q \cdot r - \int_\om u\, \dive r & = 0, && \text{ for all } r \in H(\dive, \om), \text{ and } \\ \int_\om v\, \dive q & = \int_\om f \,v, && \text{ for all } v \in L^2(\om). \end{aligned} \end{split}\]
To obtain a numerical method, starting from the above weak form, we use the DG space \(P_p(\oh)\) in place of \(L^2(\om)\) and the Raviart-Thomas finite element space (introduced earlier)
in place of \(H(\dive, \om).\)
The Raviart-Thomas (RT) mixed method of order \(p\) finds \(u_h \in P_p(\oh)\) and \(q_h \in R_{h, p}\) satisfying
\[\begin{split}\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} \begin{aligned} \int_\om \kappa^{-1} q_h \cdot r_h - \int_\om u_h\, \dive r_h & = 0, && \text{ for all } r_h \in R_{h, p}, \text{ and } \\ \int_\om v_h\, \dive q_h & = \int_\om f \,v_h, && \text{ for all } v_h \in P_p(\oh). \end{aligned} \end{split}\]
In the latter equation, if we select a test function \(\newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} v_h\) that is the indicator function of one element \(K \in \oh\), which is of course contained in \(P_p(\oh)\) for any \(p \ge 0\), then we obtain
which is equivalent to Condition (2) in our definition of a conservative flux. Of course, Condition (1) is automatically satisfied by \(q_h\), since every \(q_h \in R_{h, p}\) has ⟦ \(q_h\cdot n\) ⟧ \( = 0\). This proves that the Raviart-Thomas mixed method yields conservative fluxes.
Next, we proceed to implement this method in NGSolve. There are two ways to assemble such a mixed weak form, as we shall see in the next section.
8.5. Two ways to assemble#
8.5.1. The first way: assemble a composite form#
We can think of the system using a “big” bilinear form \(C(\cdot, \cdot)\) in the product space \( \newcommand{\R}{\mathbb{R}} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} R_{h, p} \times P_p(\oh)\), i.e., for any \((q_h, u_h), (r_h, v_h) \in R_{h, p} \times P_p(\oh)\), set
Then the previous set of two equations for the Raviart-Thomas mixed method can be described as the single equation using the composite bilinear form:
To compute with this form, we need the product space
\(\newcommand{\R}{\mathbb{R}}
\newcommand{\dive}{\mathop{\mathrm{div}}}
\newcommand{\grad}{\mathop{\mathrm{grad}}}
\newcommand{\om}{\varOmega}
\newcommand{\oh}{\varOmega_h}
R_{h, p} \times P_p(\oh)\). This space is represented by the object RW introduced below.
p = 4
R = ng.HDiv(mesh, order=0, RT=True)
W = ng.L2(mesh, order=0)
RW = R * W # Alternately: RW = ng.FESpace([R, W])
Note that trial and test functions in RW have two components. Otherwise the procedure for declaring the bilinear form is exactly the same as what we saw previously.
from ngsolve import div, dx
q, u = RW.TrialFunction()
r, v = RW.TestFunction()
C = ng.BilinearForm(RW)
C += ((1/kappa) * q*r - u*div(r) - div(q)*v) * dx
C.Assemble(); type(C.mat)
ngsolve.la.SparseMatrixd
We now have a sparse assembled matrix representing the matrix of the big bilinear form \(C(\cdot, \cdot)\) on the product finite element space.
8.5.2. The second way: assemble a block system#
In the second approach, we make individual components of the bilinear form and put them together, i.e., we make the following “small” bilinear forms:
Then the RT mixed-method equations can be rewritten using these forms as
Note that the bilinear form \(b\) uses different trial and test spaces. Implementing this requires the use of keyword arguments trialspace and testspace that we have not seen previously, but is self-explanatory. There is no need to define a product space in this approach.
q, r = R.TnT()
u, v = W.TnT()
b = ng.BilinearForm(trialspace=R, testspace=W)
b += -div(q)*v * dx
a = ng.BilinearForm(R)
a += (1/kappa) * q*r * dx
b.Assemble(); a.Assemble();
When working with such component forms, you also need a facility to place them into appropriate places of a block partitioned matrix. Let \(\mathtt A\) and \(\mathtt B\) denote the stiffness matrices of the forms \(a(\cdot, \cdot)\) and \(b(\cdot, \cdot)\) made in the code block above, respectively. Then, what we want is the matrix
where \({}^T\) denotes the transpose.
BB = ng.BlockMatrix([[a.mat, b.mat.T],
[b.mat, None ]])
type(BB)
ngsolve.la.BlockMatrix
This BB object represents the above-mentioned block matrix \(\mathbb B\).
Note that the None element in the code above is how we represent the zero block sparsely.
Obviously \(\mathbb B\) and \(C\) represent the same linear operator. Block structures are just a convenience feature to work with components (without replicating the already allocated memory of the components). Block matrices expect block vectors when asked to perform matrix multiplication.
Be warned that it’s not a good idea to directly perform a generic vector operation between a block vector object and a regular NGSolve vector. In some cases NGSolve will raise an exception and will not allow you to do so. Also, when you work with block vectors, it is a good idea to ensure that you have variable names bound to its individual component vectors in your python workspace. Here is an exercise to get some practice on working with the BlockMatrix structures and with matrix-vector products in NGSolve.
Question for discussion:
How would you write Python code to verify that the application of the block matrix \(\mathbb{B}\) and the matrix \(C\) of the composite bilinear form to the same vector gives the same result. (Suggestion: Look up the
componentattribute of an NGSolve grid function on a product space.)
8.6. Solving for the thermal flux#
We now return to solve the model problem stated earlier. To repeat the problem, translating it to a boundary value problem, we need to solve for temperature \(T\) satisfying
To solve this using the mixed method, we have to carefully incorporate the boundary conditions into the mixed formulation. Dirichlet and Neumann conditions are respectively imposed as natural and essential conditions in the mixed formulation, in exactly the opposite manner of the primal formulation. We use the subspace of the Raviart-Thomas space where the flux boundary conditions are incorporated essentially:
Also incorporating the Dirichlet boundary conditions on \(T\) naturally, we obtain the following equations of the method: find \(T_h \in P_p(\oh)\) and \(q_h \in \Rtb\) satisfying
8.6.1. Block linear system#
Adopting the second approach to assembly mentioned above, we produce a block linear system, whose left- and right-hand sides are returned by the next function. Also note that to make the right-hand side of the first equation, we need to integrate the trace of \(r_h \cdot n\) along boundary edges, which is accomplished below by r.Trace() * n and ds provided by NGSolve.
from ngsolve import dx, ds
def MakeMixed(mesh, f, Tbdr, kappa, p=4):
"""Return the block linear system (rhs and lhs) of the
RT mixed method for solving Example T"""
R = ng.HDiv(mesh, order=p, RT=True, dirichlet='top|bot')
W = ng.L2(mesh, order=p)
q, r = R.TnT()
u, v = W.TnT()
dl = ds(definedon=mesh.Boundaries('lft|rgt'))
n = ng.specialcf.normal(mesh.dim)
b = ng.BilinearForm(trialspace=R, testspace=W)
b += -div(q) * v * dx
a = ng.BilinearForm(R)
a += (1/kappa) * q * r * dx
fR = ng.LinearForm(R)
fR += -Tbdr * (r.Trace() * n) * dl
fW = ng.LinearForm(W)
fW += -f * v * dx
with ng.TaskManager():
b.Assemble()
a.Assemble()
fR.Assemble()
fW.Assemble()
BB = ng.BlockMatrix([[a.mat, b.mat.T], [b.mat, None]])
FF = ng.BlockVector([fR.vec, fW.vec])
return BB, FF, R, W
BB, FF, R, W = MakeMixed(mesh, f, Tbdr, kappa)
Next, we need to solve for a vector x satisfying the block system BB * x = FF. In previous studies, we assembled sparse matrix objects which we could hand over to a sparse solver (using the Inverse(...) method). However, the block matrix BB does not have a sparse Inverse method. Hence let us take a moment to consider an iterative solver.
8.6.1.1. Iterative solver#
Iterative solvers built using Krylov spaces do not generally require an assembled matrix, only an object that can perform the associated matrix-vector multiplication. The block matrix BB can perform the multiplication BB * x. Although Krylov space iterative solvers can be easily implemented in a few lines of Python code, let us make use of a ready-made Krylov solver implementation that NGSolve provides. MINRES is an iterative method suitable for invertible symmetric (not necessarily positive definite) linear systems, like the block system we have produced. Its implementation is available in ng.solvers.MinRes.
The rate of convergence of iterative solvers like MINRES depends on the condition number of the system. Hence they are usually effective only when we also provide it with a good preconditioner. Recall that a preconditioner for use in solving \(\mathbb B x = b\) is a linear operator \(\mathbb P\) (acting on vectors of the same size as \(\mathbb B\)) with the property that \(\mathbb P \mathbb B\) is better conditioned than \(\mathbb B\). Once a preconditioner is given, the iterative solver, instead of solving \(\mathbb B x = b\), solves the equivalent but better conditioned system \(\mathbb P\mathbb B x = \mathbb P b\) behind the scenes.
Setting \(\mathbb P\) to the identity operator on the free degrees of freedom is tantamount to no preconditioning. To witness the performance of the solver without preconditioning, let us start with that option for \(\mathbb P\). Here is how we set the needed block identity operator.
# project to range=true/false bits of mask
identityR = ng.Projector(mask=R.FreeDofs(), range=True)
identityW = ng.IdentityMatrix(W.ndof)
PI = ng.BlockMatrix([[identityR, None],
[None, identityW]])
qT = ng.BlockVector([ng.GridFunction(R).vec,
ng.GridFunction(W).vec])
ng.solvers.MinRes(mat=BB, pre=PI, rhs=FF, sol=qT, maxsteps=30);
LinearSolver iteration 1, residual = 2.9349850200882894
LinearSolver iteration 2, residual = 2.8525538570019844
LinearSolver iteration 3, residual = 2.294830356474514
LinearSolver iteration 4, residual = 2.2023632468343575
LinearSolver iteration 5, residual = 1.7265143891534573
LinearSolver iteration 6, residual = 1.591797669476658
LinearSolver iteration 7, residual = 1.406188933386344
LinearSolver iteration 8, residual = 1.3282590813510842
LinearSolver iteration 9, residual = 1.179691970585117
LinearSolver iteration 10, residual = 1.1020337549793986
LinearSolver iteration 11, residual = 0.995812038645757
LinearSolver iteration 12, residual = 0.9656993370434038
LinearSolver iteration 13, residual = 0.8865937134575593
LinearSolver iteration 14, residual = 0.8685093304602324
LinearSolver iteration 15, residual = 0.7996103532005372
LinearSolver iteration 16, residual = 0.7911405926295936
LinearSolver iteration 17, residual = 0.7396760073834668
LinearSolver iteration 18, residual = 0.7329865091545326
LinearSolver iteration 19, residual = 0.6930839079783313
LinearSolver iteration 20, residual = 0.6865590523269943
LinearSolver iteration 21, residual = 0.6575316805574284
LinearSolver iteration 22, residual = 0.6522408925896295
LinearSolver iteration 23, residual = 0.6323027561275425
LinearSolver iteration 24, residual = 0.6272771825433502
LinearSolver iteration 25, residual = 0.612805795049337
LinearSolver iteration 26, residual = 0.6070239828105076
LinearSolver iteration 27, residual = 0.5960285341344539
LinearSolver iteration 28, residual = 0.5888856793115653
LinearSolver iteration 29, residual = 0.5792974712781549
LinearSolver iteration 30, residual = 0.5731299473401457
WARNING: LinearSolver did not converge to TOL
As you can see the convergence is abysmally slow.
Question for discussion:
Can you get these iterations to reduce the error further? (Increase
maxstepsand run again to see how many iterations you might need to do before getting anywhere close to convergence.)
8.6.1.2. Preconditioner#
We need a better preconditioner. For our block matrix
one technique is to construct a preconditioner in the following form, using the same block partitioning,
where \(\mathtt M_R\) and \(\mathtt M_W\) are the stiffness matrices of the bilinear forms \( \int_\om \kappa^{-1} q\cdot r + \int_\om (\dive q ) (\dive r),\) for \(q, r \in \Rtb\) and \( \int_\om u \, v,\) for \(u, v \in P_p(\oh),\) respectively. Since \(\mathtt M_W\) is block diagonal with one block for each element (why?), it is very easy to invert. So the cost of construction of this preconditioner is essentially the cost of inversion of \(\mathtt M_R\), a smaller matrix than the original \(\mathbb B\).
def MakePrec(R, W):
q, r = R.TnT()
mR = ng.BilinearForm(((1/kappa)*q*r + div(q)*div(r))*dx)
mR.Assemble()
P = ng.BlockMatrix([[mR.mat.Inverse(R.FreeDofs()), None],
[None, W.Mass(1).Inverse()]])
return P
PP = MakePrec(R, W)
qh = ng.GridFunction(R)
Th = ng.GridFunction(W)
qT = ng.BlockVector([qh.vec, Th.vec])
ng.solvers.MinRes(mat=BB, pre=PP, rhs=FF, sol=qT, maxsteps=15);
LinearSolver iteration 1, residual = 4.655166599019508
LinearSolver iteration 2, residual = 4.650238209773678
LinearSolver iteration 3, residual = 1.6631181609423022
LinearSolver iteration 4, residual = 0.3060735234476721
LinearSolver iteration 5, residual = 0.03971097203811199
LinearSolver iteration 6, residual = 0.0020292060009608147
LinearSolver iteration 7, residual = 0.00012449831832854423
LinearSolver iteration 8, residual = 3.397171011593384e-06
LinearSolver iteration 9, residual = 7.718494250916117e-08
Question for discussion:
Why does this preconditioner work well?
8.6.1.3. The temperature and heat flux#
MinRes overwrites the solution into the same memory location as the input qT. It is in the form of a block vector, whose components occupy the same memory as qh and Th. So to plot the temperature, we can directly use the variable name Th which is bound to the memory block of the \(T\)-component.
Draw(Th, deformation=True, settings={"Colormap":{"autoscale": True, "ncolors": 20}, "camera":{"transformations" :[{"type": "rotateX", "angle": -50}, {"type": "rotateY", "angle": -20}, {"type": "rotateZ", "angle": -10}]}});
This plot of the calculated temperature \(T\) reveals lower temperature variations on the right than on the left. You also see that there is a kink at the middle interface (indicating a discontinuity in \(\grad T\)) to accommodate a flatter \(T\) profile on the right.
Next, we plot the heat flux vector field:
Draw(qh, vectors={'grid_size': 30}, settings={"Colormap":{"autoscale": True, "ncolors": 20}, "Light": {"ambient":0.1, "diffuse":0.7, "shininess": 10,"specularity": 0.9}});
It is clear from this plot that the heat flow on the right subdomain is more than on the left (as the colors represent the magnitude of the heat flux). In fact, from the direction of the fluxes, we see that the left boundary condition is helping to dissipate some of the heat generated by the source. Also observe the apparent continuity of the normal component of the heat flux despite the jump in the material properties. This is exactly as it should be for a conservative method, as otherwise nonphysical energy creation or energy loss will be predicted across the interface.
Let us also note that even though the blob of heat source was centered in the domain, it created an uncentered temperature peak to the left. Is this shift of the peak due to the difference in boundary conditions on the left and right, or due to the difference in the left and right thermal conductivity? The utility of having a simulation tool is that we can now answer questions like this easily even when the answers are not intuitive.
Question for discussion:
Is there a primary cause for the shift of the temperature peak? Could it be primarily due to the effect of the boundary conditions? Or is it because the right half of the domain, due to its higher thermal conductivity, can transmit heat through it more efficiently? (Suggestion: Experiment further, changing the boundary conditions, and changing the thermal conductivities.)
Questions for discussion, on comparison with the primal formulation:
How will you solve the same model problem using Lagrange finite elements applied to the primal weak form? In such an approach you will compute an approximate temperature, which we denote by \(\tilde T_h\).
Using your computed temperature \(\tilde T_h\) (in the Lagrange space \( V_{h, p}\)) in your approach above, how will you compute an element-wise approximate flux \(\tilde{q}_h = -\kappa \grad \tilde{T}_h\)? Can this flux be conservative?
Can you quantify the spurious heat energy loss across \(\Gamma\) in such a method? Specifically, letting \(\Gamma\) denote the middle interface, how will you compute \(\int_\Gamma\) ⟦ \(\tilde q_h\cdot n\) ⟧ ? (Suggestion: To integrate multivalued functions on element interfaces,
dx(element_boundary=True)can be useful.)
Question for discussion, on use of other spaces in the mixed formulation:
What happens if, in the mixed formulation, you replace the RT space \(\Rtb\) with its subspace of Cartesian product of Lagrange spaces of degree at most \(p\), say for \(p=1\)?
8.7. Summary#
We studied conservation as a discrete structure-preservation property.
Conservative fluxes were defined in the finite element context.
Conservation was seen as a superconvergence property.
The Raviart-Thomas mixed method produces conservative fluxes by design.
One way to implement a mixed method is via a composite bilinear form on a product space.
Another way is to make component bilinear forms and put them together in a block matrix.
Iterative solvers like MINRES can be used to solve the block mixed systems.
A block preconditioner using smaller component inverses was effective.
Conservative fluxes do not produce spurious energy loss across interfaces.
The choice of spaces in the mixed formulation is crucial.