41

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 ...


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, ...


7

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 ...


6

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


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 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

There are basically two methods: You can disrectize the angular part via grid points, or you can discretize it via basis expansion. I will focus on spherical symmetry here, the cylindrical case is quite similar. In the basis expansion approach, one applies the ansatz $$ v(x,t) = \sum_{klm} a_{klm}(t)\,R_{klm}(r)\, Y_{lm}(\theta,\phi) $$ This is inserted ...


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 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 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

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

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

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

So after browsing the paper a bit, I think that the answer is essentially what Christian Clason stated in his comment. It seems that the original question refers to the statement just above Equation (3) in the article linked by kwesi : There, the authors say that the advection-diffusion equation (Equation (1) in the paper) $\frac{\partial c}{\partial t} + \...


3

When you discretize a PDE and your discretization blows up, the correct response is not "let's ignore that result and try a different discretization". That's like trying to fix your car by randomly replacing parts. Instead, you should stop and investigate until you understand what is wrong, then correct that. I see three possibilities here: Your code has ...


3

Providing a whole detailed solution is out of the scope of this site, but asking for references is on-topic, so here is what I would suggest to get started: There are many good books on finite difference methods (for instance, "Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems" by LeVeque, or ...


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

Particles in SPH simulations mimic the motion of material points by simply updating their positions due to their velocities. In situations like heat conduction you don't need to move them, in fact, there is no advection at all. SPH is a collocation scheme like FDM (finite difference method) but without the rigid interpretation of nodes over a grid. The SPH ...


3

I would like to add that besides the spatial discretization the temporal discretization can also introduce numerical diffusion. Consider the advection equation $$ u_t + cu_x = 0$$ where I use the subscripts for partial derivatives: $ u_i \equiv \frac{\partial u}{\partial i} $. Applying a general $\theta$-scheme to discretize the temporal derivative yields ...


3

In the time-dependent nonlinear case, if you drop the diffusive term then you have a nonlinear hyperbolic problem. Solutions will naturally generate singularities (discontinuities) in finite time. To extend the solution beyond that time, one must consider weak solutions, and uniqueness is lost. To specify a unique, physically relevant solution one ...


3

You just add the diffusion along the other dimensions. This superposition from orthogonal directions makes some sense, as they are independent. So, going by wikipedia for Fick's second law of diffusion in 1D: $$ \frac{\partial \psi}{\partial t} = D \frac{\partial^2 \psi}{\partial x^2} $$ We extend it to 2d as: $$ \frac{\partial \psi}{\partial t} = D \...


3

Essentially, the time dependent Stokes equation looks like the heat equation: $$ \frac{\partial u}{\partial t} - \nu\Delta u = f-\nabla p, $$ plus the incompressibility condition $\nabla \cdot u=0$ that for the current discussion is immaterial. Thus, the same considerations for time step choice apply as for the heat equation. Consequently, using an ...


3

In case you can easily implement the roated grid, this is possibly the easist fix, so let me answer your question about the invariance of the convection diffusion equation under rotation of the coordinate system. TL;DR Your physical problem does not depend on the coordinate system that you choose. Only the definition of the Laplace operator and the ...


3

First, write out the semi-discrete equation for $u_0$, assuming it's an interior node: $ \frac{d}{dt}(u_0) = - c \frac{\left( u_1 - u_{-1} \right)}{2h} + \alpha \frac{\left(u_1 - 2u_0 + u_{-1}\right)}{h^2}$ Then, eliminate the ghost node $u_{-1}$ by using the equation you obtained from the Neumann BC: $u_{-1} = u_1$. Substituting this condition into the ...


2

Heat is a passive scalar (under assumptions of incompressibility). But there is no heat transfer from boundaries without diffusion. So you can't really do any such heat transfer problems in a purely Lagrangian framework (unless you go stochastic). Same thing holds for permeable boundaries where the scalar diffuses into the domain.


2

A lot of this will depend on the magnitude of $D$ and mesh size $h$. For pure convection, upwinding works well; however, upwinding ends up too diffusive for $D \approx h$. If $h \ll D$, then upwinding should behave less diffusively, though only will $D=0$ will the upwind scheme show minimal numerical diffusivity. As $D$ increases relative to mesh size, an ...


2

You should never use explicit method for the diffusion equation. Implicit is unconditionally stable and just as easy to implement. Also if you use an implicit method (like backward Euler or Crank-Nicolson) it will not matter how small d is. In fact you could use a dirac delta function if you wanted. As far as speed goes, doing a implicit method with a ...


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