[Octave forge] [bim package] bim1a_advection*

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

[Octave forge] [bim package] bim1a_advection*

JuanPi
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

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


Reply | Threaded
Open this post in threaded view
|

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

c.-2


Il giorno 11 gen 2020, alle ore 00:56, JuanPi <[hidden email]> 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.





Reply | Threaded
Open this post in threaded view
|

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

JuanPi
Thanks Carlo!

> 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
>
Last time I was playing with the package got the impression the doc is
quite misleading. This is, I think, yeat another instanc eof that.
The doc of bim1a_advection_diffusion states the formula

- div (ALPHA * GAMMA (ETA grad (u) - BETA u)) = f

(note the missing mutiplication operator after GAMMA, making it look
like an evaluation of a function GAMMA()). Hence I was using  ETA and
BETA to set my problem, with ALPHA = GAMMA = 1.
When I moved to the separated matrices using bim1a_laplacia and
bim1a_advection_upwind I thought the same convection was kept, hence I
did not change the parameters.
I think that, in the vein of least surprise, one should try to set a
uniform model structure for al the functions.

>  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?)

In bim1_advection_upwind help it is stated that

 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.

This, luckily is indeed uniform to all the functions (perhaps also
promoting my failure, as BETA is very similar to the BETE in
bim1a_advection_diffusion).
so to avoid playing with the number of elements I just pass PHI
instead of BETA, with PHI = BETA * x
Actually I think you code will fail with size of BETA error, I think
the correct call is

Conv     = bim1a_advection_upwind (msh{i}, (b/a) * ones (size (msh{i})
- [1,0]));

as BETA is element-wise constant... This is why I prefer to pass PHI,
which is piecewise linear on the nodes.

Are you developing bim still in OF, or you have your own repo? I would
like to apply some doc patches, and add more demos (and maybe some
encapsulation of boilerplate code, like rhs and reaction, and checks
for dimension of the parameter functions).

Cheers



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


Reply | Threaded
Open this post in threaded view
|

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

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


Reply | Threaded
Open this post in threaded view
|

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

JuanPi
But still, the following is very puzzling (not consistent with de
docs), could you explain?

a = 1; b = 4;
# These all are equivalent
ConvDiff1 = bim1a_advection_diffusion (msh, a, 1, 1, (b/a) * ones (Nx, 1));
ConvDiff2 = bim1a_advection_diffusion (msh, 1, a, 1, (b/a) * ones (Nx, 1));
ConvDiff3 = bim1a_advection_diffusion (msh, 1, 1, a, b * ones (Nx,
1));   # Preferred

# These all are equivalent
Diff1 = bim1a_advection_diffusion (msh, a, 1, 1, 0);
Diff2 = bim1a_advection_diffusion (msh, 1, a, 1, 0);
Diff3 = bim1a_advection_diffusion (msh, 1, 1, a, 0);   # Preferred way
to define Poisson equation

Conv1 = bim1a_advection_diffusion (msh, 0, 1, 1, (b/a) * ones (Nx,
1));  # 0 because 1st argument multiplies all
Conv2 = bim1a_advection_diffusion (msh, 1, 0, 1, (b/a) * ones (Nx,
1));  # 0 because 2nd argument multiplies all
Conv3 = bim1a_advection_diffusion (msh, 1, 1, 0, b * ones (Nx, 1));  #
NaN because for some unexplained reason there is a division by the 3rd
argument (etak incode)


Why is Conv3 NaN? According to the equation in the doc, one could set
ETA to zero without affecting the rest of the equation


Reply | Threaded
Open this post in threaded view
|

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

JuanPi
> Why is Conv3 NaN? According to the equation in the doc, one could set
> ETA to zero without affecting the rest of the equation

I think I pin it down, it can be seen either as a mistake in the
documentation or in the implementation.

If bim1a_advection_diffusion implements the documented equation (my
correction in brackets)

 - div ( ALPHA * GAMMA [*] ( ETA [*] grad (u) - BETA [*] u ) ) = f

Then the following should give the same result

Nx  = 150;
msh = [linspace(0, 1, Nx).'; 1.01];

a = 1e-4;
b = 100;
f = @(x)10 * exp(-(x-0.5).^2 / 2 / 0.05^2);

ConvDiff = bim1a_advection_diffusion (msh, 1, 1, a, b * ones (Nx, 1));
Diff     = bim1a_advection_diffusion (msh, 1, 1, a, 0);
Conv     = bim1a_advection_upwind (msh, b * ones (Nx, 1));
Reaction = bim1a_reaction (msh, 1, 1);
Source   = bim1a_rhs (msh, 1, f(msh));

Ta = (ConvDiff + Reaction) \ Source;
Tb = (Conv + Diff + Reaction) \ Source;

figure(1)
plot(msh(2:end-1), Ta(2:end-1), '-', msh(2:end-1), Tb(2:end-1,:), '--')

It doesn't, however the following does gives the same as Tb

ConvDiff2 = bim1a_advection_diffusion (msh, 1, 1, a, (b/a) * ones (Nx, 1));

Tc = (ConvDiff2 + Reaction) \ Source;

figure(2)
plot(msh(2:end-1), Tc(2:end-1), '-', msh(2:end-1), Tb(2:end-1,:), '--')

Hence, either:

1) the parametrization of the equationin the docs  is wrong (it seems
ETA is also multiplying BETA), e.g.

 - div ( ALPHA * GAMMA * ETA * ( grad (u) - BETA * u ) ) = f

In which case GAMMA or ETA are redundant

2) the implementation is not correct and ETA is wrongly treated.

I hope it is 2) and that the correct implementation is easy to get.

Do you agree?