function [norm2]=scheme(D,nx,nyc,nyf,central,hom); % [norm2]=scheme(0.01, 16, 8, 64, 0, 0); % Numerical solution of 2-D singular perturbation problem for stationary % convection-diffusion equation: % - dw/dx - D*Laplacian(w) = q % x = 1: Dirichlet x = 0: homogeneous or inhomogeneous Neumann % y = 0: Dirichlet y = 1: homogeneous Neumann % % Central or upwind scheme for convection % Cell-centered grid; uniform in x, locally refinement % % % Input variables: % D = Diffusion coefficient ( = 1/[Peclet number] ) % nx = Horizontal number of cells % nyc = Vertical number of cells outside refinement zone % nyf = Vertical number of cells inside refinement zone % central = Enter 1 for central scheme or something else for upwind scheme % hom = Neumann condition at outflow boundary: 1 for homogeneous, % or something else for inhomogeneous (exact) condition % Output variables: % norm2 = L2 norm of error compared with exact solution % Copyright 2001 P. Wesseling % This program and its subprograms may be freely used, modified and distributed % under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html % Theory is given in Section 4.7 of % P. Wesseling: Principles of Computational Fluid Dynamics % Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp. % See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html % This program generates tables 4.3 and 4.4 in the book % Functions called: exact_solution, fL, q %NK: Modify for large D del = min(8*sqrt(D), 0.5); % Thickness of refinement zone % Implementation of upwind scheme by artificial diffusion in x-direction dx = 1/nx; if central == 1 Da = D; % Artificially increased diffusion coefficient else Da = D + dx/2; end % Definition of volume boundaries and grid points % x = 0 % grid |---o---|---o---|---o---|---o---| % 1 1 2 2 3 J % |<---- del----->| % y=0 K1 cells | K2 cells y=1 % grid |-o-|-o-|-o-|-o-|---o---|---o---|---o---| % 1 1 2 2 3 3 J = nx; K1 = nyf; K2 = nyc; x = [1:J]'*dx - dx/2; % Horizontal coordinates of cell centers K = K1 + K2; by = zeros(K+1,1); % Vertical coordinates of cell boundaries dy = del/K1; % Vertical mesh size in refinement zone for k = 1:K1+1, by(k) = (k-1)*dy; end dy = (1-del)/K2; % Vertical mesh size in unrefined zone for k = K1+2:K+1, by(k) = del + (k-K1-1)*dy; end n=length(by); y=zeros(n-1,1); % Vertical coordinates of cell centers dy=zeros(n-1,1); % Vertical mesh sizes for k = 1:n-1 y(k) = 0.5*(by(k)+by(k+1)); dy(k) = by(k+1)-by(k); end rhs = zeros(J*K,1); % Right-hand side of differential equation for k = 1:K % Lexicographical ordering of unknowns: for j = 1:J % j: x-direction k: y-direction rhs(j + (k-1)*J) = q(x(j),y(k),D)*dx*dy(k); end end % Modification of rhs due to boundary conditions % x = 0: outflow: homogeneous or inhomogeneous Neumann, depending on % parameter hom for k = 1:K m = 1 + (k-1)*J; rhs(m) = rhs(m) + dy(k)*(Da - 0.5*dx)*fL(y(k),D,hom); end % x = 1: Dirichlet for k = 1:K m = k*J; rhs(m) = rhs(m) + dy(k)*(1 + 2*Da/dx)*exact_solution(1,y(k),D); end % y = 1: homogeneous Neumann, y = 0: Dirichlet for j = 1:J, rhs(j) = rhs(j) + 2*dx*D*exact_solution(x(j),0,D)/dy(1); end % Matrix construction by flux evaluation n = J*K; z = zeros(n,1); d = [-J -1 0 1 J]; A = spdiags([ z z z z z], d',n,n); % Declaration of sparsity pattern for k = 1:K % Loop over vertical faces m1 = (k-1)*J; m = 1 + m1; A(m,m) = A(m,m) + dy(k); a = - dy(k)/2 + dy(k)*Da/dx; b = - dy(k)/2 - dy(k)*Da/dx; for j = 2:J m = j + m1; A(m-1,m-1) = A(m-1,m-1) + a; A(m-1,m) = A(m-1,m) + b; A(m,m-1) = A(m,m-1) - a; A(m,m) = A(m,m) - b; end m = J+m1; A(m,m) = A(m,m) + 2*Da*dy(k)/dx; end for j=1:J % Loop over horizontal faces; m = j; A(m,m) = A(m,m) + 2*D*dx/dy(1); for k = 2:K m = j + (k-1)*J; a = 2*dx*D/(dy(k) + dy(k-1)); A(m-J,m-J) = A(m-J,m-J) + a; A(m-J,m) = A(m-J,m) - a; A(m,m-J) = A(m,m-J) - a; A(m,m) = A(m,m) + a; end end numerical_solution = A\rhs; exsol = zeros (K*J,1); % Exact solution for k = 1:K for j = 1:J, exsol(j + (k-1)*J) = exact_solution(x(j),y(k),D); end end % Text output error = numerical_solution - exsol; n = K*J; norm1 = norm(error,1)/n; norm2 = norm(error,2)/sqrt(n); [norminf,m] = max(abs(error)); %NK: Modify figures % Solution profile at x = x(1) hh = zeros(K,1); for k = 1:K, hh(k) = numerical_solution(1 + (k-1)*J); end figure (1), clf, subplot(2,2,1), hold on, plot(y,hh,'o') yy = 0:0.02:1; exact = zeros (length(yy)); for kk = 1:length(yy), exact(kk) = exact_solution(x(1), yy(kk), D); end plot(yy,exact);xlabel('Y');ylabel('\phi(y)'); title(['Solution profile at x = ', num2str(x(1))]);grid on; % Solution profile at y = y(kkk) vv = zeros(J,1); kkk = 1; y(kkk) for j = 1:J, vv(j) = numerical_solution(j + (kkk-1)*J); end subplot(2,2,4), plot(x,vv,'o') for kk = 1:length(yy), exact(kk) = exact_solution( yy(kk), y(kkk), D); end hold on, plot(yy,exact);xlabel('X');ylabel('\phi(x)'); grid on; title(['Solution profile at y = ', num2str(y(kkk))]); % Solution profile at x = x(J) hh = zeros(K,1); for k = 1:K, hh(k) = numerical_solution(J + (k-1)*J); end subplot(2,2,2), plot(y,hh,'o'), hold on for kk = 1:length(yy), exact(kk) = exact_solution(x(J), yy(kk), D); end plot(yy,exact) xlabel('Y');ylabel('\phi(y)'); title(['Solution profile at x = ', num2str(x(J))]);grid on; % Contour plot consol = zeros(K,J); for k = 1:K,for j = 1:J, consol(k,j) = numerical_solution((k-1)*J + j);end,end subplot(2,2,3),colormap autumn; contour(x,y,consol);xlabel('X');ylabel('Y');title('Contour plot of numerical solution');grid on; numsol = zeros(K,J); % waterfall plot for k = 1:K,for j = 1:J, numsol(k,j) = numerical_solution(j + (k-1)*J);end,end figure(2),clf, surfc(x,y,numsol,'EdgeColor','none');xlabel('X');ylabel('Y'); camlight left; lighting phong;view(200,70), colormap autumn;grid on;