Editing 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 <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem.
<!-- The complete code for this example is on [[Agora]] at this [http://agora.octave.org/snippet/1bqV link] -->.

We want to solve the equation

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

<math>  \varphi(x, y)\ =\ x + y </math>

with mixed Dirichlet / Neumann boundary conditions

<math> u(x, y) = u_d(x, y)\qquad \mbox{ on } \Gamma_D </math>

<math> -( \varepsilon\ \nabla u(x, y) - \nabla \varphi(x,y)\ u(x, y) )  \cdot \mathbf{n} = j_N(x, y)\qquad \mbox{ on } \Gamma_N</math>

<b> Create the mesh and precompute the mesh properties </b>

To define the geometry of the domain we can use [http://gmsh.geuz.org gmsh].

the following gmsh input 

<pre>
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};
</pre>

will produce the geometry below

[[File:fiume.png]]

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

{{Code|Meshing the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
[mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1);
[mesh] = bim2c_mesh_properties (mesh);
</syntaxhighlight>}}

to see the mesh you can use functions from the [[fpl_package|fpl package]]

{{Code|Visualizing the mesh for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
pdemesh (mesh.p, mesh.e, mesh.t)
view (2)
</syntaxhighlight>}} 

[[File:fiume_msh.png]]


<b> Set the coefficients for the problem:</b>

Get the node coordinates from the mesh structure

{{Code|Get mesh coordinates in the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
xu     = mesh.p(1,:).';
yu     = mesh.p(2,:).';
</syntaxhighlight>
}}


Get the number of elements and nodes in the mesh

{{Code|Get number of elements in the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
nelems = columns (mesh.t);
nnodes = columns (mesh.p);
</syntaxhighlight>}} 

{{Code|Set value of coefficients for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
epsilon = .1;
phi     = xu + yu;
</syntaxhighlight>}} 

<b> Construct the discretized operators</b>

{{Code|Discretized operators for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi);
Mass    = bim2a_reaction (mesh, 1, 1);
b       = bim2a_rhs (mesh,f,g);
A       = AdvDiff + Mass;
</syntaxhighlight>}}

<b> To Apply Boundary Conditions, partition LHS and RHS</b>

The tags of the sides are assigned by gmsh we let <math> \Gamma_D </math> be composed by sides 1 and 2
and <math> \Gamma_D </math> be the rest of the boundary

{{Code|Boundary conditions for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
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); 
</syntaxhighlight>}}


<B> Solve for the tracer density</B>

{{Code|Compute solution of the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
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);
</syntaxhighlight>}}

<b> Compute the fluxes through Dirichlet sides</b><br>

{{Code|Boundary fluxes in the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd;
</syntaxhighlight>}}


<B> Compute the gradient of the solution </B>

{{Code|Gradient of solution in the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
[gx, gy] = bim2c_pde_gradient (mesh, u);
</syntaxhighlight>}}

<B> Compute the internal Advection-Diffusion flux</B>

{{Code|Total flux for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi);
</syntaxhighlight>}}

<B> Export data to VTK format</B>

The resut can be exported to vtk format to visualize with [[http://www.paraview.org|paraview]] 
or [[https://wci.llnl.gov/codes/visit/|visit]]

{{Code|Export the solution of the 2D problem to vtk|<syntaxhighlight lang="octave" style="font-size:13px">
fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1);
</syntaxhighlight>}}

you can also plot your data directly in Octave using <code> pdesurf </code>

{{Code|Rubbersheet visualization of the solution of the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px">
pdesurf (mesh.p, mesh.t, u)
</syntaxhighlight>}}

it will look like this

[[File:fiume_sol_pdesurf.png|500px]]

=== 3D Time dependent problem ===

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

The equation being solved is

<math> \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] </math>

<math> ~\varphi = x + y - z </math>

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

The initial condition is

<math> u = \exp (- \left(\frac{x-.2}{.2}\right)^2 - \left(\frac{y-.2}{.2}\right)^2 - \left(\frac{z-.2}{.2}\right)^2)</math>

{{Code|Define the 3D problem|<syntaxhighlight lang="octave" style="font-size:13px"> 
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
</syntaxhighlight>}}

[http://youtu.be/2E6Z_AcV8CQ This is a video] showing the .3 isosurface of the solution.

== External links ==
* [http://octave.sourceforge.net/bim/index.html BIM package at Octave Forge].

[[Category:Octave-Forge]]