-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathPySHEMAT.py
4744 lines (4318 loc) · 189 KB
/
PySHEMAT.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
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
PySHEMAT is a free set of Python modules to create and process input
files for fluid and heat flow simulation with SHEMAT (http://137.226.107.10/aw/cms/website/zielgruppen/gge/research_gge/~uuv/Shemat/?lang=de)
******************************************************************************************
PySHEMAT is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
PySHEMAT is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with PyTOUGH. If not, see <http://www.gnu.org/licenses/>.
******************************************************************************************
For module, update and documentation of PySHEMAT please see:
https://github.com/flohorovicic/PySHEMAT
If you use PySHEMAT for a scientific study, please cite our publication in Computers and Geoscience:
J. Florian Wellmann, Adrian Croucher, Klaus Regenauer-Lieb, Python scripting libraries for subsurface fluid and heat flow simulations with TOUGH2 and SHEMAT, Computers & Geosciences,
Available online 21 October 2011, ISSN 0098-3004, 10.1016/j.cageo.2011.10.011.
"""
import os, sys
import re # for regular expression fit, neccessary for stupid nlo files
# from matplotlib import use
# use("Agg")
import matplotlib as m
import numpy as np
class Shemat_file:
"""Class for SHEMAT simulation input and output files
Object definition for SHEMAT simulation files. The Object
can be used for simulation input files (.nml) and output files (.nlo).
Object methods enable a direct access of all variables and parameters
defined in the input file. Further methods are available for a
simplified generation of result plots in 2-D sections and to export simulation
results in a variety of different formats.
**Arguments**:
- *new_filename* = string: filename in case an empty file is created
**Optional keywords**:
- *offscreen* = True/False: set variables for offscreen rendering, e.g. to create plots on
a remote machine via ssh
"""
# class needed? maybe good as a "containter -> might be good
# to assure compatibility with future versions...
def __init__(self, filename='', **kwds):
"""Initialization of SHEMAT object
**Arguments**:
- *new_filename* = string: filename in case an empty file is created
- *filename* = string : filename of SHEMAT file to load
**Optional keywords**:
- *offscreen* = True/False: set variables for offscreen rendering, e.g. to create plots on
a remote machine via ssh
"""
if filename == '':
print "create empty file"
if kwds.has_key('new_filename'): self.filename = kwds['new_filename']
else:
self.filelines = self.read_file(filename)
self.filename = filename
self.idim = int(self.get("IDIM"))
self.jdim = int(self.get("JDIM"))
self.kdim = int(self.get("KDIM"))
if kwds.has_key('offscreen') and kwds['offscreen']:
m.use('Agg')
def read_file(self, filename):
"""Open and read a SHEMAT .nml or .nlo file
**Arguments**:
- *filename* = string: filename
**Returns**:
List of lines in the file
"""
try:
file = open(filename, "r")
except IOError, (nr, string_err):
print "Can not open file " + filename + ": " + string_err + " Err#" + str(nr)
print "Please check file name and directory and try again"
raise IOError
# check if number of entries correct
filelines = file.readlines()
file.close()
# set local variables
return filelines
def write_file(self, filename):
"""Write SHEMAT object to file
**Arguments**:
- *filename* = string: filename of SHEMAT file (extension .nml is
automatically assigned if not given)
"""
if filename[-4:] != ".nml":
# print "Add extesion .nml to filename " + filename
filename += ".nml"
try:
file = open(filename,"w")
except IOError, (nr, string_err):
print "Can not open file" + filename + ": " + string_err + " Err#" + str(nr)
print "Please check file name and directory and try again"
exit(0)
print "Write new SHEMAT file: " + filename
file.writelines(self.filelines)
file.close()
def adjust_ctl_file(self, old_filename, new_filename, **kwds):
"""Adjust an existing SHEMAT control file with a new filename
SHEMAT control files (.ctl) are required to run a SHEMAT simulation
from the command line. This function changes the filename settings
within a control file in the current directory from old_filename to
new_filename.
**Arguments**:
- *old_filename* = string: filename of original SHEMAT file in directory
- *new_filename* = string: filename of new SHEMAT file
**Optional keywords**:
- *backup* = False/ True: create shemat_ctl.bak backup file, default true
- *ctl_filename* = "shemat.ctl" : set other ctl filename
"""
# check, if extension at filename
if new_filename[-4:] != ".nml":
# print "Add extesion .nml to filename " + filename
new_filename += ".nml"
if kwds.has_key("ctl_filename"):
ctl_file = kwds["ctl_filename"]
else:
ctl_file = "shemat.ctl"
try:
file = open(ctl_file,"r")
except IOError, (nr, string_err):
print "Can not open file " + ctl_file + ": " + string_err + " Err# " + str(nr)
from os import getcwd
print "Please check if file exists in working directory " + getcwd()
exit(0)
ctl_lines = file.readlines()
file.close()
# write backup file
if kwds.has_key("backup") and not kwds["backup"]:
print "do not write backup file"
from os import remove
remove(ctl_file)
else:
print "Rename " + ctl_file + " to backup file shemat_ctl.bak"
from os import rename, remove
try:
rename(ctl_file, "shemat_ctl.bak")
except WindowsError:
try:
remove("shemat_ctl.bak")
rename(ctl_file, "shemat_ctl.bak")
except WindowsError, (nr, string_err):
print "Error renaming file: " + string_err + " Err# " + str(nr)
exit(0)
for i,l in enumerate(ctl_lines):
if "SHEMAT Control File" in l:
ctl_lines[i] = "# SHEMAT Control File (generated by Shemat_file.py)\n"
if l[0] != "#":
if old_filename in l:
ctl_lines[i] = new_filename + "\n"
if old_filename[0:-4] + " 3" in l:
ctl_lines[i] = new_filename[0:-4] + " 3\n"
# write new file
try:
file = open(ctl_file,"w")
except IOError, (nr, string_err):
print "Can not open file " + ctl_file + ": " + string_err + " Err# " + str(nr)
from os import getcwd
print "Please check if file exists in working directory " + getcwd()
exit(0)
file.writelines(ctl_lines)
file.close()
def coupling_info(self):
"""Determine coupling settings
Read coupling information out of variable KOPLNG and
dechipher.
**Returns**:
String with clear text about coupling settings"""
coupl_str = self.get("KOPLNG")
# dechipher
str_out = ""
if coupl_str[0] == "C":
str_out = str_out + "Coupled simulation of flow and heat\n"
elif coupl_str[1] == "O":
str_out = str_out + "density = function of concentration flow and chemical reaction:\n"
if coupl_str[2] == "U":
str_out = str_out + "\t- Kozeny-Carman-Pape\n"
else:
str_out = str_out + "\t- other coupling, see book pg. 37\n"
elif coupl_str[3] == "P":
str_out = str_out + "rock thermal cond = f(T)\n"
else:
str_out = str_out + "no coupling defined\n"
return str_out
def fixed_parameter_info(self):
"""Determine status of fixed parameters and return as string
**Returns**:
String with clear text about fixed parameters
"""
s = ""
f = self.get_array("KOPPX_FLAG")
v = self.get_array("KOPPX_FIELD")
n = ['Density:\t\t', 'Thermal Conductivity:\t',
'Compressibility:\t', 'Thermal capacity:\t', 'Viscosity:\t\t']
for i,f in enumerate(f):
s += n[i]
if f == 0: s += "not fixed\n"
else: s+= "%f (fixed) \n"
return s
def get(self, var_name, line=1):
"""Get the value of a scalar variable
Determines the value of a scalar variable in the SHEMAT object
**Arguments**:
- *var_name* = string : Name of scalar variable
**Optional keywords**:
- *line* = int : Number of lines for multiline variables
**Returns**:
String with variable
"""
for (i,l) in enumerate(self.filelines):
if var_name in l:
if line == 1:
return self.filelines[i+1]
break
else:
lines = []
for j in range(line):
lines.append(self.filelines[i+1])
return lines
break
def get_array(self, var_name):
"""Get the value of an array variable
Array variables (e.g. temperature, pressure, etc.) in SHEMAT are stored
in 1-D arrays in a compressed format. With this method, the variables
are decompressed and returned as a 1-D list. The method also adjusts
special boundary condition settings which are partly implemented as
negative values of pressure, porosity and permeability.
**Arguments**:
- *var_name* = string : Name of scalar variable
**Returns**:
- 1-D list with values
"""
for (i,l) in enumerate(self.filelines):
# Include: explicit definition of var_name to prevent double-matching
if var_name[0] == "#":
if not var_name in l: continue
elif not "# " + var_name in l: continue
if True:
# problem with .nlo files: strangely other value separators
# (csv) and unfortunately also some extra newlines...
# thus: read data until next # separator (can this cause problems?)
j = 1
data_raw = ''
if 'AREAKT' in var_name: data_raw = data_raw + self.filelines[i+j]
else:
while "#" not in self.filelines[i+j]:
data_raw = data_raw + self.filelines[i+j]
# print "used %d" % j
j += 1
# data_raw = self.filelines[i+1]
# if "#" not in self.filelines[i+2]: data_raw = data_raw + self.filelines[i+2]
data_raw = re.split(', *| *', data_raw)
data = []
# print "length of data for " + var_name + " %d" % len(data_raw)
if len(data_raw) == 1: # new_line at end of string!
if data_raw[0][-1] == '\n':
data_raw[0] = data_raw[0][0:-1] # last "entry" is newline command!
else:
data_raw[0] = data_raw[0][0:-2]
else:
if data_raw[-1] == '\n':
data_raw = data_raw[0:-1] # last "entry" is newline command!
for d in data_raw:
# if multiple definition with "*": split
if "*" in d:
d1 = d.split("*")
for i in range(int(d1[0])):
data.append(float(d1[1]))
else:
if d == '' or d == ' ' or d == '\n': continue # last value is an empty string
try:
data.append(float(d))
except ValueError:
print "Possibly a problem with a value, probably empty string: " + d
print "While reading %s " % var_name
print "If all values are correct, this error could be due to an incorrect"
print "line ending, for example when opening a file created on a windows machine"
print "with PySHEMAT on Linux or Mac"
# check, if Boundary Conditions are affected
if var_name == "POR" or var_name == "PERM" or var_name == "PRES":
# check, if bcs are already read:
try:
self.diri_temp
except AttributeError:
self.get_bcs()
# now: check data itself and return as positive value!
for i,d in enumerate(data):
if d < 0:
data[i] = -d
try:
data
except UnboundLocalError:
print "Variable " + var_name + " not defined in Shemat file!"
return None
return data
def get_bcs(self):
"""Determine the boundary conditions and return as 1-D list
Dirichlet boundary conditions are stored in a .nml file as:
- conc: negative PRES values
- temp: negative POR values
- head: neagitve PERM values
With this method, boundary conditions are evaluated and stored in
object variables self.diri_temp, self.diri_conc, and self.diri_head"""
# initialize object variables
print "get BCs and store in object variables"
self.diri_conc = []
self.diri_temp = []
self.diri_head = []
for var_name in ["POR", "PRES", "PERM"]:
if var_name == "POR":
print "read temperature boundary conditions"
elif var_name == "PRES":
print "read concentration boundary conditions"
elif var_name == "PERM":
print "read hydraulic head boundary conditions"
for (i,l) in enumerate(self.filelines):
if var_name in l:
# problem with .nlo files: strangely other value separators
# (csv) and unfortunately also some extra newlines...
# thus: read data until next # separator (can this cause problems?)
j = 1
data_raw = ''
while "#" not in self.filelines[i+j]:
data_raw = data_raw + self.filelines[i+j]
# print "used %d" % j
j += 1
# data_raw = self.filelines[i+1]
# if "#" not in self.filelines[i+2]: data_raw = data_raw + self.filelines[i+2]
data_raw = re.split(', *| *', data_raw)
for d in data_raw[0:-1]: # last "entry" is newline command!
# if multiple definition with "*": split
if "*" in d:
d = d.split("*")
for i in range(int(d[0])):
if d == '' : continue
# now, if value is negative: BC!
if float(d[1]) < 0: tmp = True
else: tmp = False
if var_name == "POR":
self.diri_temp.append(tmp)
elif var_name == "PRES":
self.diri_conc.append(tmp)
elif var_name == "PERM":
self.diri_head.append(tmp)
else:
if d == '' : continue
try:
if float(d) < 0: tmp = True
except ValueError:
continue
else: tmp = False
if var_name == "POR":
self.diri_temp.append(tmp)
elif var_name == "PRES":
self.diri_conc.append(tmp)
elif var_name == "PERM":
self.diri_head.append(tmp)
return True
def update_bcs(self):
"""Update boundary conditions after changing arrays self.diri_temp, etc.
Dirichlet boundary conditions are stored in a .nml file as:
- conc: negative PRES values
- temp: negative POR values
- head: neagitve PERM values
With this method, boundary conditions are evaluated and stored in
object variables self.diri_temp, self.diri_conc, and self.diri_head
..note : Only works, if self.diri_temp, self.diri_conc, self.diri_head are already object attributes!
"""
por = self.get_array("POR")
self.set_array("POR", por)
perm = self.get_array("PERM")
self.set_array("PERM", perm)
pres = self.get_array("PRES")
self.set_array("PRES", pres)
def set(self, var_name, value, line=1):
"""Set a SHEMAT variable to a specific value
**Arguments**:
- *var_name* = string : name of the SHEMAT variable
- *value* = string or number : variable value"""
for (i,l) in enumerate(self.filelines):
if var_name in l:
self.filelines[i+line] = str(value) + "\n"
def set_array(self, var_name, value_list, **kwds):
"""Set a SHEMAT array variable to values in a 1-D list
Set variable 'var_name' with values in 'value_list' in .nml_file;
mulitpliers "*" are constructed (as used in Shemat .nml file).
The SHEMAT boundary condition definitions are considered:
- self.diri_conc: Dirichlet BC for concentration => PRES neg
- self.diri_temp: Dirichlet BC for temperature => POR neg
- self.diri_head: Dirichlet BC for hydr head => PERM neg
these object variables/ local variables are automatically
set if PERM, POR or PRES are set
**Arguments**:
- *var_name* = string : Name of SHEMAT variable
- *value_list* = 1-D list of strings or numbers : Values
**Optional keywords**:
- float_type = 'high_res', 'normal' : precision of floating point; normal: 2 digits only
"""
value = ""
# check, if Boundary Conditions are affected
if var_name == "POR" or var_name == "PERM" or var_name == "PRES":
# check, if bcs are already read:
try:
self.diri_temp
except AttributeError:
self.get_bcs()
# now: set negative values, if BCs are defined at list positions
if var_name == "POR":
for i,l in enumerate(value_list):
if self.diri_temp[i]:
value_list[i] = -l
if var_name == "PERM":
for i,l in enumerate(value_list):
if self.diri_head[i] == True:
value_list[i] = -l
if var_name == "PRES":
for i,l in enumerate(value_list):
if self.diri_conc[i] == True:
value_list[i] = -l
# search variable position in .nml file
for (i,l) in enumerate(self.filelines):
if var_name in l:
# construct variable in correct format with multiplier "*"
n = 1
for (j,val) in enumerate(value_list):
try:
if val == value_list[j+1]:
n += 1
continue
except IndexError: pass
if n == 1:
if var_name == "PERM" or var_name == "HPR":
# check if value == 0, then: assign 0 only
if val == 0:
value += "0 "
else:
value += "%.2e " % val
elif var_name == "GEOLOGY": # or var_name == "ANISOJ" or var_name == "ANISOI":
# integer value
value += "%d " % val
elif var_name == "QBASAL3D" or var_name == "# QBASAL3D":
# assign four digits
value += "%.4e " % val
elif kwds.has_key('float_type'):
if kwds['float_type'] == 'high_res':
value += "%.2e " % val
else:
value += "%.2f " % val
else:
if var_name == "PERM" or var_name == "HPR":
# check if value == 0, then: assign 0 only
if val == 0:
value += "%d*0 " % n
else:
value += "%d*%.2e " % (n,val)
elif var_name == "GEOLOGY": # or var_name == "ANISOJ" or var_name == "ANISOI":
# integer value
value += "%d*%d " % (n,val)
elif var_name == "QBASAL3D":
# assign three digits
value += "%d*%.3e " % (n,val)
elif kwds.has_key('float_type'):
if kwds['float_type'] == 'high_res':
value += "%d*%e " % (n, val)
else:
value += "%d*%.2f " % (n, val)
n = 1
# set in .nml file
self.filelines[i+1] = value + "\n"
def set_array_from_xyz_structure(self,var_name, xyz_structure_list,**kwds):
"""Set a SHEMAT variable from a 3-D list of values
Set variable 'var_name' with values in 'xyz_structure_list' in .nml_file;
xyz_structure_list can be one of those derived from get_array_as_xyz_structure
mulitpliers "*" are constructed (as used in Shemat .nml file)
The SHEMAT boundary condition definitions are considered:
- self.diri_conc: Dirichlet BC for concentration => PRES neg
- self.diri_temp: Dirichlet BC for temperature => POR neg
- self.diri_head: Dirichlet BC for hydr head => PERM neg
these object variables/ local variables are automatically
set if PERM, POR or PRES are set.
**Arguments**:
- *var_name* = string : Name of SHEMAT variable
- *value_list* = 1-D list of strings or numbers : Values
**Optional keywords**:
- float_type = 'high_res', 'normal' : precision of floating point; normal: 2 digits only
"""
value = ""
# decompose xyz-list into norma value list
idim = int(self.get("IDIM"))
jdim = int(self.get("JDIM"))
kdim = int(self.get("KDIM"))
value_list = []
for k in range(kdim):
for j in range(jdim):
for i in range(idim):
value_list.append(xyz_structure_list[i][j][k])
# check, if Boundary Conditions are affected
if var_name == "POR" or var_name == "PERM" or var_name == "PRES":
# check, if bcs are already read:
try:
self.diri_temp
except AttributeError:
self.get_bcs()
# now: set negative values, if BCs are defined at list positions
if var_name == "POR":
for i,l in enumerate(value_list):
if self.diri_temp[i]:
value_list[i] = -l
if var_name == "PERM":
for i,l in enumerate(value_list):
if self.diri_head[i] == True:
value_list[i] = -l
if var_name == "PRES":
for i,l in enumerate(value_list):
if self.diri_conc[i] == True:
value_list[i] = -l
# search variable position in .nml file
for (i,l) in enumerate(self.filelines):
if var_name in l:
# construct variable in correct format with multiplier "*"
n = 1
for (j,val) in enumerate(value_list):
try:
if val == value_list[j+1]:
n += 1
continue
except IndexError: pass
if n == 1:
if var_name == "PERM" or var_name == "HPR":
value += "%.2e " % val
elif kwds.has_key('float_type'):
if kwds['float_type'] == 'high_res':
value += "%.2e " % val
else:
value += "%.2f " % val
else:
if var_name == "PERM" or var_name == "HPR":
value += "%d*%.2e " % (n,val)
elif kwds.has_key('float_type'):
if kwds['float_type'] == 'high_res':
value += "%d*%e " % (n, val)
else:
value += "%d*%.2f " % (n, val)
n = 1
# set in .nml file
self.filelines[i+1] = value + "\n"
def create_formations_ids(self):
"""create array self.formation_ids with formation ids from geology data array"""
# property zones in SHEMAT are defined from 1 to n => sufficient to find the maximum value of geology array
geology = self.get_array("GEOLOGY")
self.formation_ids = range(1,int(max(geology))+1)
return self.formation_ids
def array_to_xyz_structure_object(self, array):
"""restructure array into x,y,z 3-D structure as data[i][j][k]
i,j,k are counters in x,y,z-direction, determined from the object itself
**Arguments**:
- *array* = list or array to be restructured
**Returns**:
Restructured 3-D list of data values
"""
from numpy import size
if not (self.idim*self.jdim*self.kdim == size(array)):
print "new dimensions not consistent with array size"
print "can not create new array!"
return
# initialize output array
data = []
# iterate over dimensions and re-sort array
n = 0
# initialize array in x,y,z-order
for i in range(self.idim):
tmp2 = []
for j in range(self.jdim):
tmp = [0 for k in range(self.kdim)]
tmp2.append(tmp)
data.append(tmp2)
# now, fill data structure in z,y,x loops for correct x,y,z order
for k in range(self.kdim):
for j in range(self.jdim):
for i in range(self.idim):
data[i][j][k] = array[n]
n += 1
return data
def array_to_xyz_structure(self,array,idim,jdim,kdim):
"""restructure array into x,y,z 3-D structure as data[i][j][k]
with i,j,k: counters in x,y,z direction
**Arguments**:
- *array*: array to be restructured
- *idim* = int : dimension of new array in x-direction (default: model dimensions)
- *jdim* = int : dimension of new array in y-dircetion
- *kdim* = int : dimension of new array in z-direction
**Returns**:
Restructured 3-D list of data values
"""
# check if array length and new dimensions are consistent
from numpy import size
if not (idim*jdim*kdim == size(array)):
print "new dimensions not consistent with array size"
print "can not create new array!"
return
# initialize output array
data = []
# iterate over dimensions and re-sort array
n = 0
# initialize array in x,y,z-order
for i in range(idim):
tmp2 = []
for j in range(jdim):
tmp = [0 for k in range(kdim)]
tmp2.append(tmp)
data.append(tmp2)
# now, fill data structure in z,y,x loops for correct x,y,z order
for k in range(kdim):
for j in range(jdim):
for i in range(idim):
data[i][j][k] = array[n]
n += 1
return data
def get_cell_centres(self):
"""Calculate centre of cells in absolute values
Cell centres are stored in object variables
self.centre_x, self.centre_y, self.centre_z
"""
# reload cell boundaries
self.get_cell_boundaries()
# calculate centres for x
self.centre_x = []
for i in range(len(self.boundaries_x[:-1])):
self.centre_x.append((self.boundaries_x[i+1]-self.boundaries_x[i])/2.+self.boundaries_x[i])
# calculate centres for y
self.centre_y = []
for i in range(len(self.boundaries_y[:-1])):
self.centre_y.append((self.boundaries_y[i+1]-self.boundaries_y[i])/2.+self.boundaries_y[i])
# calculate centres for z
self.centre_z = []
for i in range(len(self.boundaries_z[:-1])):
self.centre_z.append((self.boundaries_z[i+1]-self.boundaries_z[i])/2.+self.boundaries_z[i])
def update_model_from_geomodeller_xml_file(self, geomodeller_xml_file, **kwds):
"""Determine the model geology from a geological model created with GeoModeller
With this method, it is possible to map the distribution of geological units
from a 3-D geological model, created with GeoModeller, on the mesh used
in the SHEMAT simulation. This is here done with a direct dll access to read
the geology data into shemat object.
.. warning: directory path to GeoModeller bin directory has to be adjusted in code!
**Arguments**:
- *geomodeller_xml_file* = string : XML file of GeoModeller project
**Optional keywords**:
- *compute* = True/False: (re-) compute geological model before processing
- *lower_left_x* = float : x-position of lower-left corner (default: model range)
- *lower_left_y* = float : y-position of lower-left corner (default: model range)
- *lower_left_z* = float : z-position of lower-left corner (default: model range)
.. note: for more information on GeoModeller see the `GeoModeller page <http:www.geomodeller.com\>`_
"""
import geomodeller_api as g_api
g_api.initialise_geomodeller_api(r'C:\Geomodeller\GeoModeller_1502\bin')
g_api.load_model(geomodeller_xml_file)
if kwds.has_key('compute') and kwds['compute']==True:
g_api.compute()
(geo_dx,geo_dy,geo_dz) = g_api.get_model_extent()
(xmin,ymin,zmin,xmax,ymax,zmax) = g_api.get_model_bounds()
if kwds.has_key('lower_left_x') and kwds['lower_left_x'] != '':
xmin = kwds['lower_left_x']
if kwds.has_key('lower_left_y') and kwds['lower_left_y'] != '':
ymin = kwds['lower_left_y']
if kwds.has_key('lower_left_z') and kwds['lower_left_z'] != '':
zmin = kwds['lower_left_z']
# check model extents of GeoModeller model and SHEMAT file
self.get_cell_centres()
geo_new = [] # new geology array
n_cells = len(self.centre_z) * len(self.centre_y) * len(self.centre_x)
i = 0
for z in self.centre_z:
for y in self.centre_y:
for x in self.centre_x:
# geo_new.append(g_api.get_lithology(x,y,-geo_dz+z)) # old version
print "process cell %8d of %d, %4.1f Percent" % (i, n_cells, (float(i)/float(n_cells)*100))
geo_new.append(g_api.get_lithology(xmin + x, ymin + y, zmin + z))
i += 1
self.set_array("GEOLOGY", geo_new)
def update_properties_from_csv_list(self, csv_file_name):
"""update SHEMAT properties from csv file, based on property
distribution in GEOLOGY array;
This function can be used
- to populate SHEMAT input file with new parameters
- to update SHEMAT input file when GEOLOGY array was changed
(e.g. after update of GEOLOGY array from GeoModeller input file)
conventions for .csv file:
Header line neccessary, should include GEOLOGY column and one column
for all the properties to be updated
Header can contain a first column with Title 'NAME' for name of formation (not used,
just for better overview in .csv file)
Optional property lines in .csv-file:
PERM_FUNC: indicating that definition of permeability function in following lines
PERM_RANDOM: indicating that definition of permeability random distribution follows
WLFM0_RANDOM: indicating that thermal conductivity random distribution follows
RHOCM_RANDOM: indicating that thermal diffusivity random distribution follows
"""
geol = self.get_array("GEOLOGY")
# open csv file
csv = open(csv_file_name,'r')
csv_lines = csv.readlines()
csv.close()
csv_header = csv_lines[0].split(',')
# remove end of line character from last entry
csv_header[-1] = csv_header[-1].rstrip()
# test if file contains another column (e.g. a 'Name' column) before GEOLOGY column
if csv_header[0] == 'GEOLOGY':
start = 1
else:
start = 2 # in this case, GEOLOGY column HAS TO BE in the second column!!
csv_datalines = []
for i_csv,l in enumerate(csv_lines[1:]):
# print "line %d: %s" % (i_csv,l)
if "PERM_FUNC" in l:
print "permeability function defined"
# test correct definition
if csv_lines[i_csv+2].split(',')[0] != "Type":
print "please define type (kwd Type) in input line %d" % i_csv+2
break
if csv_lines[i_csv+3].split(',')[0] != "Lines":
print "please define number of lines (kwd Lines) in input line %d" % i_csv+3
break
# deconstruct input file and create permeability_function_dict with formation id as key
type = csv_lines[i_csv+2].split(',')[1]
lines = csv_lines[i_csv+3].split(',')[1]
perm_func = {}
perm_func['type'] = type
if type == 'linear_global':
for n_line in range(int(lines)):
line = csv_lines[i_csv+5+n_line].split(',')
func_dict = {'k_min' : float(line[start]),
'k_min_depth' : float(line[start+1]),
'k_max' : float(line[start+2]),
'k_max_depth' : float(line[start+3])}
perm_func[line[start-1]] = func_dict
if "PERM_RANDOM" in l:
print "permeability random distribution defined"
# test correct definition
if csv_lines[i_csv+2].split(',')[0] != "Type":
print "please define type (kwd Type) in input line %d" % i_csv+2
break
if csv_lines[i_csv+3].split(',')[0] != "Lines":
print "please define number of lines (kwd Lines) in input line %d" % i_csv+3
break
# deconstruct input file and create permeability_random_dict with formation id as key
type = csv_lines[i_csv+2].split(',')[1]
lines = csv_lines[i_csv+3].split(',')[1]
perm_rand = {}
perm_rand['type'] = type
if type == 'lognormal':
for n_line in range(int(lines)):
line = csv_lines[i_csv+5+n_line].split(',')
rand_dict = {'log_stdev' : line[start]}
perm_rand[line[start-1]] = rand_dict
else:
csv_datalines.append(l.split(','))
if csv_header[0] != "GEOLOGY" and csv_header[1] != "GEOLOGY":
print "first or second column should contain GEOLOGY data!"
return None
prop_array = []
for prop in csv_header[start:]:
p = self.get_array(prop)
prop_array.append(p)
# initialise new array
prop_array_new = []
for p in prop_array:
prop_array_new.append([])
# determine row number of geology units (for correct mapping below)
geol_pos = {}
for i,line in enumerate(csv_datalines):
geol_pos[int(line[1])] = i
# print geol_pos
# create new property array
for i,g in enumerate(geol):
for j,prop in enumerate(prop_array):
try:
val = csv_datalines[geol_pos[int(g)]][j+start]
prop_array_new[j].append(float(val))
except KeyError:
print "Geology: %d" % int(g-1)
print j
print start
# print "Value: %f" % csv_datalines
print 80*'*'
print "csv file not complete: no values for geology unit %d found!" % g
print "please check file and try again"
print 80*'*'
raise KeyError("Input file not complete, index for one geology unit not found!")
# update propperty arrays in Shemat file
for i,prop in enumerate(csv_header[start:]):
self.set_array(prop,prop_array_new[i])
# now: test, if permeability function defined and if so: (re-) compute permeability for thir formation
try:
perm_func
perm_func_defined = True
except NameError:
pass
# now: test, if permeability randm distribution defined and if so: (re-) compute permeability for thir formation
try:
perm_rand
perm_rand_defined = True
except NameError:
pass
try:
perm_func_defined
print "Permeability function defined, recompute permeabilities for formations"
# define peermeability function
def permeability(z, total_depth, k_min, k_min_depth, k_max, k_max_depth, **kwds):
"""define permeability as a function of depth"""
from numpy import log10
# k = 10**(log10(k_max) - z * (log10(k_max) - log10(k_min)) / (k_max_depth - k_min_depth))
k = 10**(log10(k_max) - (k_max_depth - z) * (log10(k_max) - log10(k_min)) / (k_max_depth - k_min_depth))
return k
# load relevant parameters
geol_xyz = self.get_array_as_xyz_structure("GEOLOGY")
idim = int(self.get("IDIM"))
jdim = int(self.get("JDIM"))
kdim = int(self.get("KDIM"))
dz = self.get_array("DELZ")
# calculate total depth, assuming constant layer thickness!!
total_depth = sum(dz)
perm_xyz = self.get_array_as_xyz_structure("PERM")
# formation iD
for j in range(jdim):
for i in range(idim):
for k in range(kdim):
g = str(int(geol_xyz[i][j][k]))
if perm_func.has_key(g):
# perm_xyz[i][j][k] = permeability((total_depth - sum(dz[0:k])), total_depth, 1E-15, 1E-13)
perm_xyz[i][j][k] = permeability((total_depth - sum(dz[0:k])),
total_depth,
perm_func[g]['k_min'],
perm_func[g]['k_min_depth'],
perm_func[g]['k_max'],
perm_func[g]['k_max_depth']
)
else:
perm_xyz[i][j][k]
self.set_array_from_xyz_structure("PERM", perm_xyz)
except NameError:
pass
return None
def update_property_from_dict(self,property,property_dict):
"""update the value of one property according to geology id;
properties are passed as a dictionary with assignments, e.g:
{1 : 2.9, 2 : 3, ...}
property = string : property name, according to SHEMAT variable name
property_dict = dict : assigned values; can be defined for a subset of
geological units only, all not defined are not changed
"""
geol = self.get_array("GEOLOGY")
prop_array = self.get_array(property)
prop_array_new = []
for i,g in enumerate(geol):
if property_dict.has_key(g):
prop_array_new.append(property_dict[g])
else:
prop_array_new.append(prop_array[i])
self.set_array(property, prop_array_new)
def get_model_extent(self):
"""determine model extent in x,y,z direction from self.boundaries_x,y,z
return as tuple (extent_x, extent_y, extent_z) and store in self.extent_x, self.extent_y, self.extent_z"""
# recompute boundaries (always or trigger with kwd?)
self.get_cell_boundaries()
self.extent_x = max(self.boundaries_x) - min(self.boundaries_x)
self.extent_y = max(self.boundaries_y) - min(self.boundaries_y)
self.extent_z = max(self.boundaries_z) - min(self.boundaries_z)
return (self.extent_x, self.extent_y, self.extent_z)
def get_cell_boundaries(self):
"""calculate cell boundaries from arrays DELX, DELY, DELZ and save
to self.boundaries_x, self.boundaries_y, self.boundaries_z"""
delx = self.get_array("DELX")
dely = self.get_array("DELY")
delz = self.get_array("DELZ")
# x
b = []
laststep = 0
b.append(laststep)
for step in delx:
b.append(laststep+step)
laststep += step
self.boundaries_x = b
# y
b = []
laststep = 0
b.append(laststep)
for step in dely:
b.append(laststep+step)
laststep += step
self.boundaries_y = b
# z
b = []
laststep = 0
b.append(laststep)
for step in delz:
b.append(laststep+step)
laststep += step
self.boundaries_z = b
def get_array_as_xyz_structure(self,var_name):
"""read variable and order into 3-D structure as
data[i][j][k]
with i,j,k: counters in x,y,z direction
arguments:
var_name : property from SHEMAT file (read again in this function! rewrite?)
"""
# initialize output array
data = []
# read array from SHEMAT file
ori_array = self.get_array(var_name)
# read array length from SHEMAT file
idim = int(self.get("IDIM"))
jdim = int(self.get("JDIM"))