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?