Manually changing CellVariable in a spatially-dependent advection-diffusion problem

[taoyun]

I'm able to apply a piece of old Stackoverflow code and solve for the transient solutions of a 1D Fokker-Planck equation, where the coefficients of advection and diffusion vary with space. The problem describes an animal's movement (its probabilistic location) as it moves across a landscape. The working code (based on FiPy 3.1) is as follows:

import fipy
import numpy as np
import matplotlib.pyplot as plt

nx = 100
dx = 0.5
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.getCellCenters()[0]

# set initial Gaussian location
distr = 10/(0.5*numerix.sqrt(2*numerix.pi)) * numerix.exp(-((x-5.) ** 2) / (2 * 0.5**2))
phi = CellVariable(name="solution variable", mesh=mesh, value=distr)

# zero-flux boundary conditions
BCs = (FixedFlux(faces=mesh.getFacesRight(),  value=0.),
       FixedFlux(faces=mesh.getFacesLeft(),   value=0.))

# visualization
viewer = Viewer(vars=(phi), limits={'ymin': 0., 'ymax': 6.})

# baseline coefficients
D = 10.
C = -1.
# the effect of space on baseline movement
omega = ((x-25)/50.)**2

eqX = TransientTerm() == ExplicitDiffusionTerm(CellVariable(mesh=mesh,value=D*(omega**2))) + \
                         ConvectionTerm(CellVariable(mesh=mesh,value=[C*omega]))

# time steps
timeStepDuration = 0.9 * dx**2 / (2 * 2.5)
steps = 1000
t = timeStepDuration * steps

# solve
plt.plot(x.value, D*omega, 'r-')
for step in range(steps):
    print 'step', step
    eqX.solve(var=phi, boundaryConditions=BCs, dt=timeStepDuration)
    viewer.plot()

Question 1:
If I wish to manually change all omega.value < 0.05 to 0.05, how would I do that? omega.setValue(0.05, where=omega.value < 0.05 gives an TypeError.

Quesion 2:
I noticed that the above code does not specify the variable to solve in the equation. In fact, running the line

eqX = TransientTerm() == ExplicitDiffusionTerm(CellVariable(mesh=mesh,value=D*(omega**2)), var=phi) + \
                         ConvectionTerm(CellVariable(mesh=mesh,value=[C*omega]), var=phi)

gives an ExplicitVariableError. This is problematic because I'm not sure how one would model two moving, potentially interacting, individuals, phi1 and phi2, at once. Any advice would be much appreciated!

[wd15]

The code as written doesn't work, obviously. Grid1D isn't defined for example. Also, use a later version of FiPy. I'm using 3.4.3.

Question 1: It's in the attached code and it works fine. Might have something to do with your older version of FiPy

Question 2: I just think you need to add var=phi into the TransientTerm and omit it from solve. See code below.

import fipy
import numpy as np
import matplotlib.pyplot as plt
from fipy import Grid1D, numerix, CellVariable, FixedFlux, Viewer, TransientTerm, ExplicitDiffusionTerm, ConvectionTerm

nx = 100
dx = 0.5
mesh = Grid1D(nx=nx, dx=dx)
x, = mesh.cellCenters

# set initial Gaussian location
distr = 10/(0.5*numerix.sqrt(2*numerix.pi)) * numerix.exp(-((x-5.) ** 2) / (2 * 0.5**2))
phi = CellVariable(name="solution variable", mesh=mesh, value=distr)

# zero-flux boundary conditions
BCs = (FixedFlux(faces=mesh.facesRight,  value=0.),
       FixedFlux(faces=mesh.facesLeft,   value=0.))

# visualization
viewer = Viewer(vars=(phi), limits={'ymin': 0., 'ymax': 6.})

# baseline coefficients
D = 10.
C = -1.
# the effect of space on baseline movement
omega = ((x-25)/50.)**2

eqX = TransientTerm(var=phi) == ExplicitDiffusionTerm(CellVariable(mesh=mesh,value=D*(omega**2)), var=phi) + \
                         ConvectionTerm(CellVariable(mesh=mesh,value=[C*omega]), var=phi)

# time steps
timeStepDuration = 0.9 * dx**2 / (2 * 2.5)
steps = 1000
t = timeStepDuration * steps

# solve
plt.plot(x, D*omega, 'r-')
for step in range(steps):
    print('step', step)
    phi.setValue(0.05, where=omega < 0.05)
    eqX.solve(boundaryConditions=BCs, dt=timeStepDuration)
    viewer.plot()