Helmholtz open-arc Dirichlet solver via Bruno-Lintner formulation

chunkie/+chnk/+helm2d/fmm.m

@@ -68,6 +68,10 @@
         srcuse.charges = coefs(2)*sigma(:).';
         srcuse.dipstr  = coefs(1)*sigma(:).';
         srcuse.dipvec  = srcinfo.n(1:2,:);
+    case {'sc'}
+        % c_sp*S' + c_s*S — single FMM call with charges=sigma; combine the
+        % potential and the (target-normal-dotted) gradient at the end.
+        srcuse.charges = sigma(:).';
 end
 
 if ( isstruct(targinfo) )
@@ -79,7 +83,7 @@
 pg = 0;
 pgt = min(nargout, 2);
 switch lower(type)
-    case {'sprime', 'dprime', 'cprime','sp','dp','cp'}
+    case {'sprime', 'dprime', 'cprime','sp','dp','cp','sc'}
         pgt = max(pgt, 2);
 end
 U = hfmm2d(eps, zk, srcuse, pg, targuse, pgt);
@@ -96,6 +100,17 @@
             end
             varargout{1} = ( U.gradtarg(1,:).*targinfo.n(1,:) + ...
                              U.gradtarg(2,:).*targinfo.n(2,:) ).';
+        case {'sc'}
+            % c_sp*S' + c_s*S
+            if ( ~isfield(targinfo, 'n') )
+                error('CHUNKIE:helm2d:fmm:normals', ...
+                    'Targets require normal info when evaluating Helmholtz kernel ''sc''.');
+            end
+            coefs = varargin{1};
+            spval = ( U.gradtarg(1,:).*targinfo.n(1,:) + ...
+                      U.gradtarg(2,:).*targinfo.n(2,:) ).';
+            sval  = U.pottarg.';
+            varargout{1} = coefs(1)*spval + coefs(2)*sval;
         otherwise
             error('CHUNKIE:helm2d:fmm:pot', ...
                 'Potentials not supported for Helmholtz kernel ''%s''.', type);

chunkie/+chnk/+helm2d/kern.m

@@ -202,6 +202,17 @@
   submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
   submat = coef(1)*submatd + coef(2)*submats;
 
+% S' + S combined: coef(1)*S' + coef(2)*S
+case {'sc', 'spcombined'}
+  targnorm = targinfo.n(:,:);
+  coef = ones(2,1);
+  if(nargin == 5); coef = varargin{1}; end
+  [submats,grad] = chnk.helm2d.green(zk,src,targ);
+  nx = repmat((targnorm(1,:)).',1,ns);
+  ny = repmat((targnorm(2,:)).',1,ns);
+  submatsp = grad(:,:,1).*nx + grad(:,:,2).*ny;
+  submat = coef(1)*submatsp + coef(2)*submats;
+
 % normal derivative of combined field
 case {'cp', 'cprime'}
   coef = ones(2,1);

chunkie/@kernel/helm2d.m

@@ -118,6 +118,20 @@
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs);
         obj.sing = 'hs';
 
+    case {'sc', 'spcombined'}
+        % "S' + S combined": c_sp * S' + c_s * S, with coefs = [c_sp, c_s].
+        % Used for impedance BIE on a Neumann/sound-hard boundary (e.g.
+        % photonics K_imp = 2*S' + 2*eta*S on cap edges).
+        if ( nargin < 3 )
+            warning('Missing combined coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'sc';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'sc', coefs);
+        obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'sc', sigma, coefs);
+        obj.sing = 'log';
+
     case {'all', 'trans_sys', 'tsys'}
         if ( nargin < 3 )
             warning('Missing transmission coefficients. Defaulting to [1,1;1,1].');

chunkie/chunkermat.m

@@ -77,6 +77,14 @@
 %                    adaptive quadrature for near touching panels on
 %                    different chunkers within rcip
 %           opts.eps = (1e-14) tolerance for adaptive quadrature
+%           opts.open_arc_eye = boolean (false), if true and opdims=[2,2],
+%                    overwrite the (1,1) ndim-sub-block of each per-chunker
+%                    diagonal block with -I before RCIP correction. This is
+%                    used by the Bruno-Lintner open-arc Helmholtz Dirichlet
+%                    formulation, where K = [Z, c*S; c*T, Z] and the caller
+%                    adds eye(nsys) outside; the override makes the (1,1)
+%                    diagonal vanish and the (2,2) diagonal equal I, giving
+%                    the system [0, c*S; c*T, I].
 %           opts.srhs_eval = function handle ([]), enables the singular-RHS
 %                    RCIP recursion (Helsing & Karlsson 2022, vector form).
 %                    Signature:
@@ -181,6 +189,7 @@
 adaptive_correction = false;
 rcip_adaptive_correction = false;
 sing = 'log';
+open_arc_eye = false;
 
 % get opts from struct if available
 
@@ -232,6 +241,16 @@
 if (isfield(opts,'rcip_adaptive_correction'))
     rcip_adaptive_correction = opts.rcip_adaptive_correction;
 end
+if isfield(opts,'open_arc_eye')
+    open_arc_eye = opts.open_arc_eye;
+end
+% Bruno-Lintner block-kernel sub-block size for the (1,1) override:
+%   helmos / Helmholtz block kernel:    opdims=[2,2], sub_ndim=1 (default)
+%   open-arc Stokes mobility:           opdims=[4,4], sub_ndim=2
+open_arc_eye_subdim = 1;
+if isfield(opts,'open_arc_eye_subdim') && ~isempty(opts.open_arc_eye_subdim)
+    open_arc_eye_subdim = opts.open_arc_eye_subdim;
+end
 
 % forcewlchs: kernel-split self/adjacent-panel correction for Laplace
 % S/D/Sp/Dp kernels. See doc above for accepted forms.
@@ -540,6 +559,32 @@
             end
         end
     end
+    if open_arc_eye
+        % Bruno-Lintner open-arc override: overwrite the top-left
+        % d-by-d sub-block of each diagonal point block with -I_d, so
+        % that after the caller's +eye(nsys) the block system reads
+        % [0, c*A; c*B, I_d]. Default d=1 (helmos opdims=[2,2], scalar
+        % sub-blocks); d=2 for the open-arc Stokes mobility 4x4 block
+        % kernel where each sub-block is itself a 2x2 Stokes operator.
+        d = open_arc_eye_subdim;
+        if opdims(1) ~= opdims(2) || opdims(1) ~= 2*d
+            error('chunkermat: open_arc_eye requires opdims = [2*subdim, 2*subdim] (got [%d,%d], subdim=%d)', ...
+                  opdims(1), opdims(2), d);
+        end
+        npts2 = size(sysmat_tmp, 1);
+        Npts  = npts2 / opdims(1);
+        block_off = (0:Npts-1) * opdims(1);
+        for jj = 1:d
+            for ii = 1:d
+                lin = sub2ind(size(sysmat_tmp), block_off + ii, block_off + jj);
+                if ii == jj
+                    sysmat_tmp(lin) = -1;
+                else
+                    sysmat_tmp(lin) = 0;
+                end
+            end
+        end
+    end
 
     if adaptive_correction
         flag = flagnear(chnkr,chnkr.r(:,:));

chunkie/demo/geo_corners.m

@@ -0,0 +1,26 @@
+function cg = geo_corners(ncorner, zk, ppw)
+%GEO_CORNERS chain of ncorner+1 line segments connecting alternating
+% heights, mirroring testhelmos's icase=2 with curvetype=1 (default).
+% Vertices: ncorner+2 (2 endpoints + ncorner interior corners).
+% Edges: ncorner+1 (line segments between consecutive vertices).
+
+if nargin < 1 || isempty(ncorner), ncorner = 1; end
+if nargin < 3 || isempty(ppw), ppw = 10; end
+pref = []; pref.k = 16;
+maxchunklen = pref.k/ppw/abs(zk)*2*pi;
+
+% Same vertex layout as testhelmos's curvetype=1 in ossetup_corners
+nsing  = ncorner + 2;
+ncurve = ncorner + 1;
+d = 1;
+z = zeros(nsing,1);
+z(1:2:nsing) = (-ncorner/2-1)*d - 0.5i + (1:2:nsing)*d;
+z(2:2:nsing) = (-ncorner/2-1)*d + (2:2:nsing)*d + 0.5i;
+
+verts = [real(z).'; imag(z).'];
+edge2verts = [1:ncurve; 2:ncurve+1];
+
+cparams = struct('ta',0,'tb',1,'maxchunklen',maxchunklen,'nover',2);
+cg = chunkgraph(verts, edge2verts, [], cparams, pref);
+cg = balance(cg);
+end

chunkie/demo/geo_linesegment.m

@@ -0,0 +1,20 @@
+function cg = geo_linesegment(zk, ppw, klam)
+%GEO_LINESEGMENT chunkgraph of a line segment from (-1,-0.2) to (1,-0.2)
+% (matching testhelmos's icase=0 placement).
+%
+% Inputs:
+%   zk   - wavenumber
+%   ppw  - points per wavelength (default 10)
+%   klam - desired chunk length scaling factor; if empty, default = 1.
+
+if nargin < 2 || isempty(ppw), ppw = 10; end
+if nargin < 3 || isempty(klam), klam = 1; end
+
+pref = []; pref.k = 16;
+maxchunklen = klam * pref.k/ppw/abs(zk)*2*pi;
+verts = [-1, 1; -0.2, -0.2];
+edge2verts = [1; 2];
+cparams = struct('ta',0,'tb',1,'maxchunklen',maxchunklen,'nover',2);
+cg = chunkgraph(verts, edge2verts, [], cparams, pref);
+cg = balance(cg);
+end

chunkie/demo/geo_spiral.m

@@ -0,0 +1,31 @@
+function cg = geo_spiral(zk, ppw)
+%GEO_SPIRAL chunkgraph for the Bruno-Lintner 2012 logarithmic spiral
+% gamma(t) = exp(-1-5i + (2+10i) t), t in [0,1].
+% Endpoints: gamma(0) = exp(-1-5i), gamma(1) = exp(1+5i).
+
+if nargin < 2 || isempty(ppw), ppw = 10; end
+pref = []; pref.k = 16;
+maxchunklen = pref.k/ppw/abs(zk)*2*pi;
+
+a = exp(-1 - 5i);
+b = 2 + 10i;
+fchnk = @(t) spiralcurve(t,a,b);
+
+z0 = a;            % gamma(0)
+z1 = a*exp(b);     % gamma(1)
+verts = [real(z0), real(z1); imag(z0), imag(z1)];
+edge2verts = [1;2];
+
+cparams = struct('ta',0,'tb',1,'maxchunklen',maxchunklen,'nover',2);
+cg = chunkgraph(verts, edge2verts, {fchnk}, cparams, pref);
+cg = balance(cg);
+end
+
+function [r,d,d2] = spiralcurve(t,a,b)
+% z(t)=a*exp(b*t),  z'=a*b*exp(b*t),  z''=a*b^2*exp(b*t)
+ebt = exp(b*t(:).');
+z   = a*ebt;       zp  = a*b*ebt;     zpp = a*b*b*ebt;
+r   = [real(z);  imag(z)];
+d   = [real(zp); imag(zp)];
+d2  = [real(zpp);imag(zpp)];
+end

chunkie/demo/helmos.m

@@ -0,0 +1,145 @@
+function out = helmos(geo, bctype, zk, opts)
+%HELMOS Helmholtz open-arc solver via Bruno-Lintner integral representation.
+%
+% Geometry-agnostic driver. Caller provides the chunkgraph (with all
+% endpoints/corners/junctions as vertices) and a target list. Uses chunkie's
+% RCIP machinery via chunkermat/chunkerkerneval with the open-arc identity
+% override (opts.open_arc_eye=true) and rcipsav-aware field evaluation
+% (opts.rcipsav).
+%
+% Status:
+%   Dirichlet ('d')  -- delivers 10+ digits at typical grid distances
+%                       (validated by helmos_smoke / helmos_subpanel_test).
+%   Neumann   ('n')  -- formulation in place but limited to ~7 digits at
+%                       present because chunkie's hypersingular 'dp'
+%                       kernel does not yet have a kernel-split smooth
+%                       quadrature like 'd'/'s' do; |sigma| grows with
+%                       mesh refinement. Needs follow-up (port testhelmos
+%                       Toper correction or add splitinfo to 'dp').
+%
+% Inputs:
+%   geo - struct with fields:
+%       .cg       - chunkgraph
+%       .targets  - 2 x nt target locations
+%       .ubdry_fn - function handle ubdry_fn(srcinfo, bctype) -> column,
+%                   the boundary data for the BIE on the chunkgraph nodes
+%   bctype - 'd' (Dirichlet) or 'n' (Neumann)
+%   zk - wavenumber
+%   opts - optional struct
+%       .nsub        - RCIP refinement (default 30)
+%       .tol         - GMRES tol (default 1e-12)
+%       .gmres_maxit - default 200
+%       .verbose     - print progress (default true)
+%
+% Output: struct with fields .sol, .iter, .ufield, .timings, .rcipsav,
+%         .Keval. The Keval is returned so the caller can postprocess.
+
+if nargin < 4 || isempty(opts), opts = []; end
+nsub  = getfielddef(opts,'nsub',30);
+tol   = getfielddef(opts,'tol',1e-12);
+maxit = getfielddef(opts,'gmres_maxit',200);
+vb    = getfielddef(opts,'verbose',true);
+rcip_adap = getfielddef(opts,'rcip_adaptive_correction',false);
+
+cg = geo.cg;
+S  = kernel('helm','s',  zk);
+T  = kernel('helm','dp', zk);
+D  = kernel('helm','d',  zk);
+Z  = kernel.zeros();
+c  = 2.0;
+
+% Bruno-Lintner formulation:
+%   Dirichlet:   u = c S (-c T) sigma  -- field eval kernel = c S applied to tau
+%   Neumann:     u = c D (-c S) sigma  -- field eval kernel = c D applied to tau
+% Block matrix K + open_arc_eye gives [0 cA; cB I] form, where (A,B) =
+% (S, T) for Dirichlet and (T, S) for Neumann.  After GMRES, the second
+% density component is the auxiliary tau used by the field eval.
+if bctype == 'd'
+    K            = kernel([Z, c*S; c*T, Z]);
+    Keval        = kernel([Z, c*S]);
+    eval_type    = 's';        % chnk.kernsplit type for close eval
+    fwlchs_layout = {{'zero',0}, {'s', c}; {'dp', c}, {'zero',0}};
+elseif bctype == 'n'
+    K            = kernel([Z, c*T; c*S, Z]);
+    Keval        = kernel([Z, c*D]);   % NOTE: D, not T -- field uses double layer
+    eval_type    = 'd';
+    fwlchs_layout = {{'zero',0}, {'dp', c}; {'s', c}, {'zero',0}};
+else
+    error('HELMOS: bctype must be ''d'' or ''n''');
+end
+
+src = struct('r', cg.r(:,:), 'n', cg.n(:,:));
+ubdry = geo.ubdry_fn(src, bctype, zk);
+
+npts  = cg.npt;
+nsys  = K.opdims(1)*npts;
+rhs   = zeros(nsys,1);
+rhs(1:K.opdims(1):end) = ubdry;
+
+% forcewlchs covers all singular self+adj entries for the helmos block
+% kernel, so GGQ self/adj corrections in chunkermat would be discarded.
+% Use 'native' (smooth GL) and let wlchs+open_arc_eye fill in the
+% diagonal/near-singular pieces.
+matopts = struct('rcip',true,'nsub_or_tol',nsub, ...
+                 'open_arc_eye',true,'rcip_savedepth',inf, ...
+                 'rcip_adaptive_correction',rcip_adap, ...
+                 'quad','native', ...
+                 'forcewlchs', struct('zk',zk,'layout',{fwlchs_layout}));
+matfree = getfielddef(opts,'matrix_free',false);
+if matfree
+    matopts.corrections = true;   % return only sparse cormat
+end
+if vb, fprintf('[helmos] building %s system (npts=%d, nsys=%d)...\n', ...
+        ternary(matfree,'sparse','dense'), npts, nsys); end
+t = tic; [A_or_cor, ~, rcipsav] = chunkermat(cg, K, matopts);
+t_mat = toc(t);
+if vb, fprintf('[helmos] matrix build: %.2f s\n', t_mat); end
+
+if matfree
+    chnkrmerge = merge(cg.echnks);
+    src.r  = cg.r(:,:);  src.d  = cg.d(:,:);
+    src.d2 = cg.d2(:,:); src.n  = cg.n(:,:);
+    opdims_mat = [K.opdims(1); K.opdims(2)];
+    fmm_smooth_opts = struct('forcefmm', getfielddef(opts,'forcefmm',true));
+    apply = @(x) x + A_or_cor*x + ...
+        chnk.chunkerkerneval_smooth(chnkrmerge, K, opdims_mat, x, src, [], fmm_smooth_opts);
+else
+    A = A_or_cor + eye(nsys);
+    apply = @(x) A*x;
+end
+
+if vb, fprintf('[helmos] gmres (%s)...\n', ternary(matfree,'matrix-free','dense')); end
+t = tic; [sol, flag, relres, iter] = gmres(apply, rhs, [], tol, maxit);
+t_solve = toc(t);
+if vb
+    if iscell(iter), it = iter{end}; else, it = iter(end); end
+    fprintf('[helmos] gmres: flag=%d relres=%.2e iter=%d in %.2f s\n', flag, relres, it, t_solve);
+end
+
+if vb, fprintf('[helmos] field evaluation at %d targets...\n', size(geo.targets,2)); end
+t = tic;
+% Single chunkerkerneval call: rcipsav + forcewlchs (block-kernel struct
+% form) handles both the smooth far-field and the per-vertex per-edge
+% kernel-split corner correction.  The layout selects which density
+% component is the active one (column 2 = tau) and which kernsplit
+% kernel type to apply (eval_type, with prefactor c).
+eval_layout = {{'zero',0}, {eval_type, c}};
+eval_opts = struct('rcipsav',  {rcipsav}, ...
+                   'forcewlchs', struct('zk',zk,'layout',{eval_layout}), ...
+                   'forcefmm', getfielddef(opts,'forcefmm',true));
+ufield = chunkerkerneval(cg, Keval, sol, geo.targets, eval_opts);
+t_eval = toc(t);
+if vb, fprintf('[helmos] eval: %.2f s\n', t_eval); end
+
+out = struct('sol',sol, 'iter',iter, 'ufield',ufield, ...
+             'rcipsav',{rcipsav}, 'Keval',Keval, ...
+             'timings',struct('matrix',t_mat,'solve',t_solve,'eval',t_eval));
+end
+
+function v = getfielddef(s,f,d)
+    if isfield(s,f), v = s.(f); else, v = d; end
+end
+
+function v = ternary(cond, a, b)
+    if cond, v = a; else, v = b; end
+end