Define the inflow and outflow on interior boundaries

General Discussion
1

Hello everyone, I’m trying to implement an upwind DG scheme, but this requires defining the inflow and outflow on the interior boundaries, and I don’t know how to represent them. This has been bothering me for a long time.

Define mesh Th = square(128, 128, [2*x-1, 2*y-1]); Assume that the velocity \mathbf{u}_{h}^n has already been obtained. uh=[u1old,u2old] with fespace Vh(Th, [P2,P2]); We have \sigma_{h}^{n+1},w_{h}\in W_{h} with fespace Wh(Th, P2dc);

Now I want to reproduce an upwind DG scheme in W_{h}.
\left( \frac{\sigma_{h}^{n+1}-\sigma_{h}^n}{\delta t},w_{h} \right) + \left( \sum_{e}\left(\hat{\sigma}_{h}^{n+1}\mathbf{u}_{h}^n,[[w_{h}]]_{e}\right) - (\sigma_{h}^{n+1}\mathbf{u}_{h}^n,\nabla w_{h})\right) = 0.

For any function w continuous on each K\in \sum, we define
[[w]] = (w^{-} - w^{+})\cdot \mathbf{n}_e
to be the jump of the function w on an edge e=\partial K^{-}\cap \partial K^{+}, with w^{\pm} denoting the trace of w from K^{\pm} onto e, and \mathbf{n}_{e} denoting the normal vector on e pointing towards K^{+}.

The \hat{\sigma}_{h}^{n+1} is the numerical flux, which is a single-valued function defined at the cell interfaces and in general depending on the values of the numerical solution \sigma_{h}^{n+1} from both sides of the interface. Here, the numerical flux on an edge e can be chosen according to the upwind principle, namely,

\hat{\sigma}_{h}^{n+1} = (\sigma_{h}^{n+1})^{-}\quad \text{on }e_{+}^n,\qquad \hat{\sigma}_{h}^{n+1} = (\sigma_{h}^{n+1})^{+}\quad \text{on }e_{-}^n

where e_{-}^n and e_{+}^n denote the numerical inflow and outflow parts of an edge e at time t=t_{n}, defined by
e_{-}^n = \{\mathbf{x}\in e:(\mathbf{u}_{h}^n\cdot \mathbf{n}_{e})(\mathbf{x})\leq 0\},\quad\quad e_{+}^n = \{\mathbf{x}\in e:(\mathbf{u}_{h}^n\cdot \mathbf{n}_{e})(\mathbf{x})> 0\}.

I want to know how to implement the crucial feature in FreeFem++.
\sum_{e}\left(\hat{\sigma}_{h}^{n+1}\mathbf{u}_{h}^n,[[w_{h}]]_{e}\right)

2

It seems to me that this sum can be written (at least for the interior edges) as
\displaystyle\sum_T\sum_{e \subset T} \int_e\sigma_h^{n+1} \max(\mathbf{u}_h^n\cdot \mathbf{n}_e,0)(-1)[w_h]_e
where [w]_e=w^+ - w^- and “-” means the side of T, “+” the other side, \mathbf{n}_e goes from “-” to “+” ie points outside T, and \sigma_h^{n+1} is evaluated from the side of T.
This can be coded as

load "Element_PkEdge"
verbosity=0;

border bb(t=0.,1.){x=cos(t*2.*pi);y=sin(t*2.*pi);}
mesh Th=buildmesh(bb(20));

fespace Vh(Th,[P2,P2]);
Vh [u1old,u2old];
[u1old,u2old]=[-y,x];

fespace Wh(Th,P2dc);
Wh sig,wh;

fespace Wedgeh(Th,P2edge);// functions P2 on each edge (continuous through edge)
//this space could be useful

varf inout(sig,wh)=
  intalledges(Th)(sig*max(u1old*N.x+u2old*N.y,0.)*(-1.)*jump(wh))
  ;
4

Thank you very much, Professor, I understand now. Also, what if I choose the boundary condition
\sigma(\mathbf{x},t)|_{\Gamma_{\mathbf{u}}}=g(\mathbf{x},t) with \Gamma_{\mathbf{u}}=\{\mathbf{x}\in\partial \Omega: \mathbf{u}(\mathbf{x})\cdot \mathbf{n}<0\}

macro h(t)( ...) //

Th = change(Th, flabel = ( (label < 5) && (u1old*N.x + u2old*N.y < 0) ) ? 10 : label);

solve inout(sig,wh)=
  int2d(Th)(sig*wh/dt - (sig * u1old * dx(wh) + sig * u2old * dy(wh)) )
- int2d(Th)(sigold*wh/dt)
+ intalledges(Th)(sig*max(u1old*N.x+u2old*N.y,0.)*(-1.)*jump(wh))
- int1d(Th,10)( h(t) * (u1old*N.x+u2old*N.y) * (-1) * wh );

Is this correct, Professor?

5

I would put
+int1d(Th)( h(t) * max(-(u1old*N.x+u2old*N.y),0.) * (-1) * wh )
in order to represent an edge integral which is the boundary of a “ghost triangle” outside of the domain, but I do not guarantee that this is correct…
(and do not use the line Th=change())

6

I understand what you mean, Professor. But why can’t I use Th=change()?

I’ve provided the exact solutions for \mathbf{u} and \sigma, so I think I should use \Gamma_{\mathbf{u}}=\{\mathbf{x}\in\partial \Omega: \mathbf{u}(\mathbf{x})\cdot \mathbf{n}<0\} with the exact solutions \mathbf{u} instead of max(-(u1old*N.x+u2old*N.y),0.).

Another point, intalledges(Th)(sig*max(u1old*N.x+u2old*N.y,0.)*(-1.)*jump(wh)) is it unnecessary to distinguish between internal edges and boundary edges in this command?

7

You can use the change() and/or the \Gamma_u, but I think it is more in the spirit of the method to write the boundary term similarly as the interior terms.
I think that the intalledges has to include both internal and boundary edges.

8

Thank you, Professor. I think I understand how FreeFem++ handles the DG scheme.