11
I would not insist on demanding a geometrical meaning from Galerkin methods in general. There is a connection, but it becomes less meaningful as you extend it further and further. (In a sense, it would be better to call Galerkin methods "generalized projection methods".) To really understand the connection between collocation and Galerkin methods requires ...
10
The reason is that with the exception of linear problems, if you do a Fourier (or other) decomposition in time, you end up with a significant number of problems that are coupled globally in time. In other words, you have to solve lots of problems on the entire time interval concurrently. That will typically bust your computational or memory budget. The ...
8
I see several issues:
The DFT computed with fft puts the zero mode at the beginning of the array, and if you want to compute the derivative, it is necessary to apply fftshift/ifftshift to the array N to make sure the derivative is correct. It is easy to see for yourself what the correct expression is by working it out with pen and paper, and see also the ...
8
Added after my initial answer:
It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition numbers) in the table while giving maximum absolute errors rather than norm relative errors or maximum elementwise relative errors (these are all different measures.) ...
8
"Get more RAM" may be one of your best options. :) Prices are reasonably low right now, and it's one of the best upgrades you can gift your computer anyway. 10k x 10k is borderline but still doable on modern computers: that matrix takes $10^4 \times 10^4 \times 8$ bytes, that is, 760 MiB. On my laptop that code runs without problems.
Another option is ...
6
I hope I understood the question correctly. They try to compute exactly the same thing, so they really are equivalent. I'll use Chebyshev polynomials because they are easy to analyze.
Given a function $f(x)$ on $[-1,1]$, the spectral interpolant is the truncation of
$$ \begin{aligned}
f(x) &= \sum_{n\geq0} \bar a_n T_n(x), \\
\bar a_n &= \frac{1+[n&...
6
For domains that are logically square or cubic (like your quadrilateral), you can use the tensor product (dimension-by-dimension) approach. That is, generate the 2D Gauss-Lobatto point matrix as the tensor product of your 1D Gauss-Lobatto point vectors.
An example: if your 1D Gauss-Lobatto points are $(x_1,x_2)$, then in 2D you get the following four points:...
5
Let $w_k$ be the k-th column of W (the kth eigen-vector) and $v_k$ be the k-th element of v (the kth eigen-value). Then we can write:
$$ G = \sum_k w_k v_k w_k^T $$
which, element-wise, is equivalent to:
$$ G_{i,j} = \sum_k v_k w_{i,k} w_{j,k} $$
where the indices are column-major.
Since we always get $w_k$ twice, an overall sign on a column of $W$ cancels.
...
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 ...
4
First, some notation for clarity. We're talking about solving $u_{xx}=f(x)$. So $u$ is the solution and $f$ is the right hand side.
What makes you think either of your solutions is or is not correct?
To check your results, just compute two derivatives of your solution (using finite differences or FFT, either one) and compare with $f$. In your example, $f$ ...
4
For something with a spectral flavor in time, look at deferred correction methods, starting with this paper. I would argue that they're not spectral in the usual sense of the word, but they give you a family of arbitrary-order Runge-Kutta methods, so if you think of "refining" by increasing the order (by adding more nodes), then the convergence can be ...
4
Chebyshev are orthogonal wrt to a weight function which is singular at the end-points. When you approximate a function f(x) with Chebyshev, the convergence of the approximations is not affected by the values of f or its derivatives at the end points. The convergence rate depends only on the smoothness of the function f and not on its boundary values. This is ...
4
This is probably too late to help you, but here is a compact code that will do points, weights and first derivatives for Gauss, Lobatto or either Radau.
function [x,w,A] = OCnonsymGLReig(n,meth)
% code for nonsymmetric orthogonal collocation applications on 0 < x < 1
% n - interior points
% meth = 1,2,3,4 for Gauss, Lobatto, Radau (right), ...
4
Pros:
With trigonometric basis functions your problem size is $N \text{log}(N)$ instead of $N^2$.
Stabilization techniques are easy to implement and cheap:
Filtering in the modal space.
Zero padding in the modal space.
No aliasing due to the Galerkin ansatz.
Energy/Entropy stable disctretizations, e.g. via a skew symmetric implementation, are quite easy.
...
4
Ok, here comes the answer promised in the comment section. Let's start the other way round, going from a general grid to Gaussian grids and further constructions such as spectral elements.
In grid methods, one basically selects a number of $N+1$ gridpoints $\{x_k\}_{k=0}^{N}$. As basis functions, one can use Lagrange polynomials constructed over these nodes, ...
3
The key steps are to consider the advection equation
$u_t + au_x = 0$
where $a=\omega/k$ is the advection speed. Exact solutions to this equation is of the form $u(x,t) = f(x-at)$, where $f(y)$ is an arbitrary function.
For example, discretize using a standard Galerkin method we derive the weak form
$\int_\Omega v u_t dx + \int_\Omega a u_x = 0$
Assuming ...
3
For incompressible spectral DNS, often you'll see reference to John Kim, Parviz Moin and Robert Moser (1987): Turbulence statistics in fully developed channel flow at low Reynolds number, Journal of Fluid Mechanics, 177, pp 133-166, which includes details on how continuity is enforced. You might also find the book Spectral Methods in Fluid Dynamics by ...
3
One approach to acceleration would go as follows: Assume that $A^t=A(u^t)$ is the matrix you try to solve with, i.e., you are looking to solve the linear systems
$$
A^t x^t = b^t.
$$
Let me assume for a moment that $f^t \ge 0$, then $A^t$ is a symmetric and positive definite matrix. (If my assumption should be wrong, then it is still symmetric but may no ...
3
1) One usually chooses one or two representative wall-normal distances and presents the spectra for <u'u'>, <v'v'>, etc. over the homogeneous directions. For example, see Figure 11 within Spalart 1988 (http://dx.doi.org/10.1017/s0022112088000345).
2) By "scale" do you mean a lengthscale? If so, a lengthscale quantifying what aspect of the flow?...
3
I'm assuming that you know how Chebyshev collocation methods work (but if not, let me know and I'll explain a bit more); a good introduction is Nick Trefethen's Spectral Methods in Matlab as well as his Approximation Theory and Approximation Practice (in particular Chapter 21).
(But note that this code does Legendre collocation, not Chebyshev collocation!)
...
3
First, Galerkin in his article from 1915 does not discuss any weighting other the $\omega(x) = 1$. The Galerkin method is a direct generalization of the Rayleigh-Ritz method, and a variational procedure cannot be constructed with any other weight. That approach was the only one used until the FFT came about. If you want to use FFT, then you must use ...
3
$c$ depending on time is not the issue. You will use an RK scheme which takes care of this. The issue is $c$ is discontinuous in $x$. I recommend SBP-SAT schemes for this.
(1) Derive an energy equation at PDE level. (2) Search literature for SBP-SAT schemes which enforce interface conditions via SAT penalty terms, which are designed to mimic the energy ...
3
Indeed the problem was that you were trying to calculate your nonlinear term incorrectly and you forgot that:
$$\mathcal{F}[(f(x))^{2}] \neq (\mathcal{F}[f(x)])^{2}$$
I changed your code to this:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
nu = 1
L = 100
nx = 1024
t0 = 0
tN = 200
dt = 0.05
nt = int((tN - t0) / 0....
3
There are some unknowns in what you are doing but for simplicity, suppose we want to find $u(t)$ as discrete times $t_1, t_2, \cdots, t_n$. Let $\textbf{F} = [F(t_1), F(t_2), \cdots, F(t_n)]^T$ and $\textbf{u} = [u(t_1), u(t_2), \cdots, u(t_n)]^T$ be column vectors representing $F$ and $u$ evaluated at the desired times. From your problem statement, you wish ...
2
Goofy question, are you dealiasing in some fashion or is it unnecessary in your setup? The issue sounds convective in some sense otherwise I'd expect your diffusion dominated limit to misbehave. If not, you might see if using the skew-adjoint form of the convective operator improves things somewhat (see https://github.com/RhysU/suzerain/blob/master/...
2
In the comments you mentioned that $f$ has a very special form. Let $u:[a,b]\rightarrow \mathbb{R}$. Then $f:[a,b]\rightarrow\{0,1\}$ is defined for some set $S \subset [a,b]$ as follows:
$$f(x) = \begin{cases} 1 & \text{if $x \in S$,} \\ 0 & \text{otherwise.} \end{cases}$$
As Bill Barth mentioned, you can now split the matrix $A$ into two parts as
$$...
2
You are correct that the $k=0$ mode corresponds to the mean (volume averaged) velocity in the domain. Because the equations for Stokes flow are Galilean invariant they are undetermined up to a constant (the mean flow). Usually this ambiguity is resolved by the boundary conditions, but purely periodic boundary conditions offer no such resolution. Instead, you ...
2
If you want to use nonperiodic boundary conditions, you could change the Fourier basis in the $y$-direction for your spectral method to a Chebyshev basis. Fast Chebyshev Transforms make use of Fast Fourier Transforms; both are linear transformations on function spaces, so switching between the two amounts to a change of basis, and you can reuse significant ...
2
Thanks for Kirill's detailed answer, which clarifies all the confusion in my head. According to Kirill's answer and the materials he provided, now I want to generalize it a bit to common cases.
Let us suppose $\{F_n(x)\}$ is a set of orthogonal polynomials on $[-1,1]$, i.e.,
\begin{equation}
\int_{-1}^1 F_m(x)F_n(x)w(x)dx=g_m\delta_{mn},
\end{equation}
...
2
Frequently some sort of velocity-pressure splitting is employed which gives a step in the method where the velocity and pressure are coupled directly.
Only top voted, non community-wiki answers of a minimum length are eligible