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 |
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. |
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 > 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 |
> 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 = 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? |
Free forum by Nabble | Edit this page |