14
The tradeoffs below apply equally to DG and to spectral elements (or $p$-version finite elements).
Changing the order of an element, as in $p$-adaptivity, is simpler for modal bases because the existing basis functions do not change. This is generally not relevant to performance, but some people like it anyway. Modal bases can also be filtered directly ...
11
Stability does indeed mean that small changes in the data lead to small changes in the solution. This can be shown for the linear advection equation through the energy method; in proving the stability of a discretization, we seek to mimic the energy estimate of the exact problem. Considering the PDE,
$$
\frac{\partial u}{\partial t} + a \frac{\partial u}{\...
10
For parabolic/elliptic PDE's, I highly recommend Beatrice Riviere's book: Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation.
For hyperbolic PDE's and general (i.e. nonlinear) conservation laws, I recommend Hesthaven & Warburton's book: Nodal discontinuous Galerkin methods: algorithms, analysis, and ...
9
One reason DG methods may receive more attention as a parallel method is that it is readily seen that the method is inherently local to an element. The coupling in DG methods is weak, as it only occurs through adjacent edges (or faces in 3d). So, for triangles or quads DG will communicate to three or four processers at most, respectively. Whereas CG methods ...
9
You are confused about different concepts. A mesh is really just a collection of cells defined by the vertices of the mesh and which vertices together form each cell. Consequently, a mesh is an entirely geometric object. It knows nothing about finite element spaces one may define on it -- that's for a later step in your workflow. It's also important that ...
8
From my many years writing FEM software, I believe that the statement that DG schemes are better suited to parallelization than CG schemes is apocryphal. It is frequently used in introductions of DG papers as a justification for DG methods, but I have never seen it substantiated with a reference that actually investigated the question. It is similar to every ...
8
The argument is misleading. The different spaces that are used are $P_1=\text{span}\{1,x,y\}$ and $Q_1=\text{span}\{1,x,y,xy\}$. For higher orders, they are
$P_2=\text{span}\{1,x,y,x^2,y^2\}$ and $Q_2=\text{span}\{1,x,y,x^2,y^2,xy,x^2y,xy^2,x^2y^2\}$. If you continue this, it is true that the number of basis functions for the $Q_k$ spaces grows faster than ...
6
The problem with ABC in the form of derivative operators at the boundary is that the exact ABCs are non-local in space and often in time, which makes them hard to use in numerical simulations. You can derive approximate ABCS based on series expansions, but this can be tedious and it may be difficult to get stable schemes. This approach was made famous by a ...
6
The penalty usually scales as $O(p^2/h)$ for polynomial discontinuous Galerkin methods regardless of dimension, but for a more precise constant, it's usually enough to take the penalty large enough to guarantee coercivity with respect to the DG norm. Shahbazi gives an expression for the penalty over a given face which just depends on $p$ and geometric ...
6
The restrictive CFL of DG schemes typically comes from the combination of high order accuracy and a compact stencil (see this reference for example). The CFL depends on bounding the variational form in terms of the $L^2$ norm of the solution, which depends on derivative and traces of polynomials. Bounds for each of these quantities (using Bernstein or ...
6
Of course, the Fourier transform is a linear operator. So, you have the kinetic energy defined as: $E(\mathbf{r}) = \frac{1}{2} \mathbf{u}(\mathbf{r}) \cdot \mathbf{u}(\mathbf{r})$. The Fourier transform of $\mathbf{u}$ and $E$ are:
$$\tilde{E}(\mathbf{k}) = \int_{\Omega} E(\mathbf{r}) e^{-i \mathbf{k} \cdot \mathbf{r}} d^{3} \mathbf{r}$$
$$\tilde{\mathbf{...
answered Sep 20 '19 at 20:25
Mithridates the Great
2,1642 gold badges6 silver badges25 bronze badges
5
The problem was solved recently: Lederer/Schöberl: Polynomial robust stability analysis for H(div)-conforming finite elements for the Stokes equations
5
Your observed quadratic convergence indicates that the Jacobian is likely correct. Have you looked at the solutions for your under-resolved configurations? Galerkin optimality uses the operator norm, which contains only the symmetric part, thus the solution of the discrete system could be quite different from the projection of the exact solution. This ...
5
Since you haven't had an answer yet, I'll reformulate my comment.
Saying that a method is $p$-th order accurate implies that a polynomial manufactured solution of lesser order can be captured exactly. For example, a 2nd order method will represent a linear solution up to machine precision. This has often helped me in finding implementation issues. It's ...
5
Just as it is the case with most general statements about numerical schemes, the answer depends on the exact circumstances you're looking at. First of all, the advantages of DG concerning parallelization mainly pay off in case of explicit time-integration schemes because of
The cell-local mass matrix of DG schemes (so you don't have to apply the inverse of ...
5
The concept of $p$-nonconforming meshes can also apply to continuous FEM, not just DG. For continuous FEM, continuity is enforced by enforcing a single set of degrees of freedom on a face shared between two elements. If those elements have varying polynomial degrees, the trace space on the face must be made the same. This may be done by restricting the ...
5
"Spectral methods" usually means methods which make use of global basis functions. Fourier spectral methods use sine and cosine, and are used when you have periodic boundary conditions. Chebyshev methods use Chebyshev polynomials and are useful in non-periodic cases. These two methods are used in DNS, see e.g., hit3d which uses fourier and periodic bc, and ...
5
I believe that for odd polynomial degrees you can get superconvergence by one order at the jump. For even degrees I don't think that's the case.
For the odd case, we consider a smooth function $u:\mathbb{R}\to\mathbb{R}$. Consider two adjacent intervals $I^-=[x_{i-1},x_i]$ and $I^+=[x_i,x_{i+1}]$ both of length $h$, and denote the $L^2$ projection of $u$ on ...
4
I was curious to see some answers to this question, but somehow nobody bothers to reply...
Regarding literature, I really like the book Spectral/hp Element Methods for Computational Fluid Dynamics (there's also a cheaper soft-cover version now) and also the book of Hesthaven and Warburton. These two go into quite some detail that will help you implement the ...
4
The short answer to your question is yes, you can invert a large number of small, independent matrices on a GPU, and more than likely you can do it efficiently. The best way to go about it, however, is not such a straight forward answer. I can think of three possible implementations, although I will admit right from the start I have never attempted these, ...
4
The penalty term does indeed act as a stabilization - coercivity is usually shown using norms which depend on this penalty, and you typically need a "sufficiently large penalty" to show convergence of symmetric interior penalty DG.
I believe one exception is NIPG (non-symmetric interior penalty DG, they flip the sign of the numerical flux terms involving ...
4
No. The statement means that the finite element method, without modifications, will lead to numerical solutions for conservation laws that violate many of the physically reasonable constraints (e.g., that solutions have certain monotony properties, that entropy can only grow, of that the density can not become negative -- in mathematical language, the ...
4
DG methods are a more general description, referring to finite element methods which represent the solution in a discontinuous manner (usually with coupling between elements through a numerical flux).
The specific numerical flux tends to differ between type of equation (hyperbolic vs parabolic/elliptic) and formulation (2nd order equation vs first order ...
4
Since you read Riviere's book and know how to assemble element integrals I assume you are familiar with concepts such as
transformation to reference elements,
numerical quadrature.
These techniques are applied analogously to edge integrals - with two differences:
Integrals are transformed to integrals over edges of the reference element, and then a ...
4
Basically, you have the following set of derivations:
Strong form of the PDE -> Finite differences
Integral form of the PDE -> Finite volumes
Weak (variational) form of the PDE -> Galerkin methods (including FEM and DGFEM)
4
As far is I know, a proof for the general case has not been presented so far. The most detailed analysis has probably been performed by Krivodonova & Qin 2013 in this work. Quoting from the conclusion:
We have derived a closed form expressions for the eigenvalues of the DG spatial discretization applied to the onedimensional linear advection equation ...
4
For the discontinuous Galerkin method, the choice of control points is entirely unimportant because they conceptually lie in the interior of the cell (even if they are physically on the boundary of the cell). We think of a control point that is conceptually located at a vertex to be one that all cells adjacent to that vertex share; and of one on an edge that ...
4
We are considering the one-dimensional scalar conservation law,
$$
\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x}= 0, \quad x \in \Omega, \quad t > 0,
$$
subject to appropriate initial and boundary conditions. For a DG method, we would like to seek solutions $u_h(\cdot, t)$ in the space $V_h \subset L^2(\Omega)$ containing functions that ...
4
Taking your questions in order, the definition of $D_r$ as the operator converting the nodal values of $u_h$ into derivatives follows immediately from the definition, although the notation is a little inexact,:
$$ u_h(r) = \sum_i u(r_i)l_i(r)$$
hence
$$\frac{du_h}{dr} = \sum_j u(r_j)\frac{dl_j}{dr} $$
so that
$$\left.\frac{du_h}{dr}\right|_{r_i} = \sum_j u(...
4
We can break down the code
rhsu = -a*rx.*(Dr*u) + LIFT*(Fscale.*(du));
into the two parts
-a*rx.*(Dr*u)
and
LIFT*(Fscale.*(du));
The first part is simply taking a derivative in reference space by multiplying with the matrix Dr, then transforming into physical space by multiplying with the Jacobian rx and then the coefficient a. Next we need to do the ...
Only top voted, non community-wiki answers of a minimum length are eligible