[xammm]

Dear forum users!

I'm currently implementing a reaction-advection-diffusion system on a closed 2D-surface (surface of a sphere) using C. First I started using finite elements which works just fine for the diffusion, however I could not get the advection to work. After doing some research on the internet I've came to the conclusion that it's better to use finite volume for the advection part instead.

So I've now implemented a vertex-centered finite volume routine with upwind scheme for the advection. However it's still unstable! My test is a single Gaussian pulse which is supposed to be advected around the sphere.

Here's a typical plot I get:



which is a cut at theta = pi/2. phi goes from 0 to 2*pi. The red curve is the initial condition and you can clearly see that this is unstable even though diffusion is still pretty strong. I would like to have a code, where the pulse gets advected by (e.g.) one rotation and has still the same shape. There's also a video:

http://www.xammm.com/surf.mpeg

I'm using a Crank-Nickolson/Cholesky decomposition scheme.

Is there an easy way to stabilize this?

[FMDenaro]

what about the cfl number at which your code was tested?

[cfdnewbie]

The mesh you are using looks indeed fine enough to rule out major dispersion and dissipation from the spatial scheme.... I assume it is a first order approximation (constants in cells?) .

As Filippo pointed out, the temporal scheme could be the issue. Whatever your CFL number, reduce it by a factor of 10, rerun the simulation and check the results.

Just out of interest: What flux function are you using?

[FMDenaro]

Originally Posted by cfdnewbie

The mesh you are using looks indeed fine enough to rule out major dispersion and dissipation from the spatial scheme.... I assume it is a first order approximation (constants in cells?) .

As Filippo pointed out, the temporal scheme could be the issue. Whatever your CFL number, reduce it by a factor of 10, rerun the simulation and check the results.

Just out of interest: What flux function are you using?

If I have well understood, the flux function should be piece-wise constant ...

[Rami]

Hi xammm,

1. First, your solution exhibits over- and under- shoots already at t=1. It indicates your convection treatment is inefficient. Try using a 1st order upwind as a starting point to remove that behavior. Only then enhance it to higher order.
2. I suggest you start by pure convection, so that you know what the expected solution is.
3. I also suggest to start with a rectangular domain in a Cartesian coordinate system with periodic BCs, to rule out possible bugs arising from the spherical coordinates treatment.
4. I join the recommendation to decrease the CFL No.
5. A comment: In convection-dominated problems FEM should work fine, but needs more effort. FVM is simpler to formulate (but less consistent).

Best luck,

Rami

[FMDenaro]

Originally Posted by Rami

Hi xammm,

1. First, your solution exhibits over- and under- shoots already at t=1. It indicates your convection treatment is inefficient. Try using a 1st order upwind as a starting point to remove that behavior. Only then enhance it to higher order.
2. I suggest you start by pure convection, so that you know what the expected solution is.
3. I also suggest to start with a rectangular domain in a Cartesian coordinate system with periodic BCs, to rule out possible bugs arising from the spherical coordinates treatment.
4. I join the recommendation to decrease the CFL No.
5. A comment: In convection-dominated problems FEM should work fine, but needs more effort. FVM is simpler to formulate (but less consistent).

Best luck,

Rami

What do you mean?
Perhaps, owing to the second order time integration, is not strange that the solution shows over/under-estimation even with first-order upwind ..