diff --git a/chunkie/+chnk/+axissymlap2d/g0funcall.m b/chunkie/+chnk/+axissymlap2d/g0funcall.m index fe2a3bdc..5e86aeb4 100644 --- a/chunkie/+chnk/+axissymlap2d/g0funcall.m +++ b/chunkie/+chnk/+axissymlap2d/g0funcall.m @@ -62,19 +62,7 @@ iffwd = 0; end - end - - %disp(['inside g0mall, x = ' num2str(x)]) - %if (iffwd == 0) - % disp(['running backward recursion']); - %end - - %if (iffwd == 1) - % disp(['running backward recursion']); - %end - - %return - + end gvals = zeros(maxm+1,1); @@ -177,7 +165,6 @@ dernext = 0; der = 1; - % run the downward recurrence for j = 1:(nterms-maxm+1) i = nterms-j+1; diff --git a/chunkie/+chnk/+axissymlap2d/g0funcall_vec.m b/chunkie/+chnk/+axissymlap2d/g0funcall_vec.m new file mode 100644 index 00000000..1f01f2c1 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/g0funcall_vec.m @@ -0,0 +1,57 @@ +function [gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = g0funcall_vec(r, rp, dr, z, zp, dz, m) +% +% chnk.axissymlap2d.g0funcall evaluates a collection of axisymmetric Laplace +% Green's functions, defined by the expression: +% +% gfunc(n) = pi*rp * \int_0^{2\pi} 1/|x - x'| e^(-i n t) dt +% +% it is assumed that x = (x,0,z) otherwise the integral above should pick up a +% phase factor out front of exp(i*n*phi), where phi is the azimuthal coordinate +% of x in cylindrical coordinates. +% +% The extra factor of rp (and maybe pi?) out front makes subsequent interfacing +% with RCIP slightly easier. Modes 0 through maxm are returned, with gval(1) = +% mode 0 and gval(maxm+1) = mode maxm. The function is even, so g_{-n} = g_n. +% +% The above scaling should be consistent with what is in +% chnk.axissymlap2d.gfunc, which is for merely the zero-mode +% + + r = reshape(r,[1,size(r,1),size(r,2)]); + rp = reshape(rp,[1,size(rp,1),size(rp,2)]); + dr = reshape(dr,[1,size(dr,1),size(dr,2)]); + z = reshape(z,[1,size(z,1),size(z,2)]); + zp = reshape(zp,[1,size(zp,1),size(zp,2)]); + dz = reshape(dz,[1,size(dz,1),size(dz,2)]); + + t = (dz.^2+dr.^2)./(2.*r.*rp); + chi = t+1; + + [qm, qmd, qmdd] = chnk.axissymlap2d.qleg_half_miller_vec(t,m); + + gval = 2*pi*sqrt(rp./r).*qm; + gdz = 2*pi*sqrt(rp./r).*qmd ... + ./(rp.*r).*dz; + + rfac = -r/2.*qm+(-(1+t).*r+rp).*qmd; + gdrp = 2*pi*sqrt(rp./r)./(rp.*r) ... + .*rfac; + + rfac = -rp/2.*qm+(-(1+t).*rp+r).*qmd; + gdr = 2*pi*sqrt(rp./r)./(rp.*r) ... + .*rfac; + + rfac = 1./(rp.*r).*qmd + (dz./(rp.*r)).^2.*qmdd; + gdzz = 2*pi*sqrt(rp./r).*rfac; + + rfac = -3./(2*r.^2.*rp).*qmd + (-chi./(r.^2.*rp) + 1./(r.*rp.^2)).*qmdd; + gdrz = 2*pi*sqrt(rp./r).*dz.*rfac; + + rfac = -3./(2*rp.^2.*r).*qmd + (-chi./(rp.^2.*r) + 1./(rp.*r.^2)).*qmdd; + gdrpz = 2*pi*sqrt(rp./r).*dz.*rfac; + + rfac = 1./(4*r.*rp).*qm ... + + (2*chi./(rp.*r) - 3./(2*r.^2) - 3./(2*rp.^2)).*qmd ... + + (-chi./r+1./rp).*(-chi./rp + 1./r).*qmdd; + gdrpr = 2*pi*sqrt(rp./r).*rfac; +end \ No newline at end of file diff --git a/chunkie/+chnk/+axissymlap2d/gfunc.m b/chunkie/+chnk/+axissymlap2d/gfunc.m index 1381c888..5e91c4b6 100644 --- a/chunkie/+chnk/+axissymlap2d/gfunc.m +++ b/chunkie/+chnk/+axissymlap2d/gfunc.m @@ -1,4 +1,4 @@ -function [gval, gdz, gdr, gdrp] = gfunc(r, rp, dr, z, zp, dz) +function [gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = gfunc(r, rp, dr, z, zp, dz) % % chnk.axissymlap2d.gfunc evaluates the zeroth order axisymmetric Laplace @@ -14,8 +14,8 @@ n = 3; t = (dz.^2+dr.^2)./(2.*r.*rp); - [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t); - + [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t); + gval = domega3*(rp./r).^((n-2)/2).*q0; gdz = domega3*(rp./r).^((n-2)/2).*q0d ... ./(rp.*r).*dz; @@ -27,6 +27,27 @@ rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d; gdr = domega3*(rp./r).^((n-2)/2)./(rp.*r) ... .*rfac; + + + % Amandin hessian stuff + chi = 1+t; + q0dd = ((1/4)*q0 + 2*chi.*q0d)./(-t.^2-2*t); + + rfac = 1./(rp.*r).*q0d + (dz./(rp.*r)).^2.*q0dd; + gdzz = 2*pi*sqrt(rp./r).*rfac; + + rfac = -3./(2*r.^2.*rp).*q0d + (-chi./(r.^2.*rp) + 1./(r.*rp.^2)).*q0dd; + gdrz = 2*pi*sqrt(rp./r).*dz.*rfac; + + rfac = -3./(2*rp.^2.*r).*q0d + (-chi./(rp.^2.*r) + 1./(rp.*r.^2)).*q0dd; + gdrpz = 2*pi*sqrt(rp./r).*dz.*rfac; + + rfac = 1./(4*r.*rp).*q0 ... + + (2*chi./(rp.*r) - 3./(2*r.^2) - 3./(2*rp.^2)).*q0d ... + + (-chi./r+1./rp).*(-chi./rp + 1./r).*q0dd; + gdrpr = 2*pi*sqrt(rp./r).*rfac; + + % % check the scaling here correctly diff --git a/chunkie/+chnk/+axissymlap2d/green_0th_mode.m b/chunkie/+chnk/+axissymlap2d/green_0th_mode.m new file mode 100644 index 00000000..e4545393 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/green_0th_mode.m @@ -0,0 +1,53 @@ +function [val, grad, hess] = green_0th_mode(src, targ, origin) +% +% CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace 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 +% +% Returns for gradient are: +% grad = d_{r}, d_{r'}, d_{z}, d_{z'} +% +% Returns for hess are: +% hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z} + +[~, ns] = size(src); +[~, nt] = size(targ); + +rt = repmat(targ(1,:).',1,ns); % r +rs = repmat(src(1,:),nt,1); % r' +r = (rt + origin(1)); +rp = (rs + origin(1)); +dr = rt-rs; % r - r' +z = repmat(targ(2,:).',1,ns); +zp = repmat(src(2,:),nt,1); +dz = z-zp; % z - z' + +[gs,gdzs,gdrs,gdrps,gdrpr,gdzz,gdrz,gdrpz] = chnk.axissymlap2d.gfunc(r,rp,dr,z,zp,dz); + +const = 1/(4*pi^2); + +val = gs*const; +grad = zeros(nt, ns, 4); +grad(:,:,1) = gdrs; +grad(:,:,2) = gdrps; +grad(:,:,3) = gdzs; +grad(:,:,4) = -gdzs; +grad = grad*const; + +hess = zeros(nt, ns, 4); +hess(:,:,1) = gdrpr; +hess(:,:,2) = -gdzz; +hess(:,:,3) = -gdrz; +hess(:,:,4) = gdrpz; +hess = hess*const; + +% if anynan(hess) +% error('Execution stopped: NaN value detected.'); +% end + +end diff --git a/chunkie/+chnk/+axissymlap2d/green_modal.m b/chunkie/+chnk/+axissymlap2d/green_modal.m new file mode 100644 index 00000000..09f7745d --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/green_modal.m @@ -0,0 +1,58 @@ +function [val, grad, hess] = green_modal(src, targ, origin, m) +% +% CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace modal green's function +% for the given sources and targets. +% +% Note: that the first coordinate is r, and the second z. +% The code is not efficient in the sense that val, grad, hess +% are always internally computed independent of nargout +% +% Returns for gradient are: +% grad = d_{r}, d_{r'}, d_{z}, d_{z'} +% +% Returns for hess are: +% hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z} +% +% m is the mode we return +% +% all_modes = true: means we return size [m*nt, ns] instead of [nt, ns] + +[~, ns] = size(src); +[~, nt] = size(targ); + +rt = repmat(targ(1,:).',1,ns); % r +rs = repmat(src(1,:),nt,1); % r' +r = (rt + origin(1)); +rp = (rs + origin(1)); +dr = rt-rs; % r - r' +z = repmat(targ(2,:).',1,ns); +zp = repmat(src(2,:),nt,1); +dz = z-zp; % z - z' + +[gs,gdzs,gdrs,gdrps,gdrprs,gdzzs,gdrzs,gdrpzs] = chnk.axissymlap2d.g0funcall_vec(r,rp,dr,z,zp,dz,m); + + +const = 1/(4*pi^2); + +val = gs(m,:,:); +val = reshape(val,[nt,ns]); +val = const*val; + +grad = zeros(nt, ns, 4); +grad(:,:,1) = gdrs(m,:,:); +grad(:,:,2) = gdrps(m,:,:); +grad(:,:,3) = gdzs(m,:,:); +grad(:,:,4) = -gdzs(m,:,:); +grad = const*grad; + +hess = zeros(nt, ns, 4); +hess(:,:,1) = gdrprs(m,:,:); +hess(:,:,2) = -gdzzs(m,:,:); +hess(:,:,3) = -gdrzs(m,:,:); +hess(:,:,4) = gdrpzs(m,:,:); +hess = const*hess; + + + + +end diff --git a/chunkie/+chnk/+axissymlap2d/kern.m b/chunkie/+chnk/+axissymlap2d/kern.m index 9560b421..b4bca02d 100644 --- a/chunkie/+chnk/+axissymlap2d/kern.m +++ b/chunkie/+chnk/+axissymlap2d/kern.m @@ -74,8 +74,28 @@ end if strcmpi(type, 's') +<<<<<<< HEAD + submat = chnk.axissymlap2d.green_0th_mode(src, targ, origin); +end + +if strcmpi(type, 'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + nxt = repmat((targnorm(1,:)).',1,ns); + nyt = repmat((targnorm(2,:)).',1,ns); + nxs = repmat(srcnorm(1,:), nt, 1); + nys = repmat(srcnorm(2,:), nt, 1); + + % hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z} + [~, ~, hess] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + submat = (hess(:,:,1).*nxt.*nxs + hess(:,:,3).*nys.*nxt ... + + hess(:,:,4).*nxs.*nyt + hess(:,:,2).*nyt.*nys); +end + +======= submat = chnk.axissymlap2d.green(src, targ, origin); end +>>>>>>> upstream/lap_axisym diff --git a/chunkie/+chnk/+axissymlap2d/kern_0th_mode.m b/chunkie/+chnk/+axissymlap2d/kern_0th_mode.m new file mode 100644 index 00000000..97228706 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/kern_0th_mode.m @@ -0,0 +1,119 @@ +function submat = kern_0th_mode(srcinfo, targinfo, origin, type) +%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D +% +% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targinfo,type) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). Note that the normal information is obtained +% by taking the perpendicular to the provided tangential deriviative +% info and normalizing +% +% Here the first and second components correspond to the r and z +% coordinates respectively. +% +% Kernels based on G(x,y) = \int_{0}^{\pi} 1/(d(t)) \, dt \, +% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with +% x = (r,z), and y = (r',z') +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} G(x,y) +% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) +% +% Input: +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targinfo - description of targets in ptinfo struct format, +% if info not relevant (d/d2) it doesn't need to +% be provided. sprime requires tangent info in +% targinfo.d +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'sprime', normal derivative of single +% layer S' +% +% +% Output: +% submat - the evaluation of the selected kernel for the +% provided sources and targets. the number of +% rows equals the number of targets and the +% number of columns equals the number of sources +% +% + +src = srcinfo.r; +targ = targinfo.r; + +[~, ns] = size(src); +[~, nt] = size(targ); + +if strcmpi(type, 'd') + srcnorm = srcinfo.n; + [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + nx = repmat(srcnorm(1,:), nt, 1); + ny = repmat(srcnorm(2,:), nt, 1); + % dr'*nr' + dz'*nz' + submat = -(grad(:,:,2).*nx + grad(:,:,4).*ny); +end + +if strcmpi(type, 's') + [val, ~] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + submat = val; +end + +if strcmpi(type, 'sprime') + targnorm = targinfo.n; + [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + nx = repmat((targnorm(1,:)).',1,ns); + ny = repmat((targnorm(2,:)).',1,ns); + % dr*nr + dz*nz + submat = -(grad(:,:,1).*nx + grad(:,:,3).*ny); +end + +if strcmpi(type, 'sgradr') + [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + submat = grad(:,:,1); +end + +if strcmpi(type, 'sgradz') + [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + submat = grad(:,:,3); +end + +if strcmpi(type, 'dgradr') + h = 1e-200; + targ_cr = targinfo; + targ_cr.r(1,:) = targ_cr.r(1,:) + 1i*h; + D_cr = kern_0th_mode(srcinfo, targ_cr, origin, 'd'); + submat = imag(D_cr)/h; +end + +if strcmpi(type, 'dgradz') + h = 1e-200; + targ_cz = targinfo; + targ_cz.r(2,:) = targ_cz.r(2,:) + 1i*h; + D_cz = kern_0th_mode(srcinfo, targ_cz, origin, 'd'); + submat = imag(D_cz)/h; +end + +if strcmpi(type, 'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + nxt = repmat((targnorm(1,:)).',1,ns); + nyt = repmat((targnorm(2,:)).',1,ns); + nxs = repmat(srcnorm(1,:), nt, 1); + nys = repmat(srcnorm(2,:), nt, 1); + + % hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z} + [~, ~, hess] = chnk.axissymlap2d.green_0th_mode(src, targ, origin); + submat = (hess(:,:,1).*nxt.*nxs + hess(:,:,3).*nys.*nxt ... + + hess(:,:,4).*nxs.*nyt + hess(:,:,2).*nyt.*nys); +end + +end \ No newline at end of file diff --git a/chunkie/+chnk/+axissymlap2d/kern_modal.m b/chunkie/+chnk/+axissymlap2d/kern_modal.m new file mode 100644 index 00000000..4b82db69 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/kern_modal.m @@ -0,0 +1,93 @@ +function submat = kern_modal(srcinfo, targinfo, origin, type, m) +%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D +% +% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targinfo,type) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). Note that the normal information is obtained +% by taking the perpendicular to the provided tangential deriviative +% info and normalizing +% +% Here the first and second components correspond to the r and z +% coordinates respectively. +% +% Kernels based on G(x,y) = \int_{0}^{\pi} 1/(d(t)) \, dt \, +% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with +% x = (r,z), and y = (r',z') +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} G(x,y) +% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) +% +% Input: +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targinfo - description of targets in ptinfo struct format, +% if info not relevant (d/d2) it doesn't need to +% be provided. sprime requires tangent info in +% targinfo.d +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'sprime', normal derivative of single +% layer S' +% m - int, mode +% +% +% Output: +% submat - the evaluation of the selected kernel for the +% provided sources and targets. the number of +% rows equals the number of targets and the +% number of columns equals the number of sources +% +% + +src = srcinfo.r; +targ = targinfo.r; + +[~, ns] = size(src); +[~, nt] = size(targ); + +if strcmpi(type, 'd') + srcnorm = srcinfo.n; + [~, grad] = chnk.axissymlap2d.green_modal(src, targ, origin, m); + nx = repmat(srcnorm(1,:), nt, 1); + ny = repmat(srcnorm(2,:), nt, 1); + % dr'*nr' + dz'*nz' + submat = -(grad(:,:,2).*nx + grad(:,:,4).*ny); +end + +if strcmpi(type, 's') + [val, ~] = chnk.axissymlap2d.green_modal(src, targ, origin, m); + submat = val; +end + +if strcmpi(type, 'sprime') + targnorm = targinfo.n; + [~, grad] = chnk.axissymlap2d.green_modal(src, targ, origin, m); + nx = repmat((targnorm(1,:)).',1,ns); + ny = repmat((targnorm(2,:)).',1,ns); + % dr*nr + dz*nz + submat = -(grad(:,:,1).*nx + grad(:,:,3).*ny); +end + +if strcmpi(type, 'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + nxt = repmat((targnorm(1,:)).',1,ns); + nyt = repmat((targnorm(2,:)).',1,ns); + nxs = repmat(srcnorm(1,:), nt, 1); + nys = repmat(srcnorm(2,:), nt, 1); + % hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z} + [~, ~, hess] = chnk.axissymlap2d.green_modal(src, targ, origin, m); + submat = (hess(:,:,1).*nxt.*nxs + hess(:,:,3).*nys.*nxt ... + + hess(:,:,4).*nxs.*nyt + hess(:,:,2).*nyt.*nys); +end + +end \ No newline at end of file diff --git a/chunkie/+chnk/+axissymlap2d/qleg_half_miller.m b/chunkie/+chnk/+axissymlap2d/qleg_half_miller.m new file mode 100644 index 00000000..bcc9ff42 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/qleg_half_miller.m @@ -0,0 +1,157 @@ +function [qm, qmd, qmdd] = qleg_half_miller(t, m) + + % + % Uses Miller's Algorithm to compute the half Legendre functions of the + % second kind + % + % Inputs + % t: chi - 1 + % m: all the modes we want to compute + % + % Outputs + % qm: [Q_{-1/2},...,Q_{m-1/2}] + % qmd: [dQ_{-1/2},...,dQ_{m-1/2}] + % qmdd: [d^2 Q_{-1/2},...,d^2 Q_{m-1/2}] + + % compute chi + chi = t+1; + + % initialize the seeds for the Legendre functions + [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t); + q1d = (-q0 + chi.*q1)/(2*(chi+1)*t); + + % check if we can run the forward recurrence + run_forward = 0; + if (chi < 1.005) + run_forward = 1; + end + + qm = zeros(m+1,1); + qmd = zeros(m+1,1); + qmdd = zeros(m+1,1); + + % run the forward reccurence + if run_forward == 1 + qm(1) = q0; + qm(2) = q1; + qmd(1) = q0d; + qmd(2) = q1d; + for i = 1:(m-1) + qm(2+i) = 4*(i-1)/(2*i-1)*chi*qm(1+i) - (2*i-3)/(2*i-1)*qm(i); + qmd(2+i) = (-(m-1/2)*chi.*qm(2+i) + (m+1/2)*qm(1+i))./(1-chi.^2); + end + + for i = 1:m+1 + qmdd(i) = (-(i-3/2)*(i-1/2)*qm(i) + 2*chi.*qmd(i))./(1-chi.^2); + end + + return + end + + % Millers Algorithm: + + % 1. run the forward reccurence until it has blown up + % recurrence intialization + fprev = 1; + fprevprev = 0; + + % NOTE: + % should maybe wait until second derivatives also blow up + + maxiter = 100000; % max number of extra terms + upbound = 1.0e19; % blow up tolerance + nterms = m+10; % minimum number of extra terms + for i = m:maxiter + f = 4*(i-1)/(2*i-1)*chi*fprev - (2*i-3)/(2*i-1)*fprevprev; + d = (-(m-1/2)*chi.*fprev + (m+1/2)*fprevprev)./(1-chi.^2); + + if (abs(f) >= upbound) + if (abs(d) >= upbound) + nterms = i+1; + break + end + end + + fprevprev = fprev; + fprev = f; + end + + % 2. run the backward reccurence + % recurrence intialization + fprev = 1; + fprevprev = 0; + + % run the backward reccurence from nterms to m + for j = 1:(nterms-m+1) + i = nterms-j+1 +1; + + f = 4*(i-1)/(2*i-3)*chi*fprev - (2*i-1)/(2*i-3)*fprevprev; + + fprevprev = fprev; + fprev = f; + end + + qm(m) = fprev; + qm(m+1) = fprevprev; + + % run the backward reccurence from m to 1 + for j = 1:(m-1) + i = m-1-j+1 +1; + qm(i-1) = 4*(i-1)/(2*i-3)*chi*qm(i) - (2*i-1)/(2*i-3)*qm(i+1); + end + + % 3. compute error and scale solution + ratio = q0/qm(1); + qm = qm.*ratio; + + % NOTE: + % combine these loops to speed things up ? + + % 4. compute 1st derivatives + qmd(1) = q0d; + for i = 1:m + qmd(i+1) = (-(i-1/2)*chi.*qm(i+1) + (i-1/2)*qm(i))./(1-chi.^2); + end + + % 5. compute 2nd derivatives + for i = 1:m+1 + qmdd(i) = (-(i-3/2)*(i-1/2)*qm(i) + 2*chi.*qmd(i))./(1-chi.^2); + end +end +%{ +function [mask_blowup, mask_stable, nterms] = get_masks(t, m) + % run the forward reccurence until it has blown up + % recurrence intialization + fprev = ones(size(t,1),size(t,2)); + fprevprev = zeros(size(t,1),size(t,2)); + + % NOTE: + % should maybe wait until second derivatives also blow up + + maxiter = 100000; % max number of extra terms + ubound = 1.0e20; % upper bound tolerance + lbound = 1.0e15; % lower bound tolerance + nterms = m; % minimum number of extra terms + for i = m:maxiter + f = 4*(i-1)/(2*i-1)*chi*fprev - (2*i-3)/(2*i-1)*fprevprev; + d = (-(m-1/2)*chi.*fprev + (m+1/2)*fprevprev)./(1-chi.^2); + + if (max(abs(f),[],"all") >= ubound) + if (max(abs(d),[],"all") >= ubound) + nterms = i+1; + break + end + end + + fprevprev = fprev; + fprev = f; + end + + if (nterms < 10) + nterms = 10; + end + + mask_blowup = (abs(f) >= lbound) & (abs(d) >= lbound); + mask_stable = ~mask_blowup; +end +%} \ No newline at end of file diff --git a/chunkie/+chnk/+axissymlap2d/qleg_half_miller_vec.m b/chunkie/+chnk/+axissymlap2d/qleg_half_miller_vec.m new file mode 100644 index 00000000..c5bf1f98 --- /dev/null +++ b/chunkie/+chnk/+axissymlap2d/qleg_half_miller_vec.m @@ -0,0 +1,181 @@ +function [qm, qmd, qmdd] = qleg_half_miller_vec(t, m) + + % + % Uses Miller's Algorithm to compute the half Legendre functions of the + % second kind + % + % Inputs + % t: chi - 1 + % m: all the modes we want to compute + % + % Outputs + % qm: [Q_{-1/2},...,Q_{m-1/2}] + % qmd: [dQ_{-1/2},...,dQ_{m-1/2}] + % qmdd: [d^2 Q_{-1/2},...,d^2 Q_{m-1/2}] + % + + % compute chi + chi = t+1; + + % check if we can run the forward recurrence + if (m > 12307) + mask1 = 1.00000005d0 <= (chi < 1.005); + elseif (m > 4380) + mask_forward = 1.0000005d0 <= (chi < 1.005); + elseif (m > 1438) + mask_forward = 1.000005d0 <= (chi < 1.005); + elseif (m > 503) + mask_forward = 1.00005d0 <= (chi < 1.005); + elseif (m > 163) + mask_forward = 1.0005d0 <= (chi < 1.005); + else + mask_forward = chi < 1.005; + end + mask_back = ~mask_forward; % backward reccurence mask + + qm = zeros(m+1,size(t,2),size(t,3)); + qmd = zeros(m+1,size(t,2),size(t,3)); + qmdd = zeros(m+1,size(t,2),size(t,3)); + + % run the forward reccurence + if (sum(mask_forward,'all') ~= 0) + [qm(:,mask_forward),qmd(:,mask_forward),qmdd(:,mask_forward)] = forward_reccurence(t(mask_forward),m); + end + + % run the backward reccurence + if (sum(mask_back,'all') ~= 0) + [qm(:,mask_back),qmd(:,mask_back),qmdd(:,mask_back)] = backward_reccurence(t(mask_back),m); + end + +end + +function [qm,qmd,qmdd] = forward_reccurence(t, m) + t = t(:); + chi = t+1; + + % initialize the seeds for the Legendre functions + [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t); + q1d = (-q0 + chi.*q1)./(2*(chi+1).*t); + + qm = zeros(m+1,size(t,1)); + qmd = zeros(m+1,size(t,1)); + qmdd = zeros(m+1,size(t,1)); + + qm(1,:) = q0; + qm(2,:) = q1; + qmd(1,:) = q0d'; + qmd(2,:) = q1d'; + + for i = 1:(m-1) + j = i+1; + qm(2+i,:) = 4*(j-1)/(2*j-1)*chi'.*qm(1+i,:) - (2*j-3)/(2*j-1)*qm(i,:); + qmd(2+i,:) = ((j-1/2)*qm(1+i,:) - (j-1/2)*chi'.*qm(2+i,:))./(-t'.^2-2*t'); + end + + for i = 1:m+1 + j = i-1; + qmdd(i,:) = (-(j-1/2)*(j+1/2)*qm(i,:) + 2*chi'.*qmd(i,:))./(-t'.^2-2*t'); + end +end + +function [qm,qmd,qmdd] = backward_reccurence(t, m) + t = t(:); + chi = t+1; + + % initialize the seeds for the Legendre functions + [q0,~,q0d] = chnk.axissymlap2d.qleg_half(t); + + % Millers Algorithm: + qm = zeros(m+1,size(t,1)); + qmd = zeros(m+1,size(t,1)); + qmdd = zeros(m+1,size(t,1)); + + % 1. run the forward reccurence until it has blown up + % recurrence intialization + fprev = ones(size(t)); + fprevprev = zeros(size(t)); + + % NOTE: + % should maybe wait until second derivatives also blow up + + maxiter = 1000; % max number of extra terms + upbound = 1.0e17; % blow up tolerance + nterms = zeros(size(t)); % number of extra terms + + % list of recurrences that have blown up + already_blowup = zeros(size(t)); + for i = m:maxiter + j = i+1; + f = 4*(j-1)/(2*j-1)*chi.*fprev - (2*j-3)/(2*j-1)*fprevprev; + d = ((j-1/2)*fprevprev - (j-1/2)*chi.*fprev)./(-t.^2-2*t); + + % check which values blew up + f_blowup = (abs(f) >= upbound); + + % check which derivatives blew up + d_blowup = (abs(d) >= upbound); + + % check where both blew up + both_blowup = f_blowup.*d_blowup; + + % check only the new ones that blew up + new_blowup = (both_blowup > already_blowup); + + % set the nterms for the new ones that blew up + nterms = nterms + (i+1)*new_blowup; + + % update the list of ones that have blown up + already_blowup = already_blowup + new_blowup; + + fprevprev = fprev; + fprev = f; + end + + % replace zeros in nterms with maxiter + maxedout = nterms == 0; + nterms = nterms + maxiter*maxedout; + + % 2. run the backward reccurence + % recurrence intialization + fprev = ones(size(t)); + fprevprev = zeros(size(t)); + + % run the backward reccurence from nterms to m + nterm_max = max(nterms, [], 'all'); + for i = 1:(nterm_max-m+1) + j = nterm_max-i+1+1; + + f = 4*(j-1)/(2*j-3)*chi.*fprev - (2*j-1)/(2*j-3)*fprevprev; + + % check if we should start the backwards reccurence yet + should_compute = (j <= nterms); + + % only update if we have started the reccurence for that element + fprevprev = fprev.*should_compute; + fprev = f.*should_compute + fprev.*(~should_compute); + end + + qm(m,:) = fprev'; + qm(m+1,:) = fprevprev'; + + % run the backward reccurence from m to 1 + for i = 1:(m-1) + j = m-1-i+1+1; + qm(j-1,:) = 4*(j-1)/(2*j-3)*chi'.*qm(j,:) - (2*j-1)/(2*j-3)*qm(j+1,:); + end + + % 3. compute error and scale solution + ratio = q0'./qm(1,:); + qm = qm.*ratio; + + % 4. compute 1st derivatives + qmd(1,:) = q0d'; + for i = 1:m + qmd(i+1,:) = (-(i-1/2)*chi'.*qm(i+1,:) + (i-1/2)*qm(i,:))./(-t'.^2-2*t'); + end + + % 5. compute 2nd derivatives + for i = 1:m+1 + qmdd(i,:) = (-(i-3/2)*(i-1/2)*qm(i,:) + 2*chi'.*qmd(i,:))./(-t'.^2-2*t'); + end +end \ No newline at end of file diff --git a/chunkie/@kernel/axissymlap2d.m b/chunkie/@kernel/axissymlap2d.m index c6187ce3..1fef9a1f 100644 --- a/chunkie/@kernel/axissymlap2d.m +++ b/chunkie/@kernel/axissymlap2d.m @@ -1,27 +1,5 @@ -function obj = axissymlap2d(type) -%KERNEL.AXISSYMHELM2D Construct the axissymmetric Helmholtz kernel. -% KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) -% constructs the axissymmetric single-layer Helmholtz kernel with -% wavenumber ZK. -% -% KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) -% constructs the axissymmetric double-layer Helmholtz kernel with -% wavenumber ZK. -% -% KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) -% constructs the normal derivative of the axissymmetric single-layer -% Helmholtz kernel with wavenumber ZK. -% -% KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or -% KERNEL.AXISSYMHELM2D('combined', ZK, COEFS) -% constructs the combined-layer axissymmetric Helmholtz kernel with -% wavenumber ZK and parameter COEFS, -% i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + -% COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK). -% -% NOTES: The axissymetric kernels are currently supported only for purely -% real or purely imaginary wave numbers -% +function obj = axissymlap2d(type, m) +%KERNEL.AXISSYMLAP2D Construct the axissymmetric Laplace kernel. % See also CHNK.AXISSYMHELM2D.KERN. if ( nargin < 1 ) @@ -36,25 +14,32 @@ case {'s', 'single'} obj.type = 's'; - obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 's'); - obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 's'); + obj.eval = @(s,t) chnk.axissymlap2d.kern_modal(s, t, [0,0], 's',m); + obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern_modal(s, t, o, 's',m); obj.fmm = []; obj.sing = 'log'; case {'d', 'double'} obj.type = 'd'; - obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'd'); - obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'd'); + obj.eval = @(s,t) chnk.axissymlap2d.kern_modal(s, t, [0,0], 'd',m); + obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern_modal(s, t, o, 'd',m); obj.fmm = []; - obj.sing = 'log'; + obj.sing = 'smooth'; case {'sp', 'sprime'} obj.type = 'sp'; - obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprime'); - obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprime'); + obj.eval = @(s,t) chnk.axissymlap2d.kern_modal(s, t, [0,0], 'sprime',m); + obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern_modal(s, t, o, 'sprime',m); obj.fmm = []; obj.sing = 'log'; + case {'dp', 'dprime'} + obj.type = 'dp'; + obj.eval = @(s,t) chnk.axissymlap2d.kern_modal(s, t, [0,0], 'dprime',m); + obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern_modal(s, t, o, 'dprime',m); + obj.fmm = []; + obj.sing = 'hs'; + otherwise error('Unknown axissym Laplace kernel type ''%s''.', type); diff --git a/devtools/test/lap_axisym_tests/test_axissymlap_green.m b/devtools/test/lap_axisym_tests/test_axissymlap_green.m new file mode 100644 index 00000000..5c9e9081 --- /dev/null +++ b/devtools/test/lap_axisym_tests/test_axissymlap_green.m @@ -0,0 +1,134 @@ +% Test values and gradient of Green's third identity for Laplace eq. when +% the target is on the boundary for the mth mode: +% 1/2 u = D u - S du/dn +% and +% 1/2 u' = D' u - S' du/dn + +% using chunkies hypersingular GGQ because chunkerkerneval only works for +% off/near boundary evaluation + +clearvars; +close all; +format long e; +%addpaths_loc(); + +%% geometry +[chnkr] = get_torus_geometry(); +charge = [1.0;0.0;-2.0]; % source (axisymm data) + +npts = chnkr.npt; % total number of points in discretization +src = chnkr.r(:,:); % generating curve +n_src = chnkr.n(:,:); % normals + +% plot geometry +%plot(chnkr); + +% convert to cylindrical coordinates +charge_cyl = [sqrt(charge(1)^2 + charge(2)^2);atan2(charge(2),charge(1));charge(3)]; + +p_modes = 5; % total number of positive fourier modes +n_modes = 2*p_modes + 1; % total number of fourier modes (must be odd for pos/0/neg) +n_angles = n_modes; % number of angles/rotations +m = 3; % mode to test (should be positive and < p_modes, m=1 => 0th mode) + +%% compute u and du/dn for densities +u = zeros(n_angles,npts); +dudn = zeros(n_angles,npts); +for i=1:n_angles + % get polar/cartesian coordinates + theta = (i-1)*2*pi/n_angles; + + x=src(1,:).*cos(theta); + y=src(1,:).*sin(theta); + z=src(2,:); + pts=[x;y;z]; + + n_x=n_src(1,:).*cos(theta); + n_y=n_src(1,:).*sin(theta); + n_z=n_src(2,:); + n_pts = [n_x;n_y;n_z]; + + rvec = pts - charge; + r = vecnorm(rvec); + r3 = vecnorm(rvec).^3; + + % construct densities + dudn(i,:) = sum(rvec .* n_pts, 1) ./ (4*pi*r3); + u(i,:) = 1./(4*pi*r); +end + +% Reorder FFT output to match +modes = -p_modes:p_modes; + +u_fft = fft(u, n_modes, 1) / n_angles; % FFT (normalized) +dudn_fft = fft(dudn, n_modes, 1) / n_angles; % FFT (normalized) + +u_m = fftshift(u_fft, 1); % puts negative freqs first +dudn_m = fftshift(dudn_fft, 1); % puts negative freqs first + +%% evaluate operators +opts = []; +opts.rcip = false; +opts.forcesmooth = false; +opts.l2scale = false; + +Sdudn = 0; +Spdudn = 0; +Du = 0; +Dpu = 0; +origin = [0,0]; + +% define operators +S = kernel('axissymlaplace','s',m); +Sp = kernel('axissymlaplace','sprime',m); +D = kernel('axissymlaplace','d',m); +Dp = kernel('axissymlaplace','dprime',m); + +% evaluate at target +S_m = chunkermat(chnkr, S, opts); +Sp_m = chunkermat(chnkr, Sp, opts); +D_m = chunkermat(chnkr, D, opts); +Dp_m = chunkermat(chnkr, Dp, opts); + +% the correct index for the mth mode +i = p_modes + m; + +% convolve +Sdudn_m = S_m*dudn_m(i,:)'; +Spdudn_m = Sp_m*dudn_m(i,:)'; +Du_m = D_m*u_m(i,:)'; +Dpu_m = Dp_m*u_m(i,:)'; + +% check if it satisfies greens identity +val_error = norm(0.5*u_m(i,:)' + Sdudn_m - Du_m); +der_error = norm(0.5*dudn_m(i,:)' + Spdudn_m - Dpu_m); + +% display error +disp(['mode:' num2str(modes(i))]); +disp(['green identity error = ' num2str(val_error)]); +disp(['derivative error = ' num2str(der_error)]); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [chnkobj] = get_torus_geometry() + pref = []; + pref.k = 16; % points per chunk + %pref.nchmax = 6; + + cparams = []; + cparams.eps = 1.0e-10; + %cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = 2; + %cparams.nchmin = 4; + + ctr = [3 0]; + narms = 0; + amp = 0.25; + + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); +end diff --git a/devtools/test/lap_axisym_tests/test_axissymlap_neumann.m b/devtools/test/lap_axisym_tests/test_axissymlap_neumann.m new file mode 100644 index 00000000..6dbf5f52 --- /dev/null +++ b/devtools/test/lap_axisym_tests/test_axissymlap_neumann.m @@ -0,0 +1,154 @@ +%{ +- 3D Laplace's equation +- Neumann boundary condition +- Torus boundary +- Single and Double layer potential representation +- mth mode (non-axisymmetric B.C.) +%} + +clearvars; +close all; +format long e; + +%% geometry +% target is where we evaluate the solution +[chnkr,target,charge1,charge2] = get_torus_geometry(); + +npts = chnkr.npt; % total number of points in discretization +src = chnkr.r(:,:); % generating curve +n_src = chnkr.n(:,:); % normals + +% plot geometry +%plot(chnkr); + +p_modes = 3; % number of positive fourier modes +n_modes = 2*p_modes + 1; % number of fourier modes (must be odd for pos/0/neg) +n_angles = n_modes; % number of angles/rotations +strength = 1.0; % strength of charges + +%% discretization +% compute f (boundary condition) +f = zeros(n_angles,npts); +for i=1:n_angles + % get polar/cartesian coordinates + theta = (i-1)*2*pi/n_angles; + + x=src(1,:).*cos(theta); + y=src(1,:).*sin(theta); + z=src(2,:); + pts=[x;y;z]; + + n_x=n_src(1,:).*cos(theta); + n_y=n_src(1,:).*sin(theta); + n_z=n_src(2,:); + n_pts = [n_x;n_y;n_z]; + + % set up neumann boundary data (with both charges) + rvec1 = pts - charge1; + r3_1 = vecnorm(rvec1).^3; + fn1 = sum(rvec1 .* n_pts, 1) ./ (4*pi*r3_1); + + rvec2 = pts - charge2; + r3_2 = vecnorm(rvec2).^3; + fn2 = -sum(rvec2 .* n_pts, 1) ./ (4*pi*r3_2); + + f(i,:) = strength*(fn1 + fn2); +end + +% Reorder FFT output to match +modes = -p_modes:p_modes; +f_fft = fft(f, n_modes, 1) / n_angles; % FFT (normalized) +f_m = fftshift(f_fft, 1); % puts negative freqs first + +%% solve +% solve the integral equation for each fourier mode +opts = []; +opts.rcip = false; +opts.forcesmooth = false; +opts.l2scale = false; + +sigma1_m = zeros(n_modes,npts); % single layer density +sigma2_m = zeros(n_modes,npts); % double layer density +origin = [0,0]; +for i=1:n_modes + m = abs(modes(i))+1; + Sp = kernel('axissymlaplace','sprime',m); + Dp = kernel('axissymlaplace','dprime',m); + + % Build the system matrix + Sp_m = chunkermat(chnkr, Sp, opts) - 0.5*eye(npts); + Dp_m = chunkermat(chnkr, Dp, opts); + + % Enforce zero-mean constraint for compatability condition + Sp_m = Sp_m + onesmat(chnkr); + Dp_m = Dp_m + onesmat(chnkr); + + % Solve the linear system + rhs = f_m(i,:)'; + sigma1_m(i,:) = gmres(Sp_m, rhs, [], 1e-12, npts); + sigma2_m(i,:) = gmres(Dp_m, rhs, [], 1e-12, npts); +end + +%% solution building +opts.verb = false; +opts.quadkgparams = {'RelTol', 1e-10, 'AbsTol', 1.0e-10}; + +% target in cylindrical coordinates (r,theta,z) +target_cyl = [sqrt(target(1)^2 + target(2)^2);atan2(target(2),target(1));target(3)]; +target_new = [target_cyl(1);target_cyl(3)]; + +u1_sol = 0; % single layer solution +u2_sol = 0; % double layer solution +for i=1:n_modes + m = abs(modes(i))+1; + S = kernel('axissymlaplace','s',m); + D = kernel('axissymlaplace','d',m); + + % evaluation of operators + u1_m = chunkerkerneval(chnkr, S, sigma1_m(i,:), target_new, opts); + u2_m = chunkerkerneval(chnkr, D, sigma2_m(i,:), target_new, opts); + + % fourier composition + u1_sol = u1_sol + real(u1_m * exp(1i * modes(i) * target_cyl(2))); + u2_sol = u2_sol + real(u2_m * exp(1i * modes(i) * target_cyl(2))); +end + +% compute the exact solution explicitly +% only unique up to a constant since neumann BC +r1 = norm(target - charge1); +r2 = norm(target - charge2); +u_true = strength*1.0/(4*pi)*(1/r1 - 1/r2); + +% compute the error of 1st/2nd kind integral equation +err1 = norm(u1_sol-u_true) +err2 = norm(u2_sol-u_true) + + + + +%% gemotry functions +function [chnkobj,target,charge1,charge2] = get_torus_geometry() + pref = []; + pref.k = 16; % points per chunk + %pref.nchmax = 2; + + cparams = []; + %cparams.eps = 1.0e-10; + %cparams.nover = 1; + cparams.ifclosed = true; + cparams.ta = 0; + cparams.tb = 2*pi; + cparams.maxchunklen = 2; + %cparams.nchmin = 8; + + ctr = [3 0]; + narms = 0; + amp = 0.25; + + chnkobj = chunkerfunc(@(t) starfish(t, narms, amp, ctr), cparams, pref); + chnkobj = sort(chnkobj); + + target = [3;0.0;-0.7]; + charge1 = [1.0;0.0;3.0]; + charge2 = [1.0;0.0;-3.0]; +end diff --git a/devtools/test/lap_axisym_tests/test_g0funcall_vec.m b/devtools/test/lap_axisym_tests/test_g0funcall_vec.m new file mode 100644 index 00000000..a80389ef --- /dev/null +++ b/devtools/test/lap_axisym_tests/test_g0funcall_vec.m @@ -0,0 +1,98 @@ + +% +% verify that the modeal Green's functions for Laplace are working properly +% + +h = 0.5; + +r = 1.0; +rp = r + h; +dr = -h; % r - r' + +z = 0.5; +zp = z-h; +dz = h; % z - z' + + + +maxm = 10; + +exact = zeros(maxm+1,1); +exact_gdz = zeros(maxm+1,1); +exact_gdr = zeros(maxm+1,1); +exact_gdrp = zeros(maxm+1,1); + +exact_gdzz = zeros(maxm+1,1); +exact_gdrz = zeros(maxm+1,1); +exact_gdrpz = zeros(maxm+1,1); +exact_gdrpr = zeros(maxm+1,1); + +zk = 0; +nq = 2000; +hhh = 0.0001; + +% evaluate the kernel using brute force trapezoidal integration, and then +% compare with the recurrence relation code. Note: The "exact" derivatives are +% estimated using finite difference, so it is expected that the error is O(h^2) + +for m = 1:(maxm+1) + dint = chnk.axissymlap2d.gkm_brute(r, rp, z, zp, zk, m-1, nq); + exact(m) = dint*4*pi*pi*rp; + + % calculate z derivative using finite differences + dint2 = chnk.axissymlap2d.gkm_brute(r, rp, z+hhh, zp, zk, m-1, nq); + dint1 = chnk.axissymlap2d.gkm_brute(r, rp, z-hhh, zp, zk, m-1, nq); + exact_gdz(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp; + + % gdzz + exact_gdzz(m) = (dint2 - 2*dint + dint1)/(hhh^2)*4*pi*pi*rp; + + % calculate r derivative using finite differences + dint2 = chnk.axissymlap2d.gkm_brute(r+hhh, rp, z, zp, zk, m-1, nq); + dint1 = chnk.axissymlap2d.gkm_brute(r-hhh, rp, z, zp, zk, m-1, nq); + exact_gdr(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp; + + % calculate rp derivative using finite differences + dint2 = chnk.axissymlap2d.gkm_brute(r, rp+hhh, z, zp, zk, m-1, nq); + dint1 = chnk.axissymlap2d.gkm_brute(r, rp-hhh, z, zp, zk, m-1, nq); + exact_gdrp(m) = (dint2-dint1)/(2*hhh)*4*pi*pi*rp; + + % gdrz + dpp = chnk.axissymlap2d.gkm_brute(r+hhh, rp, z+hhh, zp, zk, m-1, nq); % +r, +z + dpm = chnk.axissymlap2d.gkm_brute(r+hhh, rp, z-hhh, zp, zk, m-1, nq); % +r, -z + dmp = chnk.axissymlap2d.gkm_brute(r-hhh, rp, z+hhh, zp, zk, m-1, nq); % -r, +z + dmm = chnk.axissymlap2d.gkm_brute(r-hhh, rp, z-hhh, zp, zk, m-1, nq); % -r, -z + exact_gdrz(m) = (dpp - dpm - dmp + dmm) / (4 * hhh^2)*4*pi*pi*rp; + + % gdrpz + dpp = chnk.axissymlap2d.gkm_brute(r, rp+hhh, z+hhh, zp, zk, m-1, nq); % +r, +z + dpm = chnk.axissymlap2d.gkm_brute(r, rp+hhh, z-hhh, zp, zk, m-1, nq); % +r, -z + dmp = chnk.axissymlap2d.gkm_brute(r, rp-hhh, z+hhh, zp, zk, m-1, nq); % -r, +z + dmm = chnk.axissymlap2d.gkm_brute(r, rp-hhh, z-hhh, zp, zk, m-1, nq); % -r, -z + exact_gdrpz(m) = (dpp - dpm - dmp + dmm) / (4 * hhh^2)*4*pi*pi*rp; + + %gdrpr + dpp = chnk.axissymlap2d.gkm_brute(r+hhh, rp+hhh, z, zp, zk, m-1, nq); % +r, +z + dpm = chnk.axissymlap2d.gkm_brute(r+hhh, rp-hhh, z, zp, zk, m-1, nq); % +r, -z + dmp = chnk.axissymlap2d.gkm_brute(r-hhh, rp+hhh, z, zp, zk, m-1, nq); % -r, +z + dmm = chnk.axissymlap2d.gkm_brute(r-hhh, rp-hhh, z, zp, zk, m-1, nq); % -r, -z + exact_gdrpr(m) = (dpp - dpm - dmp + dmm) / (4 * hhh^2)*4*pi*pi*rp; +end + + +% +% compare with routine +% +[gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = chnk.axissymlap2d.g0funcall_vec(r, rp, dr, z, zp, dz, maxm); + +disp(['gval error = ' num2str(norm(gval-exact))]); +disp(['gdz error = ' num2str(norm(gdz-exact_gdz))]); +disp(['gdr error = ' num2str(norm(gdr-exact_gdr))]); +disp(['gdrp error = ' num2str(norm(gdrp-exact_gdrp))]); +disp(['gdzz error = ' num2str(norm(gdzz-exact_gdzz))]); +disp(['gdrz error = ' num2str(norm(gdrz-exact_gdrz))]); +disp(['gdrpz error = ' num2str(norm(gdrpz-exact_gdrpz))]); +disp(['gdrpr error = ' num2str(norm(gdrpr-exact_gdrpr))]); + + +