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

[c.]

Il giorno 11 gen 2020, alle ore 00:56, JuanPi <address@hidden> ha scritto:

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_<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
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

Hi JPi,

The convention being adopted in BIM is that Diffusion Advection equations be written in the form:

   - (a * (u'(x) - b u(x)))' = f

if the system you intend to solve is

  - A u'' + B u = f

with A and B constant, you have to set

  a = A

  b = B/A

Here is an example showing the correct convergence rate (order 1 for upwind, independent of the peclet number) 

on a steady state problem :

nn = 20 * 2.^ (1:5);

for ii = 1 : numel (nn)

  n = nn(ii);

  mesh      = linspace(0,1,n+1)';

  uex       = @(r) - r.^2 + 1;

  Nnodes    = numel(mesh);

  Nelements = Nnodes-1;

  D = 1;  v = 1;  sigma = 0;

  alpha  = D*ones(Nelements,1);

  gamma  = ones(Nnodes,1);

  eta    = ones(Nnodes,1);

  beta   = 1/D*v*ones(Nelements,1);

  delta  = ones(Nelements,1);

  zeta   = sigma*ones(Nnodes,1);

  f      = @(r) 2*D - 2*v.*r + sigma*uex(r);

  rhs    = bim1a_rhs(mesh, ones(Nelements,1), f(mesh));

  S = bim1a_laplacian(mesh,alpha,gamma);

  A = bim1a_advection_upwind(mesh, beta);

  R = bim1a_reaction(mesh, delta, zeta);

  S += (A+R);

  u = zeros(Nnodes,1); u([1 end]) = uex(mesh([1 end]));

  u(2:end-1) = S(2:end-1,2:end-1)\(rhs(2:end-1) - S(2:end-1,[1 end])*u([1 end]));

  err(ii) = norm (u - uex(mesh), inf);

endfor

loglog (nn, err, nn, 1./nn)

In your example, if I understand correctly what you want to do, I suggest changing the following :

  Diff     = bim1a_laplacian (msh{i}, 1, a);

 Conv     = bim1a_advection_upwind (msh{i}, b * msh{i});

to

 Diff     = bim1a_laplacian (msh{i}, a, 1);
 Conv     = bim1a_advection_upwind (msh{i}, (b/a) * ones (size (msh{i})));

(why did you multiply b by the mesh node coordinates? was that intentional?)

c.