-
Notifications
You must be signed in to change notification settings - Fork 1
/
Main.py
executable file
·145 lines (132 loc) · 3.79 KB
/
Main.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 20 14:22:11 2020
@author: el
"""
import matplotlib.pyplot as plt
from param_acquisition import Geometry,ParamMVG, ParamGPRMAX
import os
os.chdir('/home/el/Codes/Porchet-GPR')
#%% Param MVG
# Teneur en eau résiduelle
tr = 0.008
# Teneur en eau à saturation
ts = 0.35
# Teneur en eau initiale
ti = 0.07
# Perméabilité à saturation
Ks = 0.215
# param fitting retention n
n = 5
# param fitting retention alpha
alpha = 0.03
pVg=ParamMVG(tr=tr, ts=ts, ti=ti, Ks=Ks, n=n, alpha=alpha)
pVg.porosity = pVg.ts
#%% Geometrie
#def des paramètres géométriques
geometry=Geometry()
#Domaine de calcul (en cm)
#Tolerance
geometry.tol=10**(-7)
# largeur
geometry.xmin=0
geometry.xmax=40
# hauteur (elevation)
geometry.emin=0
geometry.emax = 80
# profondeur du trou en cm
geometry.dtrou = 30
# elevation du fond du trou
geometry.etrou = geometry.emax - geometry.dtrou
# rayon du trou en cm
geometry.r=2.5
# hauteur d'eau imposée au fond du trou en cm
geometry.h_eau=10#5.0
# pas de la maille en cm
#geometry.dx = 0.1
geometry.dx = 1
# profondeur sous le trou (cm) jusqu'où on souhaite un maillage affiné.
geometry.zaff= 20
#largeur horizontal de la zone affinée (cm)
geometry.waff=20
# elevation de l'affinage
geometry.eaff=geometry.etrou-geometry.zaff
# contrainte d'angle min pour mesh
geometry.quality=33
# maximum triangle size (m*²)
geometry.area=5
# tupple for mesh generation
geometry.smooth=[1,5]
#%% def de paramètres SWMS
#Temps d'infiltration où à lieu le calcul de chaque trace (minutes)
#temps=[1.00, 2.00]
temps=[0.17, 0.33, 0.50, 0.67, 0.83, 1.00, 2.00, 3.00, 4.00, 5.00, 6.00]
# Temps max de calcul SWMS2D au delà duquel on arrète le calcul (secondes)
tmax_SWMS2D = 600
#tmax_SWMS2D = 10
nT=len(temps)
#%% Definition des param gprMax
paramGPRMAX=ParamGPRMAX()
# Domaine de calcul (cm)
paramGPRMAX.xmin = geometry.xmin
paramGPRMAX.xmax = geometry.xmax
paramGPRMAX.zmin = geometry.emin
paramGPRMAX.zmax = geometry.emax
# Taille des mailles (cm)
################################################
paramGPRMAX.dx =0.5## attention ici
# Electrical conductivity of the medium
paramGPRMAX.sigma=0.0000
# Relative dielectric permittivity of water
paramGPRMAX.eps_w=80.1
# Relative dielectric permittivity of PVC
paramGPRMAX.eps_pvc=2
# Relative dielectric permittivity of pure silice
paramGPRMAX.eps_s=2.5
# Ricker signal central frequency (Hz)
paramGPRMAX.wave_freq = 1000e6
# Frequence max du signal EM (Hz)
paramGPRMAX.freq_max = 2.8 * paramGPRMAX.wave_freq
# Distance between hole middle and source (m)
paramGPRMAX.d_emet = 0.18
# Distance between hole middle and receiving antenna (m)
paramGPRMAX.d_recept = 0.22
# param qui raffine le pas spatial (par défaut 10 d'après doc gprmax)
paramGPRMAX.spatial_step = 5
# Trace time window (ns)
paramGPRMAX.time = 30e-9
#time_step_stability_factor (pas utilisé pour le moment...)
paramGPRMAX.fac_dt = 0.01
#%% Forward
from Forward import Forward
[TWT,Vol]=Forward(geometry,pVg,paramGPRMAX,temps,tmax_SWMS2D)
#import multiprocessing as mp
#print("Number of processors: ", mp.cpu_count())
#%%
#temps=[0.17, 0.33, 0.50, 0.67, 0.83, 1.00, 2.00, 3.00, 4.00, 5.00, 6.00]
XX=temps.copy()
XX.insert(0, 0)
plt.close('all')
fig,ax1=plt.subplots(1,1,figsize=(25,15))
ax1.plot(XX,TWT,'r',label='TWT(ns)')
ax1.grid()
ax1.set_xlabel('Minutes')
ax1.set_ylabel('TWT(ns)')
ax2 = ax1.twinx()
ax2.plot(XX,Vol,'b',label='Volume')
ax2.grid()
ax2.set_ylabel('Volume(mL)')
fig.legend()
#%%
from outils import rada_plot
import h5py
filepath='./OUTTESRdtrou30_rtrou3.0_tr10/ts0.35_ti0.07_tr0.008_n5_alpha0.03_Ks0.215/'
filename='radargram__merged.out'
f = h5py.File(filepath + filename, 'r')
path = '/rxs/rx1/'
data = f['%s%s' % (path, 'Ez')][:,:]
plt.figure(figsize=(10,15))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.imshow(data[:,:],aspect=0.005)