- 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" );
Environment env( _argc=argc, _argv=argv,
_desc=opts,
_about=about(_name="myadvection",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org") );
auto mesh = unitSquare();
auto Xh = Pch<1>( mesh );
auto u = Xh->element( "u" );
auto v = Xh->element( "v" );
double epsilon = option(_name="epsilon").as<double>();
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.);
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 ) ) );
auto l = form1( _test=Xh );
l+= integrate( _range=elements( mesh ), _expr=f*id( v ) );
a += on( _range=boundaryfaces( mesh ), _rhs=l, _element=u,
_expr=cst(0.) );
a.solve( _rhs=l, _solution=u );
auto e = exporter( _mesh=mesh );
e->add("u",u);
e->save();
}
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.
|
|
|
\(\epsilon=1\) | \(\epsilon=0.01\) | \(\epsilon=0.0001\) |