-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRBCMotionInUnboundedFlows.m
187 lines (162 loc) · 6.74 KB
/
RBCMotionInUnboundedFlows.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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
%% Spectral Galerkin BIE implementation for an RBC motion and deformation
%% in an unbounded domain under the ambient flow field
%%
clear all; close all; clc;
addpath(genpath('../SpectralBoundaryIntegralMethod.m'))
verbose_Plot = false;
Starttime = tic;
%% Input the model and parameters for the analysis from Models folder
LoadElasticRBC_Shear_N16
% LoadElasticRBC_Parabolic_N16
% LoadMemViscosityRBC_Shear_N16
% LoadMemViscosityRBC_Parabolic_N16
%% Create grid for representation: [N+1 x 2*N+1]
[nlat, nlon, thet, phi, weight] = GridOnSphere(N);
%% Create grid for integration points on unit sphere:
[NGthet, NGphi, eta, wg] = setup_integration_grid(N, NGSphere);
%% Create geometry of RBC: undeformed geometry of an RBC
Xi = getRBCInitialGeometry(thet,phi,InitXi,InitOrient);
%% Initial conditions:
%% spherical harmonics coefficients of the undeformed RBC coordinates
[aXi,bXi] = shagcm(Xi);
%% Compute the first and second fundamental form coefficients
[~, ~, EE, FF, GG, WW, JXibrev, ~, ~, ~, LL, MM, NN] = ...
coefficientsOfFundamentalForm(aXi, bXi, UpSampleFactor);
%% Check the geometry and position of RBC
if verbose_Plot
figure('Color','white')
hold on
VisualizeGeometry(nlat, nlon, aXi, bXi, 'r', true)
axis on
xlabel('X'); ylabel('Y'); zlabel('Z');
end
%% Create rotation matrices
RotationMatrix = RotationMatrixForSHCoefficients(N, nlat, nlon, thet, phi);
%% Masks to go between Spherepack and vector representations of SH coeff
mask_a = repmat(triu(true(N+1),0),1,1,3);
mask_b = mask_a;
mask_b(1,:,:) = false;
%% Set up output file
fidTime = fopen(['Time_',name,'.dat'],'w');
fidCoord = fopen(['Coord_',name,'.dat'],'w');
fidMemFor = fopen(['MemFor_',name,'.dat'],'w');
fidSol = fopen(['Sol_',name,'.dat'],'w');
%% Set up time integration
cu = zeros(3*(N+1)^2,1);
%% Initialization
xi = Xi; axi = aXi; bxi = bXi;
viscousStress_prev = zeros(UpSampleFactor*N+1, 2*UpSampleFactor*N+1, 4);
epsilbrev_prev = zeros(UpSampleFactor*N+1, 2*UpSampleFactor*N+1, 4);
%% Time-stepping
for nstep = 0:NSTEPS
if (nstep==0 || mod(nstep,SaveAtIncrementalSteps)==0)
%% Write the time to file
fwrite(fidTime, Time, 'double');
%% Write the position of cell to file
cxi = [axi(mask_a); bxi(mask_b)];
fwrite(fidCoord, cxi, 'double');
end
%% Ambient simple shear flow
if ShearFlow
uinf = zeros(size(Xi));
uinf(:,:,1) = ShearRate*xi(:,:,2); % Shear in y-direction
% uinf(:,:,1) = ShearRate*xi(:,:,3); % Shear in z-direction
[auinf,buinf] = shagcm(uinf);
if Relaxation && nstep > RelaxationStep
uinf = zeros(size(Xi));
auinf = zeros(size(axi));
buinf = zeros(size(bxi));
end
end
%% Ambient parabolic velocity field
if ParabolicFlow
uinf = zeros(size(Xi));
uinf(:,:,1) = aParabolic*(bParabolic-(xi(:,:,2).^2 + xi(:,:,3).^2));
[auinf,buinf] = shagcm(uinf);
if Relaxation && nstep > RelaxationStep
uinf = zeros(size(Xi));
auinf = zeros(size(axi));
buinf = zeros(size(bxi));
end
end
%% Compute membrane traction
[f, viscousStress_prev, epsilbrev_prev] = ...
MembraneForces(EE, FF, GG, WW, JXibrev, ...
LL, MM, NN, ...
axi, bxi, ...
UpSampleFactor, ...
ES, ED, EB, ...
StVenantKirchhoff, ...
NeoHookean, Skalak, ...
MembraneViscoelasticity, ...
mu_Mem, Tau, DT, ...
viscousStress_prev, ...
epsilbrev_prev);
[af,bf] = shagcm(f);
%% Solve Boundary Integral Equation
[u, au, bu, cu, FLAG, RELRES, ITER] = ...
SolveSpectralBIE(axi, bxi, xi, ...
RotationMatrix, ...
NGSphere, NGthet, NGphi, ...
eta, wg, ...
af, bf, auinf, buinf, ...
nlat, nlon, mask_a, mask_b, ...
mu, lam, N, ToleranceGMRES, ...
cu);
if (nstep==0 || mod(nstep,SaveAtIncrementalSteps)==0)
%% Write output to file
writef = f(:)';
fwrite(fidMemFor,writef,'double');
writecu = cu';
fwrite(fidSol,writecu,'double');
end
%% Display
f_norm = sqrt(f(:,:,1).^2+f(:,:,2).^2+f(:,:,3).^2);
u_norm = sqrt(u(:,:,1).^2+u(:,:,2).^2+u(:,:,3).^2);
fprintf('# of step %d;\t remaining # of steps %d; \n', ...
nstep, NSTEPS-nstep)
fprintf('normf = %g;\t normu = %g\n', max(max(f_norm)),max(max(u_norm)))
fprintf('DT: %g;\t Time %d \n', DT, Time);
fprintf('# of iteration GMRES %d\n', ITER(2))
fprintf('# of stableCounter %d\n', stableCounter)
fprintf('# of counterTime %d\n', counterTime)
fprintf('\n')
%% Update state (Forward Euler time scheme)
axi = axi + DT*au;
bxi = bxi + DT*bu;
xi = shsgcm(axi,bxi);
%% Set the center of the mass of RBC to the initial position, InitXi
if ParabolicFlow
[gxi_thet, gxi_phi] = gradgcm(axi, bxi);
normalVector = cross(gxi_thet, gxi_phi, 3);
xi_dot_normalVector = sum(xi.*normalVector,3);
xi_xi_dot_normalVector = xi.*xi_dot_normalVector;
V = sum(sum(xi_dot_normalVector.*weight))*(2*pi/nlon)/3;
Int_xi_xi_dot_normalVector_dA = ...
sum(sum(xi_xi_dot_normalVector.*weight))*(2*pi/nlon)/4;
CenterMass = reshape(Int_xi_xi_dot_normalVector_dA./V,3,1);
xi = xi - reshape((CenterMass - InitXi),1,1,3);
end
[axi,bxi] = shagcm(xi);
%% Update current time
Time = Time + DT;
%% Update time step increment
[DT, counterTime, IterationCounter, stableCounter, maxStableDT] = ...
updateDT(DT, minDT, maxDT, ...
ITER, IterationCounter, ...
counterTime, maxCounterTime, ...
stableCounter, maxStableDT, ...
nstep);
%% Return if time is larger than the final time
if Time - DT > EndTime
load handel
sound(y,Fs)
fclose('all');
TimeInHours = toc(Starttime)/60/60 % in hours
return
end
end
load handel
sound(y,Fs)
fclose('all');
RunTimeInHours = toc(Starttime)/60/60