Minimal diffusion - advection example

[domagojhack]

Hi,
I am trying to solve a simple 2D advection-diffusion equation using py-pde

The equation is:

I understand how to implement first part of the equation with:

eq = pde.PDE({"C": f"{D} * laplace(C)"})

But cannot find any description or examples that deal with vector fields to implement the divergence part of the equation.

I am asking for pointers on how to implement the advection part in py-pde, if possible. I assume I need to write my own PDEBase class but I can't find a comparable example in docs.

Best regards.

[david-zwicker]

Using vectors in pde.PDE is certainly possible, but it depends a bit on what you want to do. First, where does your vector V come from? Is it a constant value or a constant field? Also, boundary conditions might matter and you will need to specify them. Finally, advection problems are notoriously difficult to solve and you might run into instabilities – we haven't played with this very much, yet.

[domagojhack]

Thanks for the fast reply.

For sake of simplicity V is a constant field. The system is open so there should not be some special boundary conditions.

By putting the formula in the context, C is the concentration of something diffusing in the air while V is the vector field describing wind in the domain.
---edit:
In discussions, I found one of your examples (Fokker-Planck diffusion equation #166)

eq = pde.PDE({'P': f'-divergence(P * V) + {D}/2 * laplace(P)'}, consts={'V': velocity})

This is something similar to my issue. Now the next question is how can I update constant V in every time step? V is not time-dependent but I want to change it manually.

[david-zwicker]

Here is a simple example that advects a Gaussian spot:

grid = pde.CartesianGrid([[-10, 10], [-10, 10]], 32)
field = pde.ScalarField.from_expression(grid, 'exp(-x**2-y**2)')
velocity = pde.VectorField.from_expression(grid, ['1', '1'])
eq = pde.PDE({'c': 'laplace(c) - dot(velocity, gradient(c))'}, consts={'velocity': velocity})
result = eq.solve(field, t_range=5)
result.plot()

The initial fiel,
,
gets turned into

[domagojhack]

Thank you for your example,

I am adding stochasticity across the to the V VectorField.

Currently, I am doing that by running my simulation like the following example:

for t in range(sim_len):
    V = pde.VectorField(grid, GetWind(angle, speed))
    eq = pde.PDE(
    {"C": f"{D} * laplace(C) - divergence(C*V)"},
    consts={"V": V, "Cs": Cs_})

    sol = eq.solve(state, 1, dt=1e-3)
    # plotting and storing code ....

where I change my wind vector by some function when I add stochasticity.

I am getting the simulation that I want but at the cost of the performance. The reason is, I am calling the solve function in every iteration so I can update my V. Is there a way to update my V in every time step without the for loop?

[david-zwicker]

There is no mechanism to update the constants in the pde.PDE class since they are baked into the compiled rhs for speed. If you need more flexibility it might be best to implement a custom PDE class, where you have full flexibility on specifying the right hand side. However, for optimal speed you would need to implement the numba compiled right hand side by overwriting _make_pde_rhs_numba; see the documentation for initial help.

[domagojhack]

Thank you for your reply, I am going to definitely try to make a numba compatible custom PDE class.

Cheers