% Material visualization with velocity clear all close all dir1='../../'; % input data path dir2='../../'; % output figure path % variables on nodes name = 'RigSubc'; % comp name_t = 'RigSubt'; % T % number nn1 = 500;%[0 10:20:1200]; % extension of data format ext = '.txt'; % defining size of the model xsize = 5050; ysize = 1100; % defining zoomed area xbeg = 1800; xend = 3000; ybeg = 0; yend = 400; % Defining zoom: YES > 0; NO < 0 zoom = 1; % Defining melting area s1 = 0*ones(1,length(nn1)); s2 = 0*ones(1,length(nn1)); s3 = 0*ones(1,length(nn1)); s4 = 0*ones(1,length(nn1)); % main loop - insert numbers of images for idx=1:length(nn1) cd (dir1) fname = [name,num2str(nn1(idx)) ext]; %read input data fid = fopen(fname,'r'); if(fid<=0) fname disp('No such file...') break; end A = textscan(fid, '%n'); fclose(fid); Data = A{1}; Time = Data(1); nx = Data(2); ny = Data(3); Data_vec = Data(4:end); % define coordinates dx = xsize/(nx-1); dy = ysize/(ny-1); gx = 0:dx:xsize; gy = 0:dy:ysize; s0 = dx*dy; % the area of grid (km^2) if(zoom>0) % find beginning/end index % x direction indx1 = find(gx==xbeg); if(length(indx1)==0) indx1 = find(gx>xbeg,1); end indx2 = find(gx==xend); if(length(indx2)==0) indx2 = find(gx>xend,1); end % y direction indy1 = find(gy==ybeg); if(length(indy1)==0) indy1 = find(gy>ybeg,1); end indy2 = find(gy==yend); if(length(indy2)==0) indy2 = find(gy>yend,1); end % zoomed area coordinates gx = gx(indx1:indx2); gy = gy(indy1:indy2); end % initialize composition matrix ColorGrid = ones(ny,nx); % colormap for composition C_map=[1 1 1;... %0 0.50588 0.99608 0.78824;... %1 0.4754 0.251 0.2;... %2 0.68235 0.34118 0;... %3 0.75294 0.75294 0.75294;... %4 101/250 102/250 101/250;... %5 139/250 140/250 140/250;... %6 13/250 137/250 97/250;... %7 0.392157 0.6 0;... %8 184/250 206/250 185/250;... %9 227/250 240/250 221/250;... %10 120/250 178/250 177/250;... %11 12/250 117/250 134/250;... %12 25/250 25/250 112/250;... %13 0.7 0.12 0.3;... %14 0.4 0.4 0.0;... %15 0.2 0.6 0;... %16 13/250 137/250 97/250;... %17 0.3 0.5 0.6;... %18 0 0 0;... %19 0 0 0;... %20 0 0 0;... %21 0 0 0;... %22 1 1 0.31765;... %23 1 0.90196 0.18824;... %24 0.46667 0.46667 0.23529;... %25 0.50196 0.50196 0;... %26 0.92549 0.43922 0.99608;... %27 0.92549 0.43922 0.99608;... %28 1 1 0;... %29 1 1 0;... %30 0.99216 0.38824 0.30196;... %31 0 0 0;... %32 0 0 0;... %33 1 0 0;... %34 1 0 0;... %35 1 0 0;... %36 0.92549 0.43922 0.99608;... %37 0.92549 0.43922 0.99608]; %38 %Fill the colorgrid with the data of opened file num = 1; ind = 1; while num0) ColorGrid = ColorGrid(indy1:indy2,indx1:indx2); end %save memory clear Data_vec time_all(idx)=Time/1e6; % -------temperature ------- if(10) format long g fname_t = [name_t,num2str(nn1(idx)) ext]; fid = fopen (fname_t,'r'); B = textscan (fid, '%n')'; fclose(fid); Data = B{1}; nx_t = Data(2); ny_t = Data(3); Data_vec_t= Data(4:end); Data_vec_t= Data_vec_t-273; % form Kelvin to Celcius ColorGrid_t = ones(ny_t,nx_t); ColorGrid_t = reshape(Data_vec_t,ny_t,nx_t); % define coordinates for T dx_t = xsize/(nx_t-1); dy_t = ysize/(ny_t-1); gx_t = 0:dx_t:xsize; gy_t = 0:dy_t:ysize; if(zoom>0) % find beginning/end index % x direction indx1_t = find(gx_t==xbeg); if(length(indx1_t)==0) indx1_t = find(gx_t>xbeg,1); end indx2_t = find(gx_t==xend); if(length(indx2_t)==0) indx2_t = find(gx_t>xend,1); end % y direction indy1_t = find(gy_t==ybeg); if(length(indy1_t)==0) indy1_t = find(gy_t>ybeg,1); end indy2_t = find(gy_t==yend); if(length(indy2_t)==0) indy2_t = find(gy_t>yend,1); end % zoomed area coordinates gx_t = gx_t(indx1_t:indx2_t); gy_t = gy_t(indy1_t:indy2_t); end % zoom in T if(zoom>0) ColorGrid_t = ColorGrid_t(indy1_t:indy2_t,indx1_t:indx2_t); end end name2 = 'RigSub'; cmyr=100*365*24*3600; % read data icount = 0; ti = 1; for idx = 1:length(nn1) icount = icount+1 close all cd(dir1) fileprn = [name2,num2str(nn1(idx)),'.prn']; % Opening data file fdata=fopen(fileprn,'rb'); % Read sizes of variables A=fread(fdata,4,'uchar'); % Read model parameters % Grid resolution xnumx=fread(fdata,1,'int64'); ynumy=fread(fdata,1,'int64'); % Markers per cell mnumx=fread(fdata,1,'int64'); mnumy=fread(fdata,1,'int64'); % Number of markers marknum=fread(fdata,1,'int64'); % Model sizes xsize=fread(fdata,1,'float64'); ysize=fread(fdata,1,'float64'); % Pressure value pinit=fread(fdata,5,'float64'); % Gravity GXKOEF=fread(fdata,1,'float64'); GYKOEF=fread(fdata,1,'float64'); % Number of rocks rocknum=fread(fdata,1,'int32'); % Number of Boundary conditions bondnum=fread(fdata,1,'int64'); % Stage,time n1=fread(fdata,1,'int32'); timesum=fread(fdata,1,'float64'); if(timesum>1e3) time=sprintf('%3.3f',timesum/1e6); time=str2num(time); else time=timesum/1e6; end time_all(icount) = time; % Ma % Skip general information curpos0 = 4*1; % skip grid parameters curpos1 = 12*8 + 2*8 + 4 + 8 + 4 + 8; % Skip rock informations curpos2 = rocknum*(1*4 + 24*8); % Skip nodes information curpos3 = xnumx*ynumy*(3*4 + 3*8 + 16*4 + 1*8 + 5*4); %change % skip gridline positions curpos4 = (xnumx + ynumy)*4; % skip boundary conditions equations curpos5 = (bondnum-1)*(4*4 + 3*8); % skip marker infos curpos6 = (marknum+1)*(9*4+1); % total value should equals to the real prn size (in byte) curpos = curpos0 + curpos1 + curpos2 + curpos3 + curpos4 + curpos5 + curpos6; % initialitialise parameter sum_nu = 0; % mean effective viscosity sum_eii = 0; %------- ------------load gridlines --------------------------- curpos_grid = curpos0 + curpos1 + curpos2 + curpos3; fseek(fdata,curpos_grid,'bof') gx_gy = fread(fdata,xnumx+ynumy,'float32'); gx2 = gx_gy(1:xnumx)/1000; gy2 = gx_gy(xnumx+1:end)/1000; % ------------ fast: write data into a 2d matrix ------------- % find the starting position % sxx curpos_sxx = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 3*4; fseek(fdata,curpos_sxx,'bof'); % only read sxx, skip useless information skip1 = 3*4 + 3*8 + 16*4 + 1*8 + 4*4; %change % 2d matrix sxx = fread(fdata, [ynumy,xnumx],'float32',skip1); % syy curpos_syy = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 4*4; fseek(fdata,curpos_syy,'bof'); % only read sxx, skip useless information skip_syy = skip1; % 2d matrix syy = fread(fdata, [ynumy,xnumx],'float32',skip_syy); % sxy curpos_sxy = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 5*4; fseek(fdata,curpos_sxy,'bof'); skip2 = skip1; sxy = fread(fdata, [ynumy,xnumx],'float32',skip2); % exx curpos_exx = curpos0 + curpos1 + curpos2 + 3*4+3*8; fseek(fdata,curpos_exx,'bof'); skip3 = skip1; exx = fread(fdata, [ynumy,xnumx],'float32',skip3); % eyy curpos_eyy = curpos0 + curpos1 + curpos2 + 3*4+3*8 +4; fseek(fdata,curpos_eyy,'bof'); skip3 = skip1; eyy = fread(fdata, [ynumy,xnumx],'float32',skip3); % exy curpos_exy = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 2*4; fseek(fdata,curpos_exy,'bof'); skip4 = skip1; exy = fread(fdata, [ynumy,xnumx],'float32',skip4); % vx curpos_vx = curpos0 + curpos1 + curpos2 + 4; fseek(fdata,curpos_vx,'bof'); skip3 = skip1; vx = fread(fdata, [ynumy,xnumx],'float32',skip3); vx = vx*cmyr; % vy curpos_vy = curpos0 + curpos1 + curpos2 + 2*4; fseek(fdata,curpos_vy,'bof'); skip3 = skip1; vy = fread(fdata, [ynumy,xnumx],'float32',skip3); vy = vy*cmyr; % RMS velocity rms_vel(icount) = ((sum(sum(vx.^2)) + sum(sum(vy.^2)))/(2*xnumx*ynumy))^0.5; rms_vx(icount) = ((sum(sum(vx.^2)))/(xnumx*ynumy))^0.5; rms_vy(icount) = ((sum(sum(vy.^2)))/(xnumx*ynumy))^0.5; % shear heating Hs = 2*sxx.*exx + 2*sxy.*exy; aa = find(Hs<1e-30); Hs(aa) = 1e-30; % density curpos_dens = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 6*4; fseek(fdata,curpos_dens,'bof'); dens = fread(fdata, [ynumy,xnumx],'float32',skip1); % second invariant sii = (sxx.^2 + sxy.^2).^(1/2); eii = (exx.^2*0.5 + eyy.^2*0.5 + exy.^2).^(1/2); eii(find(eii<=0)) = 1e-30; % viscosity curpos_nu = curpos0 + curpos1 + curpos2 + 3*4+3*8 + 7*4; fseek(fdata,curpos_nu,'bof'); skip5 = skip1; nu = fread(fdata, [ynumy,xnumx],'float32',skip5); % pr(pressure) curpos_pr = curpos0 + curpos1 + curpos2; fseek(fdata,curpos_pr,'bof'); %only read sxx, skip useless information skip3 = skip1;%跳转到下一个节点该数据的位置进行读取 %2d matrix pr = fread(fdata,[ynumy,xnumx],'float32',skip3); % temperature curpos_T = curpos0 + curpos1 + curpos2 + 3*4+3*8+ 15*4; fseek(fdata,curpos_T,'bof'); T = fread(fdata, [ynumy,xnumx],'float32',skip1); T=T-273; % mext curpos_mext = curpos_T + 8+4*4; %change fseek(fdata,curpos_mext,'bof'); %change mext = fread(fdata, [ynumy,xnumx],'float32',skip1); %change % meltfrac curpos_meltfrac = curpos_T +8+5*4; %change fseek(fdata,curpos_meltfrac,'bof'); %change meltfract = fread(fdata, [ynumy,xnumx],'float32',skip1); %change % mean viscosity % entire box mean_nu_box(idx) = mean(mean(nu)); % mantle lithosphere nodecount1 = find(T<=1300); nodecount2 = find(dens>3000); nodecount3 = intersect(nodecount1,nodecount2); mean_nu_lith(idx) = mean(nu(nodecount3)); % asthenosphere nodecount4 = find(T>1300); mean_nu_asth(idx) = mean(nu(nodecount4)); % ----------------------- Plot ------------------------- lw = 1.2; ft = 12; % ------------------------ comp plot ------------------------ if(10) cd (dir2) figure(1) step=1; colormap(C_map); pcolor(gx,gy,ColorGrid); shading interp if(10) hold on [C,h] = contour(gx_t(1:step:end),gy_t(1:step:end),ColorGrid_t(1:step:end,1:step:end),... [500:400:1300],'color','[0.5 0.5 0.5]','LineWidth',0.2); text_handle = clabel(C,h,'color','k','labelspacing',10000,'fontsize',6); hold off end if(10) step1=15; hold on; yyid=190; scale=2; h1=quiver(gx2(70:step1:430), gy2(20:step1:yyid), vx(20:step1:yyid,70:step1:430), vy(20:step1:yyid,70:step1:430),0,'k','linewidth',0.5); hu=get(h1,'UData'); hv=get(h1,'VData'); set(h1,'UData',scale*hu,'VData',scale*hv); hold off; end caxis ([0 38]) axis ij equal tight xlim ([1800 3000]) ylim ([0 400]) set(gca,'xtick',1800:200:3000); title(['Time = ', num2str(Time/1e6),' Ma'],'fontsize',ft) set(gca,'position',[0 0 1 1]) if zoom >= 0 filename = ['Comp_z_prn',num2str(nn1(idx)),'.png']; set(gca,'layer','top','position',[0.07 0.08 0.9 0.9],'linewidth',0.01) else filename = ['Comp_prn',num2str(nn1(idx))]; set(gca,'layer','top','linewidth',lw,'fontsize',ft) end print ('-dpng', '-r800',filename); end end end