# An Advection-Diffusion-Reaction System: Chemical Reaction Problem

[1]:

try:
import dune.femdg
except ImportError:
print("This example needs 'dune.femdg' - skipping")
import sys
sys.exit(0)

import numpy
from matplotlib import pyplot
from ufl import *
from dune.ufl import DirichletBC
from dune.grid import structuredGrid
from dune.fem.plotting import plotPointData as plot
from dune.fem.space import lagrange, dgonb
from dune.fem.scheme import galerkin
from dune.femdg import femDGOperator
from dune.femdg.rk import femdgStepper

gridView = structuredGrid([0,0],[2*numpy.pi,2*numpy.pi],[20,20])


Method to compute the velocity based on a scalar elliptic equation

[2]:

def computeVelocity():
streamSpace = lagrange(gridView, order=1, dimRange=1)
Psi  = streamSpace.interpolate(0,name="streamFunction")
u,v  = TrialFunction(streamSpace), TestFunction(streamSpace)
x    = SpatialCoordinate(streamSpace)
streamScheme = galerkin([form == 0, DirichletBC(streamSpace,[0]) ])
streamScheme.solve(target=Psi)
return as_vector([-Psi[0].dx(1),Psi[0].dx(0)])


Chemical reaction model

[3]:

class Model:
transportVelocity = computeVelocity()
def S_e(t,x,U,DU):
P1 = as_vector([0.2*pi,0.2*pi]) # midpoint of first source
P2 = as_vector([1.8*pi,1.8*pi]) # midpoint of second source
f1 = conditional(dot(x-P1,x-P1) < 0.2, 1, 0)
f2 = conditional(dot(x-P2,x-P2) < 0.2, 1, 0)
f  = conditional(t<5, as_vector([f1,f2,0]), as_vector([0,0,0]))
r = 10*as_vector([U[0]*U[1], U[0]*U[1], -2*U[0]*U[1]])
return f - r
def F_c(t,x,U):
return as_matrix([ [*(Model.velocity(t,x,U)*u)] for u in U ])
def maxWaveSpeed(t,x,U,n):
return abs(dot(Model.velocity(t,x,U),n))
def velocity(t,x,U):
return Model.transportVelocity
def F_v(t,x,U,DU):
return 0.02*DU
def physical(t,x,U):
return conditional(U[0]>=0,1,0)*conditional(U[1]>=0,1,0)*\
conditional(U[2]>=0,1,0)
boundary = {range(1,5): as_vector([0,0,0])}

fig = pyplot.figure(figsize=(20,20))
plot(Model.transportVelocity, gridView=gridView, gridLines=None, vectors=[0,1],figure=fig)


Time evolution scheme

[4]:

space    = dgonb( gridView, order=3, dimRange=3)
u_h      = space.interpolate([0,0,0], name='u_h')
operator = femDGOperator(Model, space, limiter="scaling")
stepper  = femdgStepper(order=3, operator=operator)
operator.applyLimiter(u_h)

t        = 0
saveTime = 1.5

fig = pyplot.figure(figsize=(30,20))
c = 0
while t < 9:
operator.setTime(t)
t += stepper(u_h)
if t > saveTime:
saveTime += 1.5
u_h.plot(gridLines=None, figure=(fig, 231+c), colorbar=False)
c += 1

femDGOperator: Limiter = scaling


All three components

[5]:

from dune.fem.plotting import plotComponents
from matplotlib import ticker
plotComponents(u_h, gridLines=None, level=1,
colorbar={"orientation":"horizontal", "ticks":ticker.MaxNLocator(nbins=4)})


This page was generated from the notebook chemical_nb.ipynb and is part of the tutorial for the dune-fem python bindings