From 3293da094cd1890c8d01ecb9cef45e2c4d004763 Mon Sep 17 00:00:00 2001 From: pbrubeck Date: Sat, 30 Nov 2019 10:41:26 +0000 Subject: [PATCH] greedy coloring --- Stabilization/adf2.m | 8 ++-- Stabilization/bubfilt.m | 8 ++-- Stabilization/color_greedy.m | 15 ++++++ Stabilization/get_graph.m | 51 -------------------- Stabilization/my_gmres.m | 87 ++++++++++++++++++++++++++++++++++ Stabilization/toposort_loops.m | 64 +++++-------------------- 6 files changed, 120 insertions(+), 113 deletions(-) create mode 100644 Stabilization/color_greedy.m create mode 100644 Stabilization/my_gmres.m diff --git a/Stabilization/adf2.m b/Stabilization/adf2.m index 8fe0540..c449a02 100644 --- a/Stabilization/adf2.m +++ b/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 diff --git a/Stabilization/bubfilt.m b/Stabilization/bubfilt.m index 62a789d..12a612b 100644 --- a/Stabilization/bubfilt.m +++ b/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 \ No newline at end of file diff --git a/Stabilization/color_greedy.m b/Stabilization/color_greedy.m new file mode 100644 index 0000000..6732a61 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Stabilization/get_graph.m b/Stabilization/get_graph.m index 7fcc298..d3e9cab 100644 --- a/Stabilization/get_graph.m +++ b/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 diff --git a/Stabilization/my_gmres.m b/Stabilization/my_gmres.m new file mode 100644 index 0000000..04223b0 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Stabilization/toposort_loops.m b/Stabilization/toposort_loops.m index fbbbcc6..d50b4c6 100644 --- a/Stabilization/toposort_loops.m +++ b/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(nba0) - 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 \ No newline at end of file