Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@eldila
eldila committed Jul 26, 2013
1 parent 0146d41 commit ebaeab3
Show file tree
Hide file tree
Showing 31 changed files with 11,358 additions and 0 deletions.
47 changes: 47 additions & 0 deletions examples/Elliptical Membrane/EllipticalMembrane.m
View file
@@ -0,0 +1,47 @@
% In this example, we look at an elliptical membrane in a fluid that is
% initially at rest. The ellipse should oscillate and eventually settle to
% a circle with a prescribed radius.

% Add PATH reference in order to run solver
addpath('../../solver/Peskin-TwoStep');
addpath('../../solver/utils');

% The number of grid points.
N = 2*round(2.^4);
Nb = 3*N;

% Parameter values.
mu = 1; % Viscosity.
sigma = 1e4; % Spring constant.
rho = 1; % Density.
rmin = 0.2; % Length of semi-minor axis.
rmax = 0.4; % Length of semi-major axis.
req = sqrt(rmin*rmax);
L = 0; % Resting length.

% Time step and final time.
Tfinal = .04;
dt = 5e-5;
NTime = floor(Tfinal./dt)+1;
dt = Tfinal ./ NTime;

% The initial velocity is zero.
IC_U = @(X,Y)zeros(size(X));
IC_V = @(X,Y)zeros(size(Y));

% The membrane is an ellipse.
IC_ChiX = @(S)0.5 + rmax * cos(2*pi*S);
IC_ChiY = @(S)0.5 + rmin * sin(2*pi*S);

% The action function.
Action = @( dx, dy, dt, indexT, X, Y, U, V, chiX, chiY, Fx, Fy)...
PlotMembrane(X, Y, U, V, chiX,chiY,1);

% Do the IB solve.
[X, Y, S, U, V, chiX, chiY] = ...
IBSolver(mu, rho, sigma, L, IC_U, IC_V, IC_ChiX, IC_ChiY,...
N, N, 1, 1, Nb, NTime, Tfinal, Action);

% Remove PATH reference to avoid clutter
rmpath('../../solver/Peskin-TwoStep');
rmpath('../../solver/utils');
49 changes: 49 additions & 0 deletions examples/Fiber Decay Rate/FibreDecayRates.m
View file
@@ -0,0 +1,49 @@
% The purpose of this simulation is to test the decay rate of the maximum
% height of a fibre. We begin with a fibre wrapped around the periodic
% domain with an initially sinusoidal profile.
% The fibre should oscillate with a decaying amplitude. The rate of decay
% can be compute as in Chapter 3 of John Stockie's PhD thesis. Here we
% assume that the resting length of the fibre is zero (L = 0).
%

% Add PATH reference in order to run solver
addpath('../../solver/Peskin-TwoStep');
addpath('../../solver/utils');

% The number of grid points.
N = 64;
Nb = 3*N;

% Parameter values.
L = 0; % Resting length.
mu = 0.1; % Viscosity.
sigma = 1; % Spring constant.
rho = 1; % Density.
A = 0.05; % Initial height of the fibre.

% Time step and final time.
Tfinal = 1;
dt = 5e-3;
NTime = floor(Tfinal/dt);
Tfinal = dt*NTime;

% The initial velocity is zero.
IC_U = @(X,Y)zeros(size(X));
IC_V = @(X,Y)zeros(size(Y));

% The membrane is a sinusoidally perturbed line wrapping around the domain.
IC_ChiX = @(S)S;
IC_ChiY = @(S)0.5 + A*sin(2*pi*S);

% Get the maximum height of the fibre at each time.
Action = @( dx, dy, dt, indexT, X, Y, U, V, chiX, chiY, Fx, Fy)...
PlotMembrane(X, Y, U, V, chiX,chiY,1);

% Do the IB solve.
[X, Y, S, U, V, chiX, chiY] = ...
IBSolver(mu, rho, sigma, L, IC_U, IC_V, IC_ChiX, IC_ChiY, ...
N, N, 1, 1, Nb, NTime, Tfinal, Action);

% Remove PATH reference to avoid clutter
rmpath('../../solver/Peskin-TwoStep');
rmpath('../../solver/utils');
134 changes: 134 additions & 0 deletions solver/Peskin-TwoStep/CalculateForce.m
View file
@@ -0,0 +1,134 @@
function [Fx, Fy, fx, fy] = CalculateForce...
(chiX, chiY, x, y, ds, dx, dy, Nb, Nx, Ny, Lx, Ly, sigma, L)
%
% This function computes the components of the forcing coming from
% interactions between the fluid and the immersed boundary.
% The components of the force are given by
% F(x,y) = int_Gamma f(s) * delta(x - chiX(s)) * delta(y - chiY(s)) * ds
% where
% s parameterises the immersed boundary,
% f(s) is the force density on the boundary,
% chiX, chiY are the x and y positions of points on the membrane,
% delta is a delta function.
%
% The force density is computed using a simple linear spring model with
% resting length L. This leads to a force density of
% f = [\chi_s * ( 1 - L / abs(\chi_s) ) ]_s
% where the subscript s refers to a derivative with respect to s, which
% parameterises the membrane.
%
% INPUTS: chiX A vector of length Nb that gives the current
% x-coordinates of the position of the immersed bounday
% at time-step n-1/2.
% chiY A vector of length Nb that gives the current
% y-coordinates of the position of the immersed bounday
% at time-step n-1/2.
% x,y Vectors of length N that give the x and y values
% occuring in the Eulerian grid.
% ds The step size in hte Lagrangian grid.
% dx The step size in x-direction for the Eulerian grid.
% dy The step size in y-direction for the Eulerian grid.
% Nb The total number of grid points in the Lagrangian
% grid.
% Nx The total number of grid points along x-dimension
% of the Eulerian grid.
% Ny The total number of grid points along y-dimension
% of the Eulerian grid.
% Lx The length of the domain along x-direction.
% Ly The length of the domain along y-direction.
% sigma The spring constant used to compute the force density.
% L The resting length of the membrane.
%
% OUTPUTS: Fx, Fy The x and y components of the forcing resulting from
% the fluid-membrane interactions. These are N X N
% matrices that give the force at each point on the
% grid.
% fx, fy The x and y components of the force density on the
% membrane. These are Nb X 1 vectors.
%
% Authors: Jeffrey Wiens and Brittany Froese, Copyright 2011-2012
%

% Compute the derivatives of the membrane positions with respect to s.
% These are centred differences evaluated at the points (s+1/2), which is
% equivalent to apply a forward difference to chi_s.
% For the forward difference on the periodic domain, the value of s_{Nb+1}
% is equivalent to the value of s_1.
dChiX = PeriodicDisplacement(chiX([2:Nb,1]), chiX, Lx) / ds;
dChiY = PeriodicDisplacement(chiY([2:Nb,1]), chiY, Ly) / ds;

% We also want the magnitude of the derivative.
dChi = sqrt(dChiX.^2 + dChiY.^2);

% Now we compute the tension, which is also given at the points (s+1/2),
% and points in the direction of the tangent vector.
Tx = sigma * dChiX .* (1 - L ./ dChi);
Ty = sigma * dChiY .* (1 - L ./ dChi);

% The force density is computed at the points s by taking a narrow centred
% difference of the tension. This is equivalent to computing a backward
% difference of T_{s+1/2}.
% This assumes a simple linear spring model with resting length L.
% For the backward difference on the periodic domain, the value of s_0
% is equivalent to the value of s_{Nb}.
fx = (Tx - Tx([Nb,1:Nb-1])) / ds;
fy = (Ty - Ty([Nb,1:Nb-1])) / ds;

% Convert the force density into a (sparse) diagonal matrices.
fxds = spdiags(fx * ds,0,Nb,Nb);
fyds = spdiags(fy * ds,0,Nb,Nb);

% Find out how many intervals of the grid the support of the 1D delta
% function will cover.
supp = delta1D();

% Get the indices of the points x_i, y_j where the 1D delta functions
% delta(x - chiX_l) and delta(y - chiY_l) are non-zero.
% These will be Nb X supp and supp X Nb sized matrices, respectively.
indX = IndicesOfNonZeroDelta1D(chiX, Nb, Nx, dx, supp);
indY = IndicesOfNonZeroDelta1D(chiY, Nb, Ny, dy, supp)';

% We also need a matrix of indices that chi can take, which will need to
% have the same size as the other index matrices so we can do the
% subtraction at these points.
% This is Nb X supp.
indChi = repmat((1:Nb)',1,supp);

% Now we will construct sparse matrices to hold the values of the 1D delta
% functions at the points x_i - chiX_l
% where i = 1,..,Nx and l = 1,..,Nb
% stored in an Nb X Nx matrix
% and at the points y_j - chiY_l
% where j = 1,..,Ny and l = 1,..,Nb
% stored in an Ny X Nb matrix
deltaX = sparse(Nb, Nx);
deltaY = sparse(Ny, Nb);

% We only need to evaluate the distance between points in the grid and
% points on the membrane that will contribute to the delta function. There
% is no point in computing all possible values of (x_i - chiX_l) and
% (y_j - chiY_l) when most of these will not contribute.
% The distances x_i - chiX_l are stored in an Nb X supp matrix.
% The distances y_j - chiY_l are stored in a supp X Nb matrix.
distX = PeriodicDistance(x(indX),chiX(indChi),Lx);
distY = PeriodicDistance(y(indY),chiY(indChi'),Ly);

% Now update the delta functions to include any non-zero contributions.
% These still have sizes (Nb X Nx) and (Ny X Nb).
% Each delta function is evaluated only at Nb * 4 points. The other points
% will not contribute.
% These are still stored using a sparse data structure.
deltaX(sub2ind([Nb, Nx], indChi, indX)) = delta1D(distX, dx);
deltaY(sub2ind([Ny, Nb], indY, indChi')) = delta1D(distY, dy);

% Compute the components of the force via the (discrete) line integral
% F(x,y) = int_Gamma f(s) * delta(x - chiX(s)) * delta(y - chiY(s)) * ds.
% We recall that each delta function (deltaX, deltaY) will have only Nb*4
% non-zero entires, while the force density matrices (fxds, fyds) are
% diagonal and have only Nb non-zero entries.
% The multiplication is performed using the sparse matrix structure, so
% this should be much less expensive than multiplying full matrices.
% The result will be an Ny X Nx matrix.
Fx = full(deltaY * fxds * deltaX);
Fy = full(deltaY * fyds * deltaX);

135 changes: 135 additions & 0 deletions solver/Peskin-TwoStep/FluidSolve.m
View file
@@ -0,0 +1,135 @@
function [Uh, Vh, U, V] = FluidSolve(U, V, Fx, Fy, rho, mu, dx, dy, dt, Nx, Ny, Lx, Ly, MatDx, MatDy, MatLap, IndX, IndY)
%
% This function solves the incompressible Navier-Stokes (NS) equations
% rho*u_t = -rho*u*u_x + rho*v*u_y + mu*laplacian(u) - p_x + Fx
% rho*v_t = -rho*u*v_x + rho*v*v_y + mu*laplacian(v) - p_y + Fy
% u_x + v_y = 0.
%
% We are using Peskin's two-step algorithm where the advection terms
% is expressed in skew symmetric form.
%
% INPUTS:
% U, V The x and y components of the fluid velocity.
% These should be input as matrices.
% Fx,Fy The x and y components of the forcing on the
% grid. These are matrices the same size as U and
% V.
% rho The fluid density, which is a scalar.
% mu The diffusion coefficient, which is a scalar.
% dx, dy, dt The spatial step-size. The value for a FULL time-step.
% Nx, Ny The number of grid points along each side of the
% domain.
% Lx, Ly The length of the domain.
% MatDx, MatDy Centered difference approximation to the first derivative
% in the x and y direction with periodic BC.
% MatLap Centered 5-point laplacian stencil
% IndX, IndY Matrix containing indices in x- and y-direction.
% Used as the wave number in FFT solution.
%
% OUTPUTS: Uh, Vh, U, V The updated x and y components of the fluid
% velocity at a full and half time-step.
% These will be in matrix format.
%
% Authors: Jeffrey Wiens and Brittany Froese, Copyright 2011-2012
%

%
% Evolve the Fluid half a time-step
%

% Construct FFT Operator
A_hat = 1 + 4*mu*dt/rho*( (sin(pi*IndX/Nx)/dx).^2 + (sin(pi*IndY/Ny)/dy).^2 );

% Vectorize and construct Diagonal Matrix
UVec = U(:);
VVec = V(:);
UMat = spdiags(UVec, [0], length(UVec), length(UVec));
VMat = spdiags(VVec, [0], length(VVec), length(VVec));
I = speye(length(VVec));

% Construct right hand side in linear system
rhs_u = .5*dt/rho*( - .5*rho*(UMat*MatDx*UVec + VMat*MatDy*UVec) ...
- .5*rho*(MatDx*(UVec.^2) + MatDy*(VVec.*UVec)) ...
+ Fx(:) ...
) + UVec;

rhs_v = .5*dt/rho*( - .5*rho*(UMat*MatDx*VVec + VMat*MatDy*VVec) ...
- .5*rho*(MatDx*(UVec.*VVec) + MatDy*(VVec.^2)) ...
+ Fy(:) ...
) + VVec;

rhs_u = reshape(rhs_u,Ny,Nx);
rhs_v = reshape(rhs_v,Ny,Nx);

% Perform FFT
rhs_u_hat = fft2(rhs_u);
rhs_v_hat = fft2(rhs_v);

% Calculate Fluid Pressure
p_hat = (1.0i/dx*sin(2*pi*IndX/Nx).*rhs_u_hat + 1.0i/dy*sin(2*pi*IndY/Ny).*rhs_v_hat)./...
( dt*(sin(2*pi*IndX/Nx)/dx).^2 + dt*(sin(2*pi*IndY/Ny)/dy).^2 );

% Zero out modes.
p_hat(1,1) = 0;
p_hat(1,Ny/2+1) = 0;
p_hat(Nx/2+1,Ny/2+1) = 0;
p_hat(Nx/2+1,1) = 0;

% Calculate Fluid Velocity
u_hat = (rhs_u_hat + 1.0i*dt/dx*sin(2*pi*IndX/Nx).*p_hat)./A_hat;
v_hat = (rhs_v_hat + 1.0i*dt/dy*sin(2*pi*IndY/Ny).*p_hat)./A_hat;

% 2d inverse of Fourier Coefficients
Uh = real(ifft2(u_hat));
Vh = real(ifft2(v_hat));

%
% Evolve the Fluid a full time-step
%

% Construct FFT Operator
A_hat = 1 + 2*mu*dt/rho*( (sin(pi*IndX/Nx)/dx).^2 + (sin(pi*IndY/Ny)/dy).^2 );

% Vectorize and construct Diagonal Matrix
UhVec = Uh(:);
VhVec = Vh(:);
UhMat = spdiags(UhVec, [0], length(UhVec), length(UhVec));
VhMat = spdiags(VhVec, [0], length(VhVec), length(VhVec));
I = speye(length(VVec));

% Construct right hand side in linear system
rhs_u = dt/rho*( - .5*rho*(UhMat*MatDx*UhVec + VhMat*MatDy*UhVec) ...
- .5*rho*(MatDx*(UhVec.^2) + MatDy*(VhVec.*UhVec)) ...
+ Fx(:) + mu/2*MatLap*UVec ...
) + UVec;

rhs_v = dt/rho*( - .5*rho*(UhMat*MatDx*VhVec + VhMat*MatDy*VhVec) ...
- .5*rho*(MatDx*(UhVec.*VhVec) + MatDy*(VhVec.^2)) ...
+ Fy(:) + mu/2*MatLap*VVec ...
) + VVec;

rhs_u = reshape(rhs_u,Ny,Nx);
rhs_v = reshape(rhs_v,Ny,Nx);

% Perform FFT
rhs_u_hat = fft2(rhs_u);
rhs_v_hat = fft2(rhs_v);

% Calculate Fluid Pressure
p_hat = (1.0i/dx*sin(2*pi*IndX/Nx).*rhs_u_hat + 1.0i/dy*sin(2*pi*IndY/Ny).*rhs_v_hat)./...
( dt*(sin(2*pi*IndX/Nx)/dx).^2 + dt*(sin(2*pi*IndY/Ny)/dy).^2 );

% Zero out modes.
p_hat(1,1) = 0;
p_hat(1,Ny/2+1) = 0;
p_hat(Nx/2+1,Ny/2+1) = 0;
p_hat(Nx/2+1,1) = 0;

% Calculate Fluid Velocity
u_hat = (rhs_u_hat + 1.0i*dt/dx*sin(2*pi*IndX/Nx).*p_hat)./A_hat;
v_hat = (rhs_v_hat + 1.0i*dt/dy*sin(2*pi*IndY/Ny).*p_hat)./A_hat;

% 2d inverse of Fourier Coefficients
U = real(ifft2(u_hat));
V = real(ifft2(v_hat));

0 comments on commit ebaeab3

Please sign in to comment.