notes.md

Notes on this DG implemetation

We wish to solve some differential equation on a subdivided domain. Instead of breaking some domain down into its components, we start with the components, and patch them together into a whole system.

Domains

Domains are regions with their own discretization rule. One domain may be used for the entire problem, or multiple domains may be stitched together using fluxes.

Spectral Element 2D

spec_elem.py:spectral_element_2D

The standard spectral element, a domain mapping to a reference square $[-1,1]^2$.

Spectral Element Mesh

spec_elem.py:spectral_mesh_2D

A mesh of 2D spectral elements, stitched together conformally. That is, connected edges are completely flush together, with their collocation nodes lining up perfectly together. Basis functions are stitched together, meaning that nodes on shared edges have support in both elements.

Problems

Diffusion problem

$$\partial_t u=c\operatorname{lapl}(u),~~~~u\in\Omega$$ $$u|_{\Gamma_D}=g_D,~~~\left.\left(\frac{\partial u}{\partial n}\right)\right|_{\Gamma_N}=g_N$$ The weak form of the diffusion problem becomes: $$\int_{\Omega} v\partial_t u~dV=-\int_{\Omega}(\nabla u)\cdot\nabla( cv)~dV+\int_{\Gamma}cv(\nabla u)\cdot n,dS$$ $$=-\int_{\Omega}c(\nabla u)\cdot(\nabla v)+(\nabla c)\cdot(v\nabla u)~dV+\int_{\Gamma}cv(\nabla u)\cdot n,dS$$

Over a basis ${\phi_k}$, the Galerkin problem is $$\langle\phi_i,\phi_j\rangle_{L^2} \dot u^j= -(c\phi_i,\phi_j){H^1}u^j+F_i$$ If we write $c=c^k\phi_k$, and expand the gradient of $c\phi_i$ as we did in the integral above, $$(c\phi_i,\phi_j){H^1}=\int_{\Omega}(c^k\phi_k(\nabla \phi_{i})+\phi_ic^k(\nabla \phi_k))\cdot\nabla \phi_jdV$$ If we are using SEM, then we approximate the integrals as $$\langle \phi_i,\phi_j\rangle_{L^2}=\int_\Omega \phi_i\phi_jdV\approx \sum_{\ell=0}^MJ(\vec x_\ell)\phi_i(\vec x_\ell)\phi_j(\vec x_\ell)=J(\vec x_i)\delta_{i,j}$$ $$(c\phi_i,\phi_j){H^1}\approx \sum{\ell=0}^M J(\vec x_\ell)w_\ell\left(c^k\delta_k^\ell(\nabla\phi_i)(\vec x_\ell)+\delta_i^\ell c^k(\nabla\phi_k)(\vec x_\ell)\right)\cdot (\nabla \phi_j)(\vec x_\ell)$$ $$=\left(\sum_{\ell=0}^M J(\vec x_\ell)w_\ell c^\ell(\nabla\phi_i)(\vec x_\ell)\cdot (\nabla \phi_j)(\vec x_\ell)\right) +J(\vec x_i)w_i c^k(\nabla\phi_k)(\vec x_i)\cdot (\nabla \phi_j)(\vec x_i)$$ where $J$ denotes the Jacobian of the transform between real space and the reference element.

When subdividing for discontinuous Galerkin, the additional boundaries created should have a flux equivalence.

Question on Flux Equivalence

If we take the weak form and combine it over two domains, one would expect to recover the weak form on the entire domain. This would imply that the corresponding flux terms cancel on the shared boundary. Since this is over all test functions, would this imply *c partial_n u* is continuous across the boundary, or are there other weak-problem considerations that need to be made (for example, the $\nabla c$ becoming Dirac for discontinuous c, adding another boundary term)?

Boundary Condition Enforcement

Handling the boundary conditions can be done in one of to ways. Either the space of functions can be reduced, or the flux can be fixed to enforce the conditions. It may be good to study both approaches.

Space discretization gives us the system $$M\dot U +KU=F(\dots)$$ To enforce Dirichlet conditions here, we can set $(M^{-1}F)i=(M^{-1}K_1U)i+\partial_t (g_D)i$ for $i$ that has $U_i$ correspond to the Dirichlet boundary. This should be equivalent to modifying the problem to a homogeneous Dirichlet by setting $\tilde u\leftarrow u-u^0$, for $u^0$ a function agreeing with $g_D$ on the boundary for the case of a diagonal mass matrix. To see this, let $\mathring V^h\subseteq V^h$ be the finite spaces of functions for Galerkin, with ${\phi_1,\dots,\phi_M}$ as a basis for $\mathring V^h\subseteq H_0^1$, and ${\phi_1,\dots,\phi{M+B}}$ as a basis for $V^h$. Two formulations are as follows: $$\langle \phi_i,\partial_t u\rangle = - \langle \nabla \phi_i, \nabla u\rangle + f(\phi_i,\nabla u),~~~~i=1,\dots,M+B,~~~u|{\partial \Omega} = g_D$$ $$\langle \phi_i,\partial_t \tilde u\rangle = - \langle \nabla \phi_i, \nabla \tilde u\rangle - \langle \nabla \phi_i, u^0\rangle + f(\phi_i,\nabla u) - \langle \phi_i,\partial_t u^0\rangle,~~~~i=1,\dots,M$$ with some $u^0\in V^h$ having $u^0|{\partial \Omega} = \Pi_{V^h}(g_D)$ ($g_D$ projected into $V^h$). The first problem solves for $u\in V^h$, while the second solves for $\tilde u\in \mathring V^h$. Expanding $u = \sum u_j\phi_j$, noting $u_j=\tilde u_j + u^0_j$, and writing inner products as the corresponding matrices, we can write the first problem as

$$M_{i,j}\dot{u}j = - K{i,j} (\tilde u_j + u_j^0) + f(\phi_i,\nabla \phi_j)u_j,~~~i=1,\dots,M+B$$ where we cannot distribute the dot product over the expansion of $u$ yet. The second problem yields the same RHS: $$M_{i,j}\dot{\tilde u}j + M{i,j}\dot{u}^0_j=-K_{i,j}\tilde u_j-K_{i,j}u^0_j+f(\phi_i,\nabla \phi_j)u_j,~~~i=1,\dots,M$$ Locking $u|{\partial \Omega} = g_D$ is based on manipulating $\dot u$ after the mass inversion. The notation here hides the fact that the inversion is different in each problem. However, in the case that $M{i,j}=0$ for $i\le M<j$, we've isolated whatever $M^{-1}$ would have given us on the boundary from what it would have given us on the interior.

in the case of $u^0_j=0$ for $j\le M$ and $\tilde u_j=0$ for $j > M$ should be equivalent to f

We should theoretically be able to enforce Neumann conditions modifying the flux term with the substitution $(\nabla u)\cdot n \to g_N$. The other option would be at the Galerkin level, constraining $U$, which sounds pretty hard. I would wonder how different they are in terms of (enforce -> discretize) versus (discretize -> enforce).

Wave Equation (2nd order)

$$\partial_{tt} u=\nabla\cdot(c^2\nabla u),~~~~u\in\Omega$$ $$u|{\Gamma_D}=g_D,~~~\left.\left(\frac{\partial u}{\partial n}\right)\right|{\Gamma_N}=g_N$$ The weak form of the diffusion problem becomes: $$\int_{\Omega} v\partial_{tt} u~dV=-\int_{\Omega}c^2 (\nabla u)\cdot(\nabla v)~dV+\int_{\Gamma}c^2v(\nabla u)\cdot n,dS$$ Here, we keep the $c^2$ together, since that is how we will store the speed field. The terms are very similar to that of the diffusion problem above.

Boundary Condition Enforcement

For a function $u^0$ with $u^0|{\Gamma_D} = g_D$, the weak problem can be found as $$\int{\Omega} v\partial_{tt} w~dV=-\int_{\Omega}c^2 (\nabla w)\cdot(\nabla v)dV - \int_{\Omega} v\partial_{tt} u^0dV-\int_{\Omega}c^2 (\nabla u^0)\cdot(\nabla v) +\int_{\Gamma_N}c^2vg_N,dS$$ over test functions $v$ that evaluate to zero on $\Gamma_D$. The last terms have no $w$, dependence, so can be combined into a functional $f(v)$. Solving for $w$ in the same space as where $v$ varies, we can reconstruct $u = w+u^0$. For the quadrature for SEM, this is equivalent to just working with $u$ on the full space (no restricting the test functions), and enforcing dirichlet by explicitly setting the value of $u$ on the boundary.

Wave Equation (1st order system)

$$\partial_t \sigma = {c^2} \nabla u$$ $$\partial_t u = \nabla \cdot \sigma$$

$$\int_\Omega \tau \cdot \partial_t \sigma dV= \int_{\Omega }{c^2} \tau\cdot \nabla udV=-\int_{\Omega}{c^2} u\nabla\cdot \taudV + \int_{\partial\Omega}{c^2} u\tau\cdot ndS$$ $$\int_\Omega v\partial_t udV = \int_\Omega v\nabla\cdot \sigmadV = -\int_{\Omega}\sigma\cdot \nabla vdV + \int_{\partial\Omega}v\sigma\cdot ndS$$ $$\langle V,\partial_t U\rangle + a(U,V)=\int_{\partial\Omega}V^T\begin{pmatrix}0&0&{c^2} n^1\0&0&{c^2} n^2\n^1&n^2&0 \end{pmatrix}UdS$$ where $V=(\tau^1,\tau^2,v)^T,U=(\sigma^1,\sigma^2,u)^T$. This matrix, if we call it $A$, has eigenpairs $$0,\begin{pmatrix}-n^2\ n^1\ 0\end{pmatrix};~~~~-c,\begin{pmatrix}n^1\ n^2\ -1/c\end{pmatrix};~~~~c,\begin{pmatrix}n^1\ n^2\ 1/c\end{pmatrix}$$ The 1D advection equation ($\partial_t u +\gamma\partial_xu=0$) has weak form $$\int_{\Omega}v\partial_t udV=\int_{\Omega}\gamma u\partial_x v~dV-\int_{\partial\Omega}uv ,(\gamma n)~dS$$ where $-\gamma n$ tales the place of $A$ ($n$ is the outward facing "surface" normal). A dG upwind scheme would use the value of $u$ on this side if $\gamma n > 0$ (the normal and velocity are in the same direction). Motivated by this, we would use $U$ on this side for directions with negative eigenvalue. More specifically, we can write $$A = -c e_- \omega^- + c e_+ \omega^+$$ where $e_\pm$ are the corresponding eigenvectors and $\omega^\pm$ are the corresponding covectors to the eigenbasis. When evaluating the flux, we would use $U$ corresponding to the upwind values (this element plugs into $\omega^-$ and the adjacent element plugs into $\omega^+$). The corresponding covectors are $$\omega^-=\begin{pmatrix}\frac{n^1}{2}& \frac{n^2}{2} & -\frac{c}{2}\end{pmatrix},~~~~\omega^+=\begin{pmatrix}\frac{n^1}{2}& \frac{n^2}{2} & \frac{c}{2}\end{pmatrix}$$

TODO

• Clean up lagrange_deriv(), maybe move it to GLL_UTIL. tensordot will likely make it a lot better.
• general einsum optimizations?
• study the Dirichlet condition as a restriction. Is it equivalent once we use SEM's quadrature, or are we missing something?