Skip to content

Commit

Permalink
fixing stokes kernels
Browse files Browse the repository at this point in the history
  • Loading branch information
@mrachh
mrachh committed Aug 22, 2024
1 parent 696fdc9 commit 4302a62
Show file tree
Hide file tree
Showing 9 changed files with 385 additions and 39 deletions.
40 changes: 30 additions & 10 deletions chunkie/+chnk/+stok2d/fmm.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,9 @@
% note that targinfo must contain targ.n
% type = 'ctrac' traction corresponding to combined layer,
% note that targinfo must contain targ.n
% type = 'sgrad' gradient corresponding to single layer,
% type = 'dgrad' gradient corresponding to double layer,
% type = 'cgad' gradient corresponding to combined layer,
% sigma - density
% varargin{1} - coef in the combined layer formula, otherwise
% does nothing
Expand All @@ -64,12 +67,12 @@
[m, ns] = size(srcuse.sources);
sigma = reshape(sigma, [m, ns]);
switch lower(type)
case {'svel', 'strac', 'spres'}
srcuse.stoklet = 1/(2*pi)*sigma;
case {'dvel', 'dtrac', 'dpres'}
case {'svel', 'strac', 'spres', 'sgrad'}
srcuse.stoklet = 1/(2*pi*mu)*sigma;
case {'dvel', 'dtrac', 'dpres', 'dgrad'}
srcuse.strslet = -1/(2*pi)*sigma;
srcuse.strsvec = srcinfo.n(1:2,:);
case {'cvel', 'ctrac', 'cpres'}
case {'cvel', 'ctrac', 'cpres', 'cgrad'}
coef = varargin{1};
srcuse.stoklet = 1/(2*pi*mu)*coef(2)*sigma;
srcuse.strslet = -1/(2*pi)*coef(1)*sigma;
Expand All @@ -89,7 +92,7 @@
ifppregtarg = 1;
case {'spres', 'dpres', 'cpres'}
ifppregtarg = 2;
case {'strac', 'dtrac', 'ctrac'}
case {'strac', 'dtrac', 'ctrac', 'sgrad', 'dgrad', 'cgrad'}
ifppregtarg = 3;
end

Expand All @@ -112,23 +115,35 @@
varargout{2} = reshape(U.pretarg, [nt,1])*mu;
end
if nargout > 2
varargout{3} = reshape(U.gradtarg, [m, m, nt])*mu;
varargout{3} = reshape(U.gradtarg, [m, m, nt]);
end
if nargout > 3
str_err = ['CHUNKIE:stok2d:fmm:pgrad',...
'Pressure gradients not currently supported for', ...
'Stokes kernel ''%s''.'];
str_err = ['CHUNKIE:stok2d:fmm: too many output ',...
'arguments for Stokes kernel ''%s''.'];
error(str_err, type);
end
case {'spres', 'dpres', 'cpres'}
varargout{1} = U.pretarg.'*mu;
if nargout > 1
str_err = ['CHUNKIE:stok2d:fmm: too many output ',...
'arguments for Stokes kernel ''%s''.'];
error(str_err, type);
end
case {'sgrad', 'dgrad', 'cgrad'}
varargout{1} = reshape(U.gradtarg, [m*m*nt, 1]);
if nargout > 1
str_err = ['CHUNKIE:stok2d:fmm: too many output ',...
'arguments for Stokes kernel ''%s''.'];
error(str_err, type);
end

case {'strac', 'dtrac', 'ctrac'}
if ( ~isfield(targinfo, 'n') )
error('CHUNKIE:stok2d:fmm:normals', ...
'Targets require normal info when evaluating Stokes kernel ''%s''.', type);
end
du = reshape(U.gradtarg, [m, m, nt]);
dut = permute(u, [2,1,3]);
dut = permute(du, [2,1,3]);
eu = du + dut;
euxx = squeeze(eu(1,1,:));
euxy = squeeze(eu(1,2,:));
Expand All @@ -140,6 +155,11 @@
f(1:2:end) = -p.*ntx + euxx.*ntx + euxy.*nty;
f(2:2:end) = -p.*nty + euxy.*ntx + euyy.*nty;
varargout{1} = mu*f;
if nargout > 1
str_err = ['CHUNKIE:stok2d:fmm: too many output ',...
'arguments for Stokes kernel ''%s''.'];
error(str_err, type);
end
otherwise
error('CHUNKIE:stok2d:fmm:pot', ...
'Potentials not supported for Stokes kernel ''%s''.', type);
Expand Down
141 changes: 128 additions & 13 deletions chunkie/+chnk/+stok2d/kern.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,41 @@
else
varargout = {Kxx, Kxy, Kyy};
end

case 'sgrad'
r4 = r2.^2;
r4inv = 1./r4./(4*pi*mu);
Ku1x1 = rx.*(ry.^2 - rx.^2).*r4inv;
Ku1x2 = ry.*(ry.^2 - rx.^2).*r4inv;

Ku1y1 = (-3.*rx.^2.*ry - ry.^3).*r4inv;
Ku1y2 = rx.*(rx.^2 - ry.^2).*r4inv;

Ku2x1 = ry.*(ry.^2 - rx.^2).*r4inv;
Ku2x2 = (-3.*rx.*ry.^2 - rx.^3).*r4inv;

Ku2y1 = rx.*(rx.^2 - ry.^2).*r4inv;
Ku2y2 = ry.*(rx.^2 - ry.^2).*r4inv;

if ( nargout == 1 )
% Interleave
K = zeros(4*Nt, 2*Ns);
K(1:4:end, 1:2:end) = Ku1x1;
K(1:4:end, 2:2:end) = Ku1x2;
K(2:4:end, 1:2:end) = Ku1y1;
K(2:4:end, 2:2:end) = Ku1y2;

K(3:4:end, 1:2:end) = Ku2x1;
K(3:4:end, 2:2:end) = Ku2x2;
K(4:4:end, 1:2:end) = Ku2y1;
K(4:4:end, 2:2:end) = Ku2y2;

varargout = {K};
else
varargout = {Ku1x1, Ku1x2, Ku1y1, Ku1y2, Ku2x1, Ku2x2, Ku2y1, Ku2y2};
end



case {'dvel', 'd'}

Expand Down Expand Up @@ -112,15 +147,20 @@
r6 = r2.^3;
nx_s = srcinfo.n(1,:);
ny_s = srcinfo.n(2,:);

nx_t = targinfo.n(1,:).';
ny_t = targinfo.n(2,:).';

rn_s = rx.*nx_s + ry.*ny_s;

rn_t = rx.*nx_t + ry.*ny_t;

nn = nx_t.*nx_s + ny_t.*ny_s;

Kxx = mu/pi * ( -8*rx.^2.*rn_t.*rn_s ./ r6 + ...
(rx.*nx_t.*rn_s + rx.^2.*nn + rn_t.*rn_s + nx_s.*rx.*rn_t) ./ r4 + ...
nx_t.*ny_t ./ r2 );
nx_t.*nx_s ./ r2);

Kyy = mu/pi * ( -8*ry.^2.*rn_t.*rn_s ./ r6 + ...
(ry.*ny_t.*rn_s + ry.^2.*nn + rn_t.*rn_s + ny_s.*ry.*rn_t) ./ r4 + ...
ny_t.*ny_s ./ r2 );
Expand All @@ -142,15 +182,56 @@
else
varargout = {Kxx, Kxy, Kyy};
end
case 'dgrad'
r4 = r2.^2;
r6 = r2.^3;
nx = srcinfo.n(1,:);
ny = srcinfo.n(2,:);
rn = rx.*nx + ry.*ny;
overpi = 1/pi;

r4inv = 1./r4;
r6inv = 1./r6;
Ku1x1 = (rn.*rx.*r4inv + rx.*rn.*r4inv + rx.^2.*nx.*r4inv - 4*rx.^3.*rn.*r6inv)*overpi;
Ku1x2 = (rn.*ry.*r4inv + rx.*ry.*nx.*r4inv - 4*rx.^2.*ry.*rn.*r6inv)*overpi;

Ku1y1 = (rx.^2.*ny.*r4inv - 4*rx.^2.*ry.*rn.*r6inv)*overpi;
Ku1y2 = (rx.*rn.*r4inv + rx.*ry.*ny.*r4inv - 4*rx.*ry.^2.*rn.*r6inv)*overpi;

Ku2x1 = (ry.*rn.*r4inv + rx.*ry.*nx.*r4inv - 4*rx.^2.*ry.*rn.*r6inv)*overpi;
Ku2x2 = (ry.^2.*nx.*r4inv - 4*rx.*ry.^2.*rn.*r6inv)*overpi;

Ku2y1 = (rn.*rx.*r4inv + rx.*ry.*ny.*r4inv - 4*rx.*ry.^2.*rn.*r6inv)*overpi;
Ku2y2 = (rn.*ry.*r4inv + ry.*rn.*r4inv + ry.^2.*ny.*r4inv - 4*ry.^3.*rn.*r6inv)*overpi;

if ( nargout == 1 )
% Interleave
K = zeros(4*Nt, 2*Ns);
K(1:4:end, 1:2:end) = Ku1x1;
K(1:4:end, 2:2:end) = Ku1x2;
K(2:4:end, 1:2:end) = Ku1y1;
K(2:4:end, 2:2:end) = Ku1y2;

K(3:4:end, 1:2:end) = Ku2x1;
K(3:4:end, 2:2:end) = Ku2x2;
K(4:4:end, 1:2:end) = Ku2y1;
K(4:4:end, 2:2:end) = Ku2y2;


varargout = {K};
else
varargout = {Ku1x1, Ku1x2, Ku1y1, Ku1y2, Ku2x1, Ku2x2, Ku2y1, Ku2y2};
end


case 'c'

eta = varargin{1};
coefs = varargin{1};
[Sxx, Sxy, Syy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 's');
[Dxx, Dxy, Dyy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 'd');
Kxx = Dxx + eta*Sxx;
Kyy = Dyy + eta*Syy;
Kxy = Dxy + eta*Sxy;
Kxx = coefs(1)*Dxx + coefs(2)*Sxx;
Kyy = coefs(1)*Dyy + coefs(2)*Syy;
Kxy = coefs(1)*Dxy + coefs(2)*Sxy;

if ( nargout == 1 )
% Interleave
Expand All @@ -165,11 +246,11 @@
end
case 'cpres'

eta = varargin{1};
coefs = varargin{1};
[Sx, Sy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 'spres');
[Dx, Dy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 'dpres');
Kx = Dx + eta*Sx;
Ky = Dy + eta*Sy;
Kx = coefs(1)*Dx + coefs(2)*Sx;
Ky = coefs(1)*Dy + coefs(2)*Sy;

if ( nargout == 1 )
% Interleave
Expand All @@ -183,12 +264,12 @@

case 'ctrac'

eta = varargin{1};
coefs = varargin{1};
[Sxx, Sxy, Syy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 'strac');
[Dxx, Dxy, Dyy] = chnk.stok2d.kern(mu, srcinfo, targinfo, 'dtrac');
Kxx = Dxx + eta*Sxx;
Kxy = Dxy + eta*Sxy;
Kyy = Dyy + eta*Syy;
Kxx = coefs(1)*Dxx + coefs(2)*Sxx;
Kxy = coefs(1)*Dxy + coefs(2)*Sxy;
Kyy = coefs(1)*Dyy + coefs(2)*Syy;

if ( nargout == 1 )
% Interleave
Expand All @@ -201,8 +282,42 @@
else
varargout = {Kxx, Kxy, Kyy};
end
case 'cgrad'
coefs = varargin{1};
[Su1x1, Su1x2, Su1y1, Su1y2, Su2x1, Su2x2, Su2y1, Su2y2] = ...
chnk.stok2d.kern(mu, srcinfo, targinfo, 'sgrad');
[Du1x1, Du1x2, Du1y1, Du1y2, Du2x1, Du2x2, Du2y1, Du2y2] = ...
chnk.stok2d.kern(mu, srcinfo, targinfo, 'sgrad');


Ku1x1 = coefs(1)*Du1x1 + coefs(2)*Su1x1;
Ku1x2 = coefs(1)*Du1x2 + coefs(2)*Su1x2;

Ku1y1 = coefs(1)*Du1y1 + coefs(2)*Su1y1;
Ku1y2 = coefs(1)*Du1y2 + coefs(2)*Su1y2;

Ku2x1 = coefs(1)*Du2x1 + coefs(2)*Su2x1;
Ku2x2 = coefs(1)*Du2x2 + coefs(2)*Su2x2;

Ku2y1 = coefs(1)*Du2y1 + coefs(2)*Su2y1;
Ku2y2 = coefs(1)*Du2y2 + coefs(2)*Su2y2;

if ( nargout == 1 )
% Interleave
K = zeros(4*Nt, 2*Ns);
K(1:4:end, 1:2:end) = Ku1x1;
K(1:4:end, 2:2:end) = Ku1x2;
K(2:4:end, 1:2:end) = Ku1y1;
K(2:4:end, 2:2:end) = Ku1y2;

K(3:4:end, 1:2:end) = Ku2x1;
K(3:4:end, 2:2:end) = Ku2x2;
K(4:4:end, 1:2:end) = Ku2y1;
K(4:4:end, 2:2:end) = Ku2y2;

varargout = {K};
else
varargout = {Ku1x1, Ku1x2, Ku1y1, Ku1y2, Ku2x1, Ku2x2, Ku2y1, Ku2y2};
end

otherwise
error('Unknown Stokes kernel type ''%s''.', type);
Expand Down
8 changes: 6 additions & 2 deletions chunkie/@kernel/kernel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,12 @@
% 'cp'
% 'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp'
% 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt'
% 'stokes' ('stok', 's') 'svel', 'spres', 'strac',
% 'dvel', 'dpres', 'dtrac'
% 'stokes' ('stok', 's') 'svel', 'spres',
% 'strac', 'sgrad'
% 'dvel', 'dpres',
% 'dtrac', 'dgrad'
% 'cvel', 'cpres',
% 'ctrac', 'cgrad'
% 'zeros' ('zero','z')
% 'axis sym helmholtz' 's' 'd' 'sp' 'c'
% ('axissymh', 'axissymhelm')
Expand Down
40 changes: 37 additions & 3 deletions chunkie/@kernel/stok2d.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,10 @@
% constructs the single-layer Stokes kernel for traction with viscosity
% MU.
%
% KERNEL.STOK2D('sgrad', MU) or KERNEL.ELAST2D('sgradient', MU)
% constructs the single-layer Stokes kernel for velocity gradient with
% viscosity MU.
%
% KERNEL.STOK2D('dvel', MU) or KERNEL.ELAST2D('dvelocity', MU)
% constructs the double-layer Stokes kernel for velocity with viscosity
% MU. KERNEL.STOK2D('d', MU) and KERNEL.STOK2D('double', MU) are
Expand All @@ -26,7 +30,10 @@
% constructs the double-layer Stokes kernel for traction with viscosity
% MU.
%
% %
% KERNEL.STOK2D('dgrad', MU) or KERNEL.ELAST2D('dgradient', MU)
% constructs the double-layer Stokes kernel for velocity gradient with
% viscosity MU.
%
% KERNEL.STOK2D('cvel', MU) or KERNEL.ELAST2D('cvelocity', MU)
% constructs the combined-layer Stokes kernel for velocity with viscosity
% MU. KERNEL.STOK2D('c', MU) and KERNEL.STOK2D('comb', MU) are
Expand All @@ -40,6 +47,10 @@
% constructs the combined-layer Stokes kernel for traction with viscosity
% MU.
%
% KERNEL.STOK2D('cgrad', MU) or KERNEL.ELAST2D('cgradient', MU)
% constructs the combined-layer Stokes kernel for velocity gradient with
% viscosity MU.
%
%

% See also CHNK.STOK2D.KERN.
Expand Down Expand Up @@ -69,19 +80,26 @@
obj.sing = 'log';

case {'spres', 'spressure'}
obj.type = 'svel';
obj.type = 'spres';
obj.eval = @(s,t) chnk.stok2d.kern(mu, s, t, 'spres');
obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'spres', sigma);
obj.opdims = [1, 2];
obj.sing = 'pv';

case {'strac', 'straction'}
obj.type = 'svel';
obj.type = 'strac';
obj.eval = @(s,t) chnk.stok2d.kern(mu, s, t, 'strac');
obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'strac', sigma);
obj.opdims = [2, 2];
obj.sing = 'smooth';

case {'sgrad', 'sgradient'}
obj.type = 'sgrad';
obj.eval = @(s,t) chnk.stok2d.kern(mu, s, t, 'sgrad');
obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'sgrad', sigma);
obj.opdims = [4, 2];
obj.sing = 'pv';

case {'dvel', 'dvelocity', 'd', 'double'}
obj.type = 'dvel';
obj.eval = @(s,t) chnk.stok2d.kern(mu, s, t, 'dvel');
Expand All @@ -103,6 +121,14 @@
obj.opdims = [2, 2];
obj.sing = 'hs';

case {'dgrad', 'dgradient'}
obj.type = 'dgrad';
obj.eval = @(s,t) chnk.stok2d.kern(mu, s, t, 'dgrad');
obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'dgrad', sigma);
obj.opdims = [4, 2];
obj.sing = 'hs';


case {'cvel', 'cvelocity', 'c', 'comb'}
if ( nargin < 2 )
warning('Missing combined layer parameter coefs. Defaulting to 1.');
Expand Down Expand Up @@ -142,6 +168,14 @@
obj.opdims = [2, 2];
obj.sing = 'hs';

case {'cgrad', 'cgradient'}
obj.type = 'cgrad';
obj.eval = @(s,t) coefs(1)*chnk.stok2d.kern(mu, s, t, 'dgrad') + ...
coefs(2)*chnk.stok2d.kern(mu, s, t, 'sgrad');
obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'cgrad', sigma, coefs);
obj.opdims = [4, 2];
obj.sing = 'hs';

otherwise
error('Unknown Stokes kernel type ''%s''.', type);

Expand Down
Loading

0 comments on commit 4302a62

Please sign in to comment.