Evolutive advection coefficient with time in FiPy, advection - diffusion equation (2D)

Question

I am trying to develop a model that represent the stork migration from France to middle Africa. I chose to go with FiPy.

I would like to set the time as a full migration time : First, storks diffuse randomly in France, but when the time gets to a limit value, migration begins and an advection coefficient appear.

I tried first with the Diffusion coefficient and it works :

# Model parameters


D = Variable(value=(0))
ux = 100
uy = -150
q  = 0



Lx, Ly = 2500, 2500
Nx, Ny = 100, 100
dx, dy = Lx/Nx, Ly/Ny
Cx, Cy = Lx/4, 3*Ly/4
sigma  = Lx/10
T  = 10
dt = T/1000
Np = 10     # Number of snapshots of the solution (one every Nt/Np time steps)
Nt = Np*np.ceil(T/(Np*dt)).astype('int')   # Nt must be an integer divisible by Np



# Define the grid/mesh
mesh = Grid2D(nx=Nx, ny=Ny, dx=dx, dy=dy)
x, y = mesh.cellCenters[0], mesh.cellCenters[1]

# Define the model variable and set the boundary conditions
phi = CellVariable(name="numerical solution", mesh=mesh, value=numerix.exp(-((x-Cx)**2 + (y-Cy)**2)/(2*sigma**2)))
meshBnd = mesh.facesLeft | mesh.facesRight | mesh.facesBottom | mesh.facesTop





# impose zero flux on all 4 boundaries
phi.faceGrad.constrain(0., where=meshBnd)

# Define and then solve the equation
eq = TransientTerm() == DiffusionTerm(coeff=D) \
    - ExponentialConvectionTerm(coeff=(ux,uy)) \
    + ImplicitSourceTerm(coeff=q)
my_sol = np.zeros((Np,Nx*Ny))      # Matrix with Np solution snapshots
my_sol[0,:] = phi
k = 1

for step in np.arange(1,Nt):
    if step < Nt/2:
        print(step)
        D.setValue(0)
    else:
        D.setValue(5000) # setting D=5000 
    eq.solve(var=phi, dt=dt)
    if np.mod(step,Nt/Np)==0:
        print(step,k)
        my_sol[k,:] = phi
        k += 1

Here are the results :

First five snap shots :

Then the diffusion coefficient change :

but as I tried the same manip, there is no change and the ux uy coefficients stay the same.

# In model parameters :
D = Variable(value=(0))
ux = Variable(value=0)
uy = Variable(value=0)
q  = 0

# In equation resolution loop :

for step in np.arange(1,Nt):
    if step < Nt/2:
        print(step)
        D.setValue(0)
    else:
        D.setValue(5000)
        ux.setValue(100)
        uy.setValue(-150)
    eq.solve(var=phi, dt=dt)
    if np.mod(step,Nt/Np)==0:
        print(step,k)
        my_sol[k,:] = phi
        k += 1

But there is no change in the advection...

Thanks for your help

Answer 1

These:

ux = Variable(value=0)
uy = Variable(value=0)

are FiPy Variable objects, but this

(ux,uy)

is not. FiPy doesn't recognize the tuple and casts it to a NumPy array, so that it loses connection to any changes to ux or uy.

This should work:

coeff=(ux * FaceVariable(mesh=mesh, value=[1,0]) 
       + uy * FaceVariable(mesh=mesh, value=[0,1]))

Explanation:

• A Variable represents a value of arbitrary shape, but with no spatial relationship.
• A FaceVariable represents a field, of any rank, defined on the faces of a Mesh.
• Multiplying by FaceVariable(mesh=mesh, value=[1,0]) and FaceVariable(mesh=mesh, value=[0,1]) transforms the scalar values ux and uy into the x and y components of a rank-1 field.
• In principle, it should be possible to write ux * [1,0] + uy * [0,1] or ux * [[1],[0]] + uy * [[0],[1]] to create a vector Variable that isn't a field, but ExponentialConvectionTerm does not recognize changes in this quantity; I'm not sure why.