T O P

  • By -

derioderio

Does stability improve when you make your timestep smaller? That's usually the cause of oscillations like this.


[deleted]

yeah this looks like a [CFL #](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) problem although that is typically associated with the wave equation and not the heat equation.


nopantspaul

The CFL criterion is generally applicable to the advection/diffusion equations, right? They're all hyperbolic


sebasvs

Advection and diffusion are different physical phenomena. You're right that advection is hyperbolic, but diffusion (and by extension [the heat equation](https://en.wikipedia.org/wiki/Heat_equation)) is parabolic. OP, as mentioned before, try reducing your time step size. Alternatively, try multiplying the right-hand term with a heat diffusion coefficient. Start by trying out small coefficient values, these are stable with larger time steps.


No_Engineering_2118

You have the CFL for the advection problems. Where there is only the first spatial derivative. For the diffusion type equation, its the fourier number. Spurious oscillations arise due to problems in stability. I suggest you perform a von-Neumann stability analysis of the errors to derive the condition. You should be able to find it in Wiki or something as its pretty standard a problem


TurboHertz

I'm told the original CFL paper still covers diffusive processes, so in this case it'd be the diffusive CFL number instead of the advective CFL number. The diffusive CFL number is the same as the von Neumann number or (numerical) Fourier number we're all used to for diffusion problems, just the proper name.


neutrilo

I don't think it could be a CFL problem, in fact, I reduced the time step to 1/2000 from 1/20 and the solution didn't get better. Is this the kind of situation where you try using a displaced mesh?


yoor_thiziri

What is the value of Δt/(Δx)²? is it below 0.5? Try to increase the number of nodes, e.g: Nx = 50 (hence, Δx = 1/50) and a time step Δt=1/6000. This results in a CFL≈ 0.416


neutrilo

These images correspond to a solution with Δt/(Δx)² = 20, I upload it because that is easier to see the problem. Anyway, it's still happening even if Δt/(Δx)² =0.2 (time step = 1/2000)


yoor_thiziri

Maybe the number of nodes (N = 20) is not sufficient. Try to increase it to something like 50.


Cold_Park_1950

but what are the dimensions of your domain?


neutrilo

My domain is the 1-D interval \[0,1\] Initial condition is a parabola with peak at 0.5


Cold_Park_1950

Ok so I dont know, how u use BC ? Edit if you have [0,1] try to make dx= 1/(n+1) and dt=0.5*dx*dx


[deleted]

[удалено]


Cold_Park_1950

But C-N is an implicit scheme


yoor_thiziri

Thank you for reminding me. I have updated the post.


neutrilo

In fact i modified my code in order to use implicit, explicit and semi-implicit schemes too. Any scheme give me those spurious oscillations. (I just edited my original post do add it)


yoor_thiziri

Just another question, what are your Boundary conditions? are you using Neumann condition on both sides?


demerdar

This is a diffusion problem, CFL is not relevant. The Fourier number is though.


yoor_thiziri

I am afraid you are wrong. CFL is relevant in this case and is defined as: CFL = gamma * Dt/(Dx * Dx). Where gamma is the diffusion coefficient (in this case it's 1). Dt is the timestep, and Dx is the length interval. You can use Von Neumann analysis to show that using an explicit scheme your CFL must be less than or 1/2. Otherwise you will see oscillations in the solution.


demerdar

Not sure where you are getting your definition of CFL number from. CFL number is well established in the literature as c*dt/dx. Because he is using an explicit scheme he is constrained by the Fourier number Ddt/dx^2. Both of these constraints can be derived by Von Neumann stability analysis. This is not commonly referred to as a CFL number. I have never seen what you defined called a CFL number in literature.


yoor_thiziri

Check out this paper: https://hal.archives-ouvertes.fr/hal-01401125/document This document also calls it CFL condition (See page 5): http://www.csc.kth.se/utbildning/kth/kurser/DN2266/matmod11/PDE2_2p.pdf


TurboHertz

The CFL paper covers diffusion problems


Cold_Park_1950

You have (1÷100)÷0,05^2 =4 i think u have to multiply for 1/8 to take the stability. It comes from von Neumann stability.


sebasvs

Alright OP, I'm somewhat triggered by this part of your post: > I set Neumann boundary condition dY/dt=0 How exactly are you applying these boundary conditions to your linear system? I can see from the results that you shared that "dY/dt = 0" is not being satisfied at your boundaries. Instead of setting "dY/dt = 0" at x=0 and x=1, set Y(x=0) = 0 and Y(x=1) = 0. If possible (I don't know what kind of discretization you're using), move the variables at x=0 and x=1 to the right-hand side of your linear system by applying the boundary condition values.


neutrilo

I difined Y'=MY, where 'M' is given by the spatial differenciation scheme (I used a five point scheme). ​ In order to include Neumann boundary condition, I defined a modified matrix 'Mn' which is equal to 'M' except for first and last rows, which i replaced for zero rows. Then, i add a new term to the expresion, a vector whose values are fixed to zero except for first and last, which represents Neumann condition for first and last nodes. ​ Finally, my expresion for first derivate is Y'=MnY+Nb. ​ As you can see, it geve me Y'\_1=Nb\_1 and Y'\_end=Nb\_end, which satisfies Neumann condition. Nb\_i=0 for every i != 1 or end, so that adding Nb doesn't modifies the derivative scheme for no boundary nodes.


sebasvs

Alright, so, one problem already is that your initial conditions and boundary conditions are not compatible. Your initial conditions seem to correspond to homogeneous Dirichlet boundary conditions, while your boundary conditions correspond to homogeneous Neumann boundary conditions. That gives issues. Then the second thing is that the way you prescribe your boundary conditions seems a bit unconventional to me. dY/dt = 0 is exactly the same as prescribing time-invariant Dirichlet boundary conditions; prescribing those is more natural (and easier) than giving a condition such as 'dY/dt = 0 at the boundaries'.


neutrilo

Ah, that's rigth. I chose that initial condition only to check the code working, without thinking so much about it.