-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfluidSolver.m
122 lines (86 loc) · 3.23 KB
/
fluidSolver.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
function [u, p] = fluidSolver(u0, f, rho, mu, dt, h, N);
%
% [u] = fluidSolver(u0, f, rho, mu, dt, h, N);
%
% Takes one times step of the fluid solver from Dr. Peskin's notes
%
% Returns:
% u = solution at time dt after u0 solution
% p = pressure at time dt after u0 solution
%
% Input:
% u0 = current solution
% f = force at current time + dt/2
% rho = density
% mu = viscosity
% dt = time step
% h = spatial length
% N = number of spatial mesh points in each direction (NxN or NxN
% grid)
%
%
% License: This code is free to use for any purposes, provided
% any publications resulting from the use of this code
% reference the original code/author.
%
% Author: Samuel Isaacson ([email protected])
% Date: 11/2007
%
% Please notify the author of any bugs, and contribute any
% modifications or bug fixes back to the original author.
%
% Disclaimer:
% This code is provided as is. The author takes no responsibility
% for its results or effects.
persistent D0x D0y L D0xT D0yT LT D0T zeroIdxs isInit;
if( size(isInit) == 0 )
isInit = 1;
% OPERATORS:
[D0x, D0y] = D02DPeriodic(N,h);
L = lap2DPeriodic(N,h);
[D0xT, D0yT] = D02DPeriodicFT(N,h);
LT = lap2DPeriodicFT(N,h);
D0T = [D0xT D0yT];
% find zero entries in D0T when dotted, replace
zeroIdxs = find(dot(D0T,D0T,2) < 2*eps);
D0T(zeroIdxs,:) = 1;
end
M = N * N;
u = u0;
% 1 - mu*dt/(2*rho) * L
ILT = 1 - mu*dt / (2*rho) * LT;
% Calulate S(u)u
tmp = zeros(size(u));
tmp(:,1) = (u(:,1) .* (D0x * u(:,1)) + u(:,2) .* (D0y * u(:,1)) + D0x * (u(:,1) .* u(:,1)) + D0y * (u(:,2) .* u(:,1)) ) / 2;
tmp(:,2) = (u(:,1) .* (D0x * u(:,2)) + u(:,2) .* (D0y * u(:,2)) + D0x * (u(:,1) .* u(:,2)) + D0y * (u(:,2) .* u(:,2)) ) / 2;
% Calculate Transformed rhs
tmp = u - (dt/2)*tmp + (dt / (2*rho)) * f;
rhs = fft2D(tmp, N, M);
% Solve for Fourier Space Solutions at t + dt/2:
pT = 2 * rho / dt * dot( D0T, rhs, 2 ) ./ dot( D0T, D0T, 2 );
pT( zeroIdxs ) = 0;
uT = ( rhs - (dt/(2*rho)) * [ D0xT .* pT D0yT .* pT ] ) ./ [ILT ILT];
% Transform Back to Real Space
u = ifft2D(uT, N, M);
u = real(u);
% Calculate S(u)u
tmp(:,1) = (u(:,1) .* (D0x * u(:,1)) + u(:,2) .* (D0y * u(:,1)) + D0x * (u(:,1) .* u(:,1)) + D0y * (u(:,2) .* u(:,1)) ) / 2;
tmp(:,2) = (u(:,1) .* (D0x * u(:,2)) + u(:,2) .* (D0y * u(:,2)) + D0x * (u(:,1) .* u(:,2)) + D0y * (u(:,2) .* u(:,2)) ) / 2;
% Calculate Transformed rhs
tmp = u0 - dt*tmp + mu*dt / (2*rho) * (L * u0) + dt/rho * f;
rhs = fft2D(tmp, N, M);
% Solve for Fourier Space Solutions at t + dt:
pT = rho / dt * dot( D0T, rhs, 2 ) ./ dot( D0T, D0T, 2 );
pT( zeroIdxs ) = 0;
uT = ( rhs - (dt/(rho)) * [ D0xT .* pT D0yT .* pT ] ) ./ [ILT ILT];
% Transform Back to Real Space
u = ifft2D(uT, N, M);
u = real(u);
p = ifft2D1(pT, N, M);
p = real(p);
function [uT] = fft2D(u, N, M);
uT = [ reshape(fft2( reshape(u(:,1),N,N) ), M, 1) reshape(fft2( reshape(u(:,2),N,N) ), M, 1) ];
function [u] = ifft2D(uT, N, M);
u = [ reshape(ifft2( reshape(uT(:,1),N,N) ), M, 1) reshape(ifft2( reshape(uT(:,2),N,N) ), M, 1) ];
function [u] = ifft2D1(uT, N, M);
u = reshape(ifft2( reshape(uT,N,N) ), M, 1);