First implementation for stencil Diff-React

Author: Fabio Verbo

diffreactstencil.m

@@ -0,0 +1,103 @@
+function S = diffreactstencil(U, x_old, nx, ny, bnd, alpha, dxs)
+% DIFFREACTSTENCIL computes stencil as inner + bounds + corners
+
+bndN = bnd.bndN;
+bndS = bnd.bndS;
+bndE = bnd.bndE;
+bndW = bnd.bndW;
+
+S = zeros(size(U));
+
+%% Interior grid points
+for j = 2:ny-1          %(int j=1; j < jend; j++) {
+    for i = 2:nx-1      %(int i=1; i < iend; i++) {
+        S(i,j) = -(4.0 + alpha) * U(i,j) ...               % central point
+                                + U(i-1,j) + U(i+1,j) ... % east and west
+                                + U(i,j-1) + U(i,j+1) ... % north and south
+                                + alpha * x_old(i,j)  ... 
+                                + dxs * U(i,j) * (1.0 - U(i,j));
+    end
+end
+% end Interior grid points
+
+%% East boundary
+i = nx;                 %int i = nx - 1;
+for j = 2:ny-1          %(int j = 1; j < jend; j++)
+    S(i,j) = -(4.0 + alpha) * U(i,j)     ...
+                + U(i-1,j) + U(i,j-1) + U(i,j+1) ...
+                + alpha*x_old(i,j) + bndE(j)     ...
+                + dxs * U(i,j) * (1.0 - U(i,j));
+end
+% end East boundary
+
+%% West boundary
+i = 1;                  %int i = 0;
+for j = 2:ny-1 %(int j = 1; j < jend; j++)
+    S(i,j) = -(4.0 + alpha) * U(i,j)         ...
+                + U(i+1,j) + U(i,j-1) + U(i,j+1) ...
+                + alpha * x_old(i,j) + bndW(j)   ...
+                + dxs * U(i,j) * (1.0 - U(i,j));
+end
+% end West boundary
+
+%% North boundary, plus NE and NW corners
+j = ny;           %int j = ny - 1;
+
+% NW corner
+i = 1;            %int i = 0; // NW corner
+S(i,j) = -(4.0 + alpha) * U(i,j)         ...
+            + U(i+1,j) + U(i,j-1)       ...
+            + alpha * x_old(i,j) + bndW(j) + bndN(i) ...
+            + dxs * U(i,j) * (1.0 - U(i,j));
+% end NW corner
+
+% North boundary w/o corners
+for i = 2:nx-1        %(int i = 1; i < iend; i++)
+    S(i,j) = -(4.0 + alpha) * U(i,j)        ...
+                + U(i-1,j) + U(i+1,j) + U(i,j-1) ...
+                + alpha*x_old(i,j) + bndN(i)    ...
+                + dxs * U(i,j) * (1.0 - U(i,j));
+end
+% end North boundary w/o corners
+
+% NE corner
+i = nx;        % int i = nx-1; // NE corner
+S(i,j) = -(4.0 + alpha) * U(i,j)        ...
+            + U(i-1,j) + U(i,j-1)       ...
+            + alpha * x_old(i,j) + bndE(j) + bndN(i)    ...
+            + dxs * U(i,j) * (1.0 - U(i,j));
+% end NE corner
+
+%% South boundary, plus SW and SE corners
+j = 1;        % int j = 0;
+   
+% SW corner
+i = 1;        % int i = 0; // SW corner
+S(i,j) = -(4.0 + alpha) * U(i,j)    ...
+            + U(i+1,j) + U(i,j+1)   ...
+            + alpha * x_old(i,j) + bndW(j) + bndS(i) ...
+            + dxs * U(i,j) * (1.0 - U(i,j));
+% end SW corner
+
+% South boundary w/o corners
+for i = 2:nx-1      %(int i = 1; i < iend; i++)
+    S(i,j) = -(4.0 + alpha) * U(i,j)                ...
+            + U(i-1,j) + U(i+1,j) + U(i,j+1)    ...
+            + alpha * x_old(i,j) + bndS(i)      ...
+            + dxs * U(i,j) * (1.0 - U(i,j));
+end
+% end South boundary w/o corners
+
+% SE corner
+i = nx;     % int i = nx - 1; // SE corner
+S(i,j) = -(4.0 + alpha) * U(i,j)        ...
+            + U(i-1,j) + U(i,j+1)       ...
+            + alpha * x_old(i,j) + bndE(j) + bndS(i)    ...
+            + dxs * U(i,j) * (1.0 - U(i,j));
+% end SW corner
+
+end
+
+
+
+

diffusereact.m

@@ -0,0 +1,9 @@
+function B = diffusereact(X, Xold, alpha, dxs)
+% DIFFUSEREACT computes Diffusion/Reaction using Laplacian matrix.
+    
+    global A;
+    
+    aI = alpha*speye(size(A));
+    B  = (A - aI)*X + alpha*Xold + dxs * (X - X.^2);
+end
+

drstencil.m

@@ -0,0 +1,22 @@
+function S = drstencil(X, Xold, nx, ny, alpha, dxs)
+% DRSTENCIL computes the stencil for a grid including BCs
+
+    X    = reshape(X,nx,ny);
+    Xold = reshape(Xold,nx,ny);
+    S    = zeros(size(X));
+
+    for j = 2:ny-1
+        for i = 2:nx-1
+            S(i,j) = -(4.0 + alpha) * X(i,j) ...               % central point
+                                    + X(i-1,j) + X(i+1,j) ...  % east and west
+                                    + X(i,j-1) + X(i,j+1) ...  % north and south
+                                    + alpha * Xold(i,j)  ... 
+                                    + dxs * X(i,j) * (1.0 - X(i,j));
+        end
+    end
+
+    S = S(:);
+
+end
+
+

reactdiff.m

@@ -0,0 +1,152 @@
+% REACTDIFF script for testing stencil-based CG for Diffusion/Reaction problems.
+% 
+
+%% SET PARAMS
+
+nx      = 64;
+ny      = 64;
+assert(nx == ny, 'Grid must be square. Set nx = ny.');
+N       = nx * ny;
+x0      = 0;
+xmax    = 1;
+y0      = 0;
+ymax    = 1;
+tlim    = 0.01;
+tsteps  = 10;
+
+maxitCG = 100;
+tolCG   = 1e-4;
+
+maxitNT = 20;
+tolNT   = 1e-8;
+
+T  = linspace(0, tlim, tsteps);
+dt = tlim/tsteps;
+
+bndN = 0.00;
+bndS = 0.00;
+bndW = 0.00;
+bndE = 0.00;
+
+bnd.bndN = bndN;
+bnd.bndS = bndS;
+bnd.bndE = bndE;
+bnd.bndW = bndW;
+
+dx = (xmax-x0)/(nx);
+dy = (ymax-y0)/(ny);
+
+%D = 1;          % Diffusion constant
+%R = 1;
+F = 0.5;        % Initial condition
+
+alpha = dx^2/dt;
+dxs   = 100 / dx^2;
+
+
+global A;
+A = -gallery('poisson', nx);
+
+%% SET INITIAL CONDITIONS
+X = zeros(nx,ny);
+    
+% Center and radius of the circle
+c  = [xmax-1/4, y0+1/4];
+r  = min(c) / 2;
+
+for j = 1:nx
+    x = x0 + j*dx;
+    for k = 1:ny
+        p = [x, y0+k*dy];
+        %disp(p);
+        if norm(c-p) <= r
+            X(j,k) = F;
+        end
+    end
+end
+
+%% SET BOUNDARY CONDITIONS
+
+X(1,:)   = bndN;
+X(end,:) = bndS;
+X(:,1)   = bndW;
+X(:,end) = bndE;
+
+X = X(:);
+
+%% RUN TIME-S
+
+figure('Name',sprintf('%d x %d grid', nx, ny));
+
+global f;
+
+tTot = tic;
+
+ts = 1;
+for t = T
+    
+    fprintf('T-step %4d/%4d:\n', ts, tsteps);
+    
+    Xold = X;
+    imagesc(reshape(X, nx,nx)); colorbar; pause(0.5);
+
+    convNT = false;
+    for k = 1:maxitNT
+        
+        %f = @(X) drstencil(X, Xold, nx, ny, alpha, dxs);
+        f = @(X) diffusereact(X, Xold, alpha, dxs);
+        
+        X = f(Xold);
+        
+        fB = norm(X);
+        fprintf('NWT: it=%3d, ||f(B)|| = %.4e   |  ', k, fB);
+        if (fB < tolNT)
+            fprintf('Newton method has converged!\n');
+            convNT = true;
+            break;
+        end
+    
+        [dX, flg, relres, its] = stencil_pcg(f, X, tolCG, maxitCG);
+    
+        
+        fprintf('CG: ');
+        switch flg
+        case 1
+            fprintf('CG iterated %d times but did not converge.\n', maxitCG);
+        case 2
+            fprintf('CG preconditioner M was ill-conditioned.\n');
+        case 3
+            fprintf('CG stagnated.\n');
+        case 4
+            fprintf('CG: one of the scalar quantities is too small or too large.\n');
+        end
+        
+        if flg > 0
+            break;
+        end
+    
+        fprintf('residual %.4e after %3d iterations.\n', relres, its);
+        X    = Xold - dX;
+        Xold = X;
+        %spy(reshape(dX, nx,nx))
+        
+    end
+    
+    imagesc(reshape(X, nx,nx)); colorbar;
+    
+    if convNT
+        ts = ts+1;
+        continue;
+    else
+        error('Newton method FAILED to converge.');
+    end
+    
+        
+    
+end
+
+fprintf('-----------------------------------\n');
+fprintf('Simulation took %f seconds.\n', toc(tTot));
+
+
+

stencil_pcg.m

@@ -0,0 +1,67 @@
+function [x, flg, resNew, its] = stencil_pcg(f, b, tol, maxit)
+% STENCIL_PCG computes CG with 1st order approximation for f
+
+    its = 0;
+    flg = 1;
+
+    e = 1e-08;
+    
+    x0 = zeros(size(b));     % initial guess
+    x  = x0;
+    
+    
+    Fx  = f(x);
+    xx  = (1+e)*x;
+    Fxx = f(xx);
+    r = b - (1/e) * (Fxx - Fx);         % A*x \approx (1/e) * (Fxx-Fx)
+    p = r;
+    
+    resOld = r'*r;
+    resNew = resOld;
+    
+    if (sqrt(resOld) < tol)
+        flg = 0;
+        %fprintf('STENCIL_CG: converged in 0 iterations (initial guess is good enough!\n');
+        return;
+    end
+    
+    for its = 1:maxit
+        
+        %fprintf('STENCIL_CG: %3d ', its);
+       
+        Fx = f(x);
+        xx  = x + e*p;
+        Fxx = f(xx);
+        Ap  = (1/e) * (Fxx-Fx);     % A*p approx. by Ap
+        %fprintf(' * RES. FOR Ap = %.4e * ', norm(Ap-A*p)/norm(A*p));
+        
+        alpha = r'*r / (p'*Ap);      
+        
+        x    = x + alpha * p;
+        rnew = r - alpha * Ap;
+        
+        resNew = rnew' * rnew;
+        if (sqrt(resNew) < tol)
+            flg = 0;
+            break;
+        end
+
+        beta = resNew / resOld;   % res=<r(k-1),r(k-1)>, resNew=<r(k+1),r(k+1)>
+        p    = rnew + beta * p;
+ 
+        r      = rnew;
+        resOld = resNew;
+        %fprintf(' - RELRES: %.4e\n', sqrt(resNew));
+        
+    end
+    
+    if flg == 0
+        %fprintf('STENCIL_CG: converged in %d iterations.\n', its);
+    else
+        %fprintf('STENCIL_CG: could NOT converge.\n');
+    end
+    
+end
+
+
+