diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index 0f087f16..cd5b1cba 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -60,12 +60,35 @@ savedepth = max(savedepth,0); savedepth = min(savedepth,nsub); +% open_arc_eye Bruno-Lintner structure makes the (1,1) block of the local +% system matrix near-singular; use a Schur-complement-on-the-other-block +% inverse (chnk.rcip.myinv) inside SchurBana to avoid ~1e-23 RCOND warnings. +use_myinv = isfield(opts,'open_arc_eye') && opts.open_arc_eye; +% sub-block size for the myinv 2x2-block-of-(d x d) Schur structure +% (1 for helmos, 2 for open-arc Stokes mobility). +myinv_subdim = 1; +if isfield(opts,'open_arc_eye_subdim') && ~isempty(opts.open_arc_eye_subdim) + myinv_subdim = opts.open_arc_eye_subdim; +end + rcipsav.savedepth = savedepth; +rcipsav.use_myinv = use_myinv; +rcipsav.myinv_subdim = myinv_subdim; rcipsav.R = cell(nsub+1,1); rcipsav.MAT = cell(nsub,1); rcipsav.chnkrlocals = cell(nsub,1); +% Optional: cache full per-level MAT to enable srhs_recurse_only without +% rebuilding the (whole) chunkermat. Used by callers that re-run the +% singular-RHS recursion many times with different srhs_eval (e.g., +% I2I solves with multiple data_in inputs sharing the same operator). +cache_mat_for_srhs = isfield(opts,'cache_mat_for_srhs') && opts.cache_mat_for_srhs; +if cache_mat_for_srhs + rcipsav.MAT_full = cell(nsub,1); + rcipsav.chnkrlocals_array = cell(nsub,1); % un-merged (array of nedge) +end + % grab only those kernels relevant to this vertex if(size(fkern)==1) @@ -83,6 +106,25 @@ end rcipsav.fkernlocal = fkernlocal; +rcipsav.iedgechunks = iedgechunks; % needed for external srhs_recurse_only +rcipsav.vert0 = vert0; % corner vertex (global coord) + +% Singular-RHS RCIP support (Helsing & Karlsson, eq. 28). When opts.srhs_eval +% is provided, run the r_f^star vector recursion in tandem with the standard +% R recursion. The callback returns f at the 96*ndim GL nodes of the level-i +% type-b 6-panel local mesh (in global coordinates). +srhs_eval = []; +if isfield(opts,'srhs_eval') && ~isempty(opts.srhs_eval) + srhs_eval = opts.srhs_eval; +end +rcipsav.has_srhs = ~isempty(srhs_eval); + +% Debug instrumentation for the singular-RHS recursion. When set, prints +% per-level diagnostics and saves a per-level trace into rcipsav.srhs_trace. +srhs_debug = isfield(opts,'srhs_debug') && opts.srhs_debug; +if srhs_debug + rcipsav.srhs_trace = {}; +end % get coefficients of recentered edge chunks and figure out orientation @@ -190,6 +232,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % begin recursion proper +r_f_star = []; % singular-RHS recursion state (only used when srhs_eval set) + h0=ones(nedge,1); for level=1:nsub h = h0/2^(nsub-level); @@ -252,13 +296,23 @@ % test for opdims ~= [1,1] [MAT,opts] = chunkermat(chnkrlocal, fkernlocal, opts, ilistl); + % Apply cross-edge wLCHS at every recursive level (not just level==1). + % Each refinement level has its own cross-edge near interactions across + % the corner, all of which need the kernel-split correction. + if nedge > 1 + MAT = add_level1_cross_wlchs(MAT, chnkrlocal, ndim, opts, fkernlocal); + end % MAT = eye(nsys) + MAT; if level==1 % Dumb and lazy initializer for R, for now - %R=eye(nR); - R = inv(MAT(starL,starL)); + %R=eye(nR); + if use_myinv + R = chnk.rcip.myinv(MAT(starL,starL), myinv_subdim); + else + R = inv(MAT(starL,starL)); + end if level >= nsub-savedepth+1 rcipsav.R{1} = R; end @@ -266,12 +320,60 @@ if savedepth < nsub && level == nsub-savedepth+1 rcipsav.R{level} = R; end - R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS); + + % --- Singular-RHS RCIP forward recursion, vector form (Helsing 2022 eq 28) --- + % Run before the SchurBana update of R, since Rfstep needs R_{i-1}. + if ~isempty(srhs_eval) + % Translate the local 6-panel mesh to global coords for the callback. + chnkrlocal_global = chnkrlocal; + for ie = 1:nedge + chnkrlocal_global(ie) = chnkrlocal(ie) + ctr(:,ie); + end + b_ib = srhs_eval(chnkrlocal_global, vert0, iedgechunks(1,:), level, ndim); + b_ib = b_ib(:); + if numel(b_ib) ~= nsys + error(['CHNK.RCIP.RCOMPCHUNK: srhs_eval returned %d entries, ' ... + 'expected %d (3*k*nedge*ndim).'], numel(b_ib), nsys); + end + if level == 1 + % Lazy init: r_{f,0}^star = R_0 * b_{1b}(starL). + r_f_star = R * b_ib(starL); + end + r_f_star_new = chnk.rcip.Rfstep(MAT, R, r_f_star, b_ib(circL), ... + Pbc, PWbc, starL, circL, starS, circS, use_myinv); + + if srhs_debug + tr.level = level; + tr.norm_R = norm(R); + tr.norm_M = norm(MAT); + tr.norm_b_starL = norm(b_ib(starL)); + tr.norm_b_circL = norm(b_ib(circL)); + tr.norm_r_f_star = norm(r_f_star_new); + tr.b_ib = b_ib; % full level-i b vector + tr.r_f_star_new = r_f_star_new; + tr.r = arrayfun(@(c) c.r, chnkrlocal_global, 'UniformOutput', false); + rcipsav.srhs_trace{end+1} = tr; %#ok + end + end + + R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS,use_myinv,myinv_subdim); if level >= nsub-savedepth+1 rcipsav.R{level+1} = R; rcipsav.MAT{level} = MAT(starL,circL); rcipsav.chnkrlocals{level} = merge(chnkrlocal); + if cache_mat_for_srhs + rcipsav.MAT_full{level} = MAT; + rcipsav.chnkrlocals_array{level} = chnkrlocal; + end end + + if ~isempty(srhs_eval) + r_f_star = r_f_star_new; + end +end + +if ~isempty(srhs_eval) + rcipsav.r_f_star = r_f_star; end end @@ -291,3 +393,173 @@ end +function MAT = add_level1_cross_wlchs(MAT, chnkrlocal, ndim, opts, fkernlocal) +% Apply off-edge near (cross-edge) corrections at first RCIP level. +% +% Two modes supported: +% - opts.forcewlchs is struct with explicit per-(r,c) layout: same layout +% used for all edge pairs (legacy testhelmos / matrix-valued op path). +% - opts.forcewlchs == true (boolean) AND fkernlocal is a kernel cell +% array: per-(target_edge, source_edge) layout is auto-derived by +% inspecting fkernlocal(it, is) (handles 'd', 's', 'sp', 'dp', 'c'/'combined'). + +if ~isfield(opts,'forcewlchs') + return +end +fw = opts.forcewlchs; +have_struct = isstruct(fw) && isfield(fw,'layout') && isfield(fw,'zk'); +have_bool = (islogical(fw) || (isnumeric(fw) && isscalar(fw))) && fw && ... + nargin >= 5 && ~isempty(fkernlocal); +if ~have_struct && ~have_bool + return +end + +nedge = numel(chnkrlocal); +npts = arrayfun(@(c) c.npt, chnkrlocal); +offset = [0, cumsum(npts(1:end-1))*ndim]; + +if have_struct + layout = fw.layout; + zk = fw.zk; + if ~iscell(layout) || size(layout,1) ~= ndim || size(layout,2) ~= ndim + return + end +end + +for it = 1:nedge + for is = 1:nedge + if it == is + continue + end + for r = 1:ndim + for c = 1:ndim + if have_struct + entry = layout{r,c}; + if isempty(entry); continue; end + if iscell(entry) + typ = entry{1}; coef = entry{2}; + elseif isstruct(entry) + typ = entry.type; coef = entry.coef; + else + continue + end + if strcmpi(typ,'zero'); continue; end + [types, coefs] = decompose_kernel_spec(typ, coef); + else + % Boolean mode: derive from fkernlocal(it, is). When + % fkernlocal is a single (scalar) kernel object, use it + % uniformly for every (it, is) pair. + if isscalar(fkernlocal) + bk = fkernlocal; + else + bk = fkernlocal(it, is); + end + if ~isa(bk, 'kernel'); continue; end + if strcmpi(bk.name, 'zeros') || ... + (isprop(bk,'iszero') && bk.iszero); continue; end + if ~isfield(bk.params, 'zk') || isempty(bk.params.zk); continue; end + zk = bk.params.zk; + [types, coefs] = decompose_kernel_obj(bk); + if isempty(types); continue; end + end + + delta_total = []; + for kk = 1:numel(types) + typ_k = types{kk}; + coef_k = coefs(kk); + if coef_k == 0; continue; end + delta = chnk.kernsplit.helm2d_near_correction( ... + chnkrlocal(is), chnkrlocal(it), typ_k, zk); + if nnz(delta) == 0; continue; end + if isempty(delta_total) + delta_total = coef_k * full(delta); + else + delta_total = delta_total + coef_k * full(delta); + end + end + if isempty(delta_total); continue; end + + rows = offset(it) + r + (0:npts(it)-1)*ndim; + cols = offset(is) + c + (0:npts(is)-1)*ndim; + MAT(rows, cols) = MAT(rows, cols) + delta_total; + end + end + end +end +end + + +function [types, coefs] = decompose_kernel_spec(typ, coef) +% Convert a layout spec (typ, coef) into a list of scalar (type, coef) +% pairs supported by helm2d_near_correction. +typ = lower(typ); +if any(strcmp(typ, {'c','combined'})) + if numel(coef) ~= 2 + error('add_level1_cross_wlchs: c kernel needs coefs [c1, c2]'); + end + types = {'d','s'}; + coefs = [coef(1), coef(2)]; +elseif any(strcmp(typ, {'sc','spcombined'})) + if numel(coef) ~= 2 + error('add_level1_cross_wlchs: sc kernel needs coefs [c_sp, c_s]'); + end + types = {'sp','s'}; + coefs = [coef(1), coef(2)]; +else + types = {typ}; + coefs = coef; +end +end + + +function [types, coefs] = decompose_kernel_obj(kern) +% Inspect a kernel object and return list of (type, coef) pairs equivalent +% to its bare-kernel evaluation (relevant for cross-edge near correction). +% Returns empty if kernel type is unsupported. +types = {}; coefs = []; +t = lower(kern.type); +switch t + case {'s','single'} + c = wlchs_extract_scalar_local(kern, 's'); + types = {'s'}; coefs = c; + case {'d','double'} + c = wlchs_extract_scalar_local(kern, 'd'); + types = {'d'}; coefs = c; + case {'sp','sprime'} + c = wlchs_extract_scalar_local(kern, 'sp'); + types = {'sp'}; coefs = c; + case {'dp','dprime','t'} + c = wlchs_extract_scalar_local(kern, 'dp'); + types = {'dp'}; coefs = c; + case {'c','combined'} + if isfield(kern.params,'coefs') && numel(kern.params.coefs) == 2 + types = {'d','s'}; + coefs = [kern.params.coefs(1), kern.params.coefs(2)]; + end + case {'sc','spcombined'} + if isfield(kern.params,'coefs') && numel(kern.params.coefs) == 2 + types = {'sp','s'}; + coefs = [kern.params.coefs(1), kern.params.coefs(2)]; + end + otherwise + % unsupported (e.g., custom sum-kernel without recognized type) +end +end + + +function coef = wlchs_extract_scalar_local(kern, t) +% Probe the kernel eval at a fixed src/tgt pair, divide by the analytic +% bare-kernel value to recover the scalar multiplier (handles 2*dkern etc). +zk = kern.params.zk; +src_probe = struct('r',[0;0],'n',[1;0],'d',[1;0],'d2',[0;0]); +tgt_probe = struct('r',[1;0],'n',[1;0],'d',[1;0],'d2',[0;0]); +val = kern.eval(src_probe, tgt_probe); +switch lower(t) + case 's', base = 0.25i * besselh(0, zk); + case 'd', base = 0.25i * zk * besselh(1, zk); + case 'sp', base = kernel('helm','sp',zk).eval(src_probe, tgt_probe); + case 'dp', base = kernel('helm','dp',zk).eval(src_probe, tgt_probe); + otherwise, base = 1; +end +coef = val / base; +end diff --git a/chunkie/+chnk/+rcip/Rfstep.m b/chunkie/+chnk/+rcip/Rfstep.m new file mode 100644 index 00000000..aab55d26 --- /dev/null +++ b/chunkie/+chnk/+rcip/Rfstep.m @@ -0,0 +1,78 @@ +function r_f_new = Rfstep(MAT, A, r_f_prev, b_circ, Pbc, PWbc, ... + starL, circL, starS, circS, use_myinv, myinv_subdim) +%CHNK.RCIP.RFSTEP one level of the singular-RHS RCIP forward recursion. +% +% Implements the vector form of the R_f recursion (Helsing & Karlsson, +% bgkw6 / J. Comput. Phys. 2022, eq. 28) so that the inversion of the +% previous-level R cancels analytically: +% +% r_{f,i} = P_W^T * M_i^{-1} * y_i +% +% M_i = [ A^{-1} K_b(starL,circL) ] +% [ K_b(circL,starL) I + K_b(circL,circL) ] +% +% y_i(starL) = A^{-1} * r_{f,i-1} (avoided -- see derivation below) +% y_i(circL) = b_{ib}^circ +% +% Schur block-inverse on the (starL,starL) block plus cancellation gives: +% +% VA = K_b(circL,starL) * r_{f,i-1} % NB: the A^{-1}*A cancels +% AU = A * K_b(starL,circL) +% S = MAT(circL,circL) - K_b(circL,starL)*A*K_b(starL,circL) +% inner = S \ (b_{ib}^circ - VA) +% x(starL) = r_{f,i-1} - AU * inner +% x(circL) = inner +% +% r_{f,i}(starS) = PWbc' * x(starL) +% r_{f,i}(circS) = x(circL) +% +% No explicit inversion of A is needed — the user passes r_{f,i-1} directly +% rather than A^{-1}*r_{f,i-1}. +% +% Inputs: +% MAT - 3*k*ndim*nedge x ... type-b system matrix at level i (I+K_b) +% A - 2*k*ndim*nedge square: previous R, R_{i-1} +% r_f_prev - 2*k*ndim*nedge vector: previous r_f^star +% b_circ - k*ndim*nedge vector: f at the circL nodes of level-i type-b +% Pbc, PWbc - prolongation matrices from chnk.rcip.Pbcinit (level-c->b), +% sizes 2*k*ndim*nedge x k*ndim*nedge. (Pbc unused here but +% kept in the signature for symmetry with SchurBana.) +% starL, circL - bad/good index sets into MAT (size 2*k*ndim*nedge and +% k*ndim*nedge respectively) +% starS, circS - bad/good index sets into the output r_f^star (each +% k*ndim*nedge) +% use_myinv - optional flag: use chnk.rcip.myinv for the Schur complement +% (matches the open-arc case in SchurBana) +% +% Output: +% r_f_new - 2*k*ndim*nedge vector: r_{f,i}^star +% +% See also: chnk.rcip.SchurBana, chnk.rcip.Rcompchunk + +if nargin < 11 || isempty(use_myinv) + use_myinv = false; +end +if nargin < 12 || isempty(myinv_subdim) + myinv_subdim = 1; +end + +K_cs_sL = MAT(circL, starL); +K_sL_cs = MAT(starL, circL); + +VA = K_cs_sL * r_f_prev; % cancels A^{-1}: K_b(circL,starL)*r_{f,i-1} +AU = A * K_sL_cs; % R_{i-1} * K_b(starL,circL) +Mschur = MAT(circL, circL) - K_cs_sL * A * K_sL_cs; +if use_myinv + inner = chnk.rcip.myinv(Mschur, myinv_subdim) * (b_circ - VA); +else + inner = Mschur \ (b_circ - VA); +end + +x_starL = r_f_prev - AU * inner; +x_circL = inner; + +r_f_new = zeros(numel(starS) + numel(circS), 1); +r_f_new(starS) = PWbc' * x_starL; +r_f_new(circS) = x_circL; + +end diff --git a/chunkie/+chnk/+rcip/SchurBana.m b/chunkie/+chnk/+rcip/SchurBana.m index 8d02b71d..40945b43 100644 --- a/chunkie/+chnk/+rcip/SchurBana.m +++ b/chunkie/+chnk/+rcip/SchurBana.m @@ -1,4 +1,4 @@ - function A=SchurBana(P,PW,K,A,starL,circL,starS,circS) + function A=SchurBana(P,PW,K,A,starL,circL,starS,circS,use_myinv,myinv_subdim) %CHNK.RCIP.SCHURBANA block inverse used in RCIP % % uses a matrix block inversion formula to recursively compute the @@ -63,10 +63,21 @@ % author: Johan Helsing (part of the RCIP tutorial) + if nargin < 9 || isempty(use_myinv) + use_myinv = false; + end + if nargin < 10 || isempty(myinv_subdim) + myinv_subdim = 1; + end VA=K(circL,starL)*A; PTA=PW'*A; PTAU=PTA*K(starL,circL); - DVAUI=inv(K(circL,circL)-VA*K(starL,circL)); + Mschur = K(circL,circL)-VA*K(starL,circL); + if use_myinv + DVAUI=chnk.rcip.myinv(Mschur, myinv_subdim); + else + DVAUI=inv(Mschur); + end DVAUIVAP=DVAUI*(VA*P); A(starS,starS)=PTA*P+PTAU*DVAUIVAP; A(circS,circS)=DVAUI; diff --git a/chunkie/+chnk/+rcip/myinv.m b/chunkie/+chnk/+rcip/myinv.m new file mode 100644 index 00000000..a8e9ed35 --- /dev/null +++ b/chunkie/+chnk/+rcip/myinv.m @@ -0,0 +1,53 @@ +function Mi = myinv(M, sub_ndim) +% Compute the inverse of a 2x2 block matrix +% +% M = | A B | +% | C D | +% +% in the (degenerate or near-singular) case where B*C ~ -I and A is +% close to singular. The Schur trick on the OTHER block gives a +% well-conditioned formula: +% +% T = (A - B*D^{-1}*C)^{-1} +% M11 = T +% M12 = -T*B*D^{-1} +% M21 = -D^{-1}*C*T +% M22 = D^{-1} + D^{-1}*C*T*B*D^{-1} +% +% Per-node interleaving: +% sub_ndim = 1 (default, helmos style): components ordered as +% (a_1, b_1, a_2, b_2, ...) so A = M(1:2:end, 1:2:end), etc. +% sub_ndim = d (open-arc Stokes mobility): components ordered as +% (a_1..a_d, b_1..b_d, a_1..a_d, b_1..b_d, ...) so each "block" +% in the 2x2 Schur picture is itself a d x d sub-matrix per node. + +if nargin < 2 || isempty(sub_ndim), sub_ndim = 1; end + +np = size(M, 1); +opdim = 2 * sub_ndim; +Npts = np / opdim; +assert(mod(np, opdim) == 0, 'myinv: matrix size %d not divisible by 2*sub_ndim=%d', np, opdim); + +% Build top vs bottom block index sets. +% sub_ndim=1 -> top = 1:2:end, bot = 2:2:end (matches the original myinv). +block_off = (0:Npts-1) * opdim; +top_idx = reshape(block_off + (1:sub_ndim).', 1, []); +bot_idx = reshape(block_off + (sub_ndim+1:opdim).', 1, []); + +A = M(top_idx, top_idx); +B = M(top_idx, bot_idx); +C = M(bot_idx, top_idx); +D = M(bot_idx, bot_idx); + +T = inv(A - B/D*C); +M11 = T; +M12 = -T*B/D; +M21 = -D\C*T; +M22 = inv(D) + D\C*T*B/D; + +Mi = zeros(np); +Mi(top_idx, top_idx) = M11; +Mi(top_idx, bot_idx) = M12; +Mi(bot_idx, top_idx) = M21; +Mi(bot_idx, bot_idx) = M22; +end diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m index 6c544fb4..7a49b67d 100644 --- a/chunkie/+chnk/+rcip/rhohatInterp.m +++ b/chunkie/+chnk/+rcip/rhohatInterp.m @@ -101,6 +101,11 @@ end savedepth = rcipsav.savedepth; +use_myinv = isfield(rcipsav,'use_myinv') && rcipsav.use_myinv; +myinv_subdim = 1; +if isfield(rcipsav,'myinv_subdim') && ~isempty(rcipsav.myinv_subdim) + myinv_subdim = rcipsav.myinv_subdim; +end if ndepth <= savedepth @@ -123,7 +128,11 @@ for i = 1:ndepth R1 = rcipsav.R{nsub-i+1}; MAT = rcipsav.MAT{nsub-i+1}; - rhotemp = R0\rhohat0; + if use_myinv + rhotemp = chnk.rcip.myinv(R0, myinv_subdim)*rhohat0; + else + rhotemp = R0\rhohat0; + end rhohat0 = R1*(Pbc*rhotemp(starS,:) - MAT*rhohat0(circS,:)); if i == ndepth wt = weights(cl); diff --git a/chunkie/+chnk/+rcip/srhs_modify_rhs.m b/chunkie/+chnk/+rcip/srhs_modify_rhs.m new file mode 100644 index 00000000..9bdd1c45 --- /dev/null +++ b/chunkie/+chnk/+rcip/srhs_modify_rhs.m @@ -0,0 +1,42 @@ +function rhs_mod = srhs_modify_rhs(rhs_bare, rcipsav) +%CHNK.RCIP.SRHS_MODIFY_RHS apply singular-RHS correction at corner-star. +% +% In the singular-RHS RCIP framework (Helsing & Karlsson 2022), expressed +% in chunkie's hat_rho convention where sol = R*v_tilde + r_f_star, the +% only modification to the standard amat\rhs pipeline is at corner-star +% indices: the bare RHS values are replaced by inv(R) * r_f_star_corner. +% +% Mathematically (paper eq. for disc2, rewritten for hat_rho): +% amat * hat_rho = rhs +% amat(cs, cs) = inv(R) +% amat(cs, ~cs) = K(cs, ~cs) +% amat(~cs, ~cs) = I + K(~cs, ~cs) +% rhs(~cs) = b_coarse(~cs) (unchanged) +% rhs( cs) = inv(R) * r_f_star_corner (replaced) +% +% The far-from-corner entries of b_coarse are correct as-is from a bare +% chunkermat(rhs_kerns) * data_in evaluation; only the corner-star entries +% need the R_f-recursion result. +% +% Inputs: +% rhs_bare - bare RHS vector (e.g. chunkermat(rhs_kerns,opts_no_rcip) * data_in). +% rcipsav - per-corner cell array returned as the third output of +% chunkermat when opts.srhs_eval is set. +% +% Output: +% rhs_mod - RHS vector with corner-star entries replaced. +% +% See also: chnk.rcip.Rcompchunk, chnk.rcip.Rfstep, chunkermat. + +rhs_mod = rhs_bare; +for ivert = 1:numel(rcipsav) + rcs = rcipsav{ivert}; + if isempty(rcs); continue; end + if ~isfield(rcs, 'r_f_star') || isempty(rcs.r_f_star); continue; end + if ~isfield(rcs, 'starind') || isempty(rcs.starind); continue; end + + nsub = rcs.nsub; + R_final = rcs.R{nsub+1}; + rhs_mod(rcs.starind) = R_final \ rcs.r_f_star; +end +end diff --git a/chunkie/+chnk/+rcip/srhs_recurse_only.m b/chunkie/+chnk/+rcip/srhs_recurse_only.m new file mode 100644 index 00000000..122edee1 --- /dev/null +++ b/chunkie/+chnk/+rcip/srhs_recurse_only.m @@ -0,0 +1,62 @@ +function r_f_star = srhs_recurse_only(rcipsav_corner, srhs_eval) +%CHNK.RCIP.SRHS_RECURSE_ONLY Re-run the singular-RHS forward recursion at a +% single corner using cached level data, WITHOUT rebuilding any system matrix. +% +% Use this when the BIE operator is fixed (cached in rcipsav by a prior +% chunkermat call with opts.cache_mat_for_srhs=true) but the right-hand +% side data_in changes — e.g., when computing an I2I matrix by solving the +% same BIE for many input modes. Speeds up by skipping the per-iter +% rebuild of bulk + corner-level matrices. +% +% Inputs: +% rcipsav_corner - per-corner rcipsav cell entry produced by Rcompchunk +% with cache_mat_for_srhs=true. Required fields: +% .nsub, .nedge, .ndim, .Pbc, .PWbc, +% .starL, .circL, .starS, .circS, .use_myinv, +% .R{1..nsub+1}, .MAT_full{1..nsub}, +% .chnkrlocals_array{1..nsub}, .ctr, +% .vert0, .iedgechunks +% srhs_eval - function handle (same signature as in Rcompchunk): +% b_ib = srhs_eval(chnkrlocal_global, vert0, +% edge_indices, level, ndim) +% +% Output: +% r_f_star - 2*k*ndim*nedge vector at the corner. +% +% See chnk.rcip.Rfstep, chnk.rcip.Rcompchunk. + +nsub = rcipsav_corner.nsub; +nedge = rcipsav_corner.nedge; +ndim = rcipsav_corner.ndim; +Pbc = rcipsav_corner.Pbc; +PWbc = rcipsav_corner.PWbc; +starL = rcipsav_corner.starL; +circL = rcipsav_corner.circL; +starS = rcipsav_corner.starS; +circS = rcipsav_corner.circS; +use_myinv = rcipsav_corner.use_myinv; +ctr = rcipsav_corner.ctr; +vert0 = rcipsav_corner.vert0; +iedgechunks = rcipsav_corner.iedgechunks; + +r_f_star = []; +for level = 1:nsub + R_lev = rcipsav_corner.R{level}; + MAT_lev = rcipsav_corner.MAT_full{level}; + chnkrlocal_lev = rcipsav_corner.chnkrlocals_array{level}; + + % Translate local mesh to global coords for srhs_eval. + chnkrlocal_global = chnkrlocal_lev; + for ie = 1:nedge + chnkrlocal_global(ie) = chnkrlocal_lev(ie) + ctr(:,ie); + end + b_ib = srhs_eval(chnkrlocal_global, vert0, iedgechunks(1,:), level, ndim); + b_ib = b_ib(:); + + if level == 1 + r_f_star = R_lev * b_ib(starL); + end + r_f_star = chnk.rcip.Rfstep(MAT_lev, R_lev, r_f_star, b_ib(circL), ... + Pbc, PWbc, starL, circL, starS, circS, use_myinv); +end +end diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index d547ba38..a457365f 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -54,9 +54,25 @@ % opts.forcewlchs - if = true, use kernel-split (Helsing-style) % panel quadrature for close evaluation. Supports % Helmholtz (s/d/sp), Laplace (s/d/sp), and elasticity -% (s/d/strac/dalt; opdims=[2 2]) kernels -- kern must be -% a single kernel object. Delivers 10+ digits of accuracy -% at close off-curve targets. (false) +% (s/d/strac/dalt; opdims=[2 2]) kernels — kern must be +% a single kernel object. Delivers 10+ digits of +% accuracy at moderately close off-curve targets +% (d > 0.1*panel_length for elasticity; the genuinely +% singular pieces are wLCHS-corrected, bounded pieces +% degrade for very-close targets). Combine with +% opts.rcipsav for the open-arc / corner case to also +% use kernel-split close eval on the fine-mesh density +% reconstruction. (false) +% opts.rcipsav - cell array of per-vertex RCIP data returned as the +% third output of chunkermat when opts.rcip = true. If set, +% chunkerkerneval performs RCIP-aware field evaluation: +% the density at each vertex's starind is treated as the +% RCIP-transformed coarse density (only valid for far-field +% smooth quadrature) and replaced for close interactions by +% a fine-mesh density reconstructed via +% chnk.rcip.rhohatInterp. Requires chnkobj to be a +% chunkgraph. With nsub=30, accuracy stays at ~1e-12 down +% to target distances ~1e-5 from a vertex. % % output: % fints - opdims(1) x nt array of integral values where opdims is the @@ -145,97 +161,294 @@ if isfield(opts,'proxybylevel'); opts_use.proxybylevel = opts.proxybylevel; end % Kernel-split (Helsing-style) panel quadrature for accurate close -% evaluation of Laplace S, D, S' layer potentials. +% evaluation of Helmholtz / Laplace S, D, S' layer potentials. +% +% Two forms: +% opts.forcewlchs = true -- single scalar kernel +% opts.forcewlchs = struct('zk',zk,'layout', {{ {'zero',0}, {'s',c}, ... }}) +% -- 1xN block kernel (e.g. helmos +% Keval = kernel([Z, c*S])). The +% layout selects which density +% component is the active one and +% which kernsplit type to apply. fw_raw = []; if isfield(opts,'forcewlchs') && ~isempty(opts.forcewlchs) fw_raw = opts.forcewlchs; end -if (islogical(fw_raw) || isnumeric(fw_raw)) && isscalar(fw_raw) +fw_is_block = false; +fw_active_col = 1; +fw_ndim_dens = 1; +if isstruct(fw_raw) && isfield(fw_raw, 'layout') + fw_is_block = true; + L = fw_raw.layout; + if size(L,1) ~= 1 + error("CHUNKERKERNEVAL: forcewlchs struct layout must be 1xN (single output row)"); + end + use_wlchs = true; +elseif (islogical(fw_raw) || isnumeric(fw_raw)) && isscalar(fw_raw) use_wlchs = logical(fw_raw); else use_wlchs = false; end if use_wlchs - % Single kernel: scalar (helm/lap, opdims=[1 1]) or elasticity (opdims=[2 2]). - if size(kern,1) ~= 1 || size(kern,2) ~= 1 || ~isa(kern,'kernel') - error("CHUNKERKERNEVAL: forcewlchs=true requires a single kernel object"); - end - is_helm = strcmpi(kern.name, 'helmholtz'); - is_lap = strcmpi(kern.name, 'laplace'); - is_elast = strcmpi(kern.name, 'elasticity'); - if ~is_helm && ~is_lap && ~is_elast - error("CHUNKERKERNEVAL: forcewlchs supports Helmholtz, Laplace, or elasticity kernels"); - end - if is_elast - if ~ismember(lower(kern.type), {'s','single','d','double','strac','straction','dalt'}) - error("CHUNKERKERNEVAL: forcewlchs (elasticity) supports types s, d, strac, dalt"); + if fw_is_block + % --- struct form: decode the single non-zero entry in the layout + active_col = -1; t_active = ''; coef_active = 0; + for c_ = 1:numel(L) + entry = L{c_}; + if iscell(entry) + et = entry{1}; ec = entry{2}; + elseif isstruct(entry) + et = entry.type; + if isfield(entry,'coef'), ec = entry.coef; + else, ec = entry.coefs; + end + else + error("CHUNKERKERNEVAL: forcewlchs.layout entry must be cell or struct"); + end + if ~strcmpi(et, 'zero') + if active_col > 0 + error("CHUNKERKERNEVAL: forcewlchs.layout must have exactly one non-zero entry"); + end + active_col = c_; t_active = lower(et); coef_active = ec; + end + end + if active_col < 0 + error("CHUNKERKERNEVAL: forcewlchs.layout has no non-zero entries"); + end + if ~ismember(t_active, {'s','single','d','double','sp','sprime'}) + error("CHUNKERKERNEVAL: forcewlchs supports only s/d/sp layers (got %s)", t_active); + end + wlchs_type = t_active; + wlchs_coef = coef_active; + fw_active_col = active_col; + fw_ndim_dens = numel(L); + if isfield(fw_raw,'zk') && ~isempty(fw_raw.zk) + wlchs_family = 'helmholtz'; + wlchs_zk = fw_raw.zk; + else + wlchs_family = 'laplace'; + wlchs_zk = []; end else - if ~ismember(lower(kern.type), {'s','single','d','double','sp','sprime'}) - error("CHUNKERKERNEVAL: forcewlchs currently supports only s/d/sp layers"); + % --- boolean form: single scalar kernel (helm/lap, opdims=[1 1]) or + % a single elasticity kernel (opdims=[2 2]). + if size(kern,1) ~= 1 || size(kern,2) ~= 1 || ~isa(kern,'kernel') + error("CHUNKERKERNEVAL: forcewlchs=true requires a single kernel object"); end - end - wlchs_type = kern.type; - wlchs_family = kern.name; - - src_probe = struct('r',[0;0],'n',[1;0],'d',[1;0],'d2',[0;0]); - tgt_probe = struct('r',[1;0],'n',[1;0],'d',[1;0],'d2',[0;0]); - val_probe = kern.eval(src_probe, tgt_probe); - if is_helm - if ~isfield(kern.params,'zk') || isempty(kern.params.zk) - error("CHUNKERKERNEVAL: forcewlchs Helmholtz requires kern.params.zk"); + is_helm = strcmpi(kern.name, 'helmholtz'); + is_lap = strcmpi(kern.name, 'laplace'); + is_elast = strcmpi(kern.name, 'elasticity'); + if ~is_helm && ~is_lap && ~is_elast + error("CHUNKERKERNEVAL: forcewlchs supports Helmholtz, Laplace, or elasticity kernels"); end - wlchs_zk = kern.params.zk; - if strcmpi(wlchs_type,'s') || strcmpi(wlchs_type,'single') - base_probe = 0.25i * besselh(0, wlchs_zk); - elseif strcmpi(wlchs_type,'d') || strcmpi(wlchs_type,'double') - base_probe = 0.25i * wlchs_zk * besselh(1, wlchs_zk); + if is_elast + if ~ismember(lower(kern.type), {'s','single','d','double','strac','straction','dalt'}) + error("CHUNKERKERNEVAL: forcewlchs (elasticity) supports types s, d, strac, dalt"); + end else - base_probe = -0.25i * wlchs_zk * besselh(1, wlchs_zk); + if ~ismember(lower(kern.type), {'s','single','d','double','sp','sprime'}) + error("CHUNKERKERNEVAL: forcewlchs currently supports only s/d/sp layers"); + end end - wlchs_coef = val_probe / base_probe; - elseif is_lap - wlchs_zk = []; - if strcmpi(wlchs_type,'s') || strcmpi(wlchs_type,'single') + wlchs_type = kern.type; + wlchs_family = kern.name; + + src_probe = struct('r',[0;0],'n',[1;0],'d',[1;0],'d2',[0;0]); + tgt_probe = struct('r',[1;0],'n',[1;0],'d',[1;0],'d2',[0;0]); + val_probe = kern.eval(src_probe, tgt_probe); + if is_helm + if ~isfield(kern.params,'zk') || isempty(kern.params.zk) + error("CHUNKERKERNEVAL: forcewlchs Helmholtz requires kern.params.zk"); + end + wlchs_zk = kern.params.zk; + if strcmpi(wlchs_type,'s') || strcmpi(wlchs_type,'single') + base_probe = 0.25i * besselh(0, wlchs_zk); + elseif strcmpi(wlchs_type,'d') || strcmpi(wlchs_type,'double') + base_probe = 0.25i * wlchs_zk * besselh(1, wlchs_zk); + else + base_probe = -0.25i * wlchs_zk * besselh(1, wlchs_zk); + end + wlchs_coef = val_probe / base_probe; + elseif is_lap + wlchs_zk = []; + if strcmpi(wlchs_type,'s') || strcmpi(wlchs_type,'single') + tgt_probe.r = [2;0]; + val_probe = kern.eval(src_probe, tgt_probe); + base_probe = -log(2)/(2*pi); + else + base_probe = 1/(2*pi); + end + wlchs_coef = val_probe / base_probe; + else % elasticity + if ~isfield(kern.params,'lam') || isempty(kern.params.lam) || ... + ~isfield(kern.params,'mu') || isempty(kern.params.mu) + error("CHUNKERKERNEVAL: forcewlchs (elasticity) requires kern.params.lam and kern.params.mu"); + end + wlchs_lam = kern.params.lam; + wlchs_mu = kern.params.mu; + wlchs_zk = []; + % Probe far from r=0 for log/Cauchy stability tgt_probe.r = [2;0]; val_probe = kern.eval(src_probe, tgt_probe); - base_probe = -log(2)/(2*pi); - else - base_probe = 1/(2*pi); + switch lower(wlchs_type) + case {'s','single'}, et = 's'; + case {'d','double'}, et = 'd'; + case {'strac','straction'}, et = 'strac'; + case 'dalt', et = 'dalt'; + end + base_probe = chnk.elast2d.kern(wlchs_lam, wlchs_mu, src_probe, ... + tgt_probe, et); + wlchs_coef = val_probe(1,1) / base_probe(1,1); end - wlchs_coef = val_probe / base_probe; - else % elasticity - if ~isfield(kern.params,'lam') || isempty(kern.params.lam) || ... - ~isfield(kern.params,'mu') || isempty(kern.params.mu) - error("CHUNKERKERNEVAL: forcewlchs (elasticity) requires kern.params.lam and kern.params.mu"); + end +end + +% RCIP-aware path: zero out coarse-mesh density at vertex starinds, recurse +% to compute the contribution from the rest of the boundary, and then add +% per-vertex contributions evaluated against the reconstructed fine-mesh +% density. +if isfield(opts,'rcipsav') && ~isempty(opts.rcipsav) + if ~icgrph + error("CHUNKERKERNEVAL: opts.rcipsav requires a chunkgraph input"); + end + rcipsav = opts.rcipsav; + + if size(kern,1) > 1 + error("CHUNKERKERNEVAL: opts.rcipsav with multi-region (mk>1) " + ... + "kernel matrices is not yet supported"); + end + + % Density layout: single kernel -> ndim_dens=1; block kernel + % (struct-form forcewlchs) -> ndim_dens = numel(layout) and we pick + % out the active density component for the kernsplit corner step. + if use_wlchs && ~fw_is_block + ndim_dens = size(dens,1)/sum(arrayfun(@(c) c.npt, chnkr)); + if ndim_dens ~= 1 + error("CHUNKERKERNEVAL: forcewlchs=true with rcipsav and a " + ... + "block kernel requires the struct form: " + ... + "opts.forcewlchs = struct('zk',zk,'layout',{...})."); end - wlchs_lam = kern.params.lam; - wlchs_mu = kern.params.mu; - wlchs_zk = []; - % Probe far from r=0 for log/Cauchy stability - tgt_probe.r = [2;0]; - val_probe = kern.eval(src_probe, tgt_probe); - switch lower(wlchs_type) - case {'s','single'}, et = 's'; - case {'d','double'}, et = 'd'; - case {'strac','straction'}, et = 'strac'; - case 'dalt', et = 'dalt'; + end + + dens_zeroed = dens; + for ivert = 1:length(rcipsav) + if isempty(rcipsav{ivert}); continue; end + dens_zeroed(rcipsav{ivert}.starind) = 0; + end + opts_inner = rmfield(opts,'rcipsav'); + if fw_is_block + % The smooth far-field recursive call takes the block kernel + % directly; remove forcewlchs so the recursion doesn't re-enter + % the wlchs branch (which is for close eval only). + opts_inner = rmfield(opts_inner, 'forcewlchs'); + end + [fints,ids] = chunkerkerneval(chnkobj,kern,dens_zeroed,targobj,opts_inner); + + % Per-vertex per-edge fine-mesh contribution + targ_pts = local_target_points(targobj); + targ_n = local_target_normals(targobj); + needs_tgt_normals = use_wlchs && (any(strcmpi(wlchs_type, {'sp','sprime'})) || ... + (strcmpi(wlchs_family,'elasticity') && any(strcmpi(wlchs_type, {'strac','straction'})))); + if needs_tgt_normals + if isempty(targ_n) + error(['CHUNKERKERNEVAL: forcewlchs+rcipsav with sp/strac kernel ' ... + 'requires target normals; pass targobj as a struct with .n']); + end + targ_struct = struct('r', targ_pts, 'n', targ_n); + end + for ivert = 1:length(rcipsav) + if isempty(rcipsav{ivert}); continue; end + rcs = rcipsav{ivert}; + ngl = rcs.k; + nedge = rcs.nedge; + + rhohat = dens(rcs.starind); + [rho_cell,src_cell,wts_cell] = chnk.rcip.rhohatInterp(rhohat,rcs,rcs.nsub); + + for j = 1:nedge + rj = src_cell{j}.r; + nj = src_cell{j}.n; + dj = src_cell{j}.d; + d2j = src_cell{j}.d2; + wj = wts_cell{j}; + + npts_j = size(rj,2); + nch_j = npts_j/ngl; + assert(mod(npts_j,ngl)==0, ... + 'CHUNKERKERNEVAL: rcipsav fine-mesh point count not divisible by ngl'); + + pref_j = []; pref_j.k = ngl; + chnkr_j = chunker(pref_j); + chnkr_j = chnkr_j.addchunk(nch_j); + chnkr_j.r = reshape(rj, size(rj,1), ngl,nch_j); + chnkr_j.n = reshape(nj, size(nj,1), ngl,nch_j); + chnkr_j.d = reshape(dj, size(dj,1), ngl,nch_j); + chnkr_j.d2 = reshape(d2j, size(d2j,1), ngl,nch_j); + chnkr_j.wts = reshape(wj, ngl,nch_j); + + adj = zeros(2,nch_j); + adj(1,2:nch_j) = 1:nch_j-1; + adj(2,1:nch_j-1) = 2:nch_j; + chnkr_j.adj = adj; + + % rcs.ctr is [dim, nedge]: per-edge offset stored in Rcompchunk. + % Was using rcs.ctr(:,1) for all edges — correct for j=1, wrong + % for j>=2 (placed source panels at the wrong vertex coordinate, + % occasionally on top of targets, triggering downstream issues). + chnkr_j = chnkr_j + rcs.ctr(:,j); + + if use_wlchs + if fw_is_block + rho_j = rho_cell{j}(fw_active_col:fw_ndim_dens:end); + else + rho_j = rho_cell{j}(:); + end + if strcmpi(wlchs_family,'elasticity') + switch lower(wlchs_type) + case {'s','single'}, et = 's'; + case {'d','double'}, et = 'd'; + case {'strac','straction'}, et = 'strac'; + case 'dalt', et = 'dalt'; + end + if needs_tgt_normals + tg_e = targ_struct; + else + tg_e = targ_pts; + end + val_e = chnk.kernsplit.elast2d_panel_eval( ... + chnkr_j, et, wlchs_lam, wlchs_mu, rho_j, tg_e); + cor_j = wlchs_coef * val_e(:); + elseif strcmpi(wlchs_family,'laplace') + cor_j = wlchs_coef * chnk.kernsplit.lap2d_panel_eval( ... + chnkr_j, wlchs_type, rho_j, targ_pts); + elseif needs_tgt_normals + cor_j = wlchs_coef * chnk.kernsplit.helm2d_panel_eval( ... + chnkr_j, wlchs_type, wlchs_zk, rho_j, targ_struct); + else + cor_j = wlchs_coef * chnk.kernsplit.helm2d_panel_eval( ... + chnkr_j, wlchs_type, wlchs_zk, rho_j, targ_pts); + end + else + cor_j = chunkerkerneval(chnkr_j,kern,rho_cell{j}(:),targobj,opts_inner); + end + fints = fints + cor_j; end - base_probe = chnk.elast2d.kern(wlchs_lam, wlchs_mu, src_probe, ... - tgt_probe, et); - wlchs_coef = val_probe(1,1) / base_probe(1,1); end + return +end - % Non-rcipsav path: dispatch directly to kernsplit on the (possibly - % merged) chunker. Targets pass through as raw point coords or struct - % with .r and .n when normals are needed (sp; elasticity strac). +% Non-rcipsav path with forcewlchs: dispatch directly to kernsplit on the +% (possibly merged) chunker. Targets pass through as raw point coords. +if use_wlchs if length(chnkr) > 1 chnkr_use = merge(chnkr); else chnkr_use = chnkr; end targ_pts = local_target_points(targobj); + is_elast = strcmpi(wlchs_family,'elasticity'); is_strac_helm = any(strcmpi(wlchs_type, {'sp','sprime'})); is_strac_elast = is_elast && any(strcmpi(wlchs_type, {'strac','straction'})); needs_tgt_normals = is_strac_helm || is_strac_elast; @@ -247,34 +460,43 @@ 'requires target normals; pass targobj as a struct with .n']); end end + % Pass forcefmm into kernsplit's smooth-quadrature step. (Default off + % is safe — see helm2d_panel_eval; per-component rcipsav callers leave + % it unset, but a single global eval with many targets benefits from FMM.) hpe_opts = struct(); if isfield(opts,'forcefmm') && opts.forcefmm hpe_opts.forcefmm = true; end - if needs_tgt_normals - targ_for_eval = struct('r', targ_pts, 'n', targ_n); + if fw_is_block + dens_active = dens(fw_active_col:fw_ndim_dens:end); else - targ_for_eval = targ_pts; + dens_active = dens; end - switch lower(wlchs_family) - case 'helmholtz' - fints = wlchs_coef * chnk.kernsplit.helm2d_panel_eval(chnkr_use, ... - wlchs_type, wlchs_zk, dens, targ_for_eval, [], hpe_opts); - case 'laplace' - fints = wlchs_coef * chnk.kernsplit.lap2d_panel_eval(chnkr_use, ... - wlchs_type, dens, targ_for_eval, [], hpe_opts); - case 'elasticity' - switch lower(wlchs_type) - case {'s','single'}, et = 's'; - case {'d','double'}, et = 'd'; - case {'strac','straction'}, et = 'strac'; - case 'dalt', et = 'dalt'; - end - val = chnk.kernsplit.elast2d_panel_eval(chnkr_use, et, ... - wlchs_lam, wlchs_mu, dens, targ_for_eval, [], hpe_opts); - fints = wlchs_coef * val(:); - otherwise - error('CHUNKERKERNEVAL: forcewlchs family %s not supported', wlchs_family); + if is_elast + switch lower(wlchs_type) + case {'s','single'}, et = 's'; + case {'d','double'}, et = 'd'; + case {'strac','straction'}, et = 'strac'; + case 'dalt', et = 'dalt'; + end + if needs_tgt_normals + targ_for_eval = struct('r', targ_pts, 'n', targ_n); + else + targ_for_eval = targ_pts; + end + val = chnk.kernsplit.elast2d_panel_eval(chnkr_use, et, ... + wlchs_lam, wlchs_mu, dens_active, targ_for_eval, [], hpe_opts); + fints = wlchs_coef * val(:); + elseif strcmpi(wlchs_family,'laplace') + fints = wlchs_coef * chnk.kernsplit.lap2d_panel_eval(chnkr_use, ... + wlchs_type, dens_active, targ_pts, [], hpe_opts); + elseif needs_tgt_normals + targ_struct = struct('r', targ_pts, 'n', targ_n); + fints = wlchs_coef * chnk.kernsplit.helm2d_panel_eval(chnkr_use, ... + wlchs_type, wlchs_zk, dens_active, targ_struct, [], hpe_opts); + else + fints = wlchs_coef * chnk.kernsplit.helm2d_panel_eval(chnkr_use, ... + wlchs_type, wlchs_zk, dens_active, targ_pts, [], hpe_opts); end if icgrph && size(kern,1) > 1 ids = chunkgraphinregion(chnkobj, targ_pts); @@ -650,7 +872,7 @@ indji = (ji-1)*opdims(1); ind = (indji(:)).' + (1:opdims(1)).'; ind = ind(:); - fints(ind) = fints(ind) + fints1; + fints(ind) = fints(ind) + fints1; end @@ -658,7 +880,6 @@ end - function r = local_target_points(targobj) %LOCAL_TARGET_POINTS pull the dim x nt point array from any of the % supported targobj forms (chunker, chunkgraph, struct with .r, raw array). @@ -677,7 +898,7 @@ function n = local_target_normals(targobj) %LOCAL_TARGET_NORMALS pull a dim x nt target-normal array if the target % representation carries one; return [] otherwise. Used for forcewlchs -% with normal-derivative kernels (sp). +% with normal-derivative kernels (sp, dp). if isa(targobj,'chunker') || isa(targobj,'chunkgraph') n = targobj.n(:,:); elseif isstruct(targobj) && isfield(targobj,'n') && ~isempty(targobj.n) @@ -685,4 +906,4 @@ else n = []; end -end +end \ No newline at end of file diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 9741e829..fcf427c6 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -77,6 +77,23 @@ % adaptive quadrature for near touching panels on % different chunkers within rcip % opts.eps = (1e-14) tolerance for adaptive quadrature +% opts.srhs_eval = function handle ([]), enables the singular-RHS +% RCIP recursion (Helsing & Karlsson 2022, vector form). +% Signature: +% b_ib = srhs_eval(chnkrlocal_global, vert0, ... +% edge_indices, level, ndim) +% where chnkrlocal_global is a 1xnedge array of chunkers +% (in global coordinates) describing the type-b 6-panel +% local mesh at recursion level i, vert0 is the corner +% vertex coordinate, edge_indices(i) is the cgrph edge +% index that chnkrlocal_global(i) belongs to. Must return +% a vector of length 3*k*nedge*ndim with the RHS function +% evaluated at each GL node (ordered consistently with +% the level-i local system matrix). The resulting +% r_f^star vector is stored per corner in +% rcipsav{ivert}.r_f_star and is used by +% chnk.rcip.srhs_modify_rhs to correct the bare rhs +% vector at corner-star indices. % opts.forcewlchs = kernel-split self/adjacent-panel correction % for Helmholtz / Laplace S, D, S', D' kernels (and % Helmholtz combined 'c' = c1*D + c2*S, 'sc' = c1*S' @@ -568,7 +585,8 @@ if(icgrph && isrcip) - [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k); + k_rcip = chnkrs(1).k; + [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k_rcip); nch_all = horzcat(chnkobj.echnks.nch); npt_all = horzcat(chnkobj.echnks.npt); [~,nv] = size(chnkobj.verts); diff --git a/devtools/test/chunkerkerneval_rcipsav_Test.m b/devtools/test/chunkerkerneval_rcipsav_Test.m new file mode 100644 index 00000000..b36a26f2 --- /dev/null +++ b/devtools/test/chunkerkerneval_rcipsav_Test.m @@ -0,0 +1,103 @@ +% chunkerkerneval_rcipsav_Test.m +% +% Tests the opts.rcipsav code path in chunkerkerneval. Uses the same setup +% as devtools/test/rcipTest.m -- closed curve from two arcs of unit +% circle meeting at two corners, single Helmholtz D kernel BIE -- so the +% analytic field is computable from point sources outside. +% +% Compares: +% - naive chunkerkerneval (RCIP-unaware): expected to lose accuracy +% near each corner because dens(starind) is the RCIP-transformed +% density, only valid for far-field smooth quadrature. +% - chunkerkerneval with opts.rcipsav: expected to maintain ~1e-12 +% accuracy down to ~1e-5 from a corner (validated in step-2 smoke +% test using identical inner machinery). + +clearvars; close all; format long e; + +zk = 1.1; +funs = {@(t) circle(t), @(t) circle(t)}; +cparams1 = []; cparams1.ta = 0; cparams1.tb = pi/1.1; cparams1.nchmin = 1; +cparams2 = []; cparams2.ta = 0; cparams2.tb = pi/1.4; cparams2.nchmin = 8; +cg = chunkgraph([-1 1; 0 0], [1 2; 2 1], funs, {cparams1, cparams2}); +chnkrtotal = merge(cg.echnks); + +fkern = -2*kernel('helm','d',zk); +kerns = kernel('helm','s',zk); + +ns = 10; rng(8675309); +ts = 2*pi*rand(1,ns); +sources = 3.0*[cos(ts); sin(ts)]; +strengths = randn(ns,1); + +ubdry = kerns.eval(struct('r',sources), struct('r',chnkrtotal.r(:,:))) * strengths; + +opts = []; opts.rcip = true; opts.nsub_or_tol = 30; +opts.rcip_savedepth = inf; +fprintf('Building system matrix...\n'); +[sysmat, ~, rcipsav] = chunkermat(cg, fkern, opts); +sysmat = sysmat + eye(chnkrtotal.npt); + +fprintf('GMRES solve...\n'); +sol = gmres(sysmat, ubdry, [], 1e-13, 200); + +% Targets: a few far ones plus a sweep approaching corner at (-1,0) +ts = 2*pi*rand(1,3); +targets_far = 0.2*[cos(ts); sin(ts)]; +targets_far(:,1) = [0; 0.36]; +targets_far(:,2) = [-0.95; 0]; +targets_far(:,3) = [0.5; 0.2]; + +dvec = 10.^(-(1:0.5:6)); +targets_near = [-1; 0] + [cos(pi/4)*dvec; sin(pi/4)*dvec]; + +targets = [targets_far, targets_near]; +nt = size(targets,2); + +% Reference: free-space single-layer evaluation +utarg = kerns.eval(struct('r',sources), struct('r',targets)) * strengths; + +% (1) naive: chunkerkerneval on merged chunker, no rcipsav +opts_eval = []; +unum_naive = chunkerkerneval(chnkrtotal, fkern, sol, targets, opts_eval); + +% (2) rcipsav-aware: chunkerkerneval on chunkgraph, with rcipsav +opts_eval_rcip = []; opts_eval_rcip.rcipsav = rcipsav; +unum_rcip = chunkerkerneval(cg, fkern, sol, targets, opts_eval_rcip); + +err_naive = abs(unum_naive - utarg); +err_rcip = abs(unum_rcip - utarg); + +fprintf('\nTargets at fixed positions (far from corners):\n'); +fprintf(' %-22s %-14s %-14s\n', 'target', 'naive err', 'rcip-aware err'); +for i = 1:size(targets_far,2) + fprintf(' (%+.4f,%+.4f) %-14.3e %-14.3e\n', ... + targets_far(1,i), targets_far(2,i), err_naive(i), err_rcip(i)); +end +fprintf('\nTargets approaching corner at (-1,0):\n'); +fprintf(' %-12s %-14s %-14s\n', 'd from vert', 'naive err', 'rcip-aware err'); +for j = 1:size(targets_near,2) + i = size(targets_far,2) + j; + fprintf(' %-12.2e %-14.3e %-14.3e\n', dvec(j), err_naive(i), err_rcip(i)); +end + +% Far targets should be 1e-12 or better. +assert(max(err_rcip(1:size(targets_far,2))) < 1e-12, ... + 'rcip-aware should match naive at far targets'); +% Near-corner: 1e-10 down to d=3e-5; degrades gradually below that as +% targets enter the sub-deepest-panel regime where smooth quadrature +% on the fine mesh ceases to be enough. +n_far = size(targets_far,2); +err_near = err_rcip(n_far+1:end); +idx_above_1em5 = find(dvec >= 1e-5); +assert(max(err_near(idx_above_1em5)) < 1e-9, ... + 'rcip-aware should be ~1e-10 down to d=1e-5 from a corner'); +% The full sweep should stay better than 1e-9 for the deepest targets. +assert(max(err_near) < 1e-9, ... + 'rcip-aware should stay <1e-9 even at sub-panel distances'); +fprintf('\nTest passed.\n'); + +function [r,d,d2] = circle(t) + c = cos(t(:).'); s = sin(t(:).'); + r = [c; s]; d = [-s; c]; d2 = [-c; -s]; +end