33

Yes, there aren't too many resources on this for some reason. For a very long time, the standard goto was "just use BDF methods". This mantra was set in stone for few historical reasons: for one Gear's code was the first widely available stiff solver, and for another the MATLAB suite didn't/doesn't include an implicit RK method. However, this heuristic isn't ...


10

It is because we typically neglect higher order terms in error estimates. For example, we can show that $$ \|e\| \le C(u) h^2 + {\cal O}(h^3). $$ The point is that when $h$ is small, the cubic term is small and can be neglected. In fact, when $h$ is small, you can observe quadratic convergence. But whenever $h$ is not small (where "small" is relative to ...


9

I have collected some of my experience on debugging numerical codes here: deal.II FAQ: debugging. I don't know if that would have helped you in this particular case, but it may in others.


7

There is another approach called the adjoint method, which is commonly used in inverse problems for PDE and which is quite easy to generalize to other problems. This is going to be long. Your observations give you a field $u_{\text{exp}}$; you'd like to find a value of $\mu$ for which $\frac{1}{2}\iint(u - u_\text{exp})^2dx\hspace{2pt}dt$ is a minimum, ...


6

Because the advection-diffusion equation is linear, there are many exact solutions. One reference with many exact solutions (including source terms) is M. Th. van Genuchten and W. J. Alves Analytical Solutions of the One-Dimensional Convective-Dispersive Solute Transport Equation U. S. Department of Agriculture, Agricultural Research Service Technical ...


6

Add a right-hand-side/forcing function, and use the Method of Manufactured Solutions. Then you can have any solution you want. Also, Separation of Variables works on this PDE just fine, so it also has an analytical solution.


6

You will have a problem if $x=0$ is part of your domain because in that case your advection velocity $u=1/x$ becomes singular. In particular, there will be cells close to $x=0$ where the cell Peclet number is very large, and consequently you will be in the advection-dominated regime. You will need to stabilize your discretization.


6

I think this greatly depends on what kind of physics you are trying to model even though for some problems both approaches are viable. Lagrangian vs. Eulerian Framework For certain problems Lagrangian frameworks are better suited than their Eulerian counterparts. For instance, if one is interested in studying the current patterns or sedimentation problems ...


6

The algorithm you have implemented is explicit. Crank-Nicolson is an implicit method, and thus requires a solve.


6

The stationary equation you show transports information from the right to the left via the advection term; it also diffuses slightly. If you switch off the diffusion term altogether, then you only have transport from the right to the left, and you need to also drop the boundary condition at the left: because information is from the right to the left, nothing ...


6

We have the following problem: $$\frac{\partial u}{\partial t}+v\color{red}{\frac{\partial u}{\partial x}}-\nu\color{blue}{\frac{\partial^2u}{\partial x^2}}=0 \tag{*}$$ The function $u$ may represent for example the concentration that propagates at velocity $v>0$ and disperses in a medium with viscosity $\nu>0$. Since only we are discussing how terms ...


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

In the advection dominated case the problem can develop physics which is invisible to a computational mesh that is too coarse (say by having elements which contain many wavelengths). This leads to instability because neighboring elements could disagree considerably on what they believe are the right values to represent this physics (think interpolating high ...


5

There are whole books written on this, but you should investigate (search the web, really) upwind diffusion and Streamline Upwind Petrov Galerkin finite element methods first. There many more methods that extend these ideas or approach the problem from a different direction.


5

The finite element method has similar problems to FD with regards to stability in solving the advection-diffusion equation, i.e. the same restrictions on Peclet number apply. One remedy, also similar to that used in FD, is to use a formulation that includes "upwinding". A nice set of lecture notes that discusses how upwinding is added to FE formulations is ...


5

There are two different flavors of smoothing and stability. The spurious oscillations in convection-diffusion problems are not an artifact of the linear solver but an inevitable artifact of the discretization method. Any linear solver for non-symmetric matrices you use will give you the same wiggly answer, or it's wrong. Multi-grid is one such linear solver. ...


4

What I describe below is a version of the method of characteristics. Its applications are limited, but if it works in your case, it would be a fairly simple, fast, and accurate process. If you want to discretize on a mesh, it might not be that helpful. Assuming your solutions $\phi$ are smooth, you can reduce this to a family of ODE problems. Exchange the ...


4

The answer is YES, you can impose the Dirichlet boundary condition weakly on the inflow boundary $\partial \Omega^-$. Say for the following advective equation: $$\begin{cases} u_t + \nabla \cdot (\mathbf{v} u) = 0 \quad \text{in }\Omega, \\ u= g \quad \text{on }\partial \Omega^-, \end{cases} $$ where $\partial \Omega^-$ is the inflow boundary: $\mathbf{v}\...


4

Well, you can use Crank-Nicolson here but then you'll have to construct and solve a linear system for each time-step. That's easy to do but it would be much easier to use an ODE integrator that is available in MATLAB. Due to the diffusion operator in the RHS the implicit integrator ode23tb seems to be a good choice here but experimenting with different ...


4

An adaptive grid is a grid network that automatically clusters grid points in regions of high flow field gradients; it uses the solution of the flow field properties to locate the grid points in the physical plane. The adaptive grid evolves in steps of time in conjunction with a time-dependent solution of the governing flow field equations, which computes ...


4

Including an actual ghost cell for FVM discretization is more or less a matter of convenience. For instance, if the boundary condition on the wall is a Neumann type, you actually do not need to discretize in the $x^-$ direction since the flux is given by the boundary condition and is known. If the boundary condition is of Dirichlet type, what you need is to ...


4

Typically you would use a slope limiter (or artificial diffusion and just cross your fingers) which detects where the solution has gone negative and modifies the solution to restore positivity (often by modifying the gradient of the solution in order to maintain conservation, at least in conservative schemes like Discontinuous Galerkin and Finite Volume). ...


4

What is usually done when you have a matrix-like variable due to the dimensions of your problem and the indexing is to linearize the index. I.e., if you have $C_{i,j}$ you would replace that with $\hat{C}_k=C_{i,j}$ where $k=i\times N+j$ and $N$ is the number of points on one side of your $N\times N$ grid of points that you are doing finite differences on. ...


4

I'm assuming the limitation in 2D and 3D is storing the Jacobian. One option is to retain the time derivatives and use an explicit "pseudo" time-stepping to iterate to steady state. Normally the CFL number you need for diffusive and reactive systems might get prohibitively small. You could try nonlinear multigrid (it's also called Full Approximation Storage ...


4

The solution you believe to be inaccurate is actually by far the more accurate one; you've simply plotted it in a very deceptive way. For $\nu=2$, the exact solution is actually no bigger than about $10^{-35}$ everywhere -- it's zero for all intents and purposes. Therefore the numerical solution is correct to 10 digits -- far better than the accuracy of ...


3

Use a first order upwind (for the convection component) and a second order central difference (for the diffusion component). So the end result would be equivalent to discretising the equation, $$ \frac{\partial u}{\partial t} = \frac{\partial \boldsymbol{v}}{\partial x} + D\frac{\partial^2 u}{\partial x^2} $$ So using the $\theta$-method you will end up ...


3

You need to rederive the term $$ \frac{d}{(h_j^+ + h_j^-)^2} $$ because this appears to be wrong. It must lead to factor of 4 when switching to a uniform grid because you've taken the sum of two things that are $O(h)$ and squared it. Perhaps if you spell out the derivation, we can help you find the error.


3

I was (still am) looking for good answers for this. I work with multi-level adaptive grids where I use some sort of criterion for refinement. Folks doing FEM enjoy, rather cheap (computationally), rigorous error estimates that they use as refinement criterion. For us doing FDM/FVM, I have not had luck finding any such estimates. In this context, if you want ...


3

The equation is $\partial_{t} \psi = \partial_{x} F(\phi)$ where $\psi = \partial_x \phi$ The time-integration can be done, e.g., by explicit time-stepping: $\psi^{j+1}_i = \psi^j_i + \frac{\tau}{2 dx} (F^j_{i+1}-F^j_{i-1})$, which requires that $F$ is known at each grid node $i$ at the "old" time slice $j$. This can be achieved by adding a numerical "...


Only top voted, non community-wiki answers of a minimum length are eligible