waveSEM

DD/SchurNKP/annulusgrid.m

@@ -0,0 +1,6 @@
+function [z0,quads,curv] = annulusgrid(R1,R2)
+z0=kron([R1;R2],exp(1i*(0:3)'/4*2*pi));
+quads=[1,5,2,6; 2,6,3,7; 3,7,4,8; 4,8,1,5];
+curv=zeros(size(quads)); curv(:,1)=-1; curv(:,2)=-2;
+end
+

DD/SchurNKP/benchmarkNKP.m

@@ -0,0 +1,58 @@
+function [] = benchmarkNKP(ms)
+
+L=3;
+R1=1;  R2=1;
+z1=L;  z2=-L;
+m1=1;  m2=1;
+function F=binrhs(z,rho)
+    r1=hypot(rho,z-z1);
+    r2=hypot(rho,z-z2);
+    F1=-m1*(r1<=R1).*sinc(pi*r1/R1);
+    F2=-m2*(r2<=R2).*sinc(pi*r2/R2);
+    F=F1+F2;
+end
+metric = @(z,rho) abs(rho);
+RHS=@binrhs;
+
+[z0,quads,curv]=bingrid(R1,R2,L);
+
+% Mesh refinement
+ref=1;
+for j=1:ref
+    [~, adj, bnd] = meshtopo(quads);
+    [z0,quads,curv]=quadmeshrefine(z0,quads,curv,adj,bnd);
+end
+
+ngrid=zeros(size(ms));
+pcalls=zeros(size(ms));
+for i=1:numel(ms)
+    m=ms(i);
+    [lam, uu, X, Y, relres, calls] = meshSchurNKP(m, z0, quads, curv, 0, RHS, metric);
+    delete(findall(gcf,'Type','light'));
+    drawnow;
+    pcalls(i)=calls;
+    ngrid(i)=numel(uu);
+    display(relres);
+    display(calls);
+end
+
+figure(20);
+semilogx(ngrid,pcalls,'--*b','Linewidth',1);
+set(gcf,'defaultTextInterpreter','latex');
+set(gca,'TickLabelInterpreter','latex','fontsize',14);
+xlabel('Total gridpoints');
+title('Preconditioner calls by GMRES');
+ylim([0,100]);
+xlim([1E5,5E6]);
+set(gca,'XTick',ngrid);
+set(gca,'XTickLabel',num2str(ngrid(:),'%.2G'));
+set(gca,'XMinorTick','off');
+
+dp=0.15*min(pcalls);
+for i=1:numel(ms)
+    text(ngrid(i),pcalls(i)-dp,sprintf('$n=%d$',ms(i)),...
+        'FontSize',14,'HorizontalAlignment','center');
+end
+
+end
+

DD/SchurNKP/demoSchurNKP.m

@@ -1,28 +1,27 @@
 function [lam] = demoSchurNKP(m, shape, ref)
-% m = Number of gridpoint in 1D per quadrilateral element
+% m = Number of gridpoints in 1D per quadrilateral element
 % shape = 1, Equilateral triangle
 % shape = 2, Right triangle
 % shape = 3, Isoceles triangle
 % shape = 4, Scalene triangle
 % shape = 5, L shaped membrane
 % shape = 6, Unit disc
-% shape = 7, Newtonian binary stars
+% shape = 7, Annulus
+% shape = 8, Newtonian binary stars
 % ref = Level of refinement (0 = no refinement)
 
 function F=myrhs(x,y)
     F=1+0*x+0*y;
 end
-
-R1=1;
-R2=1;
-L=6;
+L=3;
+R1=1;  R2=1;
+z1=L;  z2=-L;
+m1=1;  m2=1;
 function F=binrhs(z,rho)
-    z1= L;
-    z2=-L;
     r1=hypot(rho,z-z1);
     r2=hypot(rho,z-z2);
-    F1=(r1<=R1).*sinc(pi*r1/R1);
-    F2=(r2<=R2).*sinc(pi*r2/R2);
+    F1=-m1*(r1<=R1).*sinc(pi*r1/R1);
+    F2=-m2*(r2<=R2).*sinc(pi*r2/R2);
     F=F1+F2;
 end
 
@@ -51,10 +50,14 @@
             [z0,quads,curv]=Lgrid();
         case 2
             [z0,quads,curv]=coregrid(1);
+            metric = @(z,rho) abs(rho);
         case 3
+            [z0,quads,curv]=annulusgrid(1,2);
+            metric = @(z,rho) abs(rho);
+        case 4
             [z0,quads,curv]=bingrid(R1,R2,L);
+            metric = @(z,rho) abs(rho);
             RHS=@binrhs;
-            metric = @(z,rho) abs(rho)+0*z;
     end
 end
 
@@ -64,9 +67,71 @@
     [z0,quads,curv]=quadmeshrefine(z0,quads,curv,adj,bnd);
 end
 
-lam = meshSchurNKP(m, z0, quads, curv, 0, RHS, metric);
-
-% delete(findall(gcf,'Type','light'));
+k=32;
+[lam, uu, X, Y, relres, pcalls] = meshSchurNKP(m, z0, quads, curv, k, RHS, metric);
+delete(findall(gcf,'Type','light'));
 % shading faceted;
+display(relres);
+display(pcalls);
+
+if shape==6
+    r=hypot(X,Y);
+    uex=(1-r.^2)/6;
+
+    err=norm(uu(:)-uex(:),'inf')/norm(uex(:),'inf');
+    disp(err);
+
+    figure(2);
+    for j=1:size(quads,1)
+        surf(X(:,:,j), Y(:,:,j), uu(:,:,j)-uex(:,:,j));
+        if j==1, hold on; end
+    end 
+    hold off;
+    shading interp;
+    view(2);  
+end
+
+
+if shape==7
+    r=hypot(X,Y);
+    uex=(-6+7*r-r.^3)./(6*r);
+
+    err=norm(uu(:)-uex(:),'inf')/norm(uex(:),'inf');
+    disp(err);
+
+    figure(2);
+    for j=1:size(quads,1)
+        surf(X(:,:,j), Y(:,:,j), uu(:,:,j)-uex(:,:,j));
+        if j==1, hold on; end
+    end 
+    hold off;
+    shading interp;
+    view(2);  
+end
+
+
+if shape==8
+    zet=X;
+    rho=abs(Y);
+    r1=hypot(rho,zet-z1);
+    r2=hypot(rho,zet-z2);
+    u1=m1*(R1./pi)^2*((r1<=R1).*(-1-R1*sinc(pi*r1/R1))+...
+                      (r1> R1).*(-R1./(r1+(r1<=R1))));
+    u2=m2*(R2./pi)^2*((r2<=R2).*(-1-R2*sinc(pi*r2/R2))+...
+                      (r2> R2).*(-R2./(r2+(r2<=R2))));
+    uex=u1+u2;
+    err=norm(uu(:)-uex(:),'inf')/norm(uex(:),'inf');
+    display(err);
+
+    figure(2);
+    for j=1:size(quads,1)
+        surf(X(:,:,j), Y(:,:,j), uu(:,:,j)-uex(:,:,j));
+        if j==1, hold on; end
+    end 
+    hold off;
+    shading interp;
+    view(2);
+end
+
 end
 

DD/SchurNKP/feedSchurNKP.m

@@ -1,4 +1,4 @@
-function [ischur,jschur,eschur,nkp,gf]=feedSchurNKP(GID,A1,B1,A2,B2,C1,C2)
+function [ischur,jschur,eschur,nkp,gf]=feedSchurNKP(GID,A1,B1,A2,B2)
 % Updates SchurNKP S given domain topology net (EWNS) and NKP (A1,B1,A2,B2)
 % Returns preconditioned Green's function gf
 m=size(A1,1);
@@ -10,14 +10,13 @@
 nkp=@(uu,I,J) A1(I,:)*uu*B1(J,:)'+A2(I,:)*uu*B2(J,:)';
 
 % Eigenfunctions
-[E1,V1,L1,D1]=eigenfunctions(A1, A2, C1);
-[E2,V2,L2,D2]=eigenfunctions(B2, B1, C2);
+[V1,L1,D1]=eigenfunctions(A1, A2);
+[V2,L2,D2]=eigenfunctions(B2, B1);
 
 % Green's function
 LL=1./((L1.*D1)*D2'+D1*(L2.*D2)');
 function uu=greenF(F,uu)
-    uu=E1*uu*E2';
-    rhs=F-E1(:,kd1)'*(A1*uu*B1'+A2*uu*B2')*E2(:,kd2);
+    rhs=F-(A1(kd1,:)*uu*B1(kd2,:)'+A2(kd1,:)*uu*B2(kd2,:)');
     uu=uu+V1(:,kd1)*(LL.*(V1(kd1,kd1)'*rhs*V2(kd2,kd2)))*V2(:,kd2)';
 end
 gf=@(F,uu) greenF(F,uu);
@@ -30,10 +29,10 @@
 T2=rd2(any(GID(:,rd2),1));
 
 % Building blocks
-MV1=E1'*A2*V1(:,kd1);
-KV1=E1'*A1*V1(:,kd1);
-MV2=E2'*B1*V2(:,kd2);
-KV2=E2'*B2*V2(:,kd2);
+MV1=A2*V1(:,kd1);
+KV1=A1*V1(:,kd1);
+MV2=B1*V2(:,kd2);
+KV2=B2*V2(:,kd2);
 
 % [E W] x [E W]
 [I,J]=meshgrid(T1,T1);
@@ -95,27 +94,17 @@
 end
 end
 
-function [E,V,L,D]=eigenfunctions(K,M,C)
+function [V,L,D]=eigenfunctions(K,M,C)
 % Computes stiffness SK and mass SM matrices in the constrained basis E and
 % constrainded eigenfunctions V and eigenvalues L
 m=size(K,1);
 kd=2:m-1;
-rd=[1,m];
-
-% Constrained basis
-E=eye(m);
-E(rd,kd)=-C(:,kd);
-E(rd,:)=C(:,rd)\E(rd,:);
-SK=E'*K*E;
-SM=E'*M*E;
-
 % Eigenfunctions
 V=zeros(m);
-[V(kd,kd),L]=eig(SK(kd,kd), SM(kd,kd), 'vector');
+[V(kd,kd),L]=eig(K(kd,kd), M(kd,kd), 'vector');
 % Diagonalized mass matrix
 % D=diag(V(kd,kd)'*SM(kd,kd)*V(kd,kd));
-D=sum(V(kd,kd).*(SM(kd,kd)*V(kd,kd)), 1).';
-V(rd,kd)=E(rd,kd)*V(kd,kd);
+D=sum(V(kd,kd).*(M(kd,kd)*V(kd,kd)), 1).';
 V(kd,kd)=bsxfun(@rdivide, V(kd,kd), sqrt(abs(D)).');
 V(kd,kd)=bsxfun(@rdivide, V(kd,kd), sign(V(kd(1),kd)));
 D=sign(D);