# 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,
_desc=makeOptions(),
_email="your.email@domain.com"));

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),