1 |
clear path |
2 |
|
3 |
global Nx Ny Nz |
4 |
global lat long dz dm mdep |
5 |
global delt_su su_its t_su delt |
6 |
global descriptor this_path |
7 |
global f deltaf Q beta r_expt r_heat H |
8 |
global time rots it |
9 |
global g Cp rho_bar alpha |
10 |
global u v t w |
11 |
global iterations |
12 |
|
13 |
|
14 |
param_file_name = ... |
15 |
input(' Please enter the name of the m-file with the parameters for this run : ','s') ; |
16 |
feval(param_file_name) ; |
17 |
|
18 |
% iterations |
19 |
|
20 |
itstart = input(' Please enter start iteration : ','s') |
21 |
itend = input(' Please enter end iteration : ','s') |
22 |
|
23 |
|
24 |
sizeit=size(iterations); |
25 |
for i=1:sizeit(1) |
26 |
iter(i)=eval(iterations(i,1:10)); |
27 |
end |
28 |
nitstart=find(iter==eval(itstart)) |
29 |
nitend=find(iter==eval(itend)) |
30 |
|
31 |
path = this_path |
32 |
cmdstr=['cd ' path ]; |
33 |
eval(cmdstr); |
34 |
path=pwd |
35 |
|
36 |
sumtheta=zeros(Nx,Ny,Nz); |
37 |
sumu=zeros(Nx,Ny,Nz); |
38 |
sumv=zeros(Nx,Ny,Nz); |
39 |
sumw=zeros(Nx,Ny,Nz); |
40 |
counter=0; |
41 |
|
42 |
for i=nitstart:nitend |
43 |
i,0 |
44 |
tfilename=(['T.' iterations((i),1:10) ]) ; |
45 |
t=rdmds(tfilename,'b'); |
46 |
|
47 |
ufilename=(['U.' iterations((i),1:10) ]) ; |
48 |
u=rdmds(ufilename,'b'); |
49 |
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:)); |
50 |
|
51 |
vfilename=(['V.' iterations((i),1:10) ]) ; |
52 |
v=rdmds(vfilename,'b'); |
53 |
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:)); |
54 |
|
55 |
wfilename=(['W.' iterations((i),1:10) ]) ; |
56 |
w=rdmds(wfilename,'b'); |
57 |
|
58 |
sumtheta=sumtheta+t; |
59 |
sumu=sumu+u; |
60 |
sumv=sumv+v; |
61 |
sumw=sumw+w; |
62 |
counter=counter+1; |
63 |
end |
64 |
|
65 |
meantheta=sumtheta/counter; |
66 |
meanu=sumu/counter; |
67 |
meanv=sumv/counter; |
68 |
meanw=sumw/counter; |
69 |
|
70 |
sumut=zeros(Nx,Ny,Nz); |
71 |
sumvt=zeros(Nx,Ny,Nz); |
72 |
sumwt=zeros(Nx,Ny,Nz); |
73 |
|
74 |
meanu=smooth3(meanu); |
75 |
meanv=smooth3(meanv); |
76 |
meanw=smooth3(meanw); |
77 |
meantheta=smooth3(meantheta); |
78 |
|
79 |
for i=nitstart:nitend |
80 |
i |
81 |
tfilename=(['T.' iterations((i),1:10) ]) ; |
82 |
t=rdmds(tfilename,'b'); |
83 |
|
84 |
ufilename=(['U.' iterations((i),1:10) ]) ; |
85 |
u=rdmds(ufilename,'b'); |
86 |
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:)); |
87 |
|
88 |
vfilename=(['V.' iterations((i),1:10) ]) ; |
89 |
v=rdmds(vfilename,'b'); |
90 |
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:)); |
91 |
|
92 |
wfilename=(['W.' iterations((i),1:10) ]) ; |
93 |
w=rdmds(wfilename,'b'); |
94 |
|
95 |
t=t-meantheta; |
96 |
u=u-meanu; |
97 |
v=v-meanv; |
98 |
w=w-meanw; |
99 |
|
100 |
sumut=sumut+u.*t; |
101 |
sumvt=sumvt+v.*t; |
102 |
sumwt=sumwt+w.*t; |
103 |
|
104 |
end |
105 |
|
106 |
sumut=sumut/counter; |
107 |
sumvt=sumvt/counter; |
108 |
sumwt=sumwt/counter; |
109 |
|
110 |
|
111 |
D=zeros(Nx,Ny,Nz); |
112 |
|
113 |
dx=dm; |
114 |
dy=dm; |
115 |
dz=-dz; |
116 |
|
117 |
D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ... |
118 |
-(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ... |
119 |
-(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz); |
120 |
|
121 |
Dw=zeros(Nx,Ny,Nz); |
122 |
ff=zeros(Nx,Ny,Nz); |
123 |
|
124 |
|
125 |
for j=1:Ny |
126 |
ff(:,j,:)=f+(j-1)*dy*beta; |
127 |
end |
128 |
|
129 |
Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)*2*dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2)); |
130 |
|
131 |
i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2)); |
132 |
Dw(i)=0; |
133 |
|
134 |
|
135 |
t0=input('Enter temperature : ') |
136 |
|
137 |
% Interpolate Dw on T=T0 surface |
138 |
|
139 |
h=zeros(Nx,Ny); |
140 |
Dh=zeros(Nx,Ny); |
141 |
dz=-dz |
142 |
|
143 |
meantheta(:,:,Nz)=0.; |
144 |
for i=1:Nx |
145 |
for j=1:Ny |
146 |
kk=find(meantheta(i,j,:)<t0); |
147 |
|
148 |
if kk(1)>1 |
149 |
h(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); |
150 |
|
151 |
Dh(i,j)=Dw(i,j,kk(1)-1) ... |
152 |
+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ... |
153 |
*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); |
154 |
else |
155 |
h(i,j)=0; |
156 |
Dh(i,j)=0; |
157 |
end |
158 |
|
159 |
end |
160 |
end |
161 |
|
162 |
hmax=max(max(h)) |
163 |
|
164 |
figure |
165 |
|
166 |
pcolor(Dh.');shading flat; colorbar; axis square |
167 |
title('mean Wstar'); |
168 |
% title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); |
169 |
|
170 |
figure |
171 |
|
172 |
pcolor(h.');shading flat; colorbar; axis square |
173 |
%pcolor(flipud(h.'));shading flat; colorbar; axis square |
174 |
%title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); |
175 |
|
176 |
|
177 |
we=0; |
178 |
ws=0; |
179 |
for i=1:Nx |
180 |
for j=1:Ny |
181 |
|
182 |
r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx; |
183 |
|
184 |
%if r<0.44 |
185 |
|
186 |
we=meanw(i,j,3)*dx*dx+we; |
187 |
ws=Dh(i,j)*dx*dx+ws; |
188 |
%end |
189 |
end |
190 |
end |
191 |
|
192 |
we |
193 |
ws |
194 |
|
195 |
|
196 |
|
197 |
|
198 |
|
199 |
|