Editing Bim package
The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.
Latest revision | Your text | ||
Line 1: | Line 1: | ||
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 == | == Tutorials == | ||
Line 5: | Line 5: | ||
This is a short example on how to use <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem. | This is a short example on how to use <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem. | ||
The coplete code for this example can is on [[Agora]] at this [http://agora.octave.org/snippet/1bqV link]. | |||
We want to solve the equation | We want to solve the equation | ||
Line 65: | Line 65: | ||
to see the mesh you can use functions from the [[fpl_package|fpl package]] | to see the mesh you can use functions from the [[fpl_package|fpl package]] | ||
<pre> | |||
pdemesh (mesh.p, mesh.e, mesh.t) | pdemesh (mesh.p, mesh.e, mesh.t) | ||
view (2) | view (2) | ||
</ | </pre> | ||
[[File:fiume_msh.png]] | [[File:fiume_msh.png]] | ||
Line 77: | Line 77: | ||
Get the node coordinates from the mesh structure | Get the node coordinates from the mesh structure | ||
<pre> | |||
xu = mesh.p(1,:).'; | xu = mesh.p(1,:).'; | ||
yu = mesh.p(2,:).'; | yu = mesh.p(2,:).'; | ||
</ | </pre> | ||
Get the number of elements and nodes in the mesh | Get the number of elements and nodes in the mesh | ||
<pre> | |||
nelems = columns (mesh.t); | nelems = columns (mesh.t); | ||
nnodes = columns (mesh.p); | nnodes = columns (mesh.p); | ||
</ | </pre> | ||
<pre> | |||
epsilon = .1; | epsilon = .1; | ||
phi = xu + yu; | phi = xu + yu; | ||
</ | </pre> | ||
<b> Construct the discretized operators</b> | <b> Construct the discretized operators</b> | ||
<pre> | |||
AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi); | AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi); | ||
Mass = bim2a_reaction (mesh, 1, 1); | Mass = bim2a_reaction (mesh, 1, 1); | ||
b = bim2a_rhs (mesh,f,g); | b = bim2a_rhs (mesh,f,g); | ||
A = AdvDiff + Mass; | A = AdvDiff + Mass; | ||
</ | </pre> | ||
<b> To Apply Boundary Conditions, partition LHS and RHS</b> | <b> To Apply Boundary Conditions, partition LHS and RHS</b> | ||
Line 110: | Line 109: | ||
and <math> \Gamma_D </math> be the rest of the boundary | and <math> \Gamma_D </math> be the rest of the boundary | ||
<pre> | |||
GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST | GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST | ||
GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST | GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST | ||
Line 118: | Line 117: | ||
ud = 3*xu; ## DIRICHLET DATUM | ud = 3*xu; ## DIRICHLET DATUM | ||
Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST | Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST | ||
</pre> | |||
<pre> | |||
Add = A(GammaD, GammaD); | Add = A(GammaD, GammaD); | ||
Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! | Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! | ||
Line 134: | Line 136: | ||
bn = b(GammaN); | bn = b(GammaN); | ||
bi = b(Omega); | bi = b(Omega); | ||
</ | </pre> | ||
<B> Solve for the displacements</B> | |||
< | <pre> | ||
temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; | temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; | ||
u = ud; | u = ud; | ||
u(GammaN) = temp(1:numel (GammaN)); | u(GammaN) = temp(1:numel (GammaN)); | ||
u(Omega) = temp(length(GammaN)+1:end); | u(Omega) = temp(length(GammaN)+1:end); | ||
</ | </pre> | ||
<b> Compute the fluxes through Dirichlet sides</b><br> | <b> Compute the fluxes through Dirichlet sides</b><br> | ||
<pre> | |||
jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd; | jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd; | ||
</ | </pre> | ||
<B> Compute the gradient of the solution </B> | <B> Compute the gradient of the solution </B> | ||
<pre> | |||
[gx, gy] = bim2c_pde_gradient (mesh, u); | [gx, gy] = bim2c_pde_gradient (mesh, u); | ||
</ | </pre> | ||
<B> Compute the internal Advection-Diffusion flux</B> | <B> Compute the internal Advection-Diffusion flux</B> | ||
<pre> | |||
[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi); | [jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi); | ||
</ | </pre> | ||
<B> Export data to VTK format</B> | <B> Export data to VTK format</B> | ||
Line 170: | Line 171: | ||
or [[https://wci.llnl.gov/codes/visit/|visit]] | or [[https://wci.llnl.gov/codes/visit/|visit]] | ||
<pre> | |||
fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1); | fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1); | ||
</ | </pre> | ||
you can also plot your data directly in Octave using <code> pdesurf </code> | you can also plot your data directly in Octave using <code> pdesurf </code> | ||
<pre> | |||
pdesurf (mesh.p, mesh.t, u) | pdesurf (mesh.p, mesh.t, u) | ||
</ | </pre> | ||
it will look like this | it will look like this | ||
Line 202: | Line 203: | ||
<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> | <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> | ||
<pre> | |||
pkg load bim | pkg load bim | ||
Line 237: | Line 238: | ||
fpl_vtk_write_field (name, msh, {U(ii,:)', 'u'}, {}, 1); | fpl_vtk_write_field (name, msh, {U(ii,:)', 'u'}, {}, 1); | ||
endfor | endfor | ||
</ | </pre> | ||
[http://youtu.be/2E6Z_AcV8CQ This is a video] showing the .3 isosurface of the solution. | [http://youtu.be/2E6Z_AcV8CQ This is a video] showing the .3 isosurface of the solution. | ||
== External links == | == External links == | ||
* [http://octave.sourceforge.net/bim/index.html BIM package at Octave Forge]. | |||
[[Category:OctaveForge]] | |||
[[Category:Packages]] | |||
[[Category: |