Axissym kernels

.gitignore

@@ -19,3 +19,6 @@ chunkie/FLAM
 chunkie/chunkergraph/test_chnkgraph_local.m
 chunkie/chunkergraph/local_tests/test_chnkgraph_local.m
 chunkie/chunkergraph/local_tests/test_chnkgraph_local2.m
+
+*k2test*
+*ktest*

chunkie/+chnk/+axissymhelm2d/dalph.m

@@ -0,0 +1,48 @@
+function [dout] = dalph(r,dr,z)
+
+    dout = {};
+
+    den = 2*r.^2+2*r.*dr+dr.^2+z.^2;
+    num = 2*r.*dr+dr.^2+z.^2;
+    val = num./den.^2;
+    ar = -2*(r+dr).*val;
+    num = -2*r.*dr-dr.^2+z.^2;
+    val = num./den.^2;
+    adr= -2*r.*val; 
+    az = 4*(r+dr).*r.*z./den.^2;
+
+    dout.ar = ar;
+    dout.adr= adr;
+    dout.az = az;
+
+    ardr = dr.^4+4*dr.^3.*r-8*dr.*r.^3-4*r.^4-z.^4;
+    ardr = 2*ardr./den.^3;
+
+    arz = 4*(r+dr).*z.*(dr.^2+2*dr.*r-2*r.^2+z.^2)./den.^3;
+    adrz= -4*r.*z.*(3*dr.^2+6*dr.*r+2*r.^2-z.^2)./den.^3;
+    azz = 4*r.*(r+dr).*(dr.^2+2*dr.*r+2*r.^2-3*z.^2)./den.^3;
+
+    dout.ardr = ardr;
+    dout.arz  = arz;
+    dout.adrz = adrz;
+    dout.azz  = azz;
+
+    r0   = sqrt(r.^2+(r+dr).^2+z.^2);
+    kr   = r./r0;
+    kdr  = (r+dr)./r0;
+    kz   = z./r0;
+    krdr = -r.*(r+dr)./r0.^3;
+    krz  = -r.*z./r0.^3;
+    kdrz = -(r+dr).*z./r0.^3;
+    kzz  = (2*r.^2+2*r.*dr+dr.^2)./r0.^3;
+
+    dout.kr  = kr;
+    dout.kdr = kdr;
+    dout.kz  = kz;
+    dout.krdr= krdr;
+    dout.krz = krz;
+    dout.kdrz= kdrz;
+    dout.kzz = kzz;
+
+end
+

chunkie/+chnk/+axissymhelm2d/der_ak_to_grad.m

@@ -0,0 +1,39 @@
+function [dout] = der_ak_to_grad(r,q,z,i,ida,idk,...
+    idaa,idak,idkk)
+
+    dout = {};
+    dout.int = i;
+
+    ds = chnk.axissymhelm2d.dalph(r,q,z);
+
+    intdr = ds.ar.*ida +ds.kr.*idk;
+    intdq = ds.adr.*ida+ds.kdr.*idk;
+    intdz = ds.az.*ida +ds.kz.*idk;
+
+    intdrq= ds.ardr.*ida  +  ds.krdr.*idk  +  ...
+            ds.ar.*ds.adr.*idaa +  ds.kr.*ds.kdr.*idkk +  ...
+        ds.ar.*ds.kdr.*idak+ds.kr.*ds.adr.*idak;
+
+    intdrz= ds.arz.*ida  +  ds.krz.*idk  +  ...
+            ds.ar.*ds.az.*idaa +  ds.kr.*ds.kz.*idkk +  ...
+        ds.ar.*ds.kz.*idak+ds.kr.*ds.az.*idak;    
+
+    intdqz= ds.adrz.*ida  +  ds.kdrz.*idk  +  ...
+            ds.adr.*ds.az.*idaa +  ds.kdr.*ds.kz.*idkk +  ...
+        ds.adr.*ds.kz.*idak+ds.kdr.*ds.az.*idak;     
+
+    intdzz= ds.azz.*ida  +  ds.kzz.*idk  +  ...
+            ds.az.*ds.az.*idaa +  ds.kz.*ds.kz.*idkk +  ...
+        ds.az.*ds.kz.*idak+ds.kz.*ds.az.*idak;    
+
+    dout.intdr = intdr;
+    dout.intdq = intdq;
+    dout.intdz = intdz;
+
+    dout.intdrq = intdrq;
+    dout.intdrz = intdrz;
+    dout.intdqz = intdqz;
+    dout.intdzz = intdzz;
+
+end
+

chunkie/+chnk/+axissymhelm2d/div_by_kap.m

@@ -0,0 +1,36 @@
+function [dout] = div_by_kap(r,q,z,dfunc)
+
+    r0   = sqrt(r.^2+(r+q).^2+z.^2);
+    i0   =  1./r0;
+    i0da =  0;
+    i0dk = -1./r0.^2;
+    i0daa=  0;
+    i0dak=  0;
+    i0dkk=  2./r0.^3;
+    [dout] = chnk.axissymhelm2d.der_ak_to_grad(r,q,z,i0,i0da,i0dk,...
+    i0daa,i0dak,i0dkk);
+
+    i   = dfunc.int.*i0;
+    idr = dfunc.int.*dout.intdr + dfunc.intdr.*dout.int;
+    idq = dfunc.int.*dout.intdq + dfunc.intdq.*dout.int;
+    idz = dfunc.int.*dout.intdz + dfunc.intdz.*dout.int;
+    idqr = dfunc.intdrq.*dout.int + dfunc.int.*dout.intdrq + ...
+        dfunc.intdr.*dout.intdq + dfunc.intdq.*dout.intdr;
+    idrz = dfunc.intdrz.*dout.int + dfunc.int.*dout.intdrz + ...
+        dfunc.intdr.*dout.intdz + dfunc.intdz.*dout.intdr;
+    idqz = dfunc.intdqz.*dout.int + dfunc.int.*dout.intdqz + ...
+        dfunc.intdq.*dout.intdz + dfunc.intdz.*dout.intdq;
+    idzz = dfunc.intdzz.*dout.int + dfunc.int.*dout.intdzz + ...
+        2*dfunc.intdz.*dout.intdz;
+
+    dout.int   = i;
+    dout.intdr = idr;
+    dout.intdq = idq;
+    dout.intdz = idz;
+    dout.intdrq = idqr;
+    dout.intdrz = idrz;
+    dout.intdqz = idqz;
+    dout.intdzz = idzz;
+
+end
+

chunkie/+chnk/+axissymhelm2d/green.m

@@ -0,0 +1,107 @@
+function [val, grad, hess] = green(k, src, targ, origin, ifdiff)
+%CHNK.AXISSYMHELM2D.GREEN evaluate the helmholtz green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code relies on precomputed tables and hence loops are required for 
+% computing various pairwise interactions.
+% Finally, the code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Since the kernels are not translationally invariant in r, the size
+% of the gradient is 3, for d/dr, d/dr', and d/dz
+%
+% Similarly the hessian is of size 6 and ordered as 
+%   d_{rr}, d_{r'r'}, d_{zz}, d_{rr'}, d_{rz}, d_{r'z}
+
+
+zr = real(k); zi = imag(k);
+if abs(zr*zi) > eps
+    error('AXISSYMHELM2D.green: Axissymmetric greens function not supported for arbitrary complex', ...
+           'wavenumbers, please input purely real or imaginary wave numbers\n');
+end
+
+if nargin == 4
+    ifdiff = 0;
+end
+
+
+% Load precomputed tables
+persistent asym_tables   
+if isempty(asym_tables)
+    asym_tables = chnk.axissymhelm2d.load_asym_tables();
+end
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+vtmp = zeros(nt, ns);
+gtmp = zeros(nt, ns, 3);
+htmp = zeros(nt, ns, 6);
+
+% Determine whether to use 
+ifun = 1;
+if abs(zr) < eps
+    ifun = 2;
+end
+
+if ifdiff
+    if abs(zi) > eps
+        error('AXISSYMHELM2D.green: Difference of greens function only supported for real wave numbers\n');
+    end
+    ifun = 3;
+end
+
+over2pi = 1/2/pi;
+kabs = abs(k);
+
+rt = repmat(targ(1,:).',1,ns);
+rs = repmat(src(1,:),nt,1);
+dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns);
+r  = (rt + origin(1))*kabs;
+dz = dz*kabs;
+dr = (rs-rt)*kabs;
+dout = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun);
+int = dout.int;
+intdz = dout.intdz;
+intdq = dout.intdq;
+intdr = dout.intdr;
+intdzz = dout.intdzz;
+intdrq = dout.intdrq;
+intdrz = dout.intdrz;
+intdqz = dout.intdqz;
+for j = 1:ns
+    darea = src(1,j) + origin(1);
+    for i = 1:nt
+        if (i >1 && size(int,1)<=1)
+            disp("catastrophic error")
+        end
+
+        vtmp(i,j) = int(i,j) * over2pi * kabs * darea;
+
+        gtmp(i,j,1) = intdr(i,j) * over2pi * kabs.^2 * darea;
+        gtmp(i,j,2) = intdq(i,j) * over2pi * kabs.^2 * darea;  
+        gtmp(i,j,3) = -intdz(i,j) * over2pi * kabs.^2 * darea;
+
+% Fix hessian scalings
+% No drr, dr'r', or dr r' currently available
+        %htmp(i,j,1) = 0;
+        %htmp(i,j,2) = 0;
+        htmp(i,j,3) = intdzz(i,j) * over2pi * kabs.^3 * darea;
+        htmp(i,j,4) = intdrq(i,j) * over2pi * kabs.^3 * darea;
+        htmp(i,j,5) = intdrz(i,j) * over2pi * kabs.^3 * darea;
+        htmp(i,j,6) = intdqz(i,j) * over2pi * kabs.^3 * darea;
+    end
+end
+
+if nargout > 0
+    val = vtmp;
+end
+
+if nargout > 1
+    grad = gtmp;
+end
+
+if nargout > 2
+    hess = htmp;
+end

chunkie/+chnk/+axissymhelm2d/green_neu_all.m

@@ -0,0 +1,87 @@
+function [valk, gradk, hessk, valik, gradik, ...
+     hessik, valdiff, graddiff, hessdiff] = green_neu_all(k, src, targ, origin)
+%CHNK.AXISSYMHELM2D.GREEN evaluate the helmholtz green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code relies on precomputed tables and hence loops are required for 
+% computing various pairwise interactions.
+% Finally, the code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Since the kernels are not translationally invariant in r, the size
+% of the gradient is 3, for d/dr, d/dr', and d/dz
+%
+% Similarly the hessian is of size 6 and ordered as 
+%   d_{rr}, d_{r'r'}, d_{zz}, d_{rr'}, d_{rz}, d_{r'z}
+
+% Load precomputed tables
+persistent asym_tables   
+if isempty(asym_tables)
+    asym_tables = chnk.axissymhelm2d.load_asym_tables();
+end
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+over2pi = 1/2/pi;
+kabs = abs(k);
+
+rt = repmat(targ(1,:).',1,ns);
+rs = repmat(src(1,:),nt,1);
+dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns);
+r  = (rt + origin(1))*kabs;
+dz = dz*kabs;
+dr = (rs-rt)*kabs;
+ifun = 4;
+[doutk, doutik, doutdiff] = chnk.axissymhelm2d.helm_axi_all(r, dr, dz, asym_tables,ifun);
+
+pfac = over2pi * kabs;
+gfac = pfac * kabs;
+hfac = gfac * kabs;
+[valk, gradk, hessk] = dout_to_vgh(doutk, ns, nt, src, origin, pfac, gfac, hfac);
+[valik, gradik, hessik] = dout_to_vgh(doutik, ns, nt, src, origin, pfac, gfac, hfac);
+[valdiff, graddiff, hessdiff] = dout_to_vgh(doutdiff, ns, nt, src, origin, pfac, gfac, hfac);
+
+end
+
+function [vtmp, gtmp, htmp] = dout_to_vgh(dout, ns, nt, src, origin, pfac, gfac, hfac)
+
+    vtmp = zeros(nt, ns);
+    gtmp = zeros(nt, ns, 3);
+    htmp = zeros(nt, ns, 6);
+
+    int = dout.int;
+    intdz = dout.intdz;
+    intdq = dout.intdq;
+    intdr = dout.intdr;
+    intdzz = dout.intdzz;
+    intdrq = dout.intdrq;
+    intdrz = dout.intdrz;
+    intdqz = dout.intdqz;
+    for j = 1:ns
+        darea = src(1,j) + origin(1);
+        pfacsc = pfac * darea;
+        gfacsc = gfac * darea;
+        hfacsc = hfac * darea;
+        for i = 1:nt
+            if (i >1 && size(int,1)<=1)
+                disp("catastrophic error")
+            end
+
+            vtmp(i,j) = int(i,j) * pfacsc;
+
+            gtmp(i,j,1) = intdr(i,j) * gfacsc;
+            gtmp(i,j,2) = intdq(i,j) * gfacsc;  
+            gtmp(i,j,3) = -intdz(i,j) * gfacsc;
+
+% Fix hessian scalings
+% No drr, dr'r', or dr r' currently available
+
+            htmp(i,j,3) = intdzz(i,j) * hfacsc;
+            htmp(i,j,4) = intdrq(i,j) * hfacsc;
+            htmp(i,j,5) = intdrz(i,j) * hfacsc;
+            htmp(i,j,6) = intdqz(i,j) * hfacsc;
+        end
+    end
+end