Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Pages
Diffusion Advection Reaction Problem

Table of Contents

Author
Feel++ Consortium
Date
2013-02-25



The diffusion advection reaction equation is a classical partial differential equation which can be found in many processes for example in chemistry or biology. This can be described by an equation containing a diffusion, an advection and a reaction term as follows,

\( \left\{ \begin{aligned} -\epsilon\Delta u + \bbeta \cdot \nabla u + \mu u & = f & \text{on}\; \Omega \;, \\ u & = 0 & \text{on}\; \partial\Omega \;, \\ \end{aligned} \right \)


We use here homogeneous Dirichlet boundary conditions.

Theory

To establish the variationnal formulation, as always we mutiply the first equation by a test function \(v\in H_0^1(\Omega)\) such that,

\[ H_0^1(\Omega) = \{ v\in H^1(\Omega),\; v=0 \; \text{on} \; \partial\Omega \} \;. \]

Then we integrate on the domain \(\Omega\),

\( \begin{aligned} - \int_\Omega \epsilon \Delta u\ v + \int_\Omega \bbeta \cdot \nabla u\ v + \int_\Omega \mu\ u\ v = \int_\Omega f\ v \;. \end{aligned} \)


We establish the variationnal formulation from the previous equation and using the Green formula, find \(u \in \in H_0^1(\Omega)\)

\( \begin{aligned} \int_\Omega \epsilon \nabla u \cdot \nabla v - \underbrace{\int_{\partial\Omega} \epsilon (\nabla u \cdot \mathbf n)\ v}_{=0} + \int_\Omega (\beta \cdot \nabla u)\ v + \int_\Omega \mu\ u\ v = \int_\Omega f v \; \quad \forall v \in H_0^1(\Omega), \end{aligned} \)


where \(\mathbf n\) is a unit outward normal vector. We can rewrite the problem, find \(u \in \in H_0^1(\Omega)\)

\[ a(u,v) = l(v) \quad \forall v \in H_0^1(\Omega), \]

where \(a\) is a bilinear form, continuous, coercive and \(l\) is a linear form.

Implementation

We choose for our example \(\mu = 1 \), \(\epsilon = 1\), \(f=1\), and \(\bbeta=(1,1)^T\).

int
main( int argc, char** argv )
{
po::options_description opts ( "Advection diffusion reaction options ");
opts.add_options()
( "epsilon", po::value<double>()->default_value( 1 ), "diffusion term coefficient" )
( "betax", po::value<double>()->default_value( 1 ), "convection term coefficient in x-direction" )
( "betay", po::value<double>()->default_value( 1 ), "convection term coefficient in y-direction" )
( "mu", po::value<double>()->default_value( 1 ), "reaction term coefficient" );
// Initialize Feel++ Environment
Environment env( _argc=argc, _argv=argv,
_desc=opts,
_about=about(_name="myadvection",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org") );
// create mesh
auto mesh = unitSquare();
// function space
auto Xh = Pch<1>( mesh );
auto u = Xh->element( "u" );
auto v = Xh->element( "v" );
// diffusion coeff.
double epsilon = option(_name="epsilon").as<double>();
// reaction coeff.
double mu = option(_name="mu").as<double>();
auto beta = vec( cst(option(_name="betax").as<double>()),
cst(option(_name="betay").as<double>()) );
auto f = cst(1.);
// left hand side
auto a = form2( _test=Xh, _trial=Xh );
a += integrate( _range=elements( mesh ),
_expr=( epsilon*gradt( u )*trans( grad( v ) )
+ ( gradt( u )*beta )*id(v)
+ mu*idt( u )*id( v ) ) );
// right hand side
auto l = form1( _test=Xh );
l+= integrate( _range=elements( mesh ), _expr=f*id( v ) );
// boundary condition
a += on( _range=boundaryfaces( mesh ), _rhs=l, _element=u,
_expr=cst(0.) );
// solve the system
a.solve( _rhs=l, _solution=u );
// export results
auto e = exporter( _mesh=mesh );
e->add("u",u);
e->save();
} // end main

Again the implementation is close to the mathematical formulation.
Here again, we create the mesh for an unit square geometry.
Then we define the function space \(X_h\) we choose as order 1 Lagrange basis function using Pch<Order>(). Note that here, the function space is the same for "trial" and "test" functions.
We declare the left and the right hand side integrals expressions for the equation (Advection_Math ).
Finally we add the Dirichlet boundary condition and we use the default solver to solve. We export the solution \(u\) for post processing.

Results

There are various solutions of this problem displayed with Paraview for \(\epsilon= 1, 0.01, 0.0001\).
Notice how the solution gets unstable for \(\epsilon = 0.0001\), this is classical and requires stabilisation methods to handle this issue.

sol-dar-1.png
sol-dar-2.png
sol-dar-3.png
\(\epsilon=1\)
\(\epsilon=0.01\)
\(\epsilon=0.0001\)