Initial Value Problem using Finite Element

Question

I am trying to implement a FEM solver for the following initial value problem \begin{align} \frac{\partial u}{\partial t} - \nabla^2 u &= f\quad \text{ in } \Omega\times (0,T)\\ u &= g\quad \text{ on }\partial\Omega\times(0,T)\\ u(x,y,0) &= u_0(x,y) \end{align} where $\Omega\subset\mathbb{R}^2$ is just a rectangular domain. The functions $f$ and $g$ can vary with time. I have discretized this IVP into $$\mathbf{B}\frac{\partial}{\partial t}\mathbf{u} + \mathbf{A}\mathbf{u} = \mathbf{l} $$ where \begin{align} \text{B}_{ij} &= \int_\Omega \phi_i\phi_j\\ \text{A}_{ij} &= \int_\Omega \nabla\phi_i\nabla\phi_j\\ l_i &= \int_\Omega \phi_i f \end{align} and I can assemble these matrices and vector. I have also discretized time to get the time-stepping scheme $$(\mathbf{B}+k\theta\mathbf{A})\mathbf{u}^{n+1} = (\mathbf{B}-k(1-\theta)\mathbf{A})\mathbf{u}^n + k(\theta \mathbf{l}^n + (1-\theta)\mathbf{l}^{n+1}), $$ where $\theta\in[0,1]$ and $k$ is the time-step.

Here comes my problem. How exactly do I impose the boundary and initial condition, implementation-wise? What changes should I make to the matrices $\mathbf{A},\mathbf{B}$, etc? Should I do anything at each time-step? I feel like it shouldn't be too hard, however when the time is introduced it suddenly confuses me.

Thanks

Answer 1

There are several equivalent ways of implementing Dirichlet boundary conditions with the finite element method. I'll give a brief overview, but you'll probably need more details, which can be found in the book Understanding and Implementing the Finite Element Method.

It may be easier to analyze this for the Poisson problem $-\nabla^2u = f$ than for the heat equation; the techniques are the same in both cases, and if you understand them for the steady-state problem then it shouldn't be a big leap to adapt them to the time-dependent one. For the initial value, that amounts to setting $u^0_i = \mathbf{u}_0(x_i)$. You'll then use this as part of the right-hand side in the linear system you've written down when you go to solve for the next time-step $\mathbf{u}^1$. The boundary conditions on $\partial\Omega$ are trickier.

One way to implement Dirichlet boundary conditions is to alter the matrices $A$ and $B$ and the load vector in order to guarantee that the boundary conditions are satisfied exactly. Suppose that $i$ is some unknown of the linear system corresponding to a Dirichlet boundary point; the modification you have to make is to zero out row $i$ of $B$, set row $i$ of $A$ to be the corresponding row in the identity matrix and set $\mathbf{l}_i$ equal to $g(x_i)$.

To illustrate this approach, suppose you had a $3 \times 3$ chunk of the linear system that looked like this:

$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 \\ l_2 \\ l_3\end{matrix}\right]$

and that $u_2$ is a degree of freedom corresponding to a boundary node. The modified system would look like this:

$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\ 0 & 1 & 0 \\ a_{31} & a_{32} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 \\ g_2 \\ l_3\end{matrix}\right]$

You can see how the modified system isn't symmetric anymore. Nonetheless, $u_1$ and $u_3$ are still coupled to $u_2$. The asymmetry of the system isn't usually an issue in practice provided that your initial guess for $u$ exactly satisfies the Dirichlet boundary conditions; try throwing the problem into a naive CG solver and see what happens. (Can you come up with a reason why this might be so?)

Another approach is to reduce the system size. Instead of solving for all of the coefficients $u_i$ in the expansion $\mathbf{u} = \sum_iu_i\phi_i$, you'll only solve for the coefficients corresponding to interior nodes. This means extracting sub-matrices out of $A$ and $B$ for those unknowns, and adding an extra component to the load vector $\mathbf{l}$ consisting of the effective forcing on the interior from boundary conditions. Referring back to the example matrices from before, we're left with the following reduced system for $u_1$ and $u_3$:

$\left[\begin{matrix} a_{11} & a_{13} \\ a_{31} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 - a_{12}g_2 \\ l_3 - a_{32}g_2\end{matrix}\right]$

Changing the linear system or deflating it are, from what I've observed, the most common approaches. You can also add on boundary conditions later as a constraint, or you can use penalty methods, but these are less common. Of the two, I tend to favor the first approach.

Answer 2

I guess you are also missing a discretized weak form for the boundary condition $u = g$, i.e. $\mathbf{G} \mathbf{u} = \mathbf{m}$, where $\mathbf{G}$ is assembled from the boundary elements (the basis function set for elements of $\partial \Omega$). I found the treatment in Volakis, et al., "Finite Element Method for Electromagnetics" to be quite helpful regarding this basis set.

Anyway, if both $f$ and $g$ change in time, you need to recompute the respective right-hand sides before the time-stepping iteration.

Just a guess: what about unifying the volume balance equation and the boundary flux equation in one? As in

$\mathbf{B}_{\Omega} \, \partial_t \mathbf{u}_{\Omega} + \mathbf{A}_\Omega \mathbf{u}_\Omega + \mathbf{G}_{\partial \Omega} \mathbf{u}_{\partial \Omega} - \mathbf{l}_{\Omega} -\mathbf{m}_{\partial \Omega} = 0$,

where the $\partial \Omega$ and $\Omega$ subscripts should indicate the restrictions to the respective subdomains (the first method in Daniel's comment above)