Drift-Diffusion PDE Model Benchmark Specification

1. Introduction

The drift-diffusion partial differential equation (PDE) describes the distribution of particles under the influence of diffusion and a drift velocity field. It often appears in physics, for instance in semiconductor device modeling or population dynamics.

2. Mathematical Formulation

The drift-diffusion equation in 2D and 3D is given by:

\[\frac{\partial \rho}{\partial t} = \nabla \cdot (D \nabla \rho - \rho \mathbf{v}),\]

where \(\rho(\mathbf{x}, t)\) is the density of the particles at position \(\mathbf{x}\) and time \(t\), \(D\) is the diffusion coefficient, and \(\mathbf{v}(\mathbf{x}, t)\) is the drift velocity.

3. Problem Specification

  • Geometry: A unit square for 2D or a unit cube for 3D.

  • Initial condition: A Gaussian distribution centered at the origin, \(\rho(\mathbf{x}, 0) = \exp(-||\mathbf{x}||^2)\).

  • Drift velocity: \(\mathbf{v}(\mathbf{x}, t) = -\frac{\mathbf{x}}{||\mathbf{x}||^2} \alpha(t)\), where \(\alpha(t) = \sin(t)\).

  • Boundary conditions: Dirichlet boundary conditions, \(\rho = 0\) on the boundary of the square/cube.

  • Diffusion coefficient: \(D = 1\).

4. Numerical Scheme

A finite difference method is used to discretize the PDE.

We discretize our domain into a grid of size \(\Delta x\), and we choose a time step \(\Delta t\). At each time step, we update the density at each grid point \(x_{i,j}\) (in 2D) or \(x_{i,j,k}\) (in 3D).

5. Input Data

For our benchmark, we will consider the following:

  • Diffusion coefficient: \(D = 1\)

  • Domain: a unit square or cube with grid size \(\Delta x = 0.01\).

  • Time: from \(t = 0\) to \(t = 1\) with a time step of \(\Delta t = 0.001\).

  • Initial density: a Gaussian, \(\rho(\mathbf{x}, 0) = \exp(-||\mathbf{x}||^2)\).

6. Code Snippet

Below is a basic Feel++ code snippet for the 2D case.

#include <feel/feel.hpp>

using namespace Feel;

int main(int argc, char**argv )
    Environment env( _argc=argc, _argv=argv,
                                  _author="Your Name",

    auto mesh = unitSquare();
    auto Vh = Pch<1>( mesh );
    auto u = Vh->element();
    auto v = Vh->element();
    auto D = 1.0;
    auto rho = exp(-norm2(P()));
    auto alpha = sin(t);
    auto drift_velocity = -P()/norm2(P())*alpha;

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh), _expr= id(v)*rho );

    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=D*gradt(u)*trans(grad(v)) - id(v)*dot(grad(u), drift_velocity) );

    a.solve(_rhs=l, _solution=u);

    // Proceed with the time stepping

Please note that the code above is a simplified example and it doesn’t cover all the steps of a proper simulation such as the time-stepping loop, updating the right-hand side, and boundary conditions.