Fipy: Beer-Lambert law with multiple reflections

multipass_beer_lambert.py

# http://matforge.org/wd15/blog
from __future__ import division
from fipy import *

# Grid things
N = 1000
D = 1
dx = D / N
m = Grid1D(nx=N, dx=dx)
I_right = CellVariable(mesh=m, value=1., name='I_right') # Intensity along the positive direction
I_left  = CellVariable(mesh=m, value=1., name='I_left')  # Intensity along the positive direction

# Constants
A = 0.
B = 1.
C = (1 - A)
alpha = CellVariable(mesh=m, value=5., name="absorption coefficient")

# Constraints --- problem here!
I_right.constrain( A * I_left[0] + C, where=m.getFacesLeft())  # x=0 BC
I_left.constrain(  B * I_right[-1],   where=m.getFacesRight()) # x=D BC

# Equations to solve
eq1 =  ExponentialConvectionTerm(var=I_right) + ImplicitSourceTerm(alpha, var=I_right) == 0
eq2 = -ExponentialConvectionTerm(var=I_left)  + ImplicitSourceTerm(alpha, var=I_left)  == 0

# Beer-lambert (holds when A and B are small)
beer_lambert  = CellVariable(mesh=m, name='beer_lambert', value=numerix.exp(-alpha * m.x))

# Analytical expression for the differential equation above
analytical = CellVariable(mesh=m, name='analytical', value=-( B * numerix.exp( 2 * alpha * m.x ) + numerix.exp( 2 * alpha * D )  ) * C * numerix.exp( -alpha * m.x ) / ( A * B - numerix.exp( 2 * alpha * D ) ))

eq = eq1 & eq2
step = 0
while step < 10:
    print step
    res= eq.sweep()
    step += 1

vi = Viewer(vars=(beer_lambert, analytical, I_left + I_right))
#vi = Viewer(vars=(I_left + I_right))
vi.plot(filename='light_A=0,B=1.png')
raw_input()

There is an issue with setting boundary conditions correctly in this example -- don't use!