[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Octave forge] [bim package] bim1a_advection*
From: |
JuanPi |
Subject: |
Re: [Octave forge] [bim package] bim1a_advection* |
Date: |
Sat, 11 Jan 2020 20:58:32 +0100 |
> 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?