kinnala / scikit-fem

Simple finite element assemblers
https://scikit-fem.readthedocs.io
BSD 3-Clause "New" or "Revised" License
495 stars 79 forks source link

Transport Equation using Discontinuous Galerkin Method #924

Closed jpwank closed 1 year ago

jpwank commented 1 year ago

Hi, I'm trying to solve the following linear transport equation using the discontinuous galerkin method (see e.g. 10.1002/zamm.200310088) $u_t + \nabla \cdot (\mathbf{a}u)=0, \qquad u(t=0)=u0$. The weak formulation of this problem (only spatial discretization) yields the term $\int{\partial\Omega} \widehat{\mathbf{a}u}\cdot\mathbf{n}vds$. With $v$ being the testfunction and $\widehat{\mathbf{a}u}$ the numerical flux. A stable formulation of the numerical flux is given by $\widehat{\mathbf{a}u}=\mathbf{a} \lbrace u \rbrace +\mathbf{C} \cdot [u]$. From this follows next to the term $\sum_{E\in\varepsilon_h} \intE [u][v] ds$ the term $\sum{E\in\varepsilon_h} \intE \lbrace u \rbrace[v] ds$. With the average $\lbrace u \rbrace =0.5(u^1+u^2) $. Since assembling the first term is implemented in jump() using asm() I was wondering if it would be possible to calculate the second term adapting the jump function. Following the section "Assembling jump terms" in the How-to guides the term can be split into $\sum{E\in\varepsilon_h}( 0.5\int_E u_1v_1ds - 0.5\int_E u_1v_2 ds+ 0.5\int_E u_2v_1ds - 0.5\int_E u_2v_2ds)$. The adapted code looks as following:

def average(w: FormExtraParams, *args):
out = []
for i, arg in enumerate(args):
    if i == 0:
        out.append(0.5 * arg)
    if i == 1:
        out.append((-1.0) ** w.idx[1] * arg)
return out[0] if len(out) == 1 else tuple(out)

Would this be a correct approach? I'm also interested to know if there are any stabilization techniques in scikit-fem for convection dominated problems using continuous elements.

kinnala commented 1 year ago

Yes, I think your approach is correct. I suggest you try it against some simple test case. If it works then we can add it to the library.

I once tried SUPG stabilization perhaps 5 years ago but the library was looking very different back then. AFAIK there is no such code available anywhere but it should be possible to implement.

kinnala commented 1 year ago

Closing due to inactivity. Feel free to open another thread under "Discussions".