Setting a time-dependent advection coefficient with FiPy

[Zan-IM]

Hi,

I am currently trying to solve a simple advection-diffusion equation on a 1D grid and I would like the advection coefficient to change over time.

My equation is: with no-flux boundary conditions.

To do so, I defined the time t as a FiPy Variable and I then write the advection coefficient u as a function of t. Right now, I am trying to make u exponentially decrease with time, so that, after some timesteps, I would expect the process to become almost solely diffusive. However, my FiPy code solves the equation ignoring the time dependence of u and seems to retain its initial value for all timesteps.

My cleaned code right now is:

import fipy as fi 
from fipy import *
import numpy as np

#### 1. Define the 1D Grid
nx = 100
dx = 1.
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters[0]              

#### 2. Define time as a FiPy Variable
t = Variable(value = 0.)

#### 3. Diffusion coefficient
D = 0.1

#### 4. Set the timestep 
timeStepDuration = 0.9 * dx**2 / (2 * D)

####  5. Set a time dependent advection coefficient 
u = 0.1*(np.exp(-t/(timeStepDuration)))

#### 6.  Set initial condition for the variable phi 
phi = CellVariable(name="solution variable",
                   mesh=mesh,
                   value=0., hasOld = True)
phi.setValue(10., where =  (35>x) & (25<x)) 

#### 7. Set number of iteration steps 
steps = 500

####  8. Write the equation to solve
eqX = TransientTerm(var = phi) == ExplicitDiffusionTerm(coeff=D, var = phi) -  ConvectionTerm(coeff = [u, ], var = phi )

#### 9. Solve the equation
for step in range(steps):
    phi.updateOld()       # update phi value 
    eqX.sweep(dt = timeStepDuration, var = phi)  #  solve the equation for each timestep 
    t.setValue(t + timeStepDuration)             # update the time 

If I save the values of the advection coefficient at each timestep, I do see that the value of u is decreasing in time, but it seems that the value which equation eqX is using remains the same.

Thank you to anybody who's willing to help! I am struggling to understand the error.

Best regards, Medea

[guyer]

I thought I knew what was wrong here (and I wasn't altogether incorrect), but this turns out to be subtle. Thanks for reporting it.

It turns out that, when you pass a scalar Variable to a ConvectionTerm, FiPy turns this into a FaceVariable, but throws away the "variableness". You can work around this by doing

####  5. Set a time dependent advection coefficient 
dir = FaceVariable(mesh=mesh, value=[[1.]], rank=1)
u = 0.1*(numerix.exp(-t/(timeStepDuration))) * dir

Beyond that, I recommend using:

• a DiffusionTerm instead of ExplicitDiffusionTerm. You can run larger time steps without worrying about stability; you're already on the edge of stability when the problem starts.
• VanLeerConvectionTerm instead of ConvectionTerm. This will do a better job of preserving the shape of the pulse.

[Zan-IM]

Thank you for solving the issue and for the recommendations. Now, everything seems to work perfectly!

[guyer]

I'm glad it works. I'm going to leave this issue open, because I think FiPy should handle this more elegantly.

[guyer]

Oh, and I forgot to mention, you should call

ConvectionTerm(coeff = u, var = phi )

[u,] is a list containing a Variable, which FiPy does not recognize as a rank-1 Variable. It may work in 1D, but it is not generally what you want. This is originally why I thought that the coefficient wasn't updating.

[Zan-IM]

Somehow I missed this last comment. I already realized that using [u,] didn't work, but I didn't know why. Thank you very much for the explanation!