Difference between revisions of "Bim package"
Line 1:  Line 1:  
−  +  = 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.  This is a short example on how to use <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem.  
Line 55:  Line 55:  
on the mesh structure  on the mesh structure  
−  <  +  {{CodeMeshing the 2D problem<syntaxhighlight lang="octave" style="fontsize:13px"> 
[mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1);  [mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1);  
[mesh] = bim2c_mesh_properties (mesh);  [mesh] = bim2c_mesh_properties (mesh);  
−  </  +  </syntaxhighlight>}} 
to see the mesh you can use functions from the [[fpl_packagefpl package]]  to see the mesh you can use functions from the [[fpl_packagefpl package]]  
Line 237:  Line 237:  
</pre>  </pre>  
−  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. 
+  
−  
[[Category:OctaveForge]]  [[Category:OctaveForge]]  
[[Category:Packages]]  [[Category:Packages]] 
Revision as of 07:55, 20 July 2012
2D Diffusion Advection Reaction example
This is a short example on how to use bim to solve a 2D Diffusion Advection Reaction problem. The coplete code for this example can is on Agora at this link.
We want to solve the equation
with mixed Dirichlet / Neumann boundary conditions
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
pdemesh (mesh.p, mesh.e, mesh.t) view (2)
Set the coefficients for the problem:
Get the node coordinates from the mesh structure
xu = mesh.p(1,:).'; yu = mesh.p(2,:).';
Get the number of elements and nodes in the mesh
nelems = columns (mesh.t); nnodes = columns (mesh.p);
epsilon = .1; phi = xu + yu;
Construct the discretized operators
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 be composed by sides 1 and 2 and be the rest of the boundary
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 displacements
temp = [Ann Ani ; Ain Aii ] \ [ jn+bnAnd*ud(GammaD) ; biAid*ud(GammaD)]; u = ud; u(GammaN) = temp(1:numel (GammaN)); u(Omega) = temp(length(GammaN)+1:end);
Compute the fluxes through Dirichlet sides
jd = [Add Adi Adn] * u([GammaD; Omega; GammaN])  bd;
Compute the gradient of the solution
[gx, gy] = bim2c_pde_gradient (mesh, u);
Compute the internal AdvectionDiffusion flux
[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]]
fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1);
you can also plot your data directly in Octave using pdesurf
pdesurf (mesh.p, mesh.t, u)
it will look like this
3D Time dependent problem
Here is an example of a 3D timedependent AdvectionDiffusion equation that uses lsode
for timestepping.
The equation being solved is
The initial condition is
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 ( ((xx0)/(2*sx)) .^2  ((yy0)/(2*sy)) .^2  ((zz0)/(2*sz)) .^2); A = bim3a_advection_diffusion (msh, .01*ones(ne, 1), 100*(x+yz)); 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.