From 591cb50b4eb27fb74de534a59fd94828025c19d9 Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Tue, 28 May 2024 20:26:35 -0400 Subject: [PATCH] method added to chunkgraph, small tweaks to tests --- chunkie/@chunkgraph/chunkgraph.m | 11 ++++-- devtools/test/chunkermatapplyTest.m | 13 ++++--- devtools/test/chunkgrphrcipTransmissionTest.m | 2 +- devtools/test/mixedbcTest.m | 35 +++++++++++++------ 4 files changed, 40 insertions(+), 21 deletions(-) diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m index 1d6ca2d..b7be770 100644 --- a/chunkie/@chunkgraph/chunkgraph.m +++ b/chunkie/@chunkgraph/chunkgraph.m @@ -295,12 +295,19 @@ sourceinfo.d2= d2s; sourceinfo.w = ws; end - + function inds = edgeinds(obj,edgelist) + ladr = cumsum([1,obj.echnks.npt]); + inds = []; + for j = 1:length(edgelist) + ej = edgelist(j); + inds = [inds,ladr(ej):(ladr(ej+1)-1)]; + end + end + % defined in other files spmat = build_v2emat(obj) end methods(Static) - end end diff --git a/devtools/test/chunkermatapplyTest.m b/devtools/test/chunkermatapplyTest.m index c488754..211828f 100644 --- a/devtools/test/chunkermatapplyTest.m +++ b/devtools/test/chunkermatapplyTest.m @@ -13,6 +13,7 @@ cparams = []; cparams.eps = 1.0e-4; +cparams.nover = 3; pref = []; pref.k = 16; narms = 5; @@ -38,7 +39,8 @@ dens = fkernsrc.fmm(1e-12,srcinfo,targinfo,strengths); udense = sys*dens; -cormat = chunkermat(chnkr,fkern,struct("corrections",true)); +start = tic; cormat = chunkermat(chnkr,fkern,struct("corrections",true)); +toc(start) sysapply = @(sigma) chunkermatapply(chnkr,fkern,-0.5,sigma,cormat); start = tic; u = sysapply(dens); t1 = toc(start); @@ -47,9 +49,9 @@ fprintf('relative apply error %5.2e\n',relerr); assert(relerr < 1e-13) -sol1 = sys\dens; +start = tic; sol1 = sys\dens; toc(start) -sol2 = gmres(sysapply, dens, [], 1e-14, 100); +start = tic; sol2 = gmres(sysapply, dens, [], 1e-14, 100); toc(start) relerr = norm(sol1-sol2)/norm(sol1); fprintf('relative solve error %5.2e\n',relerr); assert(relerr < 1e-13) @@ -87,7 +89,7 @@ udense = sys*bdry_data; cormat = chunkermat(chnkr,kerns,struct("corrections",true)); -sysapply = @(sigma) chunkermatapply(chnkr,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(chnkr,kerns,eye(2),sigma,cormat); start = tic; u = sysapply(bdry_data); t1 = toc(start); fprintf('%5.2e s : time for matrix free apply\n',t1) @@ -215,9 +217,6 @@ assert(relerr < 1e-10) - - - function [r,d,d2] = sinearc(t,amp,frq) xs = t; ys = amp*sin(frq*t); diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m index 4b7f6d0..1459439 100644 --- a/devtools/test/chunkgrphrcipTransmissionTest.m +++ b/devtools/test/chunkgrphrcipTransmissionTest.m @@ -7,7 +7,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %clear all - +clearvars addpaths_loc(); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/devtools/test/mixedbcTest.m b/devtools/test/mixedbcTest.m index 89a7b3f..11b13eb 100644 --- a/devtools/test/mixedbcTest.m +++ b/devtools/test/mixedbcTest.m @@ -20,17 +20,17 @@ jind = [jind jind + 1]; jind(jind>nverts) = 1; svals = [-ones(1,nverts) ones(1,nverts)]; -edge2verts = sparse(iind,jind,svals,nverts,nverts); -edge2verts = [edge2verts, 0*edge2verts; 0*edge2verts, edge2verts]; + +edgesendverts = [1 2 3 4 5 6 7 8; 2 3 4 1 8 5 6 7]; edir = 1:4; % indices of edges with Dirichlet conditions eneu = 5:8; % indices of edges with Neumann conditions -fchnks = cell(1,size(edge2verts,1)); +fchnks = cell(1,size(edgesendverts,2)); cparams = []; cparams.nover = 2; -[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); +[cgrph] = chunkgraph(verts,edgesendverts,fchnks,cparams); vstruc = procverts(cgrph); rgns = findregions(cgrph); @@ -39,18 +39,23 @@ % dirichlet and neumann test zk = 30; +% scale matrices so diagonal is identity (using RCIP here so required) kerns(length([edir,eneu]),length([edir,eneu])) = kernel(); kerns(edir,edir) = -2*kernel('helm', 'd', zk); % Dirichlet conditions kerns(eneu,edir) = -2*kernel('helm', 'dp', zk); % Neumann conditions -kerns(edir,eneu) = -2*kernel('helm', 's', zk); -kerns(eneu,eneu) = -2*kernel('helm', 'sp', zk); +kerns(edir,eneu) = 2*kernel('helm', 's', zk); +kerns(eneu,eneu) = 2*kernel('helm', 'sp', zk); start = tic; sysmat = chunkermat(cgrph,kerns); t1 = toc(start); fprintf('%5.2e s : time to assemble matrix\n',t1) -sys = eye(size(sysmat,1)) + sysmat; +indsdir = cgrph.edgeinds(edir); +indsneu = cgrph.edgeinds(eneu); +dval = zeros(cgrph.npt,1); +dval(indsdir) = 1; dval(indsneu) = 1; +sys = diag(dval) + sysmat; fkernsrc = kernel('helm','s',zk); sources = [1;1]; @@ -65,19 +70,27 @@ targinfo.d = merge(cgrph.echnks(eneu)).d(:,:); targinfo.n = merge(cgrph.echnks(eneu)).n(:,:); bdrydatan = fkernsrc.fmm(1e-12,srcinfo,targinfo,charges); -bdrydata = [bdrydatad;bdrydatan]; +bdrydata = -[bdrydatad;bdrydatan]; sol1 = sys\bdrydata; cormat = chunkermat(cgrph,kerns,struct("corrections",true)); -sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat); +sysapply = @(sigma) chunkermatapply(cgrph,kerns,1,sigma,cormat,struct("forcefmm",true)); +sysapply2 = @(sigma) chunkermatapply(cgrph,kerns,dval,sigma,cormat); + +sol2 = gmres(sysapply,bdrydata,100,1e-6,100); + +norm(sol1-sol2) xapply1 = sys*bdrydata; xapply2 = sysapply(bdrydata); +xapply3 = sysapply2(bdrydata); relerr = norm(xapply1-xapply2)/norm(xapply1); fprintf('relative matrix free apply error %5.2e\n',relerr); assert(relerr < 1e-10) - +relerr = norm(xapply1-xapply3)/norm(xapply1); +fprintf('relative matrix free apply error %5.2e\n',relerr); +assert(relerr < 1e-10) rmin = min(cgrph.r(:,:)'); rmax = max(cgrph.r(:,:)'); xl = rmax(1)-rmin(1); @@ -100,7 +113,7 @@ % compute layer potential based on oversample boundary start = tic; -fkernd = 2*kernel('helm', 'd', zk); +fkernd = -2*kernel('helm', 'd', zk); iddir = 1:merge(cgrph.echnks(edir)).npt; uscat = chunkerkerneval(merge(cgrph.echnks(edir)),fkernd,sol1(iddir), ... targets(:,in));