Diffusion equation

[navidcy]

Let's start with the diffusion equation in 1D

∂φ/∂t = κ ∂²φ/∂x²

where κ>0 is a constant and φ(x, t) a tracer concentration.

With no-flux boundary conditions at each end of the domain, i.e., ∂φ/∂x = 0 @ x=0, Lx we can prove the identity that the total φ is conserved. That is, ∫ φ(x, t) dx is constant with time.

With a centered difference scheme for both space and time, things should be accurate up to second order in dt, dx.

• [ ] Implement the boundary conditions with halo points so that they hold at the same accuracy as your spatial scheme in the interior of the domain.
• [ ] Repeat using your boundary condition implementation.
• [ ] Start with φ(x, t=0) which is a Gaussian centered somewhere round the middle of the domain. Make sure φ falls down -> 0 far before it reaches the boundary.
• [ ] Integrate forward up to certain time that φ develops non-zero amplitude at the boundaries.
• [ ] Track ∫ φ(x, t) dx and check how well is conserved. How does the error goes down with dt and dx with the different boundary condition implementations?

(This goes back to #5.)

[dhruvbhagtani]

This is completed for the 1D and the 2D case. Here are the notebooks: 1D and 2D. The 2D notebook isn't commented yet, so check that at your own risk! I'll comment it properly in the next couple of days.

[navidcy]

hey Dhruv, I'm all about taking risks but continuous test should help everybody avoiding taking risks :) Let's chat about how to do that. I know the ball was partly on my court and I have been busy...

I had a quick glance at the 2D notebook. Perhaps it would be useful if you plotted also the absolute differences from the analytical and numerical solutions?

[navidcy]

Also in the 1D notebook. I see you do the computation for how much the amplitude of the sine should be after the integrated time and then you plot something and say "look they match"? But how much accuracy can you get by eye-balling the figure?

The semiolog error plot vs number of points is, of course, quantitative.

Btw,

nx_inv = np.zeros(len(nx))
nx2_inv = np.zeros(len(nx))
for i,nxi in enumerate(nx):
    nx_inv[i] = 1/nxi
    nx2_inv[i] = 1/nxi**2

could be simplified to

nx_inv = 1 /nx
nx2_inv = 1 / nx**2

or you don't even need to define nx_inv and nx2_inv arrays but rather just use 1/nx or 1/nx**2 inside plot().

(This should be a remark to be made in a PR.)