greedy coloring

Stabilization/adf2.m

@@ -14,8 +14,8 @@
 ifplot=true;
 
 % GMRES settings
-maxit=100;
-tol=1e-10;
+maxit=10;
+tol=1e-9;
 %maxit=0;
 
 % Problem settings
@@ -29,7 +29,7 @@
 uy=@(x,y) 1-0*x.*(1-y.^2);
 
 % Stabilization CFL
-CFL=1.0E-2;
+CFL=3E-1/2;
 %CFL=inf;
 
 function [u]=vel(x,y,ijac)
@@ -338,7 +338,7 @@
         visit(ie)=visit(ie)+1;
         im(1:nex,1:ney)=visit;
         set(hv,'CData',im');
-        drawnow;
+%        drawnow;
     end
 end
 

Stabilization/bubfilt.m

@@ -9,15 +9,13 @@
 
 if(nargin==1)
 sigma = ones(N,1);
-p = 0.0;
+p = 1/16; p=0;
 cut = max(1,ceil(N*p));
-alpha = 1E-2;
+alpha = (5E-1/2)*max(1/2,1-p*cut);
 for k=1:cut
-    w = alpha*(((cut+1-k)/cut)^2);
     j = N+1-k;
-    sigma(j) = 1-w;
+    sigma(j) = 1-alpha*(((cut+1-k)/cut)^2);
 end
 end
-
 F = V*diag(sigma)/V;
 end

Stabilization/color_greedy.m

@@ -0,0 +1,15 @@
+function [icolor] = color_greedy(itopo,elems)
+nel=size(itopo,2);
+if(nargin==1)
+    elems=1:nel;
+end
+elems=reshape(elems,1,[]);
+icolor=zeros(nel,1);
+for e=elems
+    k=1;
+    while any(k==icolor(itopo(:,e)))
+        k=k+1;
+    end
+    icolor(e)=k;
+end
+end

Stabilization/get_graph.m

@@ -1,5 +1,4 @@
 function [iflux] = get_graph(itopo,vx,vy,wx,wy)
-
     tol=1e-6;
     nfaces=size(itopo,1);
     nel=size(itopo,2);    
@@ -50,55 +49,5 @@
         end
     end
     end
-
-
-
-    return
-    heat=zeros(nel,1);
-    esrc=1;
-    if(esrc==0)
-        heat(all(iflux>=0,2))=1;
-    else
-        heat(esrc)=1;
-    end
-
-    itemp=-ones(nfaces,nel);
-    isweep=zeros(nel,1);
-    src=find(heat>0);
-    q=length(src);
-    isweep(1:q)=src;
-    p1=1;
-    p2=q;
-    qlag=0;
-    while (q>qlag && q<=nel)
-        qlag=q;
-        for p=p1:p2
-            e=isweep(p);
-            heat(e)=2;
-            for f=1:nfaces
-                ee=itopo(f,e);
-                if(heat(ee)~=2)
-                    itemp(f,e)=1;
-                    if(heat(ee)==0)
-                        q=q+1;
-                        isweep(q)=ee;
-                        heat(ee)=1;
-                    end
-                end
-            end
-        end
-        p1=p2+1;
-        p2=q;
-    end
-
-    for e=1:nel
-        for f=1:nfaces
-            if(itopo(f,e)==e)
-                iflux(f,e)=0;
-            elseif(iflux(f,e)==0)
-                iflux(f,e)=itemp(f,e);
-            end
-        end
-    end
 end
 

Stabilization/my_gmres.m

@@ -0,0 +1,87 @@
+function [x, e] = my_gmres( A, b, x, max_iterations, threshold)
+  n = numel(b);
+  m = max_iterations;
+
+  %use x as the initial vector
+  r=b-A(x);
+
+  b_norm = norm(b);
+  error = norm(r)/b_norm;
+
+  %initialize the 1D vectors
+  sn = zeros(m,1);
+  cs = zeros(m,1);
+  e1 = zeros(n,1);
+  e1(1) = 1;
+  e=[error];
+  r_norm=norm(r);
+  Q(:,1) = r/r_norm;
+  beta = r_norm*e1;
+  for k = 1:m                                   
+
+    %run arnoldi
+    [H(1:k+1,k) Q(:,k+1)] = arnoldi(A, Q, k);
+
+    %eliminate the last element in H ith row and update the rotation matrix
+    [H(1:k+1,k) cs(k) sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);
+
+    %update the residual vector
+    beta(k+1) = -sn(k)*beta(k);
+    beta(k)   = cs(k)*beta(k);
+    error  = abs(beta(k+1)) / b_norm;
+
+    %save the error
+    e=[e; error];
+    if ( error <= threshold)
+      break;
+    end
+  end
+
+  %calculate the result
+  y = H(1:k,1:k) \ beta(1:k);
+  x = x + Q(:,1:k)*y; 
+end
+
+%----------------------------------------------------%
+%                  Arnoldi Function                  %
+%----------------------------------------------------%
+function [h, q] = arnoldi(A, Q, k)
+  q = A(Q(:,k));
+  for i = 1:k
+    h(i)= q'*Q(:,i);
+    q = q - h(i)*Q(:,i);
+  end
+  h(k+1) = norm(q);
+  q = q / h(k+1);
+end
+
+%---------------------------------------------------------------------%
+%                  Applying Givens Rotation to H col                  %
+%---------------------------------------------------------------------%
+function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
+  %apply for ith column
+  for i = 1:k-1                              
+    temp   =  cs(i)*h(i) + sn(i)*h(i+1);
+    h(i+1) = -sn(i)*h(i) + cs(i)*h(i+1);
+    h(i)   = temp;
+  end
+
+  %update the next sin cos values for rotation
+  [cs_k sn_k] = givens_rotation(h(k), h(k+1));
+
+  %eliminate H(i+1,i)
+  h(k) = cs_k*h(k) + sn_k*h(k+1);
+  h(k+1) = 0.0;
+end
+
+%%----Calculate the Given rotation matrix----%%
+function [cs, sn] = givens_rotation(v1, v2)
+  if (v1==0)
+    cs = 0;
+    sn = 1;
+  else
+    t=sqrt(v1^2+v2^2);
+    cs = abs(v1) / t;
+    sn = cs * v2 / v1;
+  end
+end

Stabilization/toposort_loops.m

@@ -9,6 +9,7 @@
 
 k=1;
 q=0;
+% Detect root
 for e=1:nel
     if(all(iflux(:,e)>=0))
         q=q+1;
@@ -24,6 +25,14 @@
     q=q+1;
     isweep(q)=e;
     icolor(e)=k;
+else
+    % apply greedy coloring of root (only needed for DG)
+    iroot=find(icolor==1); 
+    icolor=color_greedy(itopo,iroot);
+    iroot=find(icolor==1);
+    q=length(iroot);
+    isweep(1:q)=iroot;
+    icolor(icolor>1)=0;
 end
 
 t=0;
@@ -51,7 +60,7 @@
     e=isweep(p);   
     for f=1:nfaces
         ee=itopo(f,e);
-        if(iflux(f,e)>0 && ee~=e)
+        if(iflux(f,e)>=0 && ee~=e)
             ff=find(itopo(:,ee)==e);
             itopo(ff,ee)=ee;
             iflux(ff,ee)=0;
@@ -74,55 +83,4 @@
     p2=q;
 end
 [~,isweep]=sort(icolor);
-end
-
-
-
-function [itopo,iflux]=cut_loops(itopo,iflux,i1,i2)
-    if(i2>i1)
-        % TODO check if DAG
-        j1=floor((i1+i2-1)/2);
-        j2=j1+1;
-        nfaces=size(itopo,1);
-        nab=0;
-        for e=i1:j1
-            for f=1:nfaces
-                ee=itopo(f,e);
-                if(j2<=ee && ee<=i2 && iflux(f,e)>0)
-                    nab=nab+1;
-                end
-            end
-        end
-        nba=0;
-        for e=j2:i2
-            for f=1:nfaces
-                ee=itopo(f,e);
-                if(i1<=ee && ee<=j1 && iflux(f,e)>0)
-                    nba=nba+1;
-                end
-            end
-        end
-        if(nba<nab)
-            n1=j2; n2=i2;
-            m1=i1; m2=j1;
-        else    
-            m1=j2; m2=i2;
-            n1=i1; n2=j1;
-        end
-        for e=m1:m2
-            for f=1:nfaces
-                ee=itopo(f,e);
-                if(n1<=ee && ee<=n2 && iflux(f,e)>0)
-                    ff=(itopo(:,ee)==e);
-                    iflux(ff,ee)=0;
-                    itopo(ff,ee)=ee;
-                    iflux(f,e)=0;
-                    itopo(f,e)=e;
-                end
-            end
-        end
-        [itopo,iflux]=cut_loops(itopo,iflux,m1,m2);
-        [itopo,iflux]=cut_loops(itopo,iflux,n1,n2);
-    end
-end
-
+end