-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmodelisation.py
243 lines (185 loc) · 8.64 KB
/
modelisation.py
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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
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
#import math
#from scipy.interpolate import griddata
import numpy as np
import os
import time
#import matplotlib.pyplot as plt
#import sys, h5py, binascii
#import tempfile
#import pygimli as pg
#from pygimli.meshtools import appendTriangleBoundary, merge2Meshes, mesh
import warnings
#Pour aider pyflakes à analyser il vaut mieux importer les fonctions explicitement
from maillage_SWMS2D_EL import maillage_SWMS2D_EL
from maillage_SWMS2D import maillage_SWMS2D
from initial_conditions import initial_conditions
from initial_conditions_EL import initial_conditions_EL
from ecriture_fichiers_SWMS2D import ecriture_Selector_in, ecriture_Grid_in
from ecriture_fichiers_SWMS2D_EL import ecriture_Selector_in, ecriture_Grid_in_EL
from ecriture_fichiers_GPRMAX import ecriture_fichiers_GPRMAX
from ecriture_fichiers_GPRMAX_EL import ecriture_fichiers_GPRMAX_EL
from maillage_GPRMAX_EL import maillage_GPRMAX_EL
from outils import showQuality
from pygimli.meshtools import quality
#from picking_radargramme import picking
from param_acquisition import longueur_d_onde
#from joblib import Parallel, delayed
#import subprocess
def fxn():
warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
fxn()
#Def des paramètres de géométrie du modéle
# class Geometry :
# def __repr__(self):
# return "dtrou{dtrou}_rtrou{r}_tr{h_eau}".format(dtrou=self.dtrou,r=self.r,h_eau=self.h_eau)
# # Definition des paramètres MVG
# class ParamMVG(namedtuple("ParamMVG", ["ts", "ti", "tr", "n", "alpha", "Ks"])):
# #Recalcul des h0 à partir des param MVG
# def h0(self):
# return -(1/self.alpha)*(((self.ti-self.tr)/(self.ts-self.tr))**(self.n/(1-self.n))-1)**(1/self.n)
# def __repr__(self):
# return "ts{ts}_ti{ti}_tr{tr}_n{n}_alpha{alpha}_Ks{Ks}".format(ts=self.ts,ti=self.ti,tr=self.tr,n=self.n,alpha=self.alpha,Ks=self.Ks)
# # Definition des param gprMax
# class ParamGPRMAX :
# pass
# def longueur_d_onde(theta, freq, paramMVG, paramGPRMAX):
# """Compute the wavelength for a given water content and frequency (Hz)"""
# eps = CRIM(theta,paramMVG,paramGPRMAX)
# vitesse = 0.3 / math.sqrt(eps) # m/ns
# return vitesse*(10**9)/freq
def run(geometry,paramMVG,paramGPRMAX,temps,tmax_SWMS2D):
"""
Run a simulation on the model with the given parameters
The input and output files are created in a directory
myDirName whose name depends on the parameters.
If the simulation was already computed, return immediately.
Return the directory containing the results of the simulation
"""
# try: c'est pour pouvoir prendre la main si ça foire.
#try:
myDirName = "OUTProf"+repr(geometry)+"/"+repr(paramMVG)
nom='radargram'
filename = nom + '__merged.out'
dir = os.getcwd() #répertoire où on a lancer le script
if os.path.isfile(myDirName+"/"+filename) :
print('already simulated')
return myDirName
os.makedirs(myDirName,exist_ok=True)
os.chdir(myDirName)
fparam=open("Parameters","w")
fparam.write("""ParamMVG(tr={tr},ti={ti},ts={ts},n={n},alpha={alpha},Ks={Ks})\n""".format(ts=paramMVG.ts,ti=paramMVG.ti,tr=paramMVG.tr,n=paramMVG.n,alpha=paramMVG.alpha,Ks=paramMVG.Ks))
fparam.close()
os.system("cp "+dir+"/gprMax .")
os.system("cp "+dir+"/gprMaxMerge .")
# ce serait mieux de copier l'executable H2D dans le dossier où l'on travaille mais il est gros!
#os.system("cp "+dir+"/HD2/H2D .")
os.makedirs("SWMS_2D.IN",exist_ok=True)
os.makedirs("SWMS_2D.OUT",exist_ok=True)
#Définition du maillage triangulaire pour SWMS2D
#[mesh, pg_pos, mesh_pos, mesh_cells]=maillage_SWMS2D(geometry)
[mesh, pg_pos, mesh_pos, mesh_cells]=maillage_SWMS2D_EL(geometry)
#from pygimli.meshtools import mesh
#from pygimli.mplviewer import drawMesh
#from pygimli.viewer import showMesh
#showMesh(mesh,markers=True)
#matplotlib.use('Agg')
#figi=showQuality(mesh, quality(mesh))
#figi.savefig('mesh.png',format='png')
#Calcul des charges initiales en chaque noeud du maillage
p=initial_conditions_EL(mesh_pos, geometry, paramMVG)
#Création des fichiers .in nécessaires pour SMWS_2D
ecriture_Selector_in(mesh, paramMVG, temps, p)
ecriture_Grid_in_EL(mesh, p)
os.system("mv Grid.in SWMS_2D.IN")
os.system("mv Selector.in SWMS_2D.IN")
#Lancement SWMS_2D
start_running = time.time()
end_running=[]
error=os.system("timeout {} {}/HD2/H2D".format(tmax_SWMS2D,dir))
if error : #reagira seulement si error est différent de 0
print("HD2 fut trop long")
os.chdir(dir)
return myDirName
end_running.append(time.time()-start_running)
#os.popen("rm -rf *.in")
print('Running Hydrus time:'+str(end_running))
#Fichier contenant les thetas
f_thetas = "SWMS_2D.OUT/th.out"
#Number of traces to be calculated
nT=len(temps)
#Interpolation du maillage triangulaire sur une grille rectangulaire pour gprMax
#On crée un maillage rectangulaire avec les dimensions du modèle
[xv, yv, mx, my, mesh2, grid, grid_mat, eps_mat, sigma_grid_mat] = maillage_GPRMAX_EL(paramGPRMAX, paramMVG, mesh, mesh_pos[:,:2], f_thetas, nT)
# Lancement gprMax
# pas de calcul spatial le plus grossier possible (m)
dl = longueur_d_onde(paramMVG.ts,paramGPRMAX.freq_max,paramMVG,paramGPRMAX)/paramGPRMAX.spatial_step
materiaux = {}
A_tab={}
def materiau(x, sigma): # TODO: A regarder car je n'y comprend rien!
if x in materiaux: # si un materiau existe déjà pour x, il le renvoie, donc il n'en crée pas un nouveau.
return materiaux[x]
valeur1 = "sand{}".format(len(materiaux))
valeur2 = sigma
materiaux[x] = valeur1, valeur2
return valeur1, valeur2
xreg = np.arange(paramGPRMAX.xmin, paramGPRMAX.xmax + paramGPRMAX.dx, paramGPRMAX.dx, 'float')
zreg = np.arange(paramGPRMAX.zmin, paramGPRMAX.zmax + paramGPRMAX.dx, paramGPRMAX.dx, 'float')
end_running=[]
start_running = time.time()
# =============================================================================
# def rungprmax(ite,xreg,zreg,grid_mat,sigma_grid_mat,xv,yv,nom, paramMVG, paramGPRMAX, geometry, dl,materiaux):
# from ecriture_fichiers_GPRMAX import ecriture_fichiers_GPRMAX
#
# for j in range(0,len(zreg)):
# for k in range(0,len(xreg)):
# materiau(grid_mat[ite][j,k], sigma_grid_mat[ite][j,k])
#
# A = ecriture_fichiers_GPRMAX(xv.T*0.01, yv.T*0.01, grid_mat[ite], ite, nom, paramMVG, paramGPRMAX, geometry, dl, materiaux)
# fichier=nom+'_'+str(ite+1)+'.in'
#
# command="./gprMax "+fichier
# os.popen(command).readlines()
# #subprocess.run(command)
# Parallel(n_jobs=7,verbose=10)(delayed(rungprmax)(ite,xreg,zreg,grid_mat,sigma_grid_mat,xv,yv,nom, paramMVG, paramGPRMAX, geometry, dl,materiaux) for ite in range(0,nT+1))
#
# =============================================================================
for i in range(0,nT+1):
#start_ecriture = time.time()
for j in range(0,len(zreg)):
for k in range(0,len(xreg)):
materiau(grid_mat[i][j,k], sigma_grid_mat[i][j,k])
# Writing a file.in for running GPRMAX
A = ecriture_fichiers_GPRMAX_EL(xv.T*0.01, yv.T*0.01, grid_mat[i], i, nom, paramMVG, paramGPRMAX, geometry, dl, materiaux)
# plt.close('all')
# fig = plt.figure()
# ax = fig.add_subplot(111)
# gni=ax.scatter(A[:,0],A[:,1],s=30,c=A[:,2])
# cbar=fig.colorbar(gni,label='Eps' )
# cbar.minorticks_on()
# fig.savefig('Gnard'+str(i)+'.png',format='png')
#A_tab[i]=A
#end_ecriture.append(time.time()-start_ecriture)
#Lancement calcul gprMax
fichier=nom+'_'+str(i+1)+'.in'
command="./gprMax "+fichier
os.popen(command).readlines()
# Concatenate all nT traces
# fig = plt.figure()
# ax1 = fig.add_subplot(111)
# ax2 = ax1.twinx()
# ax1.plot(list(range(0,nT+1)),end_running,c='r')
# ax2.plot(list(range(0,nT+1)),end_ecriture,c='b')
# ax1.grid()
# ax2.grid()
command2="./gprMaxMerge "+ nom + "_"
os.popen(command2).readlines()
end_running.append(time.time()-start_running)
#os.popen("rm -rf *.in")
print('Running GPR time:'+str(end_running))
for i in range(0,nT+1) :
os.popen("rm -rf "+ nom + "_" + str(i+1) + ".out")
os.chdir(dir)
return myDirName
#run(0.3,10,0.024,0.097,0.401,) #Ks, n, alpha, tr, ts, ti