# [Octave forge] [bim package] bim1a_advection*

6 messages
Open this post in threaded view
|

## [Octave forge] [bim package] bim1a_advection*

 Hi again, I keep playing with the bim package. I am failing to understand the behavior (or likely the correct usage) of bim1a_advection_upwind. I noticed that for finer meshes the advection effect decreases, and converges to a solution with no advection at all. How can doe sone create correct (almost mesh independent, for resoanable meshes) advection-diffusion matrices? ** Example of the problem Given this equation Y_t - a Y_xx + b Y_x + Y = f(x) Y(0, t) = Yo Y_t(1,t) = -b Y_x(1,t)   # open boundary approx O(a^2) with Y_ is the derivative of Y w.r.t. Compare the solutions generated with this code a  = 1e-4; b  = 100; Yo = 1; # Initial condition Y0 = @(x) Yo * ones (size (x)); # Sources c = 10; f = @(x)Yo + (c - Yo) * exp(-(x-0.5).^2 / 2 / 0.05^2); # Dynamic equations function dY = dYdt(Y, t, b, A, msh, f)  dY      = f - A * Y;  dY(1)   = 0;  dY(end) = - b * (Y(end) - Y(end-1)) / (msh(end) - msh(end-1)); endfunction lsode_options ("integration method", "adams"); # Meshes Nt     = 10; t      = linspace (0, 10, Nt); msh{1} = linspace (0, 1, 50).'; msh{2} = linspace (0, 1, 100).'; # Matrices & solution Nmsh = numel(msh); for i = 1:Nmsh   Diff     = bim1a_laplacian (msh{i}, 1, a);   Conv     = bim1a_advection_upwind (msh{i}, b * msh{i});   Reaction = bim1a_reaction (msh{i}, 1, 1);   Source   = bim1a_rhs (msh{i}, 1, f(msh{i}));   A = Diff + Conv + Reaction;   printf ('Solving %d\n', i); fflush (stdout)   Y{i} = lsode (@(Y, t) dYdt(Y, t, b, A, msh{i}, Source), Y0(msh{i}), t); endfor # over meshes figure (1) col = flipud (gray (Nt+3))(3:end,:); # colors for time for i=1:Nmsh   subplot (2, 1, i);   h_ = plot (msh{i}, Y{i}.');   for j=1:Nt     set(h_(j), 'color', col(j,:));   endfor # over colors   ylabel (sprintf('Y%d',i)) endfor # over solutions xlabel ('x') Thanks -- JuanPi Carbajal https://goo.gl/ayiJzi----- “An article about computational result is advertising, not scholarship. The actual scholarship is the full software environment, code and data, that produced  the  result.” - Buckheit and Donoho
Open this post in threaded view
|

## Re: [Octave forge] [bim package] bim1a_advection*

Open this post in threaded view
|

## Re: [Octave forge] [bim package] bim1a_advection*

Open this post in threaded view
|

## Re: [Octave forge] [bim package] bim1a_advection*

 >  Instead of passing the vector field BETA directly one can pass a >      piecewise linear conforming scalar function PHI as the last input. >      In such case BETA = grad PHI is assumed. Looking at the source code I see that actually it is BETA = diff PHI, not grad PHI what is used (I see now why you were puzzled) Another misleading documentation point, if I understand correctly. Why not using the function gradient, instead of diff to compute these steps? ----- “An article about computational result is advertising, not scholarship. The actual scholarship is the full software environment, code and data, that produced  the  result.” - Buckheit and Donoho