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