Oscillatory stokes kernel and test changes

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)];

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

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);

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

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);