Commit 0befa827 authored by Lena Noack's avatar Lena Noack
Browse files

Update plot_snap combined with python routine

parent 44403c63
% Sara Vulpius, 2019
%
%
function [deg_CO2,deg_H2O] = get_sol(X_H2O_mantle,X_CO2_mantle,p_lid,T_lid)
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = X_H2O_mantle*10^-4; %ppm->weight percent of water in the melt
X_CO2 = X_CO2_mantle; %weight ppm of CO2 in the melt
% if (X_H2O<4) % personal communication with Fabrice Gaillard
hydrous=0;
% else
% hydrous=1;
% end
M_H = 1.0079;
M_C = 12.0107;
M_O = 15.9994;
mass_CO2 = M_C + 2*M_O;
mass_H2O = 2*M_H + M_O;
H2O_fl_mole_pred = 0;
H2O_dissolved = 0;
CO2_dissolved = 0;
H2O_dissolved_old = 1;
CO2_dissolved_old = 1;
dX_CO2 = 0;
dX_H2O = 0;
dX_H2O_mole = 0;
dX_CO2_mole = 0;
tol = 1e-10;
tol_fl_mole = 1e-4;
i_max = 10000; %max anzahl von iterationen
no_water = 0; % if water pressure goes to zero, set no_water to one, otherwise wrong oscillating solution
i=1;
fac=1;
fac_mult=0.9;
redo=1;
outer_steps=10;
j=1;
%do
H2O_fl_mole = 0.5; % = P_H2O/P_lid start value, exact one determined by iteration below
H2O_fl_mole_min = 0;
H2O_fl_mole_max = 1;
H2O_fl_mole_last = 0;
%fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',0,...
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,0)
dX_H2O = max(0,X_H2O-H2O_dissolved); %dX_H2O_lid = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved); %dX_CO2_lid = max(0,X_CO2-CO2_dissolved);
%while (redo>0 && j<outer_steps)
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol || abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
H2O_dissolved_old = H2O_dissolved;
CO2_dissolved_old = CO2_dissolved;
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,X_H2O,hydrous);
dX_H2O = max(0,X_H2O-H2O_dissolved); % what would go into gas bubble in wt-%
dX_CO2 = max(0,X_CO2-CO2_dissolved); % what would go into gas bubble in wt-ppm
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=X_CO2;, end %CO2_dissolved_old;, end
if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=X_H2O;, end %H2O_dissolved_old;, end
if (dX_H2O+dX_CO2>0)&&(no_water==0)
dX_H2O_mole = dX_H2O/mass_H2O/10^2 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
dX_CO2_mole = dX_CO2/mass_CO2/10^6 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
fac=fac*fac_mult; % update only by 10% to avoid oscillations between 0 and 100%
if (dX_H2O_mole+dX_CO2_mole==0)
H2O_fl_mole_pred = 0;
else
H2O_fl_mole_pred = dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole);
end
H2O_fl_mole = (1-fac)*H2O_fl_mole+fac*(H2O_fl_mole_pred);
else
no_water=1;
H2O_fl_mole = 0;% neither CO2 nor H2O go out, so set to smaller value to see if at least some CO2 goes out?
end
% if H2O_fl_mole>H2O_fl_mole_max || H2O_fl_mole<H2O_fl_mole_min
% H2O_fl_mole = 0.5*(H2O_fl_mole_max+H2O_fl_mole_min);
% end
% if H2O_fl_mole>H2O_fl_mole_last
% H2O_fl_mole_min = H2O_fl_mole_last;
% else
% H2O_fl_mole_max = H2O_fl_mole_last;
% end
%fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',i,...
%H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole))
dX_H2O = max(0,X_H2O-H2O_dissolved); %dX_H2O_lid = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved); %dX_CO2_lid = max(0,X_CO2-CO2_dissolved);
i=i+1;
%until ((abs(H2O_dissolved-H2O_dissolved_old)<tol && abs(CO2_dissolved-CO2_dissolved_old)<tol) || (i>i_max))
end
% fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',i,...
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,H2O_fl_mole_pred)
%
% if abs(H2O_fl_mole-(H2O_fl_mole_pred)) > tol_fl_mole % above did actually not converge, only time steps became so small that it looks like it converged
% fac_mult = 1-(1-fac_mult)*0.5;
% H2O_fl_mole_pred = 0;
% H2O_dissolved = 0;
% CO2_dissolved = 0;
% H2O_dissolved_old = 1;
% CO2_dissolved_old = 1;
% dX_CO2 = 0;
% dX_H2O = 0;
% i=1;
% fac=1;
% H2O_fl_mole = 0.001;
% else
% redo = 0; % converged
% end
% j=j+1;
%end
%
%if j==outer_steps || i==i_max
% fprintf('Error - did not converge, check get_sol.m, use more outer/inner iterations or change tolerance value\n')
%end
% if (i>i_max)
% %fprintf('Warning, iteration stopped after %d iterations at p=%f bar!\n',i,p_lid)
% end
deg_CO2 = dX_CO2; % wt-ppm
deg_H2O = dX_H2O*10^4; % wt-% -> wt-ppm
end
\ No newline at end of file
function make_movie(start=0,fpath=' ')
%function make_movie(start=0,fpath=' ',vals=' ') for normal plots
function make_movie(start=0,fpath=' ',vals='TWHDvRM')
if fpath==' '
[fname, fpath, fltidx] = uigetfile("data_char_val_ts.res");
endif
filenames = sprintf('%s%s', fpath, '/data_snap*')
filenames = sprintf('%s%s', fpath, 'data_snap*')
files = glob(filenames);
files = glob(filenames)
for i=1:numel(files)
[~, name] = fileparts(files{i});
nr=str2num(name(14:end));
% "File:",files(i)
file_name = strrep(files{i}, '\', '/'); % needed under Windows
if nr>=start
eval(sprintf('plot_snap(30,4500,-1,0,3,"%s", "-ascii");', files{i}));
eval(sprintf('plot_snap(30,0,-1,0,1,"%s","%s", "-ascii");', file_name, vals));
# eval(sprintf('plot_snap(30,3000,-1,0,1,"%s", "-ascii");', file_name));%files{i}));
endif
%% eval(sprintf('%s = load("%s", "-ascii");', name, files{i}));
endfor
......
This diff is collapsed.
This diff is collapsed.
No preview for this file type
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import csv
import math
import matplotlib.cm as cm
import matplotlib.colors as colors
import os, sys
# Below use call create_plots:
#create_plots(path) # call without file -> takes all PlotData* files
#create_plots(path,file) # call with file -> takes only this one file
```
%% Cell type:code id: tags:
``` python
def plot2D(xx,yy,data,title,time,colmap,vmin=None,vmax=None,nr_cont=10,name=None) :
y_scale = 5
fig = plt.figure(figsize=(2*y_scale+4.5,y_scale+1)) #4.5 for colorbar, 1 for figure title
ax = fig.add_subplot(111)
# plt.clf()
cmap = plt.get_cmap(colmap)
m = plt.cm.ScalarMappable(cmap=cmap)
cp = ax.contourf(xx, yy, data,nr_cont,cmap=cmap, vmin=vmin, vmax=vmax)
m.set_array(data)
m.set_clim(vmin, vmax)
plt.colorbar(m)
plt.axis('off')
# ax.set_title(title+' in '+str(time)+' Myr', fontsize=20)
time_str = str("{:2.1f}".format(float(time))) # time already as string but needs to be cut to lower precision
ax.set_title(title+' at '+time_str+' Myr', fontsize=20)
# plt.show()
if name!=None:
plt.savefig(name, format='png', bbox_inches='tight',dpi=100)
plt.close(fig)
#######################################################################
#######################################################################
def make_plot_snap(filename):
letters_number = filename[9:-4]
#print(letters_number)
letters = letters_number[:-7]
number = letters_number[-7:] # stays string for plot names
#print(letters)
#print(int(number))
file1 = open(path+filename, 'r')
time = file1.readline() # time in Myr
xx_dim = file1.readline() # (x,y) dim of xx-coordinates
xx_val = file1.readline() # xx-coordinates
yy_dim = file1.readline() # (x,y) dim of yy-coordinates
yy_val = file1.readline() # yy-coordinates
xx_dim = xx_dim.split() # get two strings with integers each
xd = int(xx_dim[0])
yd = int(xx_dim[1])
xx_val = np.reshape(np.array(xx_val.split()).astype(float),(yd,xd))
yy_val = np.reshape(np.array(yy_val.split()).astype(float),(yd,xd))
def set_plot_values(title,cm,vmin=None,vmax=None,nr_cont=10,nameF=None):
f_dim = file1.readline() # dimension of field
f_val = file1.readline() # values as list
f_val = np.reshape(np.array(f_val.split()).astype(float),(yd,xd))
plot2D(xx_val,yy_val,f_val,title,time,cm,vmin=vmin,vmax=vmax,nr_cont=nr_cont,name=plot_path+nameF+number+'.png')
for letter in letters:
if letter == 'T':
set_plot_values('Temperature in K','magma',vmin=0,vmax=3500,nr_cont=20,nameF='T2D')
elif letter == 'D':
set_plot_values('Depletion','inferno',vmin=0,vmax=30,nr_cont=50,nameF='D2D')
elif letter == 'H':
set_plot_values('Heat sources','rainbow',vmin=0,vmax=80,nr_cont=100,nameF='H2D')
elif letter == 'W':
set_plot_values('Water fraction in ppm','RdBu',vmin=0,vmax=140,nr_cont=100,nameF='W2D')
elif letter == 'v':
set_plot_values('Velocity in cm/yr','rainbow',vmin=0,vmax=25,nr_cont=100,nameF='v2D')
elif letter == 'R':
# set_plot_values('Material parameters','rainbow',vmin=3,vmax=7,nr_cont=4,nameF='R2D')
# set_plot_values('Material parameters','Set2',vmin=3,vmax=7,nr_cont=4,nameF='R2D')
set_plot_values('Material parameters','Paired',vmin=3,vmax=7,nr_cont=4,nameF='R2D')
elif letter == 'M':
set_plot_values(r'Radial material exchange in $kg/(m^2 yr)$','inferno',nr_cont=50,nameF='M2D')
elif letter == 't':
set_plot_values('Temperature variations in K','RdBu',vmin=-100,vmax=100,nr_cont=20,nameF='dT2D')
else:
print('Unused letter found (',letter,')')
file1.close()
```
%% Cell type:code id: tags:
``` python
def create_plots(path,file=None,start=0,end=10000000000):
plot_path = path+'plot/'
if os.path.isdir(plot_path)==False:
os.mkdir( plot_path );
os.chdir(path)
if file==None:
arr = [x for x in os.listdir() if x.startswith("PlotData")]
# print(arr)
for filename in arr:
letters_number = filename[9:-4]
number = int(letters_number[-7:])
if (number>=start and number<=end):
print(filename)
make_plot_snap(filename)
else:
make_plot_snap(file)
```
%% Cell type:code id: tags:
``` python
path = 'C:/Arbeit/Video_Venus_MOprofs/'
#path = 'C:/Arbeit/Video_Venus_MOprofs_PhaseTrans/'
#path = 'C:/Arbeit/Video_Venus_MOprofs_PhaseTrans_NoChemBuoy/'
plot_path = path+'plot/'
#file1.close()
create_plots(path)
#create_plots(path,start=3292,end=4331)
#create_plots(path,start=54616) # call without file -> takes all PlotData* files
#create_plots(path,'PlotData_TWHDvRMt0011564.txt') # call with file -> takes only this one file
```
%% Output
PlotData_TWHDvRMt0000001.txt
PlotData_TWHDvRMt0000132.txt
PlotData_TWHDvRMt0000267.txt
PlotData_TWHDvRMt0000426.txt
PlotData_TWHDvRMt0000721.txt
PlotData_TWHDvRMt0001510.txt
PlotData_TWHDvRMt0002266.txt
PlotData_TWHDvRMt0002701.txt
PlotData_TWHDvRMt0003031.txt
PlotData_TWHDvRMt0003292.txt
PlotData_TWHDvRMt0003503.txt
PlotData_TWHDvRMt0003689.txt
PlotData_TWHDvRMt0003853.txt
PlotData_TWHDvRMt0004035.txt
PlotData_TWHDvRMt0004185.txt
PlotData_TWHDvRMt0004332.txt
PlotData_TWHDvRMt0004484.txt
PlotData_TWHDvRMt0004642.txt
PlotData_TWHDvRMt0004807.txt
PlotData_TWHDvRMt0004996.txt
PlotData_TWHDvRMt0005177.txt
PlotData_TWHDvRMt0005333.txt
PlotData_TWHDvRMt0005468.txt
PlotData_TWHDvRMt0005593.txt
PlotData_TWHDvRMt0005744.txt
PlotData_TWHDvRMt0005898.txt
PlotData_TWHDvRMt0006052.txt
PlotData_TWHDvRMt0006197.txt
PlotData_TWHDvRMt0006338.txt
PlotData_TWHDvRMt0006481.txt
PlotData_TWHDvRMt0006622.txt
PlotData_TWHDvRMt0006756.txt
PlotData_TWHDvRMt0006893.txt
PlotData_TWHDvRMt0007022.txt
PlotData_TWHDvRMt0007154.txt
PlotData_TWHDvRMt0007290.txt
PlotData_TWHDvRMt0007439.txt
PlotData_TWHDvRMt0007589.txt
PlotData_TWHDvRMt0007746.txt
PlotData_TWHDvRMt0007885.txt
PlotData_TWHDvRMt0008016.txt
PlotData_TWHDvRMt0008149.txt
PlotData_TWHDvRMt0008287.txt
PlotData_TWHDvRMt0008430.txt
PlotData_TWHDvRMt0008573.txt
PlotData_TWHDvRMt0008712.txt
PlotData_TWHDvRMt0008845.txt
PlotData_TWHDvRMt0008981.txt
PlotData_TWHDvRMt0009108.txt