diff --git a/chunkie/+chnk/+ostok2d/green.m b/chunkie/+chnk/+ostok2d/green.m new file mode 100644 index 00000000..e06a619b --- /dev/null +++ b/chunkie/+chnk/+ostok2d/green.m @@ -0,0 +1,9 @@ +%%%%%%% OSCILLATORY STOKES GREEN FUNCTION %%%%%%%%% +function val = green(k,src,targ) + +% Green function corresponding to = −I∆G_BH + ∇ ⊗ ∇G_BH + +% [~,~, hess] = chnk.obihar2d.green(k,src,targ); +[~,~, hess] = chnk.flex2d.helmdiffgreen(k,src,targ); +hess = hess/(k*k); +val = [-hess(:,:,3) hess(:,:,2); hess(:,:,2) -hess(:,:,1)]; diff --git a/chunkie/+chnk/+ostok2d/kern.m b/chunkie/+chnk/+ostok2d/kern.m new file mode 100644 index 00000000..af571776 --- /dev/null +++ b/chunkie/+chnk/+ostok2d/kern.m @@ -0,0 +1,101 @@ +%%%%%%%%% KERNEL FOR LAYER POTENTIALS FOR OSCILLATORY STOKES %%%%%% +function varargout = kern(k,srcinfo,targinfo,type,varargin) +% for type 's' submat gives the usual GREEN's function for oscillatory +% STOKES operator +% [] +% for type 'd' submat gives the Double Layer for oscillatory +% STOKES operator +% [] +% for type 'c' submat gives the Combined field for oscillatory +% STOKES operator + + +src = srcinfo.r(:,:); +targ = targinfo.r(:,:); +[~, nt] = size(targ); +[~, ns] = size(src); + +% Single Layer + +if strcmpi(type,'s') + + [~,~, hess] = chnk.flex2d.helmdiffgreen(k,src,targ); + hess = hess/(k*k); + Kxx = -hess(:,:,3); + Kxy = hess(:,:,2); + Kyy = -hess(:,:,1); + if nargout == 1 + K = zeros(2*nt, 2*ns); + K(1:2:end, 1:2:end) = Kxx; + K(1:2:end, 2:2:end) = Kxy; + K(2:2:end, 1:2:end) = Kxy; + K(2:2:end, 2:2:end) = Kyy; + varargout = {K}; + else + varargout = {Kxx, Kxy, Kyy}; + end +end + +% double layer + +if strcmpi(type,'d') + +[~,~,~,der3] = chnk.flex2d.helmdiffgreen(k,src,targ); +der3 = der3/(k*k); +[~, nt] = size(targinfo.r(:,:)); +[~, ns] = size(srcinfo.r(:,:)); + + +nx = repmat(srcinfo.n(1,:),nt,1); +ny = repmat(srcinfo.n(2,:),nt,1); + +[~, grad] = chnk.lap2d.green(src, targ); + +Kxx = -((grad(:,:,1) - 2*der3(:,:,3)).*nx ... + + (der3(:,:,2) - der3(:,:,4)).*ny); + +Kxy = -((grad(:,:,2) + 2*der3(:,:,2)).*nx ... + + (der3(:,:,3) - der3(:,:,1)).*ny); + +Kyx = -((der3(:,:,2) - der3(:,:,4)).*nx ... + + (grad(:,:,1) + 2*der3(:,:,3)).*ny); + +Kyy = -((der3(:,:,3) - der3(:,:,1)).*nx ... + + (grad(:,:,2) - 2*der3(:,:,2)).*ny); + + + if nargout == 1 + K = zeros(2*nt, 2*ns); + K(1:2:end, 1:2:end) = Kxx; + K(1:2:end, 2:2:end) = Kyx; + K(2:2:end, 1:2:end) = Kxy; + K(2:2:end, 2:2:end) = Kyy; + varargout = {K}; + else + varargout = {Kxx, Kxy, Kyx, Kyy}; + end +end + +if strcmpi(type,'c') + coefs = varargin{1}; + [Sxx, Sxy, Syy] = chnk.ostok2d.kern(k, srcinfo, targinfo, 's'); + [Dxx, Dxy, Dyx, Dyy] = chnk.ostok2d.kern(k, srcinfo,targinfo, 'd'); + Kxx = coefs(1)*Dxx + coefs(2)*Sxx; + Kyy = coefs(1)*Dyy + coefs(2)*Syy; + Kxy = coefs(1)*Dxy + coefs(2)*Sxy; + Kyx = coefs(1)*Dyx + coefs(2)*Sxy; + + if ( nargout == 1 ) + % Interleave + K = zeros(2*nt, 2*ns); + K(1:2:end, 1:2:end) = Kxx; + K(1:2:end, 2:2:end) = Kyx; + K(2:2:end, 1:2:end) = Kxy; + K(2:2:end, 2:2:end) = Kyy; + varargout = {K}; + else + varargout = {Kxx, Kxy, Kyx, Kyy}; + end +end + +end \ No newline at end of file diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index 39efb846..db1cb5d9 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -9,10 +9,9 @@ % ---- ---- % 'laplace' ('lap', 'l') 's', 'd', 'sp', 'c' % 'helmholtz' ('helm', 'h') 's', 'd', 'sp', 'dp', 'c' -% 'cp' -% 'helmholtz1d' ('helm1d', 'h1d') 's' -% 'helmholtz difference' 's', 'd', 'sp', 'dp' -% ('helmdiff', 'hdiff', 'helm_diff') +% 'cp', 'fd_s', 'fd_d' +% 'fd_sprime', 'fd_dprime' +% 'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp' % 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt' % 'stokes' ('stok', 's') 'svel', 'spres', % 'strac', 'sgrad' @@ -24,9 +23,10 @@ % 'axis sym helmholtz' 's' 'd' 'sp' 'c' % ('axissymh', 'axissymhelm') % 'axis sym helmholtz difference' 's' 'd' 'sp' 'dp' -% ('axissymhdiff', 'axissymhelmdiff', 'axissymhelm_diff') +% ('axissymhdiff', 'axissymhelmdiff') % 'quasiperiodic helmholtz' 's', 'd', 'sp', 'dp', 'c' % ('helmquas', 'hq') 'cp' +% 'ostokes' ('ostok', 'os') 's', 'd', 'c' % The types may also be written in longer form, e.g. 'single', 'double', % 'sprime', 'combined', 'svelocity', 'spressure', 'straction', % 'dvelocity', 'dpressure', 'dtraction'. @@ -101,12 +101,12 @@ obj = kernel.lap2d(varargin{:}); case {'helmholtz', 'helm', 'h'} obj = kernel.helm2d(varargin{:}); - case {'helmholtz1d', 'helm1d', 'h1d'} - obj = kernel.helm1d(varargin{:}); - case {'helmholtz difference', 'helmdiff', 'hdiff', 'helm_diff'} + case {'helmholtz difference', 'helmdiff', 'hdiff'} obj = kernel.helm2ddiff(varargin{:}); case {'stokes', 'stok', 's'} obj = kernel.stok2d(varargin{:}); + case {'ostokes', 'ostok', 'os'} + obj = kernel.ostok2d(varargin{:}); case {'elasticity', 'elast', 'e'} obj = kernel.elast2d(varargin{:}); case {'zeros', 'zero', 'z'} @@ -116,7 +116,7 @@ case {'axis sym helmholtz', 'axissymh', 'axissymhelm'} obj = kernel.axissymhelm2d(varargin{:}); case {'axis sym helmholtz difference', 'axissymhdiff' ... - 'axissymhelmdiff', 'axissymhelm_diff'} + 'axissymhelmdiff'} obj = kernel.axissymhelm2ddiff(varargin{:}); case {'quasiperiodic helmholtz', 'helmquas', 'hq'} obj = kernel.helm2dquas(varargin{:}); @@ -160,9 +160,9 @@ obj = lap2d(varargin); obj = helm2d(varargin); - obj = helm1d(varargin); obj = helm2ddiff(varargin); obj = stok2d(varargin); + obj = ostok2d(varargin); obj = elast2d(varargin); obj = axissymhelm2d(varargin); obj = axissymhelm2ddiff(varargin); diff --git a/chunkie/@kernel/ostok2d.m b/chunkie/@kernel/ostok2d.m new file mode 100644 index 00000000..728f1883 --- /dev/null +++ b/chunkie/@kernel/ostok2d.m @@ -0,0 +1,67 @@ +function obj = ostok2d(type, zk, coefs) +%KERNEL.OSTOK2D Construct the Stokes kernel. +% KERNEL.OSTOK2D('dvel', zk) +% constructs the single-layer Oscillatory Stokes kernel for velocity with +% viscosity zk. KERNEL.OSTOK2D('s', zk) and KERNEL.OSTOK2D('single', zk) +% are equivalent. +% KERNEL.OSTOK2D('dvel', zk) +% constructs the Double-layer Oscillatory Stokes kernel for velocity with +% viscosity zk. KERNEL.OSTOK2D('d', zk) and KERNEL.OSTOK2D('double', zk) +% are equivalent. +% KERNEL.OSTOK2D('cvel', zk, coefs) +% constructs the Combined field Oscillatory Stokes kernel for velocity +% with viscosity zk and coefs. KERNEL.OSTOK2D('c', zk, coefs) and +% KERNEL.OSTOK2D('comb', zk, coefs) are equivalent. +% +% See also CHNK.OSTOK2D.KERN. + +% author: Kshitij Sinha + +if ( nargin < 1 ) + error('Missing Oscillatory Stokes kernel type.'); +end + +if ( nargin < 2 ) + error('Missing Oscillatory Stokes wave number k.'); +end + +obj = kernel(); +obj.name = 'ostokes'; +obj.params.zk = zk; + +switch lower(type) + case {'svel', 'svelocity', 's', 'single'} + obj.type = 'svel'; + obj.eval = @(s,t) chnk.ostok2d.kern(zk, s, t, 's'); + obj.fmm = []; + obj.opdims = [2, 2]; + obj.sing = 'log'; + case {'dvel', 'dvelocity', 'd', 'double'} + obj.type = 'dvel'; + obj.eval = @(s,t) chnk.ostok2d.kern(zk, s, t, 'd'); + obj.fmm = []; + obj.opdims = [2, 2]; + obj.sing = 'log'; + case {'cvel', 'cvelocity', 'c', 'combined'} + if ( nargin < 2 ) + warning(['Missing combined layer parameter coefs. ' ... + 'Defaulting to [1 1].']); + coefs = ones(2,1); + end + obj.type = 'cvel'; + obj.params.coefs = coefs; + obj.eval = @(s,t) coefs(1)*chnk.ostok2d.kern(zk, s, t, 'd') + ... + coefs(2)*chnk.ostok2d.kern(zk, s, t, 's'); + obj.fmm = []; + obj.opdims = [2, 2]; + obj.sing = 'log'; + otherwise + error('Unknown Oscillatory Stokes kernel type ''%s''.', type); +end + +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + +end diff --git a/devtools/test/chunkermat_ostok2dTest.m b/devtools/test/chunkermat_ostok2dTest.m new file mode 100644 index 00000000..0ca527ef --- /dev/null +++ b/devtools/test/chunkermat_ostok2dTest.m @@ -0,0 +1,152 @@ +% chunkermat_ostok2dTest0(); + +% function chunkermat_ostok2dTest0() + +% CHUNKERMAT_OSTOK2DTEST +% +% test the matrix builder and do a basic solve +clear; + +iseed = 8675309; +rng(iseed); + +cparams = []; +cparams.eps = 1.0e-10; +cparams.nover = 1; +pref = []; +pref.k = 20; +narms = 3; +amp = 0.25; +start = tic; +chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); +t1 = toc(start); + +fprintf('%5.2e s : time to build geo\n',t1); + +% sources + +ns = 10; +ts = 0.0+2*pi*rand(ns,1); +sources = starfish(ts,narms,amp); +sources = 3.0*sources; +strengths = randn(2*ns,1); +sources_n = rand(2,ns); + +% targets + +nt = 100; +ts = 0.0+2*pi*rand(nt,1); +targets = starfish(ts,narms,amp); +targets = targets.*repmat(rand(1,nt),2,1)*0.8; + +plot(chnkr, 'r.'); hold on; +plot(targets(1,:), targets(2,:), 'kx') +hold on; +plot(sources(1,:), sources(2,:), 'bo') +axis equal + + + + +targs = chnkr.r; targs = reshape(targs,2,chnkr.k*chnkr.nch); +targstau = tangents(chnkr); +targstau = reshape(targstau,2,chnkr.k*chnkr.nch); + +plot(chnkr, 'r.'); hold on; +plot(targets(1,:), targets(2,:), 'kx'); +plot(sources(1,:), sources(2,:), 'bo');hold on; + +t1 = toc(start); + +fprintf('%5.2e s : time to build geo\n',t1) + + +zk = 0.3; +kerns = kernel('ostok', 's', zk); +kernd = kernel('ostok', 'd', zk); + +% eval u on bdry + + +srcinfo = []; +srcinfo.r = sources; +srcinfo.n = sources_n; +kernmats = kerns.eval(srcinfo, chnkr); +ubdry = kernmats*strengths; + +% eval u at targets + +targinfo = []; +targinfo.r = targets; +targets_n = rand(2, nt); +targets_n = targets_n./sqrt(targets_n(1,:).^2+targets_n(2,:).^2); +targinfo.n = targets_n; +kernmatstargs = kerns.eval(srcinfo, targinfo); + +utarg = kernmatstargs*strengths; + + +% solve + +fkernd = kernel('ostok', 'd', zk); +fkerns = kernel('ostok', 's', zk); +coefs = rand(1,2); +fkernc = kernel('ostok', 'c', zk, coefs); + + +start = tic; +S = chunkermat(chnkr, fkerns); +D = chunkermat(chnkr, fkernd); +C = chunkermat(chnkr, fkernc); +t1 = toc(start); + +fprintf('%5.2e s : time to assemble matrix\n',t1) + +sysd = -0.5*eye(size(D,1)) + D; +sysd = sysd + normonesmat(chnkr)/sum(chnkr.wts(:)); + +syss = S; + +sysc = -0.5*coefs(1)*eye(size(D,1)) + C; +sysc = sysc + coefs(1)*normonesmat(chnkr)/sum(chnkr.wts(:)); + + +rhs = ubdry; +rhs = rhs(:); + +start = tic; +sold = gmres(sysd,rhs,[],1e-12,1000); +sols = gmres(syss,rhs,[],1e-12,1000); +solc = gmres(sysc,rhs,[],1e-12,1000); + + +t1 = toc(start); + +fprintf('%5.2e s : time for dense gmres\n',t1); + +% evaluate at targets and compare + +opts.usesmooth=false; +opts.verb=false; + +% SINGLE LAYER TEST + +fkerns = kernel('ostok', 's', zk); +Ssol = chunkerkerneval(chnkr, fkerns, sols, targets, opts); +relerr = norm(utarg-Ssol,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); +fprintf('S_ relative frobenius error %5.2e\n', relerr); + + +% DOUBLE LAYER TEST + +fkernd = kernel('ostok', 'd', zk); +Dsol = chunkerkerneval(chnkr, fkernd, sold, targets, opts); +relerr = norm(utarg-Dsol,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); +fprintf('D_relative frobenius error %5.2e\n', relerr); + +% COMBINED FIELD TEST + +fkernc = kernel('ostok', 'c', zk, coefs); +Dsol = chunkerkerneval(chnkr, fkernc, solc, targets, opts); +relerr = norm(utarg-Dsol,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); +fprintf('C_relative frobenius error %5.2e\n', relerr); \ No newline at end of file