Hi folks,

Can I just get a sanity check on this bit of code to check I'm not doing
something insane (I probably am) or stupid (equally likely).

The problem is a basic diffusion equation where a pile of sediment is
eroded and transport by diffusion only. The un-erodable topographic
surface (land) + sediment (phi) make up the land surface and the slope
over which the sediment is moved. The attached figure might show it more
clearly.

phi is therefore moved according to the 2nd order (phi+land). However,
I've also had to "turn off" diffusion by setting the diffusion
coefficient (nu) to zero where there is no sediment or volume is
definitely not conserved.

The equation is set up using:

tiny = 1e-10
surface.interpolate((((land+phi) + land)+abs((land+phi)-land)) / 2)
thickness.interpolate(surface-land)
limit.interpolate((surface - land)/(surface-land+tiny))

F = (inner((phi - phi_)/timestep, v)
      + limit*nu*inner(grad(phi+land), grad(v)))*dx

And then simply solved:

     solve(F == 0, phi)

in the time loop, updating limit and surface as we go.

Hopefully, the above is understandable, but if anyone has any opinions
on my (in)sanity, then I'm happy to hear them!

Cheers,
J.


--
I'm currently in Australia and replies will not be in standard UK
business hours.
Please do not feel you have to respond outside your normal work hours.

Dr Jon Hill
Senior Lecturer in Physical Geography
Department of Environment and Geography
University of York
M: +44(0)7748254812
P: +44(0)1904 324480
Twitter: @jonxhill
Web: https://jonxhill.wordpress.com/
Web: https://envmodellinggroup.github.io/