Partial Differential Equations

Storing various PDEs that can be will be solved in this course. This includes:

  • Diffusive 1D

    \[u_{t} - \mu u_{xx} = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]
  • Advective 1D

    \[u_{t} + c {u}_{x} = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]
  • Diffusive-Advective 1D

    \[u_{t} + c {u}_{x} - \mu u_{xx} = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]
  • Poisson in 1D

    \[u_{xx} = f(x) \qquad \forall \, x \in \Omega = [0, L]\]

The goal is to implement the code in python and not rely on existing methods.

Bram Lagerweij COHMAS Mechanical Engineering KAUST 2021

pde.advective(x, connect, c, num_q, order)

Time derivative of the PDE for advective diffusive problems.

\[u_{t} + c u_x = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]

Which, is converted into a weak form through:

\[\int_\Omega (\tilde{u}_t + c \tilde{u}_x ) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Where \(\tilde{u}\) is our approximation:

\[\tilde{u}(x) = \sum_{n = 1}^N \bar{u}_n \phi_n(x)\]

in which \(\bar{u}_n\) are the unknowns, socalled degrees of freedom and \(\phi_n(x)\) are the basisfunctions of FE approximation space \(V_h\). Substituting this approximation leads to:

\[\int_\Omega (\partial_t\sum_{j = 1}^N \bar{u}_j \phi_j(x) + c \partial_x \sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Which we split into different integrals:

\[\begin{split}&\int_\Omega \tilde{u}_t dV = \int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV = \sum_{j = 1}^N \partial_t \bar{u}_j \int_\Omega \phi_j(x) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ &\qquad \Rightarrow \qquad M \bar{u}\end{split}\]

For the first integral we notice that the basis functions are constant through time, only the degrees of freedom \(\bar{u}_j\) vary through time. Similarly these degrees of freedom does not affect the integral over space \(\int_\Omega dV\). Thus we can write:

\[\begin{split}\int_\Omega \tilde{u}_t dV & = \int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ & = \sum_{j = 1}^N \partial_t \bar{u}_j \int_\Omega \phi_j(x) \phi_i(x) dV \\ &\qquad \Rightarrow \qquad M \bar{u}\end{split}\]

where \(M\) is the mass matrix which combines the integral for all different basis functions. For the second term we aknoledge that the degrees of freedom have no spatial and temporal effects thus we can take them out of the integral and derivatives.

\[\begin{split}\int_\Omega u_x dV & = \int_\Omega (\partial_x\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ & = \sum_{j = 1}^N \bar{u}_j \int_\Omega \partial_x\phi_j(x) \phi_i(x) dV \\ & \Rightarrow \qquad T \bar{u}\end{split}\]

where \(T\) is the socalled transport matrix, wich only depends on the basis functions.

Now we can write our PDE in terms of linear algabra objects:

\[M \bar{u}_t + c T \bar{u} = 0\]

which we modify to be in the format is expected by the temporal solvers:

\[M \bar{u}_t = K \bar{u}\]
Parameters
  • x (array_like(float)) – Global coordinates of all degrees of freedom.

  • connect (array_like(int), shape((num_ele, dofe/ele))) – Element to degree of freedom connectivety map.

  • c (float) – Advective constant.

  • num_q (int) – Number of Gausian quadrature points.

  • order (int) – Order of the polynomial used by our element.

Returns

  • M (matrix, (sparse csr format)) – The mass matrix.

  • K (matrix, (sparse csr format)) – The combination of stiffness and transport matrix matrix scaled with the approprate constants \(K = - c T\).

  • b (vector, (dense array)) – The right hand side, because we consider a homogeneous PDE with diriclet conditions it is a zero vector.

pde.advectivediffusive(mesh, c, mu)

Time derivative of the PDE for advective diffusive problems.

\[u_{t} + c u_x - \mu u_{xx} = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]

Which, is converted into a weak form through:

\[\int_\Omega (\tilde{u}_t + c \tilde{u}_x - \mu \tilde{u}_{xx} ) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Where \(\tilde{u}\) is our approximation:

\[\tilde{u}(x) = \sum_{n = 1}^N \bar{u}_n \phi_n(x)\]

in which \(\bar{u}_n\) are the unknowns, socalled degrees of freedom and \(\phi_n(x)\) are the basisfunctions of FE approximation space \(V_h\). Substituting this approximation leads to:

\[\int_\Omega (\partial_t\sum_{j = 1}^N \bar{u}_j \phi_j(x) + c \partial_x \sum_{j = 1}^N \bar{u}_j \phi_j(x) - \mu \partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x) ) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Which we split into different integrals:

\[\int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x))\phi_i(x) dV + \int_\Omega (c\partial_x\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV - \int_\Omega (\mu\partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

For the first integral we notice that the basis functions are constant through time, only the degrees of freedom \(\bar{u}_j\) vary through time. Similarly these degrees of freedom does not affect the integral over space \(\int_\Omega dV\). Thus we can write:

\[\begin{split}\int_\Omega \tilde{u}_t dV & = \int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ & = \sum_{j = 1}^N \partial_t \bar{u}_j \int_\Omega \phi_j(x) \phi_i(x) dV \\ &\qquad \Rightarrow \qquad M \bar{u}\end{split}\]

where \(M\) is the mass matrix which combines the integral for all different basis functions. For the second term we aknoledge that the degrees of freedom have no spatial and temporal effects thus we can take them out of the integral and derivatives.

\[\begin{split}\int_\Omega u_x dV & = \int_\Omega (\partial_x\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ & = \sum_{j = 1}^N \bar{u}_j \int_\Omega \partial_x\phi_j(x) \phi_i(x) dV \\ & \Rightarrow \qquad T \bar{u}\end{split}\]

where \(T\) is the socalled transport matrix, wich only depends on the basis functions. For the thrird part we apply integration by parts, while assuming that Neumann boundary conditions:

\[\begin{split}\int_\Omega \tilde{u}_{xx} dV &\quad = \int_\Omega (\partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ &\quad = -\int_\Omega (\partial_x\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \partial_x\phi_i(x) dV \\ &\quad = \sum_{j=1}^N \bar{u}_j \int_\Omega -\partial_x \phi_j(x) \partial_x\phi_i(x) dV \\ &\quad \Rightarrow S \bar{u}\end{split}\]

where \(S\) is the stiffness matrix, which can be computed independently from the actual unknowns. Now we can write our PDE in terms of linear algabra objects:

\[M \bar{u}_t + c T \bar{u} - \mu S \bar{u} = 0\]

which we modify to be in the format is expected by the temporal solvers:

\[M \bar{u}_t = K \bar{u}\]
Parameters
  • x (array_like(float)) – Global coordinates of all degrees of freedom.

  • connect (array_like(int), shape((num_ele, dofe/ele))) – Element to degree of freedom connectivety map.

  • c (float) – Advective constant.

  • mu (float) – Diffusive constant.

  • num_q (int) – Number of Gausian quadrature points.

  • order (int) – Order of the polynomial used by our element.

Returns

  • M (matrix, (sparse csr format)) – The mass matrix.

  • K (matrix, (sparse csr format)) – The combination of stiffness and transport matrix matrix scaled with the approprate constants \(K = \mu S - c T\).

  • b (vector, (dense array)) – The right hand side, because we consider a homogeneous PDE with diriclet conditions it is a zero vector.

pde.diffusive(x, connect, mu, num_q, order)

Time derivative of the PDE for diffusion problems.

\[u_{t} - \mu u_{xx} = 0 \qquad \forall \, x \in \Omega = [0, L] \quad \& \quad t>0\]

Which, is converted into a weak form through:

\[\int_\Omega (\tilde{u}_t - \mu \tilde{u}_{xx} ) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Where \(\tilde{u}\) is our approximation:

\[\tilde{u}(x) = \sum_{n = 1}^N \bar{u}_n \phi_n(x)\]

in which \(\bar{u}_n\) are the unknowns and \(\phi_n(x)\) are the basisfunctions of FE approximation space \(V_h\). Substituting this approximation leads to:

\[\int_\Omega (\partial_t\sum_{j = 1}^N \bar{u}_j \phi_j(x) - \mu \partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x) ) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

Which we split into different integrals:

\[\int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x))\phi_i(x) dV - \int_\Omega (\mu\partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV = 0 \quad \forall \quad \phi_i \in V_h\]

For the first integral we notice that the basis functions are constant through time, only the degrees of freedom \(\bar{u}_j\) vary through time. Similarly these degrees of freedom does not affect the integral over space \(\int_\Omega dV\). Thus we can write:

\[\begin{split}\int_\Omega \tilde{u}_t dV & = \int_\Omega (\partial_t \sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ & = \sum_{j = 1}^N \partial_t \bar{u}_j \int_\Omega \phi_j(x) \phi_i(x) dV \\ &\qquad \Rightarrow \qquad M \bar{u}\end{split}\]

where \(M\) is the mass matrix which combines the integral for all different basis functions. For the second term we apply integration by parts, while assuming that Neumann boundary conditions:

\[\begin{split}\int_\Omega \tilde{u}_{xx} dV &\quad = \int_\Omega (\partial_{xx}\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \phi_i(x) dV \quad \forall \quad \phi_i \in V_h \\ &\quad = -\int_\Omega (\partial_x\sum_{j = 1}^N \bar{u}_j \phi_j(x)) \partial_x\phi_i(x) dV \\ &\quad = \sum_{j=1}^N \bar{u}_j \int_\Omega -\partial_x \phi_j(x) \partial_x\phi_i(x) dV \\ &\quad \Rightarrow S \bar{u}\end{split}\]

where \(S\) is the stiffness matrix, which can be computed independently from the actual unknowns. Now we can write our PDE in terms of linear algabra objects:

\[M \bar{u}_t - \mu S \bar{u} = 0\]

which we modify to be in the format is expected by the temporal solvers:

\[M \bar{u}_t = K \bar{u}\]
Parameters
  • x (array_like(float)) – Global coordinates of all degrees of freedom.

  • connect (array_like(int), shape((num_ele, dofe/ele))) – Element to degree of freedom connectivety map.

  • mu (float) – Diffusive constant.

  • num_q (int) – Number of Gausian quadrature points.

  • order (int) – Order of the polynomial used by our element.

Returns

  • M (matrix, (sparse csr format)) – The mass matrix.

  • K (matrix, (sparse csr format)) – The stiffeness matrix scaled with the diffusivity constant \(K = \mu S\).

  • b (vector, (dense array)) – The right hand side, because we consider a homogeneous PDE with diriclet conditions it is a zero vector.

pde.projection(mesh, fun)

Projecting a 1D function \(f(x)\) on a finite element basis.

Lets create our approximation function,

\[f_h(x) = \sum_{n=0}^N \bar{u}_n \phi_n(x)\]

as a weighted summation of the basisfunctions of approximation. Where \(phi_n\) are the basisfunctions of our FE space \(V_h\). The unknows here are the weights \(\bar{u}_n\), these we call degrees of freedom. To find these DOFs we formulate a weak form:

\[\int_\Omega (f_h(x) - f(x)) \, \phi_i(x)\, dV = 0 \quad \forall \quad \phi_i \in V_h\]

in which we substitute our approximation function and separate the knowns from the unknowns. We find:

\[\int_\Omega \phi_i(x) \sum_{j=0}^N \bar{u}_j \phi_j(x)\, dV = \int_\Omega \phi_i(x) f(x) \, dV \quad \forall \quad \phi_i \in V_h\]

As the weights \(\bar{u}_n\) are independent of location, we can take them out of the integral:

\[\sum_{j=0}^N \bar{u}_j \int_\Omega \phi_i(x) \phi_j(x)\, dV = \int_\Omega \phi_i(x) f(x)\, dV \quad \forall \quad \phi_i \in V_h\]

Which can be rewritten as a system of linear equations, which is:

\[M \, \bar{u} = b\]

Where \(M\) is a matrix and \(\bar{u}\) and \(b\) are vectors.

Parameters
  • mesh (Mesh) – The mesh object which specifies all discretization.

  • fun (callable) – Function that acts as our right hand side (nonhomogeneous term).

Returns

  • M (matrix, (sparse csr format)) – The mass matrix.

  • b (vector, (dense array)) – The right hand side, caused by the non-homogeneous behavior.