I am trying to solve the 1-dimensional diffusion problem numerically using method of lines:

$$ \frac{\partial c}{\partial t} =D \frac{\partial^2 c}{\partial z^2},$$

where the right hand side is discretised and the left hand side is handled with a standard ode solver.

I have two boundary conditions. One boundary condition dynamically changes with every timestep (dcdt[0]), while the other is always fixed at a constant value (dcdt[-1]). There is also one initial condition (concentration profile) at time zero, which is the line in dark blue in the picture shown.

Unfortunately, what happens is that scipy.integrate.ode ignores both boundary conditions!


def diffusion(t,c, h):
    dcdt = zeros(Nz)
    dcdt[0] = funcA(h)
    dcdt[-1] = B # a constant    
    for i in range(1,Nz-1):
        dcdt[i] = D * (c[i+1] - 2*c[i] + c[i-1]) / (dz*dz) 

    return dcdt

I use scipy.integrate.ode to integrate the discretisation equation. I have defined the diffusion function to set the boundary condition every time the ode(diffusion) is called.

ode = integrate.ode(diffusion).set_integrator('vode')
ode.set_initial_value(c_initial, 0).set_f_params(h0)
k = 0

while ode.successful() and k < Nt:
    dcdz = (y[1] - y[0]) / dz
    h[k] = h[k-1] + Q*(h[k-1]**2) / func(h[k-1]) * dt * dcdz
    c_out[:,k] = ode.y
    t[k] = ode.t

    ode.set_f_params(h[k]) # Change the b.c.
    ode.integrate(ode.t + dt)
    k += 1

I have checked the following things:

  • By inserting print h into the diffusion function I have checked that set_f_params has successfully passed h into the solver.

  • There are no numerical instabilities. I have also checked other solvers dopri5, stiff/non-stiff options, etc.

  • I have switched the order of functions in the while loop around. It seems that the problem is internal, i.e. ode.integrate(ode.t + dt) solves the equation and ignores the boundary conditions I want to impose before returning the output to me. I'm puzzled why it would do that.

  • One thing I have thought about is that I could use set_initial_value to set the boundary condition. But as far as I know, I have to provide a full solution vector $c(z,t)$ at a certain time, and ode.integrate is already ignoring my boundary condition at time zero.

  • 1
    $\begingroup$ Have you tried only solving for the interior nodes (i.e., don't solve the ode at $z=0$ and $z=1$ and use the boundary conditions as explicit values in the finite difference)? $\endgroup$ – Christian Clason May 19 '15 at 18:04
  • 2
    $\begingroup$ I may be missing something obvious, but wouldn't $\partial c(t,1)/\partial t=B$ lead to $c(t,1) = c(0,1) + Bt$, increasing linearly, just as in the plot? Is $B=0$ or $\partial B/\partial t=0$? Do you get the right boundary conditions if you set $D=0$, in effect solving two uncoupled ODEs? $\endgroup$ – Kirill May 19 '15 at 18:20
  • 1
    $\begingroup$ I would check running your code, but there is not enough pieces of it to run. $\endgroup$ – nicoguaro May 19 '15 at 18:43
  • 3
    $\begingroup$ Please write down the precise boundary conditions. $\endgroup$ – David Ketcheson May 20 '15 at 6:01