plot_snap.m 16.4 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
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
%
% Main function to produce field plots and profiles for one time step
%
function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D


  [data,input,program_dir,fname,fpath] = read_snap(fin); % read data from snapshot
%  cd(program_dir) % return to octave folder
  if (part==1)
%    fpath
%    fname
    pname = strrep(fname,'data_snapshot','data_particles');
%    [A,input,program_dir,fpath] = read_particles;
    A = read_particles(fpath,pname);
    A(:,1:2) = A(:,1:2) * data.D / 1000.0;
  else
    A=0;
  endif  
  cd(program_dir) % return to octave folder
  [data.p_lat,data.p_rad,data.p_xx,data.p_zz] = read_mesh(fpath,input,data.nl,data.ny,data.nr); % read mesh data (position of cells)
  cd(program_dir) % return to octave folder
  data = get_data(data,input); % get fields for plotting and dimensionalize if needed
  cd(fpath); 
  mkdir "plots";
  graphics_toolkit gnuplot

28
29
30
  
  size(data.depl)
  
Lena Noack's avatar
Lena Noack committed
31
  res=1200; %600;
Lena Noack's avatar
Lena Noack committed
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  ar=(max(max(data.p_xx_dim))-min(min(data.p_xx_dim)))/(max(max(data.p_zz_dim))-min(min(data.p_zz_dim)));

  time = data.time/data.t_yr/1000000.0
  iter = data.iter

  if (reduc==0) 
  % plot all relevant radial profiles (min/mean/max profiles)
  t_profile = plot_x_profile(data.T,data,input,"PrT","Temperature","southwest"," [K]");
  if (input.use_melt==1), t_profile = plot_x_profile(data.T,data,input,"PrTM","Temperature and Melting temp.","southwest"," [K]",log=0,melt=1);, endif
  v_profile = plot_x_profile(data.vel,data,input,"PrVel","Velocity","northeast"," [cm/yr]");
  vi_profile = plot_x_profile(data.visc,data,input,"PrVisc","Viscosity","northeast"," [log_{10} Pa\ s]",log=1);
  if(size(data.depl,1)>1), d_profile = plot_x_profile(data.depl,data,input,"PrD","Depletion","southeast","");, endif
  endif

  % plot fields
  if(max_depl==0),max_depl=max(max(data.depl));, endif
  min_T=min(min(data.T));
  if(max_T==0),max_T=max(max(data.T));, endif
  
51
  if (reduc<=2) 
Lena Noack's avatar
Lena Noack committed
52
53
54
55
  % useful colormaps are for example hot, copper, jet (also viridis, but load extra?)
  produce_plot('T2D_',"Temperature in K",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV)
  produce_plot('Vi2D_',"Viscosity in log_{10} Pa\ s",log10(data.visc),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
  produce_plot('V2D_',"Velocity in cm/yr",data.vel,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
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
  
  % data.cr_age only surface values, need to be set everywhere where there is crust
  % data.oc_vol_tot is total volume per lateral/perpendicular cell, need to get depth of crust
  nl = size(data.T,1);
  nr = size(data.T,3);
  cr_depth = zeros(nl,1);
  cr_vol = data.oc_vol_tot;
  dl = data.p_lat(2)-data.p_lat(1)
  dr = data.p_rad(2)-data.p_rad(1)
  for i=1:nl % in 3D would need do add for-loop for size(data.T,2)
    for k=nr:-1:1
      vol = 2*pi*dl*((data.p_rad(k)+dr/2)^2-(data.p_rad(k)-dr/2)^2); % vol of cell in 2D
%      vol = 4/3*pi*dl*((data.p_rad(k)+dr/2)^3-(data.p_rad(k)-dr/2)^3); % vol of cell
      if (cr_vol(i)>0);%vol)
        cr_vol(i) = cr_vol(i)-vol;
        cr_depth(i) = cr_depth(i) + 1; % number of cells from top downwards filled as crust
      else
        break
      end
    end  
  end
  fprintf('Crust is between %d and %d cells thick\n',min(cr_depth),max(cr_depth))
  crust = zeros(size(data.T));
  for i=1:size(data.T,1) % in 3D would need do add for-loop for size(data.T,2)
    crust(i,1,nr-cr_depth(i):nr)=max(data.cr_age)-data.cr_age(i);
  end
  
  produce_plot('Cr2D_',sprintf(strcat("Crust age at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),crust,input,data,res,ar,50,'bgry',min_z=-1,max_z=-1,full=fullV)
Lena Noack's avatar
Lena Noack committed
84
85
  endif
  
86
  if (reduc==0) 
Lena Noack's avatar
Lena Noack committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
  if(part==1)
    produce_plot('OP2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,21))
%%    produce_plot('PD2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
%%    produce_plot('TD2D_',"Tracer",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
  elseif(size(data.output_tracer,1)>1)
    produce_plot('O2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
  endif

  %min_z=-max(abs(min(min(data.vl))),max(max(data.vl)));
  %max_z=max(abs(min(min(data.vl))),max(max(data.vl)));
  produce_plot('Vlat2D_',"Lateral velocity in cm/yr",data.vl,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
  %min_z=-max(abs(min(min(data.vr))),max(max(data.vr)));
  %max_z=max(abs(min(min(data.vr))),max(max(data.vr)));
  produce_plot('Vrad2D_',"Radial velocity in cm/yr",data.vr,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
  if(min(min(data.strR))>0), produce_plot('StR2D_',"Strain rate in log_{10} 1/s",log10(data.strR),input,data,res,ar,50,'autumn',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.c,1)>1), produce_plot('C2D_',"Composition",data.c,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  endif

  if(size(data.ref_r,1)>1)
    if (reduc==0) 
    produce_plot('De2D_',"Density",data.ref_r,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
    produce_plot('Mrad2D_',"Radial mass transport in kg/(m^2 yr)",data.vr.*data.ref_r/100,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
109
110
    %elseif (reduc==1)
    %produce_plot('Mradp2D_',"Total radial mass transport in kg/(m^2 yr)",abs(data.vr.*data.ref_r/100),input,data,res,ar,10,'hot',min_z=-1,max_z=-1,full=fullV)
Lena Noack's avatar
Lena Noack committed
111
112
113
    endif
  endif

114
115
116
117
%  if(size(data.depl,1)>1), produce_plot('D2Dc_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),-data.depl,input,data,res,ar,50,'jet',min_z=-max_depl,max_z=0,full=fullV), endif
%  if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
%
%  if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
Lena Noack's avatar
Lena Noack committed
118
119

  if (reduc==0) 
120
  if(size(data.DeltaF,1)>1), produce_plot('DF2D_',"Melt fraction",data.DeltaF,input,data,res,ar,20,'jet',min_z=-1,max_z=-1,full=fullV), endif
121
%  if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-100,max_z=0,full=fullV), endif %input.H2O), endif
122

Lena Noack's avatar
Lena Noack committed
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
241
242
243
  if(size(data.ref_a,1)>1), produce_plot('Al2D_',"Heat expansion",data.ref_a,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif %input.nr_ph
%  if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=3,max_z=7,full=fullV), endif %input.nr_ph
  if(size(data.Cp,1)>1), produce_plot('Cp2D_',"Heat capacity",data.Cp,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.cond,1)>1), produce_plot('K2D_',"Conductivity",data.cond,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.heat,1)>1), produce_plot('H2D_',"Heat souce redistribution",data.heat,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  %if(size(data.visc_diff,1)>1), produce_plot('ViDiff2D_',"Viscosity Diff in log",log10(data.visc_diff),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  %if(size(data.CO2,1)>1), produce_plot('CD2D_',"CO2 content",data.CO2,input,data,res,ar,50,'jet',min_z=0,max_z=500,full=fullV), endif %input.CO2), endif
  if(size(data.U238,1)>1), produce_plot('HS_U2382D_',"U238 in ppm",data.U238,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.U235,1)>1), produce_plot('HS_U2352D_',"U235 in ppm",data.U235,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.Th232,1)>1), produce_plot('HS_Th2322D_',"Th232 in ppm",data.Th232,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.K40,1)>1), produce_plot('HS_K402D_',"K40 in ppm",data.K40,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  if(size(data.entropy,1)>1), produce_plot('En2D_',"Entropy",data.entropy,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  
  % lateral temperature variations
  data.T_var = data.T; %initialization size
  %[tprof,tprof] = meshgrid(t_profile(:,2)); % ToDo: speed-up if t_profile matrix set up, and then just subract matrix from matrix (no loop needed)? 2D/3D - how to do that?

  k=1;
   do
    data.T_var(:,:,k) = data.T_var(:,:,k) - t_profile(k,2); 
    k=k+1;
   until (k>data.nr)

%  min_z=-max(abs(min(min(data.T_var))),max(max(data.T_var)));
%  max_z=max(abs(min(min(data.T_var))),max(max(data.T_var)));
  produce_plot('Tvar2D_',"Lateral temperature variation in K",data.T_var,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)%,min_z,max_z)
  
  %cd(program_dir) % return to octave folder
  %colormap('default') % re-set color map
  endif %reduced output
  
end

% plot min/mean/max profiles for field x
function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dimen,log=0,melt=0)

  nl=data.nl; % for profiles independent of periodic or reflective boundary conditions
  ny=data.ny;
  nr=data.nr;

  x_profile = zeros(nr,3);
  m_profile = zeros(nr,2);
  r = data.p_rad_dim; % radius vector

  k=1;
  do
    x_profile(k,1) = min(field(1:nl,1,k));
    x_profile(k,2) = mean(field(1:nl,1,k)); % ToDo: include dVi to get right mean value for non-uniform lateral cell distribution
    x_profile(k,3) = max(field(1:nl,1,k));
    if (melt==1)
     p = (max(r)-r(k))*input.rho*input.g/10.0^6; %km -> GPa
     if (p<input.mz)
      m_profile(k,1) = input.m1s0 + input.m1s1*p + input.m1s2*p^2 + input.m1s3*p^3;
      m_profile(k,2) = input.m1l0 + input.m1l1*p + input.m1l2*p^2 + input.m1l3*p^3;
     else
      m_profile(k,1) = input.m2s0 + input.m2s1*p + input.m2s2*p^2 + input.m2s3*p^3 + input.m2s4*p^4;
      m_profile(k,2) = input.m2l0 + input.m2l1*p + input.m2l2*p^2 + input.m2l3*p^3 + input.m2l4*p^4;
     endif
    endif
    k=k+1;
  until k>nr  

  if (log==1), x_profile = log10(x_profile);, endif
  
  f = figure('visible','off');
  
  h=plot(x_profile(:,1),r,"b");
  hold on
  h2=plot(x_profile(:,2),r,"k");
  h3=plot(x_profile(:,3),r,"r");
  if (melt==1)
    h4=plot(m_profile(:,1),r,":g");
    h5=plot(m_profile(:,2),r,"--g");
  endif
  hold off
  set(h,'LineWidth',5)
  set(h2,'LineWidth',10)
  set(h3,'LineWidth',5)
  if (melt==1)
    set(h4,'LineWidth',5)
    set(h5,'LineWidth',5)
  endif
%  legend(sprintf(strcat(name,' min')),sprintf(strcat(name," mean")),sprintf(strcat(name," max")), "location", orient);
%  legend("boxoff");
  
  FN = findall(f,'-property','FontName');
  set(FN,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf');
  FS = findall(f,'-property','FontSize');
  set(FS,'FontSize',20);

%  dimen="";
%  if (data.D>1.0), dimen=" [K]";, endif
  dimenR="";
  if (data.D>1.0), dimenR=" [km]";, endif

  xlabel (sprintf(strcat(title_plot,dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',24);
  ylabel (sprintf(strcat("Radius",dimenR)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',24);
  title (sprintf(strcat(title_plot," profile")),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Bold.ttf','FontSize',30);
  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);                                                                          

%  graph=strcat("-S", num2str((ar+0.4)*res),",", num2str(res));
  print(filename,'-dpng','-S800,720');%graph); 

  clear f;
  
end

% help function for plot_snap to produce png plots for selected fields
function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,min_z=-1,max_z=-1,full=-1,partX=0,partY=0,partV=0)

%  min_z
%  max_z

  if (min_z==-1)
    min_z = min(min(field));
  endif
  if (max_z==-1)
    max_z = max(max(field));
  endif

244
  font_title=round(20*res/600);%(30*res/600);
Lena Noack's avatar
Lena Noack committed
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
  font_axes_title=round(24*res/600);
  font_axes=round(20*res/600);
  
  f = figure('visible','off');
  
  if (data.ny>1)
    ny_use=data.ny/2;
  else
    ny_use=1;
  endif
  dummy(:,:)=field(:,ny_use,:);
%  contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour); % works both for spherical and cartesian plots

  %[size(data.p_xx_dim),  size(data.p_zz_dim),  size(dummy), size(field)]
%  size(partX,1)
  if (size(partX,1)>1) 
%    fprintf("Particle plot")
%    size(partX,1)
    N=length(partX);
    colorMap = zeros(N, 3);
%    for k = 1:N
%     if (input.set_tracers==1)
%     elseif (input.set_tracers==2)
%%       colorMap(k,:) = (partV(k)-min(partV))/(max(partV)-min(partV))*[1,0,0]+(1-(partV(k)-min(partV))/(max(partV)-min(partV)))*[0,0,1]; % shifts data points between red and blue depending on value
%       colorMap(k,:) = partV(k)/max(partV)*[1,0,0]+(1-partV(k)/max(partV))*[0,0,1]; % shifts data points between red and blue depending on value
%     endif
%    end
    
    if (input.set_tracers==1)
     %colorMap(partV<0.25)
     for k = 1:N
      if (partV(k)<0.25)
        colorMap(k,:) = [1,0,0];
      elseif (partV(k)>0.75)
        colorMap(k,:) = [1,1,1];
      else
        colorMap(k,:) = [0,0,1];
      endif
     end        
    elseif (input.set_tracers==2)
      if (max(partV) != min(partV))
        colorMap(:,1) = (partV-min(partV))/(max(partV)-min(partV));
        colorMap(:,2) = 1 - (partV-min(partV))/(max(partV)-min(partV));
      else
        colorMap(:,1) = partV/max(partV);
        colorMap(:,3) = 1 - partV/max(partV);
      endif
    end
    
%    dummy(:,:)=0.4;
%    contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
%    hold on
    %scatter(partX,partY,100,colorMap,'filled')
    scatter(partX,partY,100,partV,'filled')
    %partX
    %plot(partX,partY)
%    hold off
  else  
%    contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
    [C,h] = contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
    set(h,'LineColor','none')
306
    set(gca,'visible','off')
Lena Noack's avatar
Lena Noack committed
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
  endif

%  surface(data.p_xx_dim,data.p_zz_dim,dummy,'EdgeColor','none'); % works both for spherical and cartesian plots

%  surfc(data.p_xx_dim,data.p_zz_dim,dummy,'EdgeColor','none'); % works both for spherical and cartesian plots
%  shading interp;
%  view(0,90)


% levelplot
  if (strcmp(cmap,'bgry'))
    mycolormap % call mycolormap function to get colormap (mymap) black-blue-green-red-yellow
    colormap(mymap)
  else
    colormap(cmap)
  endif
Claire Guimond's avatar
Claire Guimond committed
323
  %caxis([min_z max_z]);
Lena Noack's avatar
Lena Noack committed
324
325
326
327
328
329
330
331
332
333
334
335
  colorbar

  FN = findall(f,'-property','FontName');
  set(FN,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf');%dejavu/DejaVuSans.ttf');
  FS = findall(f,'-property','FontSize');
  set(FS,'FontSize',font_axes);

%  axis ("xy", "square");  
%  axis ([0, ar, 1, 2]);
  axis ([min(min(data.p_xx_dim)), max(max(data.p_xx_dim)), min(min(data.p_zz_dim)), max(max(data.p_zz_dim))]); % min_z max_z dazu?
  dimen="";
  if (data.D>1.0), dimen=" [km]";, endif
336
337
%  xlabel (sprintf(strcat("Lateral direction",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
%  ylabel (sprintf(strcat("Radius",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
Lena Noack's avatar
Lena Noack committed
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
  title (title_plot,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Bold.ttf','FontSize',font_title);%dejavu/DejaVuSans-Bold.ttf','FontSize',font_title);

  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);
%  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.eps'),f);

  if (full==-1)
%    fprintf("Regional Sphere")
    graph=strcat("-S", num2str((ar+0.6)*res),",", num2str(res)); % suitable for regional sphere
  else
%    fprintf("Full Sphere")
    graph=strcat("-S", num2str((ar+0.37)*res),",", num2str(res)); % suitable for full sphere
  endif
%  print(filename,'-deps',graph);   

  print(filename,'-dpng',graph);

  clear f;
  
Lena Noack's avatar
Lena Noack committed
356
end