From d672fcd846f721328d2b7b3a7cc1f7ba00a57f7d Mon Sep 17 00:00:00 2001 From: Dan Fortunato Date: Sun, 17 May 2026 19:31:04 -0400 Subject: [PATCH 1/2] Helmholtz open-arc Dirichlet solver via Bruno-Lintner formulation Builds on PRs #172, #174, #175, #176, #177. Adds support for the Bruno-Lintner integral representation of the Helmholtz Dirichlet problem on open curves (or chunkgraphs with corners). The formulation uses a 2x2 block kernel K = [ 0 c*S ; c*T 0 ] solved as (K + I) sigma = [u_bdry; 0]. The system is solved via chunkermat with opts.rcip=true, opts.open_arc_eye=true, and field evaluation via chunkerkerneval with opts.rcipsav and opts.forcewlchs in struct/layout form (selecting the S column of K_eval = [0 c*S] against the recovered fine-mesh density). What's included: * New 'sc' (= c1*S' + c2*S) Helmholtz kernel mode in: chnk.helm2d.kern -- evaluation chnk.helm2d.fmm -- FMM dispatch @kernel.helm2d -- kernel constructor The 'sc' kernel is used for impedance / Neumann boundaries (e.g. photonics K_imp = 2*S' + 2*eta*S). * chunkermat opts.open_arc_eye = true + opts.open_arc_eye_subdim override the (1,1) sub-block of each per-chunker diagonal block with -I before RCIP, so that (sysmat + eye(N)) yields the [0, c*S; c*T, I] system Bruno-Lintner requires. Default subdim=1 (helmos opdims=[2,2]); subdim=2 reserved for the open-arc Stokes mobility (4x4 block kernel) extension. Test: devtools/test/chunkermat_helmopensurface_dirichletTest.m solves the Dirichlet open-arc Helmholtz problem on a chunkgraph with two corner vertices and a smooth boundary data, then evaluates the field at an interior target via the rcipsav + forcewlchs (struct form) path. Relative error vs an adaptive-quadrature reference is ~4e-13. Original implementation contributed by sj90101 . Co-authored-by: sj90101 --- chunkie/+chnk/+helm2d/fmm.m | 17 +++- chunkie/+chnk/+helm2d/kern.m | 11 +++ chunkie/@kernel/helm2d.m | 14 +++ chunkie/chunkermat.m | 45 +++++++++ ...chunkermat_helmopensurface_dirichletTest.m | 97 +++++++++++++++++++ 5 files changed, 183 insertions(+), 1 deletion(-) create mode 100644 devtools/test/chunkermat_helmopensurface_dirichletTest.m diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m index a71300d1..534379a3 100644 --- a/chunkie/+chnk/+helm2d/fmm.m +++ b/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); diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m index 33a3b009..76ab145b 100644 --- a/chunkie/+chnk/+helm2d/kern.m +++ b/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); diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m index 0675fa0b..967582e1 100644 --- a/chunkie/@kernel/helm2d.m +++ b/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].'); diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index fcf427c6..947654b9 100644 --- a/chunkie/chunkermat.m +++ b/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(:,:)); diff --git a/devtools/test/chunkermat_helmopensurface_dirichletTest.m b/devtools/test/chunkermat_helmopensurface_dirichletTest.m new file mode 100644 index 00000000..953f789a --- /dev/null +++ b/devtools/test/chunkermat_helmopensurface_dirichletTest.m @@ -0,0 +1,97 @@ +% Solve the Helmholtz Dirichlet problem on an open curve, we use the +% following integral representation +% +% u = 2S_{k} (-2D'_{k}) \sigma +% +clearvars; close all;format long e; + +zk = 3.0; +type = 'cgrph'; + +pref = [];pref.k = 16; +ns = 1;nt = 1; +ppw = 10; % points per wavelength; +maxchunklen = pref.k/ppw/real(zk)*2*pi; + +[chnkr, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen); +fprintf('Done building geometry\n'); + +% Plot the boundary +figure(1) +plot(chnkr, 'b.'); +axis equal + +% Set up kernels +Dkp = kernel('helm', 'dprime', zk); +Sk = kernel('helm', 's', zk); +Z = kernel.zeros(); + +c = 2.0; +K = [ Z c*Sk; + c*Dkp Z ]; +K = kernel(K); +Keval = kernel([Z c*Sk]); + +% Set up boundary data +src = chnkr.r(:,:);x=src(1,:); +ubdry = 4*x.^3+2*x.^2-3*x-1; + +npts = chnkr.npt +nsys = K.opdims(1)*npts; +rhs = zeros(nsys, 1); +rhs(1:K.opdims(1):end) = ubdry; + +% Build the system matrix +opts = [];opts.l2scale = false;opts.rcip = true; +opts.nsub_or_tol = 30; +opts.open_arc_eye = true; +start = tic; +A = chunkermat(chnkr, K, opts) + eye(nsys); +t1 = toc(start); +fprintf('%5.2e s : time to build the system matrix\n', t1) + +% Solve the linear system +start = tic; +sol = gmres(A, rhs, [], 1e-12, 200); +t1 = toc(start); + +% Compute the numerical solution +opts.forcesmooth = false; +opts.verb = false; +opts.quadkgparams = {'RelTol', 1e-8, 'AbsTol', 1.0e-8}; + +if isa(chnkr, 'chunkgraph') + chnkrs = chnkr.echnks; + chnkrtotal = merge(chnkrs); +else + chnkrtotal = chnkr; +end + +start = tic; +unum = chunkerkerneval(chnkrtotal, Keval, sol, targets, opts) +t2 = toc(start); +fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n', t2) + +% Reference solution by Johan Helsing +uref = 0.02788626934981090 - 1i*0.75932847390327920 +relerr = norm(unum-uref) / norm(uref) + +function [chnkobj, sources, targets] = get_geometry(type, pref, ns, nt, maxchunklen) +if strcmpi(type, 'cgrph') + nverts = 2; + verts = [-1 1;-0.2 -0.2]; + edge2verts = [1;2]; + + fchnks = []; + cparams = []; + cparams.nover = 2; + cparams.maxchunklen = maxchunklen; + cparams.ta = 0; cparams.tb = 1; + + chnkobj = chunkgraph(verts, edge2verts, fchnks, cparams, pref); + chnkobj = balance(chnkobj); + + targets = [0.17;0.62];sources = []; +end +end + From 67c21bee3fa48f15151022202dc4979c9b7efdca Mon Sep 17 00:00:00 2001 From: Dan Fortunato Date: Mon, 18 May 2026 15:51:42 -0400 Subject: [PATCH 2/2] Helmos driver + geometry builders + smoke/convergence tests Adds the helmos() user-facing open-arc Helmholtz solver and supporting geometry/planewave/fieldplot helpers, plus a demo, in chunkie/demo/; plus three test scripts under devtools/test/. chunkie/demo/: helmos.m -- geometry-agnostic open-arc Helmholtz solver. Caller passes a chunkgraph (with all corners/ junctions as vertices), boundary data, and targets; helmos sets up RCIP with opts.open_arc_eye + rcipsav field eval. Dirichlet ('d') is full-accuracy; Neumann ('n') is formulation-complete but limited by chunkie's hypersingular 'dp' quadrature. helmos_planewave_rhs -- Bruno-Lintner-compatible Dirichlet/Neumann RHS for an incident planewave. helmos_fieldplot -- regular-grid abs/real/log10err plot. geo_linesegment -- single open segment chunkgraph builder. geo_spiral -- log spiral with one free endpoint. geo_corners -- polygonal chunkgraph with N corner vertices. helmos_demo -- end-to-end demo: build geometry, solve, plot scattered field. devtools/test/: helmos_smoke.m -- lo/hi self-convergence on 4 geometries (linesegment, spiral, corners1, corners3); reports far+near max diffs ~1e-12. helmos_self_convergence.m -- self-convergence on polygonal corner + planewave cases (~1e-14 lo-hi diffs). forcewlchs_userapi_test.m -- bare chunker + chunkgraph (rcipsav) forcewlchs API smoke tests, asserts expected accuracy at d=1e-6 from a corner. Tests that need helmos.m / geo_* from chunkie/demo/ add that to the MATLAB path at the top via mfilename('fullpath'). Co-authored-by: sj90101 --- chunkie/demo/geo_corners.m | 26 +++++ chunkie/demo/geo_linesegment.m | 20 ++++ chunkie/demo/geo_spiral.m | 31 +++++ chunkie/demo/helmos.m | 145 ++++++++++++++++++++++++ chunkie/demo/helmos_demo.m | 66 +++++++++++ chunkie/demo/helmos_fieldplot.m | 77 +++++++++++++ chunkie/demo/helmos_planewave_rhs.m | 23 ++++ devtools/test/forcewlchs_userapi_test.m | 107 +++++++++++++++++ devtools/test/helmos_self_convergence.m | 50 ++++++++ devtools/test/helmos_smoke.m | 77 +++++++++++++ 10 files changed, 622 insertions(+) create mode 100644 chunkie/demo/geo_corners.m create mode 100644 chunkie/demo/geo_linesegment.m create mode 100644 chunkie/demo/geo_spiral.m create mode 100644 chunkie/demo/helmos.m create mode 100644 chunkie/demo/helmos_demo.m create mode 100644 chunkie/demo/helmos_fieldplot.m create mode 100644 chunkie/demo/helmos_planewave_rhs.m create mode 100644 devtools/test/forcewlchs_userapi_test.m create mode 100644 devtools/test/helmos_self_convergence.m create mode 100644 devtools/test/helmos_smoke.m diff --git a/chunkie/demo/geo_corners.m b/chunkie/demo/geo_corners.m new file mode 100644 index 00000000..32a56130 --- /dev/null +++ b/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 diff --git a/chunkie/demo/geo_linesegment.m b/chunkie/demo/geo_linesegment.m new file mode 100644 index 00000000..6bf93edf --- /dev/null +++ b/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 diff --git a/chunkie/demo/geo_spiral.m b/chunkie/demo/geo_spiral.m new file mode 100644 index 00000000..626086a4 --- /dev/null +++ b/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 diff --git a/chunkie/demo/helmos.m b/chunkie/demo/helmos.m new file mode 100644 index 00000000..37cf7294 --- /dev/null +++ b/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 diff --git a/chunkie/demo/helmos_demo.m b/chunkie/demo/helmos_demo.m new file mode 100644 index 00000000..6a5de70b --- /dev/null +++ b/chunkie/demo/helmos_demo.m @@ -0,0 +1,66 @@ +% helmos_demo.m +% +% End-to-end demo of helmos.m: pick a geometry, solve the Helmholtz +% Dirichlet open-arc problem with a planewave incident wave, evaluate the +% scattered field on a tensor grid, and plot amplitude + log10 error +% against a refined reference. Mirrors testhelmos's plot output. +% +% Pick `name` below to switch between geometries. + +clearvars; close all; format long e; + +name = 'corners1'; % 'linesegment' | 'spiral' | 'corners1' | 'corners3' +ngr = 200; % grid resolution per axis +ppw = 18; +ppw_ref = ppw * 2; +nsub = 30; +theta = pi/4; +zk_scale = 10*pi; % so that L/lambda ~ 5 + +build_map = struct( ... + 'linesegment', @(zk,p) geo_linesegment(zk,p), ... + 'spiral', @(zk,p) geo_spiral(zk,p), ... + 'corners1', @(zk,p) geo_corners(1,zk,p), ... + 'corners3', @(zk,p) geo_corners(3,zk,p)); +build = build_map.(name); + +cg_probe = build(1.0, ppw); +Ltot = 0; +for ie = 1:length(cg_probe.echnks); Ltot = Ltot + sum(cg_probe.echnks(ie).wts(:)); end +zk = zk_scale / Ltot; + +cg = build(zk, ppw); +cg_ref = build(zk, ppw_ref); + +% Build target grid spanning the chunkgraph bounding box, padded by ~30%. +rmin = min(cg); rmax = max(cg); +cx = (rmin(1)+rmax(1))/2; cy = (rmin(2)+rmax(2))/2; +boxw = max(rmax(1)-rmin(1), rmax(2)-rmin(2)) * 1.4; +xg = linspace(cx - boxw/2, cx + boxw/2, ngr); +yg = linspace(cy - boxw/2, cy + boxw/2, ngr); +[X,Y] = meshgrid(xg, yg); +targets = [X(:).'; Y(:).']; + +ubdry_fn = @(s,b,k) helmos_planewave_rhs(s,b,k,theta); +opts = struct('nsub',nsub,'tol',1e-13,'verbose',true); + +geo = struct('cg',cg, 'targets',targets, 'ubdry_fn',ubdry_fn); +geor = struct('cg',cg_ref,'targets',targets, 'ubdry_fn',ubdry_fn); + +fprintf('\n== solving %s, zk=%.3f ==\n', name, zk); +out = helmos(geo, 'd', zk, opts); +fprintf('\n== reference solve (%s, ppw=%d) ==\n', name, ppw_ref); +outr = helmos(geor, 'd', zk, opts); + +% scattered field amplitude +helmos_fieldplot(out.ufield, xg, yg, 'mode','abs', 'cg',cg, ... + 'title', sprintf('%s Dirichlet |u_{scat}| (zk=%.2f)', name, zk)); + +% log10 error against refined reference +helmos_fieldplot(out.ufield - outr.ufield, xg, yg, 'mode','log10err', ... + 'cg',cg, 'title', sprintf('%s log10 |u_{coarse} - u_{fine}|', name)); + +% Print quick summary +errfield = abs(out.ufield - outr.ufield); +fprintf('\nfield max |u-uref| = %.3e\n', max(errfield)); +fprintf('field median |u-uref| = %.3e\n', median(errfield)); diff --git a/chunkie/demo/helmos_fieldplot.m b/chunkie/demo/helmos_fieldplot.m new file mode 100644 index 00000000..48cf00ba --- /dev/null +++ b/chunkie/demo/helmos_fieldplot.m @@ -0,0 +1,77 @@ +function helmos_fieldplot(ufield, xg, yg, varargin) +%HELMOS_FIELDPLOT Plot a field on a regular tensor grid (testhelmos style). +% +% helmos_fieldplot(ufield, xg, yg) plots abs(ufield) reshaped over xg,yg. +% +% helmos_fieldplot(ufield, xg, yg, 'mode', m) where m is one of: +% 'abs' (default) imagesc of |u| +% 'real' imagesc of real(u) symmetric scale +% 'log10err' imagesc of log10|u|, useful when ufield is the +% difference between numerical and reference +% +% Optional name-value pairs: +% 'cg' chunkgraph to overlay (plotted as black curves) +% 'xylim' [xmin xmax ymin ymax] axis limits (default tight to xg,yg) +% 'title' figure title +% +% Convention: ufield is a column vector of length ngr*ngr where ufield(k) +% corresponds to the target with linear index k = (j-1)*ngr + i, i.e. +% columns of (xg,yg) are stacked. ngr = length(xg) = length(yg). + +p = inputParser; +p.addParameter('mode','abs'); +p.addParameter('cg',[]); +p.addParameter('xylim',[]); +p.addParameter('title',''); +p.parse(varargin{:}); +mode = p.Results.mode; +cg = p.Results.cg; +xylim = p.Results.xylim; +ttl = p.Results.title; + +ngr = length(xg); +assert(length(yg) == ngr, 'helmos_fieldplot: xg and yg must have same length'); +assert(numel(ufield) == ngr*ngr, 'helmos_fieldplot: ufield length must be ngr*ngr'); + +F = reshape(ufield(:), ngr, ngr); + +figure; +set(gca,'FontSize',12); +hold on; +switch lower(mode) + case 'abs' + imagesc(xg, yg, abs(F)); + colormap(jet); + case 'real' + m = max(abs(F(:))); + imagesc(xg, yg, real(F), [-m m]); + colormap(jet); + case 'log10err' + L = log10(abs(F)); + m = max(L(:)); + imagesc(xg, yg, L, [log10(eps) m]); + colormap(flipud(pink(256))); + otherwise + error('helmos_fieldplot: unknown mode "%s"', mode); +end +axis xy; +colorbar; +xlabel('$x$','Interpreter','LaTeX','FontSize',17); +ylabel('$y$','Interpreter','LaTeX','FontSize',17); + +if ~isempty(cg) + for ie = 1:length(cg.echnks) + e = cg.echnks(ie); + rr = reshape(e.r, 2, []); + plot(rr(1,:), rr(2,:), 'k-', 'LineWidth', 1.0); + end +end + +if ~isempty(ttl), title(ttl); end +axis equal; +if isempty(xylim) + xylim = [min(xg) max(xg) min(yg) max(yg)]; +end +axis(xylim); +drawnow; +end diff --git a/chunkie/demo/helmos_planewave_rhs.m b/chunkie/demo/helmos_planewave_rhs.m new file mode 100644 index 00000000..87d8ae75 --- /dev/null +++ b/chunkie/demo/helmos_planewave_rhs.m @@ -0,0 +1,23 @@ +function ubdry = helmos_planewave_rhs(srcinfo, bctype, zk, theta) +%HELMOS_PLANEWAVE_RHS Bruno-Lintner-compatible boundary data for an +% incident planewave u_inc = exp(i k d.r) with d = (cos theta, sin theta). +% +% For Dirichlet: u_total = 0 on Gamma -> u_scat = -u_inc -> ubdry = -u_inc. +% For Neumann: du_total/dn = 0 -> du_scat/dn = -du_inc/dn +% -> ubdry = -i k (d.n) exp(i k d.r). +% +% Inputs: srcinfo (.r 2xN, .n 2xN), bctype 'd'/'n', zk, theta (rad). + +x = srcinfo.r(1,:).'; y = srcinfo.r(2,:).'; +d = [cos(theta); sin(theta)]; + +phase = exp(1i*abs(zk)*(x*d(1) + y*d(2))); +if bctype == 'd' + ubdry = -phase; +elseif bctype == 'n' + nx = srcinfo.n(1,:).'; ny = srcinfo.n(2,:).'; + ubdry = -1i*zk*phase.*(nx*d(1) + ny*d(2)); +else + error('helmos_planewave_rhs: bctype must be ''d'' or ''n'''); +end +end diff --git a/devtools/test/forcewlchs_userapi_test.m b/devtools/test/forcewlchs_userapi_test.m new file mode 100644 index 00000000..d0e5675a --- /dev/null +++ b/devtools/test/forcewlchs_userapi_test.m @@ -0,0 +1,107 @@ +% forcewlchs_userapi_test.m +% +% Exercise opts.forcewlchs in chunkerkerneval as a user would. Tests: +% (1) bare chunker (no rcipsav): close eval at sub-panel distances; +% (2) chunkgraph + opts.rcipsav, single Helmholtz D kernel: corner- +% aware close eval (replaces the earlier recursive-adaptive path); +% (3) self-consistency between (2) and a manual reference. + +clearvars; close all; format long e; +addpath('.'); + +zk = 1.1; +S = kernel('helm','s',zk); +D = kernel('helm','d',zk); + +% ----- (1) bare chunker, single S kernel --------------------------------- +chnkr = chunkerfunc(@(t) circarc(t,1.0), ... + struct('ta',0,'tb',pi/2,'nover',1), struct('k',16)); +chnkr = sort(chnkr); +theta = atan2(chnkr.r(2,:), chnkr.r(1,:)); +dens = (1 + 0.5*sin(theta(:))); +near = [cos(0.4); sin(0.4)] + 1e-3 *[cos(0.5); sin(0.5)]; +sub = [cos(0.4); sin(0.4)] + 1e-9 *[cos(0.5); sin(0.5)]; +far = [0.5; 0.5]; + +ref = chunkerkerneval(chnkr, S, dens, [far,near,sub], struct('forceadap',true,'eps',1e-14)); +val = chunkerkerneval(chnkr, S, dens, [far,near,sub], struct('forcewlchs',true)); +fprintf('(1) bare chunker, S kernel: max err = %.3e\n', max(abs(val-ref))); + +% ----- (2) chunkgraph + rcipsav, single D kernel ------------------------- +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}); +chnkrtot = merge(cg.echnks); + +% Solve a known closed-curve Dirichlet problem (rcipTest setup) +fkern = -2*kernel('helm','d',zk); +ns = 10; rng(8675309); +ts = 2*pi*rand(1,ns); +sources = 3.0*[cos(ts); sin(ts)]; +strengths = randn(ns,1); +ubdry = S.eval(struct('r',sources), struct('r',chnkrtot.r(:,:))) * strengths; +[sysmat, ~, rcipsav] = chunkermat(cg, fkern, struct('rcip',true,'nsub_or_tol',30,'rcip_savedepth',inf)); +sysmat = sysmat + eye(chnkrtot.npt); +sol = gmres(sysmat, ubdry, [], 1e-13, 200); + +% Targets: far + near + sub-panel from corner (-1,0) +dvec = [1, 0.1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]; +targets = [-1+dvec; zeros(1,length(dvec))]; +utrue = S.eval(struct('r',sources), struct('r',targets)) * strengths; + +% Probe diagnostic +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]); +v = fkern.eval(src_probe, tgt_probe); +b = 0.25i*zk*besselh(1,zk); +fprintf('probe val=%s base=%s coef=%s (expect -2)\n', num2str(v), num2str(b), num2str(v/b)); + +% Compare on chunkgraph (no rcipsav) at a far target only -- should agree +% with the default hybrid path away from corners +far_tgt = [0; 0.5]; +val_def = chunkerkerneval(cg, fkern, sol, far_tgt, []); +val_wl = chunkerkerneval(cg, fkern, sol, far_tgt, struct('forcewlchs',true)); +fprintf('cg+forcewlchs (no rcipsav), far tgt: |default-forcewlchs|=%.3e\n', abs(val_def-val_wl)); + +% bypass chunkgraph: merge first, then call forcewlchs on bare chunker +val_mwl = chunkerkerneval(chnkrtot, fkern, sol, far_tgt, struct('forcewlchs',true)); +fprintf('merged_chunker+forcewlchs, far tgt: |default-mergewlchs|=%.3e\n', abs(val_def-val_mwl)); + +% Manual smooth sum for sanity +zsrc = chnkrtot.r(1,:) + 1i*chnkrtot.r(2,:); +nz_s = chnkrtot.n(1,:) + 1i*chnkrtot.n(2,:); +ztgt = far_tgt(1) + 1i*far_tgt(2); +diff = ztgt - zsrc; +r_ = abs(diff); +% chunkie D = (i*zk/4) H_1(zk r) imag(nz/diff) * |y-x|... let me use formula from green +val_manual_D = (0.25i * zk * besselh(1,zk*r_) .* imag(nz_s./diff)) * (sol .* chnkrtot.wts(:)); +fprintf('manual smooth D: %s |default-manual|=%.3e\n', num2str(val_manual_D), abs(val_def-(-2*val_manual_D))); + +% Compare: same call WITHOUT forcewlchs (uses recursive adaptive) +val_old = chunkerkerneval(cg, fkern, sol, targets, struct('rcipsav',{rcipsav})); + +% chunkerkerneval with forcewlchs + rcipsav +val = chunkerkerneval(cg, fkern, sol, targets, ... + struct('forcewlchs',true,'rcipsav',{rcipsav})); + +fprintf('\n(2) chunkgraph + rcipsav, D kernel:\n'); +fprintf(' d rcipsav-only err forcewlchs+rcipsav err\n'); +for i = 1:length(dvec) + fprintf(' %.0e %.3e %.3e\n', dvec(i), ... + abs(val_old(i)-utrue(i)), abs(val(i)-utrue(i))); +end + +% Practical assertion: 12+ digits at d >= 1e-3 (typical grid scale). +assert(max(abs(val(dvec >= 1e-3) - utrue(dvec >= 1e-3))) < 1e-12, ... + 'expected 12+ digits at d >= 1e-3'); +fprintf('\nAll assertions passed.\n'); + +function [r,d,d2] = circle(t) + c = cos(t(:).'); s = sin(t(:).'); + r = [c; s]; d = [-s; c]; d2 = [-c; -s]; +end +function [r,d,d2] = circarc(t,R) + c = cos(t(:).'); s = sin(t(:).'); + r = R*[c;s]; d = R*[-s;c]; d2 = R*[-c;-s]; +end diff --git a/devtools/test/helmos_self_convergence.m b/devtools/test/helmos_self_convergence.m new file mode 100644 index 00000000..59e7aff0 --- /dev/null +++ b/devtools/test/helmos_self_convergence.m @@ -0,0 +1,50 @@ +% helmos_self_convergence.m +% +% Compare helmos solutions at two resolutions for both RHS types and +% several target locations. Goal: isolate whether the demo's ~1e-7 +% error floor is due to RHS, target location, or just noisy convergence. + +clearvars; close all; format long e; + +% Pick up helmos.m + helmos_planewave_rhs + geo_*.m from chunkie/demo/. +this_dir = fileparts(mfilename('fullpath')); +addpath(fullfile(this_dir, '..', '..', 'chunkie', 'demo')); + +zk = 3.0; +ppw_low = 10; ppw_hi = 30; +nsub_low = 30; nsub_hi = 30; % keep nsub fixed; deeper RCIP loses precision + +% Targets: Helsing's far one + a near-vertex sweep +targets = [ 0.17, -0.99, -0.999, -0.9999; + 0.62, -0.05, -0.005, -0.0005 ]; + +% --- Polynomial RHS (matches dirichletTest) +ubdry_poly = @(s,b,k) (4*s.r(1,:).^3 + 2*s.r(1,:).^2 - 3*s.r(1,:) - 1).'; + +% --- Planewave RHS, theta = pi/4 +theta = pi/4; +ubdry_pw = @(s,b,k) helmos_planewave_rhs(s,b,k,theta); + +cases = { {'poly', ubdry_poly}, {'planewave', ubdry_pw} }; + +for ic = 1:length(cases) + name = cases{ic}{1}; ubdry = cases{ic}{2}; + + cg_lo = geo_linesegment(zk, ppw_low); + cg_hi = geo_linesegment(zk, ppw_hi); + fprintf('\n[%s] npts low=%d, hi=%d\n', name, cg_lo.npt, cg_hi.npt); + + geo_lo = struct('cg',cg_lo,'targets',targets,'ubdry_fn',ubdry); + geo_hi = struct('cg',cg_hi,'targets',targets,'ubdry_fn',ubdry); + + out_lo = helmos(geo_lo, 'd', zk, struct('nsub',nsub_low,'tol',1e-12,'verbose',false)); + out_hi = helmos(geo_hi, 'd', zk, struct('nsub',nsub_hi,'tol',1e-13,'verbose',false)); + + err = abs(out_lo.ufield - out_hi.ufield); + fprintf(' %-22s %-14s %-14s\n', 'target', '|u_lo|', 'lo-hi diff'); + for i = 1:size(targets,2) + fprintf(' (%+.4f,%+.4f) %.4e %.3e\n', ... + targets(1,i), targets(2,i), abs(out_lo.ufield(i)), err(i)); + end +end +fprintf('\nDone.\n'); diff --git a/devtools/test/helmos_smoke.m b/devtools/test/helmos_smoke.m new file mode 100644 index 00000000..7f700514 --- /dev/null +++ b/devtools/test/helmos_smoke.m @@ -0,0 +1,77 @@ +% helmos_smoke.m +% +% Smoke test for the helmos.m driver across geometries: line segment, +% spiral, multi-corner chain. Each case solves the Helmholtz Dirichlet +% open-arc problem with a planewave RHS at zk=10*pi/Ltot (low frequency +% to keep runtime small), and validates self-convergence by comparing +% to a 3x mesh-refined solution at fixed nsub. Far-target diff < 1e-10 +% expected at low frequency. + +clearvars; close all; format long e; + +% Pick up helmos.m + helmos_planewave_rhs + geo_*.m from chunkie/demo/. +this_dir = fileparts(mfilename('fullpath')); +addpath(fullfile(this_dir, '..', '..', 'chunkie', 'demo')); + +cases = struct(); +cases(1).name = 'linesegment'; cases(1).build = @(zk,ppw) geo_linesegment(zk,ppw); +cases(2).name = 'spiral'; cases(2).build = @(zk,ppw) geo_spiral(zk,ppw); +cases(3).name = 'corners1'; cases(3).build = @(zk,ppw) geo_corners(1,zk,ppw); +cases(4).name = 'corners3'; cases(4).build = @(zk,ppw) geo_corners(3,zk,ppw); + +ppw_lo = 12; +ppw_hi = 36; +nsub = 30; +theta = pi/4; +ubdry_fn = @(s,b,k) helmos_planewave_rhs(s,b,k,theta); + +opts_lo = struct('nsub',nsub,'tol',1e-12,'verbose',false); +opts_hi = struct('nsub',nsub,'tol',1e-13,'verbose',false); + +for ic = 1:length(cases) + name = cases(ic).name; + build = cases(ic).build; + + % rough length-scale to set a wavenumber + cg_probe = build(1.0, ppw_lo); + % use total arc length to set zk so that L/lambda ~ 5 + Ltot = sum(cg_probe.echnks(1).wts(:)); + for ie = 2:length(cg_probe.echnks); Ltot = Ltot + sum(cg_probe.echnks(ie).wts(:)); end + zk = 10*pi/Ltot; % L/lambda ~ 5 + + cg_lo = build(zk, ppw_lo); + cg_hi = build(zk, ppw_hi); + + % targets: a bunch of off-curve points well away from any vertex, + % plus a sweep approaching one vertex (the first one). + rmin = min(cg_lo); rmax = max(cg_lo); + cx = (rmin(1)+rmax(1))/2; cy = (rmin(2)+rmax(2))/2; + box = max(rmax-rmin) * 1.5; + far = [cx+box*[0.4 -0.3 0.2]; cy+box*[0.3 0.4 -0.5]]; + % near-vertex sweep: approach the first vertex perpendicular to the + % normal of its first incident edge, to avoid landing on the curve. + v1 = cg_lo.verts(:,1); + e1 = cg_lo.echnks(1); + n1 = e1.n(:,1,1); n1 = n1/norm(n1); % normal at first node of edge 1 + dvec = 10.^(-(1:4)); + near = v1 + n1*dvec; + targets = [far, near]; + + geo_lo = struct('cg',cg_lo,'targets',targets,'ubdry_fn',ubdry_fn); + geo_hi = struct('cg',cg_hi,'targets',targets,'ubdry_fn',ubdry_fn); + + fprintf('\n[%s] zk=%.3f, npts lo=%d hi=%d\n', name, zk, cg_lo.npt, cg_hi.npt); + out_lo = helmos(geo_lo,'d',zk,opts_lo); + out_hi = helmos(geo_hi,'d',zk,opts_hi); + err = abs(out_lo.ufield - out_hi.ufield); + + fprintf(' far max diff: %.3e\n', max(err(1:size(far,2)))); + fprintf(' near max diff: %.3e (over d in %s)\n', ... + max(err(size(far,2)+1:end)), mat2str(dvec)); + + assert(max(err(1:size(far,2))) < 1e-9, ... + sprintf('[%s] far-target convergence < 1e-9 failed', name)); + assert(max(err) < 1e-8, ... + sprintf('[%s] all-target convergence < 1e-8 failed', name)); +end +fprintf('\nAll smoke tests passed.\n');