pythonFLAC/magma_try.py at main · chingchenn/pythonFLAC · GitHub

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
#!/usr/bin/env python
import math
import flac
import os,sys
import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt
model = str(sys.argv[1])
#path = '/home/jiching/geoflac/'+model+'/'
path = '/Users/ji-chingchen/Desktop/model/'+model+'/'
#path = model
print(model)
os.chdir(path)
fl = flac.FlacFromVTK();end = fl.nrec

def get_magma(start_vts=1,model_steps=end-1):
    melt=np.zeros(end)
    magma=np.zeros(end)
    yymelt=np.zeros(end)
    yychamber=np.zeros(end)
    arc_vol=np.zeros(end)
    for i in range(1,end):
        x,z=fl.read_mesh(i)
        phase = fl.read_phase(i)
        mm=fl.read_fmelt(i)
        chamber=fl.read_fmagma(i)
        melt[i] = np.max(mm)
        magma[i] = np.max(chamber)
        arc_vol[i]=np.sum(fl.read_area(i)[phase ==14])/1e6
        yymelt[i]=(fl.read_fmelt(i)*fl.read_area(i)/1e6).sum()
        yychamber[i]=(fl.read_fmagma(i)*fl.read_area(i)/1e6).sum()
    return melt,magma,yymelt,yychamber,arc_vol
def read_time(start_vts,model_steps):
    timestep=[0]
    for step in range(start_vts,model_steps):
        timestep.append(fl.time[step])
    # timestep=np.array(timestep)
    return timestep
time =fl.time
#----------------------------------------------------------------------------
#lam0 = 1e-12
#lam_tdep = 4e-2
#delT=1440 #temperature diffence from surface to depth
#lam = lam0*(1+np.exp(lam_tdep*delT))
#prod = 6e-13
#total_magma = prod*dt - fmagma*np.exp(-lam*dt)
#----------------------------------------------------------------------------
melt,magma,yymelt,yychamber,arc_vol=get_magma(1,end)
fig, (ax,ax2,ax3,ax4,ax5) = plt.subplots(5,1,figsize=(15,17))
ax.bar(time,yymelt,width=0.1,color='tomato')
ax2.plot(time,yychamber,color='orange')
ax3.bar(time,melt,width=0.1,color='tomato',label='fmelt')
ax4.plot(time,magma,color='orange')
ax5.plot(time,arc_vol,color='orange',label='magma')
#ax.set_xlabel('Time (Myr)',fontsize=20)
#ax2.set_xlabel('Time (Myr)',fontsize=20)
#ax3.set_xlabel('Time (Myr)',fontsize=20)
ax5.set_xlabel('Time (Myr)',fontsize=20)
ax.set_ylabel('melt * area',fontsize=20)
ax2.set_ylabel('chamber *area',fontsize=20)
ax3.set_ylabel('max fmelt',fontsize=20)
ax4.set_ylabel('max fmagma',fontsize=20)
ax5.set_ylabel('arc area',fontsize=20)
#ax.set_ylim(0,0.8)
#ax2.set_ylim(0,10*1e-3)
#ax3.set_ylim(0,10*1e-3)
#ax4.set_ylim(0,3*1e-5)
ax.set_xlim(0,24)
ax2.set_xlim(0,24)
ax3.set_xlim(0,24)
ax4.set_xlim(0,24)
ax5.set_xlim(0,24)
ax.grid()
ax2.grid()
ax3.grid()
ax4.grid()
ax5.grid()
ax.tick_params(axis='x', labelsize=16 )
ax2.tick_params(axis='x', labelsize=16 )
ax3.tick_params(axis='x', labelsize=16 )
ax4.tick_params(axis='x', labelsize=16 )
ax5.tick_params(axis='x', labelsize=16 )
ax.tick_params(axis='y', labelsize=16 )
ax2.tick_params(axis='y', labelsize=16 )
ax3.tick_params(axis='y', labelsize=16 )
ax4.tick_params(axis='y', labelsize=16 )
ax5.tick_params(axis='y', labelsize=16 )
ax.set_title('Model : '+model,fontsize=25)
#fig.savefig('/home/jiching/geoflac/figure/'+model+'_arc+magma.png')
fig.savefig('/Users/ji-chingchen/Desktop/figure/'+model+'_arc+magma.png')