43
edits
Bim package: Difference between revisions
Jump to navigation
Jump to search
→3D Time dependent problem
| Line 184: | Line 184: | ||
== 3D Time dependent problem == | == 3D Time dependent problem == | ||
Here is an example of a 3D time-dependent Advection-Diffusion equation that uses <code> lsode </code> | 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 | The equation being solved is | ||
| Line 191: | Line 191: | ||
= 0 \qquad \mbox{ in } \Omega \times [0, T] =[0, 1]^3 \times [0, 1] </math> | = 0 \qquad \mbox{ in } \Omega \times [0, T] =[0, 1]^3 \times [0, 1] </math> | ||
<math> \varphi = x + y - z </math> | <math> ~\varphi = x + y - z </math> | ||
<math> - \left(.01 \nabla u - u \nabla \varphi \right) \cdot \mathbf{n} = 0 </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> | |||
<pre> | |||
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 | |||
</pre> | |||
This is a video showing the .3 isosurface of the solution | |||
[[Media:bim.avi]] | |||
[[Category:OctaveForge]] | [[Category:OctaveForge]] | ||
[[Category:Packages]] | [[Category:Packages]] | ||