diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index 0d2de0f3..d3af4f76 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -53,19 +53,18 @@ opts = []; end -savedeep = false; -if isfield(opts,'savedeep') - savedeep = opts.savedeep; +savedepth = 10; +if isfield(opts,'rcip_savedepth') + savedepth = opts.rcip_savedepth; end +savedepth = max(savedepth,0); +savedepth = min(savedepth,nsub); -rcipsav.savedeep = savedeep; - -if savedeep - rcipsav.R = cell(nsub+1,1); - rcipsav.MAT = cell(nsub,1); - rcipsav.chnkrlocals = cell(nsub,1); -end +rcipsav.savedepth = savedepth; +rcipsav.R = cell(nsub+1,1); +rcipsav.MAT = cell(nsub,1); +rcipsav.chnkrlocals = cell(nsub,1); % grab only those kernels relevant to this vertex @@ -260,12 +259,15 @@ if level==1 % Dumb and lazy initializer for R, for now %R=eye(nR); R = inv(MAT(starL,starL)); - if savedeep + if level >= nsub-savedepth+1 rcipsav.R{1} = R; end end + if savedepth < nsub && level == nsub-savedepth+1 + rcipsav.R{level} = R; + end R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS); - if savedeep + if level >= nsub-savedepth+1 rcipsav.R{level+1} = R; rcipsav.MAT{level} = MAT(starL,circL); rcipsav.chnkrlocals{level} = merge(chnkrlocal); diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m index 1ed8b26f..6c544fb4 100644 --- a/chunkie/+chnk/+rcip/rhohatInterp.m +++ b/chunkie/+chnk/+rcip/rhohatInterp.m @@ -3,9 +3,9 @@ % density rhohat to the requested depth using the backward recursion % % When using an RCIP preconditioner (in the style of eq 34 of the RCIP -% tutorial), the resulting coarse level density is -% accurate for interactions separated from the vertex -% *at the coarse scale*. +% tutorial), the resulting coarse level density is accurate for +% interactions separated from the vertex *at the coarse scale*. +% % By running the recursion for the RCIP preconditioner backwards, an % equivalent density can be reconstructed on the fine mesh which is % appropriate for closer interactions (at a distance well-separated from @@ -53,26 +53,42 @@ % wts - a set of weights for integrating functions sampled at these % points - % author: Travis Askham +% author: Travis Askham -rhohatinterp = []; -srcinfo = []; -wts = []; - -k = rcipsav.k; -ndim = rcipsav.ndim; nedge = rcipsav.nedge; Pbc = rcipsav.Pbc; -PWbc = rcipsav.PWbc; -starL = rcipsav.starL; -starL1 = rcipsav.starL1; -starS = rcipsav.starS; -circL = rcipsav.circL; -circL1 = rcipsav.circL1; -circS = rcipsav.circS; -ilist = rcipsav.ilist; + +starL1 = sort(rcipsav.starL1); +starS = sort(rcipsav.starS); + +circL1 = sort(rcipsav.circL1); +circS = sort(rcipsav.circS); + +nrho = numel([circS,starS]); +ndens = numel(rhohat)/nrho; +rhohat = reshape(rhohat,[nrho, ndens]); + nsub = rcipsav.nsub; +rhohatinterp = cell(nedge,1); +srcinfo = cell(nedge,1); +wts = cell(nedge,1); + +ileftright = rcipsav.ileftright; + +% figure out which edge the indices belong to +% THIS ASSUMES WE HAVE THE SAME ORDER ON EDGES +circSedge = cell(nedge,1); ncS = numel(circS)/nedge; +circL1edge = cell(nedge,1); ncL1 = numel(circL1)/nedge; +starSedge = cell(nedge,1); nsS = numel(starS)/nedge; +starL1edge = cell(nedge,1); nsL1 = numel(starL1)/nedge; +for j = 1:nedge + circSedge{j} = circS((j-1)*ncS+1:j*ncS); + circL1edge{j} = circL1((j-1)*ncL1+1:j*ncL1); + starSedge{j} = starS((j-1)*nsS+1:j*nsS); + starL1edge{j} = starL1((j-1)*nsL1+1:j*nsL1); +end + if nargin < 3 ndepth = nsub; end @@ -84,58 +100,58 @@ ndepth = nsub; end -savedeep = rcipsav.savedeep; +savedepth = rcipsav.savedepth; -if savedeep +if ndepth <= savedepth % all relevant quantities are stored, just run backward recursion rhohat0 = rhohat; - rhohatinterp = [rhohatinterp; rhohat0(circS)]; - r = []; - d = []; - d2 = []; - n = []; - h = []; cl = rcipsav.chnkrlocals{nsub}; wt = weights(cl); - r = [r, cl.rstor(:,circL1)]; - d = [d, cl.dstor(:,circL1)]; - d2 = [d2, cl.d2stor(:,circL1)]; - n = [n, cl.nstor(:,circL1)]; - wts = [wts; wt(circL1(:))]; + + for j = 1:nedge + rhohatinterp{j} = rhohat0(circSedge{j},:); + srcinfo{j}.r = cl.rstor(:,circL1edge{j}); + srcinfo{j}.d = cl.dstor(:,circL1edge{j}); + srcinfo{j}.d2 = cl.d2stor(:,circL1edge{j}); + srcinfo{j}.n = cl.nstor(:,circL1edge{j}); + wts{j} = wt(circL1edge{j}); + end R0 = rcipsav.R{nsub+1}; for i = 1:ndepth R1 = rcipsav.R{nsub-i+1}; MAT = rcipsav.MAT{nsub-i+1}; rhotemp = R0\rhohat0; - rhohat0 = R1*(Pbc*rhotemp(starS) - MAT*rhohat0(circS)); + rhohat0 = R1*(Pbc*rhotemp(starS,:) - MAT*rhohat0(circS,:)); if i == ndepth - rhohatinterp = [rhohatinterp; rhohat0]; - r = [r, cl.rstor(:,starL1)]; - d = [d, cl.dstor(:,starL1)]; - d2 = [d2, cl.d2stor(:,starL1)]; - n = [n, cl.nstor(:,starL1)]; wt = weights(cl); - - wts = [wts; wt(starL1(:))]; + for j = 1:nedge + if ileftright(j) == 1 + rhohatinterp{j} = [rhohatinterp{j}; rhohat0([circSedge{j} starSedge{j}],:)]; + else + rhohatinterp{j} = [rhohatinterp{j}; rhohat0([starSedge{j} circSedge{j}],:)]; + end + srcinfo{j}.r = [srcinfo{j}.r, cl.rstor(:,starL1edge{j})]; + srcinfo{j}.d = [srcinfo{j}.d, cl.dstor(:,starL1edge{j})]; + srcinfo{j}.d2 = [srcinfo{j}.d2, cl.d2stor(:,starL1edge{j})]; + srcinfo{j}.n = [srcinfo{j}.n, cl.nstor(:,starL1edge{j})]; + wts{j} = [wts{j}, wt(starL1edge{j})]; + end else cl = rcipsav.chnkrlocals{nsub-i}; - rhohatinterp = [rhohatinterp; rhohat0(circS)]; - r = [r, cl.rstor(:,circL1)]; - d = [d, cl.dstor(:,circL1)]; - d2 = [d2, cl.d2stor(:,circL1)]; - n = [n, cl.nstor(:,circL1)]; wt = weights(cl); - wts = [wts; wt(circL1(:))]; + for j = 1:nedge + rhohatinterp{j} = [rhohatinterp{j}; rhohat0(circSedge{j},:)]; + srcinfo{j}.r = [srcinfo{j}.r, cl.rstor(:,circL1edge{j})]; + srcinfo{j}.d = [srcinfo{j}.d, cl.dstor(:,circL1edge{j})]; + srcinfo{j}.d2 = [srcinfo{j}.d2, cl.d2stor(:,circL1edge{j})]; + srcinfo{j}.n = [srcinfo{j}.n, cl.nstor(:,circL1edge{j})]; + wts{j} = [wts{j}, wt(circL1edge{j})]; + end end R0 = R1; end - srcinfo.r = r; - srcinfo.d = d; - srcinfo.d2 = d2; - srcinfo.n = n; - end diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 5a695487..1c48c6f3 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -65,6 +65,7 @@ % corrections for near corners if input chnkobj is % of type chunkergraph % opts.rcip_ignore = [], list of vertices to ignore in rcip +% opts.rcip_savedepth = (10), depth to save rcip info % opts.nsub_or_tol = (40) specify the level of refinements in rcip % or a tolerance where the number of levels is given by % ceiling(log_{2}(1/tol^2)); @@ -83,6 +84,9 @@ % Optional output % opts - with the updated opts structure which stores the relevant % quantities in opts.auxquads. +% rcipsav - precomputed structure of rcip data at corners +% for subsequent postprocessing of the solution at targets close +% to the corner % % Examples: % sysmat = chunkermat(chnkr,kern); % standard options @@ -135,6 +139,7 @@ isrcip = true; rcip_ignore = []; nsub = 40; +rcip_savedepth = 10; adaptive_correction = false; sing = 'log'; @@ -162,6 +167,9 @@ if(isfield(opts,'rcip')) isrcip = opts.rcip; end +if(isfield(opts,'rcip_savedepth')) + rcip_savedepth = opts.rcip_savedepth; +end if(isfield(opts,'rcip_ignore')) rcip_ignore = opts.rcip_ignore; if ~isrcip && isempty(rcip_ignore) @@ -466,6 +474,8 @@ npt_all = horzcat(chnkobj.echnks.npt); [~,nv] = size(chnkobj.verts); ngl = chnkrs(1).k; + + rcipsav = cell(nv,1); for ivert=setdiff(1:nv,rcip_ignore) clist = chnkobj.vstruc{ivert}{1}; @@ -512,11 +522,14 @@ optsrcip = opts; optsrcip.nonsmoothonly = false; optsrcip.corrections = false; + optsrcip.rcip_savedepth = rcip_savedepth; - R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ... + [R,rcipsav{ivert}] = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ... Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... sbclmat,sbcrmat,lvmat,rvmat,u,optsrcip); - + + rcipsav{ivert}.starind = starind; + sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim); if (~nonsmoothonly) @@ -574,6 +587,11 @@ vsysmat = [vsysmat;sysmat_tmp(:)]; end end + + if nargout > 2 + varargout{2} = rcipsav; + end + end diff --git a/devtools/test/rcipTest.m b/devtools/test/rcipTest.m index 4592c21d..d1f484b1 100644 --- a/devtools/test/rcipTest.m +++ b/devtools/test/rcipTest.m @@ -1,8 +1,6 @@ rcipTest0(); - function rcipTest0() - %RCIPTEST % % This file tests the rcip routines for solving the exterior dirichlet @@ -54,12 +52,22 @@ function rcipTest0() fcurve = cell(1,ncurve); figure(1) clf; hold on; +chnkrcell = cell(ncurve,1); for icurve = 1:ncurve fcurve{icurve} = @(t) circulararc(t,cpars{icurve}); chnkr(icurve) = chunkerfuncuni(fcurve{icurve},nch(icurve),cparams{icurve},pref); plot(chnkr(icurve)); quiver(chnkr(icurve)); + chnkrcell{icurve} = chnkr(icurve); end +cends = chunkends(chnkrcell{icurve}); +vert = [cends(:,1,1), cends(:,2,end)]; + +cg = chunkgraph(vert,[2 1; 1 2],chnkrcell); % equivalent chunkgraph + +%% + + iedgechunks = [1 2; 1 chnkr(2).nch]; isstart = [1 0]; nedge = 2; @@ -85,7 +93,7 @@ function rcipTest0() ts = 0.0+2*pi*rand(1,nt); targets = [cos(ts);sin(ts)]; targets = 0.2*targets; -targets(:,1) = [0.95;0]; +targets(:,1) = [0.999999;0]; targets(:,2) = [0,0.36]; targets(:,3) = [-0.95;0]; targets(:,2) = [0,0.36]; @@ -154,8 +162,6 @@ function rcipTest0() nsub = 40; -opts = []; - for icorner=1:ncorner clist = corners{icorner}.clist; @@ -180,11 +186,11 @@ function rcipTest0() [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = chnk.rcip.setup(ngl,ndim, ... nedge,isstart); - opts_use = []; iedgechunks = corners{icorner}.iedgechunks; optsrcip = []; - optsrcip.savedeep = true; + optsrcip.nsub = nsub; + optsrcip.rcip_savedepth = Inf; tic; [R{icorner},rcipsav{icorner}] = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ... Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,..., [],[],[],[],[],optsrcip); @@ -194,18 +200,48 @@ function rcipTest0() end sysmat = sysmat + eye(np); +opts = []; +opts.nsub = nsub; +opts.rcip_savedepth = Inf; +[sysmat2,~,rcipsav] = chunkermat(cg,fkern,opts); + +sysmat2 = sysmat2 + eye(np); + [sol,flag,relres,iter] = gmres(sysmat,ubdry,np,eps*20,np); +[sol2,flag,relres,iter] = gmres(sysmat2,ubdry,np,eps*20,np); % % interpolate to fine grid -ndepth = 10; +ndepth = 20; cor = cell(1,ncorner); for icorner = 1:ncorner - solhat = sol(starinds{icorner}); - [solhatinterp,srcinfo,wts] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth); + starindtmp = rcipsav{icorner}.starind; + solhat = sol(starindtmp); + [solhatinterpcell,srcinfocell,wtscell] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth); + + rtmp = []; + dtmp = []; + d2tmp = []; + ntmp = []; + datatmp = []; + solhatinterp = []; + wts = []; + for j = 1:length(srcinfocell) + rtmp = [rtmp, srcinfocell{j}.r]; + dtmp = [dtmp, srcinfocell{j}.d]; + d2tmp = [d2tmp, srcinfocell{j}.d2]; + ntmp = [ntmp, srcinfocell{j}.n]; + if isfield(srcinfocell{j},'data') + datatmp = [datatmp, srcinfocell{j}.data]; + end + solhatinterp = [solhatinterp; solhatinterpcell{j}(:)]; + wts = [wts; wtscell{j}(:)]; + + end + srcinfo = struct('r',rtmp,'d',dtmp,'d2',d2tmp,'n',ntmp,'data',datatmp); targtemp = targets(:,:) - rcipsav{icorner}.ctr(:,1); targinfo = []; @@ -242,6 +278,13 @@ function rcipTest0() fprintf('relative frobenius error %5.2e\n',relerr); fprintf('relative l_inf/l_1 error %5.2e\n',relerr2); +assert(relerr < 1e-10) +assert(relerr2 < 1e-10) + + + +end + %% @@ -255,9 +298,6 @@ function rcipTest0() -end - - function [r,d,d2] = circulararc(t,cpars) %%circulararc % return position, first and second derivatives of a circular arc