Boundary conditions

This ex­am­ple shows so­lu­tions to time-de­pen­dent dif­fu­sion equa­tions for­mu­lat­ed with dif­fer­ent bound­ary con­di­tions. Iso­lat­ed, Dirich­let (pre­scribed val­ue) and pe­ri­od­ic bound­ary con­di­tions are test­ed.

%matplotlib inline

import matplotlib.pylab as plt
from oedes import *
init_notebook()
from matplotlib import ticker

Be­low trans­port mod­el al­low­ing on­ly dif­fu­sion is im­ple­ment­ed. Dif­fu­sion co­ef­fi­cient is tak­en from pa­ram­e­ters, where it is spec­i­fied un­der key *.D (* de­notes the equa­tion)

def v_D_diffusion_only(ctx, eq):
    return 0., ctx.param(eq, 'D')

The func­tion be­low cre­ates and solves a mod­el. It calls sup­plied setup_bc to ini­tial­ize bound­ary con­di­tions.

def run(species_D, species_ic, setup_bc, L=1., t=10, dt=1e-9):
    model = models.BaseModel()
    mesh = fvm.mesh1d(L)

    # BaseModel requires presence of Poisson's equation
    # The species are assumed to be uncharged, therefore it is decoupled from the system
    # and has no effect
    model.poisson = Poisson(mesh)
    model.poisson.bc = [models.AppliedVoltage(boundary) for boundary in mesh.boundaries]

    # Create equations
    params=dict()
    for i,D in enumerate(species_D):
        k='species%d'%i
        species = models.AdvectionDiffusion(mesh, k, z=0, v_D=v_D_diffusion_only) # uncharged: z=0
        model.species.append(species)
        params[k+'.D']=D
    setup_bc(model)
    model.setUp()

    # These parameters are irrelevant for this test, but still are required to
    # evaluate the model
    params.update({'T': 300., 'electrode0.voltage': 0, 'electrode1.voltage': 0, 'electrode0.workfunction': 0,
              'electrode1.workfunction': 0, 'epsilon_r': 3.})

    # Create initial conditions
    xinit = model.X.copy()
    for species,ic in zip(model.species,species_ic):
        xinit[species.idx] = ic(mesh.cells['center'] / L)

    # Run the simulation and return the results
    c = context(model, x=xinit)
    c.transient(params,t,dt)
    return c

This func­tion makes plots of con­cen­tra­tions at dif­fer­ent times:

def plot_times(c,times,shrink=[0.5,0.5]):
    time_formatter = ticker.EngFormatter('s')
    times=np.asarray(times)
    figsize = np.asarray(plt.rcParams['figure.figsize'])*np.asarray(times.shape[::-1])*shrink
    fig,axes=plt.subplots(nrows=times.shape[0],ncols=times.shape[1],sharex=True,figsize=figsize)
    axes=np.asarray(axes).ravel()
    times=times.ravel()
    for ax,t in zip(axes,times):
        ct=c.attime(t)
        mpl=ct.mpl(fig,ax)
        mpl.allspecies()
        ax.set_yscale('linear')
        ax.set_title('t = %s'%time_formatter(t))
        # testing support
        for species in ct.model.species:
            testing.store(ct.output()[species.name+'.c'],rtol=1e-5)
    fig.tight_layout()

De­fault dif­fu­sion co­ef­fi­cients and ini­tial con­di­tions:

test_mu = [ 1,10,100]
test_ic = [
    lambda u: np.abs(u - 0.5) < 0.1,
    lambda u: np.abs(u - 0.2) < 0.1,
    lambda u: np.where(u > 0.999, 1e2, 0)
]

No-flux (isolated) boundary conditions

def setup_bc_isolated(model):
    pass
c=run(test_mu,test_ic,setup_bc_isolated)
plot_times(c,[[1e-7,1e-5],[1e-3,1e-1],[1,1e1]])