help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## Re: [Octave forge] [bim package] bim1a_advection*

 From: c. Subject: Re: [Octave forge] [bim package] bim1a_advection* Date: Sat, 11 Jan 2020 09:28:40 +0100

 Il giorno 11 gen 2020, alle ore 00:56, JuanPi 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, andconverges to a solution with no advection at all.How can doe sone create correct (almost mesh independent, forresoanable meshes) advection-diffusion matrices?** Example of the problemGiven this equationY_t - a Y_xx + b Y_x + Y = f(x)Y(0, t) = YoY_t(1,t) = -b Y_x(1,t)   # open boundary approx O(a^2)with Y_ is the derivative of Y w.r.t. Compare the solutions generated with this codea  = 1e-4;b  = 100;Yo = 1;# Initial conditionY0 = @(x) Yo * ones (size (x));# Sourcesc = 10;f = @(x)Yo + (c - Yo) * exp(-(x-0.5).^2 / 2 / 0.05^2);# Dynamic equationsfunction 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));endfunctionlsode_options ("integration method", "adams");# MeshesNt     = 10;t      = linspace (0, 10, Nt);msh{1} = linspace (0, 1, 50).';msh{2} = linspace (0, 1, 100).';# Matrices & solutionNmsh = 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 meshesfigure (1)col = flipud (gray (Nt+3))(3:end,:); # colors for timefor 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 solutionsxlabel ('x')ThanksHi JPi,The convention being adopted in BIM is that Diffusion Advection equations be written in the form:   - (a * (u'(x) - b u(x)))' = fif the system you intend to solve is  - A u'' + B u = fwith A and B constant, you have to set  a = A  b = B/AHere 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);endforloglog (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.

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]