help-octave
[Top][All Lists]

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

 From: JuanPi Subject: [Octave forge] [bim package] bim1a_advection* Date: Sat, 11 Jan 2020 00:56:27 +0100

```Hi again,

I keep playing with the bim package.
I am failing to understand the behavior (or likely the correct usage)
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

** 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_<var> is the derivative of Y w.r.t. <var>

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

# 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

-----