

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) advectiondiffusion 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 = 1e4;
b = 100;
Yo = 1;
# Initial condition
Y0 = @(x) Yo * ones (size (x));
# Sources
c = 10;
f = @(x)Yo + (c  Yo) * exp((x0.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(end1)) / (msh(end)  msh(end1));
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


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) advectiondiffusion 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 = 1e4; b = 100; Yo = 1;
# Initial condition Y0 = @(x) Yo * ones (size (x));
# Sources c = 10; f = @(x)Yo + (c  Yo) * exp((x0.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(end1)) / (msh(end)  msh(end1)); 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 = Nnodes1; 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:end1) = S(2:end1,2:end1)\(rhs(2:end1)  S(2:end1,[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.


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


> 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


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


> 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 = 1e4;
b = 100;
f = @(x)10 * exp((x0.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:end1), Ta(2:end1), '', msh(2:end1), Tb(2:end1,:), '')
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:end1), Tc(2:end1), '', msh(2:end1), Tb(2:end1,:), '')
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?

