# Dirichlet BCs - alternative implementation methods

I am having problems with solving a hyperbolic wave problem with Dirichlet BCs. I have tried reducing the time step sizes, which does not affect the results, and notices increasing the number of nodes makes the results worse. I have concluded the problem is with the boundary or certainly spatial.

My code is 1D, and the grid is structured. I have used FVM with central differencing. At the boundaries I have simply set the term at the boundary equal to the Dirichlet condition (for instance: velocity at east face $v_e = A$, where $A$ is the Dirichlet BC).

I was wondering if there are any ways using which I could get rid of the oscillations in my results by means of changing the way I am implementing the BC or any other quick fixes.

I have also tried upwinding, and did not even obtain convergence.

• What equation(s) are you solving? Can you write down the PDE? – Bill Barth Sep 10 '15 at 15:04
• standard inviscid euler equations. – melody Sep 10 '15 at 16:49
• Have you tried Lax-Wendroff with flux limiters? That's the favored approach in LeVeque's finite volume book. Also, are you sure you are setting the inflow BC and not the outflow BC? – Geoff Oxberry Sep 18 '15 at 0:14
• There's no way for us to answer your question, since you haven't specified completely the discretization you're using. – David Ketcheson Jan 17 '16 at 10:51

My hunch is that you may be prescribing boundary data at both ends of the domain (which leads to oscillations) or incorrectly specifying the boundary data. For a hyperbolic problem, you need to respect the flow of information in your numerical method. For example, consider the linear hyperbolic equation

$$\mathbf u_t = \mathbf A \mathbf u_x$$

with Dirichlet boundary conditions at the endpoints $x = x_L$ and $x = x_R$. For simplicity, we assume that $\mathbf A$ is a constant coefficient matrix with non-degenerate eigenvalues and linearly independent eigenvectors. By the spectral theorem we can write

$$\mathbf A = \mathbf P \mathbf D \mathbf P^{-1}$$ where $\mathbf D$ is a diagonal matrix of eigenvalues and $\mathbf P$ is the respective eigenvector matrix. Using the decomposition, we can rewrite our PDE as

$$\mathbf v = \mathbf P^{-1} \mathbf u$$ $$\mathbf v_t = \mathbf D \mathbf v_x$$

The system is now decoupled and the sign of the eigenvalue will tell you the direction that the characteristics $\xi = x + \lambda t$ are propagating. Characteristics are lines where the solution is constant. For example, for $\lambda_k > 0$, characteristics propagate from left to right. Consider the $(x, t)$-plane of the characteristics. At the initial time, $t = 0$, the initial condition supplies the starting point for the emanating characteristics. For $t > 0$, the characteristics travel from left to right. If we trace any characteristic from the initial time, we can see that the characteristic will eventually exit the domain at $x_R$. If we try to trace any characteristic emanating from $x_L$, we run into the issue that there is no information provided at this boundary point, thus we must prescribe the boundary condition for $v_k$ at $x_L$. If we were to prescribe a boundary condition at $x_R$ as well, this will lead to numerical instability for we are over constraining our system by telling it how the characteristics should behave at the outflow boundary.

If you are implementing the uncoupled system, then the boundary data that you should impose should still be in terms of the characteristic variables. For example, when solving

$$u_t = v_x$$ $$v_t = u_x$$

for $-1 \leq x \leq 1$, and Dirichlet boundary conditions $u(-1,t) = f_1(t)$, $v(-1,t) = f_2(t)$, $u(1,t) = f_3(t)$, $u(1,t) = f_4(t)$, it is easy to see that the characteristic variables are $(u + v)$ and $(u - v)$:

$$(u+v)_t = (u+v)_x$$ $$(u-v)_t = -(u-v)_x$$

and that we need to impose an incoming characteristic for $(u+v)$ at $x=-1$:

$$u(-1,t) + v(-1,t) = f_1(t) + f_2(t)$$

and an incoming characteristic for $(u-v)$ at $x=1$:

$$u(1,t) - v(1,t) = f_3(t) + f_4(t)$$

The first thing to do is to construct a simple test case, whose values you know analytically. So you might try a case with 7 nodes, where the solution is a simple function.

Once you have a system that can solve the simple case, then it is worth considering what would cause an oscillation in the results; from the wording it seems like you have a dynamics problem, but if this is a static solution then there are convergence acceleration approaches which can help, for example the Anderson mixing scheme for self-consistent problems.

Is there anything in the nature of the problem you are solving which would make it tend to overshoot, like a Coulomb repulsion?

Sometimes the stability of the solution can depend on scaling of the variables; if you have a large offset or datum then the differences that an FD scheme depends on can be truncated by their representation.

Finally also consider the precision of the variables in use. If you have a fine mesh and the floating point type cannot hold enough precision, then the results are just skipping about in a sparse solution space. So ensure you are using at least double precision and that there is no large datum offset in your BCs.

Without more detail, or results to show the problem, it is hard to offer more specific suggestions.