-
Notifications
You must be signed in to change notification settings - Fork 1
/
fluxODE.m
240 lines (238 loc) · 7.24 KB
/
fluxODE.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
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
% fluxODE uses matlabs boundary value solver to find the steady state solution
% of the RD equation for various parameter configurations by solving the ODE ( dv/dt = 0 ).
% Using this, it calculated the RHS flux and makes nice plots. Loops over
% nu, KonBt, Koff.
% fluxODE( plotMapFlux, plotSteady, saveMe, dirname )
%
% Inputs:
% plotFlag: structure of plot flags
% storeFlag: structure of store flags
% saveMe: save plots and outputs
% dirname: directory name for saving
% paramFile: initParam file
%
% plotFlag with fields
% plotFlag.plotMapFlux: surface plot jmax vs koff and konbt
% plotFlag.plotSteady: plot concentration profiles
%
% storeFlag with fields
% storeFlag.storeStdy: store steady state solution flag
%
% Outputs: fluxSummary with fields
% jMax: matrix of steady state flux vs koff and konbt
% jNorm: jMax ./ jDiff
% AconcStdy: matrix of A steady state profile vs koff and konbt
% CconcStdy: matrix of C steady state profile vs koff and konbt
% params: parameters of runs
%
% example run
%
% plotFlag.plotSteady = 1;
% plotFlag.plotMapFlux = 1;
% storeFlag.storeStdy = 1;
% saveMe = 1;
% dirname = 'blah';
%
% [fluxSummary] = fluxODE( plotFlag, storeFlag, saveMe, dirname, paramFile );
%
function [ fluxSummary ] = ...
fluxODE( plotFlag, storeFlag, saveMe, dirname, paramFile )
try
% Latex font
set(0,'defaulttextinterpreter','latex')
% Make up a dirname if one wasn't given
totalInput = 5;
if nargin < totalInput-1
if saveMe == 1
dirname = ['fluxODE_' num2str( randi( 100 ) )];
else
dirname = ['tempFluxODE_' num2str( randi( 100 ) ) ];
end
end
if nargin < totalInput
paramFile = 'initParams.m';
end
% move input structure fields to variables
plotMapFlux = plotFlag.plotMapFlux;
plotSteady = plotFlag.plotSteady;
storeStdy = storeFlag.storeStdy;
% can't plot steady if not storing
plotSteady = storeStdy * plotSteady;
% Add paths and output dir
addpath( genpath('./src') );
if ~exist('./steadyfiles','dir'); mkdir('steadyfiles'); end
if ~exist('./steadyfiles/ODE','dir'); mkdir('steadyfiles/ODE'); end
% print start time
Time = datestr(now);
fprintf('Starting fluxODE: %s\n', Time)
% Initparams
fprintf('Initiating parameters\n');
if exist( paramFile,'file')
fprintf('Init file: %s\n', paramFile);
run( paramFile );
elseif exist( 'initParams.m', 'file')
fprintf('Could not find init file: %s. Running initParams\n', ...
paramFile);
run( 'initParams.m');
else
fprintf('Could not find init file: %s or initParams. Copying and running template\n', ...
paramFile);
cpParams
initParams
end
% Copy master parameters input object
paramObj = paramMaster;
flagsObj = flags;
% Code can only handle one value of Bt currently
if length( paramObj.Bt ) > 1
paramObj.Bt = paramObj.Bt(1);
end
% set-up params
pfixed = paramObj.Bt;
BtFixed = paramObj.Bt;
pfixedStr = '$$ B_t $$';
[paramObj, kinParams] = paramInputMaster( paramObj, koffVary );
% Run the loops
paramNuLlp = kinParams.nuLlp;
paramKonBt = kinParams.konBt;
paramKoffInds = kinParams.koffInds;
numRuns = kinParams.numRuns;
% warn about low N
if paramObj.Nx < 1000
fprintf('Warning, very low number of grid points\n')
end
% save names
saveStrFM = 'flxss'; %flux map
saveStrSS = 'profileSS'; % steady state
saveStrMat = 'fluxSummary.mat'; % matlab files
if saveMe
dirname = [dirname '_nl' num2str( flagsObj.NLcoup )];
mkdir( dirname );
end
if plotMapFlux
% set colormap
randI = randi(100000);
figure(randI)
colormap( viridis );
close(randI)
% labels
ylab = kinParams.kinVar1strTex; % rows
xlab = kinParams.kinVar2strTex; % columns
end
if plotSteady
p1name = kinParams.p1nameTex;
p2name = kinParams.kinVar1strTex;
p3name = kinParams.kinVar2strTex;
end
% Specify necessary parameters for parfor
nlEqn = flags.NLcoup;
Da = paramObj.Da; AL = paramObj.AL; AR = paramObj.AR;
Bt = paramObj.Bt; Nx = paramObj.Nx; Lbox = paramObj.Lbox;
if strcmp( paramObj.A_BC,'Dir' ) && strcmp( paramObj.C_BC, 'Vn' )
BCstr = 'DirVn';
elseif strcmp( paramObj.A_BC,'Dir' ) && strcmp( paramObj.C_BC, 'Vn' )
BCstr = 'Dir';
elseif strcmp( paramObj.A_BC,'Vn' ) && strcmp( paramObj.C_BC, 'Vn' )
BCstr = 'Vn';
else
fprintf( 'I cannot handle those BC, doing A = Dir, C = Vn \n')
BCstr = 'DirVn';
end
% Flux matrix
jMax = zeros( 1, numRuns);
% Store steady state solutions;
AconcStdy = cell( numRuns, 1 );
CconcStdy = cell( numRuns, 1 );
% Calculated things
x = linspace(0, Lbox, Nx) ;
dx = x(2) - x(1);
if numRuns > 1 && flags.ParforFlag
recObj = 0;
parobj = gcp;
numWorkers = parobj.NumWorkers;
fprintf('I have hired %d workers\n',parobj.NumWorkers);
else
fprintf('Not using parfor\n')
numWorkers = 0;
end
% set bound diffusion or not
nuCell = cell(1, numRuns);
for ii = 1:numRuns
nuCell{ii} = { paramObj.DbParam{1}, paramNuLlp(ii) };
end
% set up koff cell
koffCell = cell( 1, numRuns );
for ii = 1:numRuns
koffCell{ii} = paramObj.koffObj.InfoCell{ paramKoffInds(ii) };
end
NxSave = Nx;
parfor (ii=1:numRuns, numWorkers)
% set params
nuCellTemp = nuCell{ii};
KonBt = paramKonBt(ii);
koffCellTemp = koffCell{ ii };
Kon = KonBt ./ BtFixed;
[AnlOde,CnlOde,~] = RdSsSolverMatBvFunc(...
Da, Kon, koffCellTemp, nuCellTemp, AL, AR, BtFixed, Lbox, BCstr, Nx, ...
nlEqn );
% calc flux
flux = - Da * ( AnlOde(end) - AnlOde(end-1) ) / dx;
% record
if storeStdy
AconcStdy{ii} = AnlOde;
CconcStdy{ii} = CnlOde;
else
AconcStdy{ii} = 0;
CconcStdy{ii} = 0;
end
jMax(ii) = flux;
end
% reshape to more intutive size---> Mat( p1, p2, p3, : )
numP1 = kinParams.numP1;
numP2 = kinParams.numP2;
numP3 = kinParams.numP3;
AconcStdy = reshape( AconcStdy, [numP1, numP2, numP3] );
CconcStdy = reshape( CconcStdy, [numP1, numP2, numP3] );
jMax = reshape( jMax, [numP1, numP2, numP3] );
% Get flux diff and normalize it
jDiff = Da * ( AL - AR ) / Lbox;
jNorm = jMax ./ jDiff;
% Steady states
if plotSteady
concSteadyPlotMultParams( AconcStdy, CconcStdy, x, ...
kinParams.p1Vec, kinParams.kinVar1, kinParams.kinVar2, ...
p1name, p2name, p3name, ...
pfixed, pfixedStr, saveMe, saveStrSS )
end
% Surface plot
if plotMapFlux
titstr = ['$$ j_{max} / j_{diff} $$; $$ B_t = $$ ' num2str(BtFixed)...
'; ' kinParams.p1nameTex ' = '];
surfLoopPlotter( jNorm, kinParams.p1Vec, kinParams.kinVar1, kinParams.kinVar2,...
xlab, ylab, titstr, saveMe, saveStrFM )
end
% store everything
fluxSummary.jMax = jMax;
fluxSummary.jNorm = jNorm;
fluxSummary.jDiff = jDiff;
fluxSummary.aConcStdy = AconcStdy;
fluxSummary.cConcStdy = CconcStdy;
fluxSummary.paramObj = paramObj;
fluxSummary.kinParams = kinParams;
% save data
if saveMe
save(saveStrMat, 'fluxSummary');
% make dirs and move
if plotSteady || plotMapFlux
movefile('*.fig', dirname);
movefile('*.jpg', dirname);
end
movefile(saveStrMat, dirname);
movefile(dirname, './steadyfiles/ODE' )
end
mytime = datestr(now);
fprintf('Finished %s: %s\n', mfilename, mytime)
catch err
fprintf('try/catch err in %s \n', mfilename('fullpath') )
fprintf('%s',err.getReport('extended') );
end