Bim package

Package for solving Diffusion Advection Reaction (DAR) Partial Differential Equations based on the Finite Volume Scharfetter-Gummel (FVSG) method a.k.a Box Integration Method (BIM).

Tutorials

2D Diffusion Advection Reaction example

This is a short example on how to use bim to solve a 2D Diffusion Advection Reaction problem. .

We want to solve the equation

$-\mathrm {div} \ (\varepsilon \ \nabla u(x,y)-\nabla \varphi (x,y)\ u(x,y)))+u(x,y)=1\qquad {\mbox{ in }}\Omega$

$\varphi (x,y)\ =\ x+y$

with mixed Dirichlet / Neumann boundary conditions

$u(x,y)=u_{d}(x,y)\qquad {\mbox{ on }}\Gamma _{D}$

$-(\varepsilon \ \nabla u(x,y)-\nabla \varphi (x,y)\ u(x,y))\cdot \mathbf {n} =j_{N}(x,y)\qquad {\mbox{ on }}\Gamma _{N}$

Create the mesh and precompute the mesh properties

To define the geometry of the domain we can use gmsh.

the following gmsh input

Point (1)  = {0, 0, 0, 0.1};
Point (2)  = {1, 1, 0, 0.1};
Point (3)  = {1, 0.9, 0, 0.1};
Point (4)  = {0, 0.1, 0, 0.1};
Point (5) = {0.3,0.1,-0,0.1};
Point (6) = {0.4,0.4,-0,0.1};
Point (7) = {0.5,0.6,0,0.1};
Point (8) = {0.6,0.9,0,0.1};
Point (9) = {0.8,0.8,0,0.1};
Point (10) = {0.2,0.2,-0,0.1};
Point (11) = {0.3,0.5,0,0.1};
Point (12) = {0.4,0.7,0,0.1};
Point (13) = {0.5,1,0,0.1};
Point (14) = {0.8,0.9,0,0.1};

Line (1)  = {3, 2};
Line (2) = {4, 1};

CatmullRom(3) = {1,5,6,7,8,9,3};
CatmullRom(4) = {4,10,11,12,13,14,2};
Line Loop(15) = {3,1,-4,2};
Plane Surface(16) = {15};

will produce the geometry below

we need to load the mesh into Octave and precompute mesh properties check out the tutorial for the msh package for info on the mesh structure

Code: Meshing the 2D problem

[mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1); [mesh] = bim2c_mesh_properties (mesh);

to see the mesh you can use functions from the fpl package

Code: Visualizing the mesh for the 2D problem

pdemesh (mesh.p, mesh.e, mesh.t) view (2)

Set the coefficients for the problem:

Get the node coordinates from the mesh structure

Code: Get mesh coordinates in the 2D problem

xu = mesh.p(1,:).'; yu = mesh.p(2,:).';

Get the number of elements and nodes in the mesh

Code: Get number of elements in the 2D problem

nelems = columns (mesh.t); nnodes = columns (mesh.p);

Code: Set value of coefficients for the 2D problem

epsilon = .1; phi = xu + yu;

Construct the discretized operators

Code: Discretized operators for the 2D problem

AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi); Mass = bim2a_reaction (mesh, 1, 1); b = bim2a_rhs (mesh,f,g); A = AdvDiff + Mass;

To Apply Boundary Conditions, partition LHS and RHS

The tags of the sides are assigned by gmsh we let $\Gamma _{D}$ be composed by sides 1 and 2 and $\Gamma _{D}$ be the rest of the boundary

Code: Boundary conditions for the 2D problem

GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST GammaN = setdiff (GammaN, GammaD); jn = zeros (length (GammaN),1); ## PRESCRIBED NEUMANN FLUXES ud = 3*xu; ## DIRICHLET DATUM Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST Add = A(GammaD, GammaD); Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! Adi = A(GammaD, Omega); And = A(GammaN, GammaD); ## shoud be all zeros hopefully!! Ann = A(GammaN, GammaN); Ani = A(GammaN, Omega); Aid = A(Omega, GammaD); Ain = A(Omega, GammaN); Aii = A(Omega, Omega); bd = b(GammaD); bn = b(GammaN); bi = b(Omega);

Solve for the tracer density

Code: Compute solution of the 2D problem

temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; u = ud; u(GammaN) = temp(1:numel (GammaN)); u(Omega) = temp(length(GammaN)+1:end);

Compute the fluxes through Dirichlet sides

Code: Boundary fluxes in the 2D problem

jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd;

Compute the gradient of the solution

Code: Gradient of solution in the 2D problem

[gx, gy] = bim2c_pde_gradient (mesh, u);

Compute the internal Advection-Diffusion flux

Code: Total flux for the 2D problem

[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi);

Export data to VTK format

The resut can be exported to vtk format to visualize with [[1]] or [[2]]

Code: Export the solution of the 2D problem to vtk

fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1);

you can also plot your data directly in Octave using pdesurf

Code: Rubbersheet visualization of the solution of the 2D problem

pdesurf (mesh.p, mesh.t, u)

it will look like this

3D Time dependent problem

Here is an example of a 3D time-dependent Advection-Diffusion equation that uses lsode for time-stepping.

The equation being solved is

${\frac {\partial u}{\partial t}}-\mathrm {div} \left(.01\nabla u-u\nabla \varphi \right)=0\qquad {\mbox{ in }}\Omega \times [0,T]=[0,1]^{3}\times [0,1]$

$~\varphi =x+y-z$

$-\left(.01\nabla u-u\nabla \varphi \right)\cdot \mathbf {n} =0\qquad {\mbox{ on }}\partial \Omega$

The initial condition is

$u=\exp(-\left({\frac {x-.2}{.2}}\right)^{2}-\left({\frac {y-.2}{.2}}\right)^{2}-\left({\frac {z-.2}{.2}}\right)^{2})$

Code: Define the 3D problem

pkg load bim x = linspace (0, 1, 40); y = z = linspace (0, 1, 20); msh = bim3c_mesh_properties (msh3m_structured_mesh (x, y, z, 1, 1:6)); nn = columns (msh.p); ne = columns (msh.t); x = msh.p(1, :).'; y = msh.p(2, :).'; z = msh.p(3, :).'; x0 = .2; sx = .1; y0 = .2; sy = .1; z0 = .8; sz = .1; u = exp (- ((x-x0)/(2*sx)) .^2 - ((y-y0)/(2*sy)) .^2 - ((z-z0)/(2*sz)) .^2); A = bim3a_advection_diffusion (msh, .01*ones(ne, 1), 100*(x+y-z)); M = bim3a_reaction (msh, 1, 1); function du = f (u, t, A, M) du = - M \ (A * u); endfunction time = linspace (0, 1, 30); lsode_options ("integration method", "adams"); U = lsode (@(u, t) f(u, t, A, M), u, time); for ii = 1:1:numel (time) name = sprintf ("u_%3.3d", ii); delete ([name ".vtu"]); fpl_vtk_write_field (name, msh, {U(ii,:)', 'u'}, {}, 1); endfor

This is a video showing the .3 isosurface of the solution.

External links

• BIM package at Octave Forge.

Scientific papers using BIM

• de Falco, C., Gatti, E., Lacaita, A. L., & Sacco, R. (2005). Quantum-corrected drift-diffusion models for transport in semiconductor devices. Journal of Computational Physics, 204(2), 533-561.

• de Falco, C., Jerome, J.W. and Sacco, R., 2009. Quantum-corrected drift-diffusion models: Solution fixed point map and finite element approximation. Journal of Computational Physics, 228(5), pp.1770-1789.

• de Falco, C. and O'riordan, E., 2010. INTERIOR LAYERS IN A REACTION DIFFUSION EQUATION WITH A DISCONTINUOUS DIFFUSION COEFFICIENT. International Journal of Numerical Analysis & Modeling, 7(3).

• de Falco, C., Sacco, R. and Verri, M., 2010. Analytical and numerical study of photocurrent transients in organic polymer solar cells. Computer Methods in Applied Mechanics and Engineering, 199(25), pp.1722-1732.

• de Falco, C., Porro, M., Sacco, R. and Verri, M., 2012. Multiscale modeling and simulation of organic solar cells. Computer Methods in Applied Mechanics and Engineering, 245, pp.102-116.

• de Falco, C., Denk, G. and Schultz, R., 2007. A demonstrator platform for coupled multiscale simulation. In Scientific Computing in Electrical Engineering (pp. 63-71). Springer Berlin Heidelberg.

• Maddalena, F., de Falco, C., Caironi, M. and Natali, D., 2015. Assessing the width of Gaussian density of states in organic semiconductors. Organic Electronics, 17, pp.304-318.

• Alì, G., Bartel, A., Culpo, M. and de Falco, C., 2010. Analysis of a PDE thermal element model for electrothermal circuit simulation. In Scientific Computing in Electrical Engineering SCEE 2008 (pp. 273-280). Springer Berlin Heidelberg.

• Culpo, M. and De Falco, C., 2008. Dynamical iteration schemes for coupled simulation in nanoelectronics. PAMM, 8(1), pp.10065-10068.

• Porro, M., de Falco, C., Verri, M., Lanzani, G. and Sacco, R., 2014. Multiscale simulation of organic heterojunction light harvesting devices. COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 33(4), pp.1107-1122.

• Ciucci, F., de Falco, C., Guzman, M.I., Lee, S. and Honda, T., 2012. Chemisorption on semiconductors: The role of quantum corrections on the space charge regions in multiple dimensions. Applied Physics Letters, 100(18), p.183106.

• Culpo, M., de Falco, C., Denk, G. and Voigtmann, S., 2010. Automatic thermal network extraction and multiscale electro-thermal simulation. In Scientific Computing in Electrical Engineering SCEE 2008 (pp. 281-288). Springer Berlin Heidelberg.

• Culpo, M., de Falco, C. and O’Riordan, E., 2010. Patches of finite elements for singularly-perturbed diffusion reaction equations with discontinuous coefficients. In Progress in Industrial Mathematics at ECMI 2008 (pp. 235-240). Springer Berlin Heidelberg.