Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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.add( feel_options() ),
_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$

Generated on Sun Dec 22 2013 13:11:16 for Feel++ by doxygen 1.8.5