Skip to content

Commit cb666e4

Browse files
Victor-M-GomesVictor-M-Gomes
committed
Fix to wave prop modeling in p4 and p3. Also fix to allow use of python < 3.9
1 parent 5972c30 commit cb666e4

File tree

9 files changed

+131
-89
lines changed

9 files changed

+131
-89
lines changed

‎geos-utils/src/geos/utils/vtk/helpers.py

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
import logging
22
from copy import deepcopy
33
from numpy import argsort, array
4-
from typing import Iterator, Optional
4+
from typing import Iterator, Optional, List
55
from vtkmodules.util.numpy_support import vtk_to_numpy
66
from vtkmodules.vtkCommonCore import vtkDataArray, vtkIdList
77
from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkFieldData
@@ -30,7 +30,7 @@ def vtk_iter( l ) -> Iterator[ any ]:
3030
yield l.GetCellType( i )
3131

3232

33-
def has_invalid_field( mesh: vtkUnstructuredGrid, invalid_fields: list[ str ] ) -> bool:
33+
def has_invalid_field( mesh: vtkUnstructuredGrid, invalid_fields: List[ str ] ) -> bool:
3434
"""Checks if a mesh contains at least a data arrays within its cell, field or point data
3535
having a certain name. If so, returns True, else False.
3636
@@ -73,7 +73,7 @@ def getFieldType( data: vtkFieldData ) -> str:
7373
return "vtkFieldData"
7474

7575

76-
def getArrayNames( data: vtkFieldData ) -> list[ str ]:
76+
def getArrayNames( data: vtkFieldData ) -> List[ str ]:
7777
if not data.IsA( "vtkFieldData" ):
7878
raise ValueError( f"data '{data}' entered is not a vtkFieldData object." )
7979
return [ data.GetArrayName( i ) for i in range( data.GetNumberOfArrays() ) ]
@@ -91,7 +91,7 @@ def getCopyArrayByName( data: vtkFieldData, name: str ) -> Optional[ vtkDataArra
9191

9292

9393
def getGlobalIdsArray( data: vtkFieldData ) -> Optional[ vtkDataArray ]:
94-
array_names: list[ str ] = getArrayNames( data )
94+
array_names: List[ str ] = getArrayNames( data )
9595
for name in array_names:
9696
if name.startswith("Global") and name.endswith("Ids"):
9797
return getCopyArrayByName( data, name )
@@ -114,11 +114,11 @@ def getNumpyArrayByName( data: vtkFieldData, name: str, sorted: bool=False ) ->
114114
arr: array = vtk_to_numpy( getArrayByName( data, name ) )
115115
if arr is not None:
116116
if sorted:
117-
array_names: list[ str ] = getArrayNames( data )
117+
array_names: List[ str ] = getArrayNames( data )
118118
sortArrayByGlobalIds( data, arr, array_names )
119119
return arr
120120
return None
121121

122122

123123
def getCopyNumpyArrayByName( data: vtkFieldData, name: str, sorted: bool=False ) -> Optional[ array ]:
124-
return deepcopy( getNumpyArrayByName( data, name, sorted=sorted ) )
124+
return deepcopy( getNumpyArrayByName( data, name, sorted=sorted ) )

‎pygeos-tools/pyproject.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ dependencies = [
2222
"numpy",
2323
"scipy",
2424
"mpi4py",
25-
"vtk",
25+
# "vtk", # no need, should use the one with geos
2626
"pyevtk",
2727
"xmltodict",
2828
"h5py",

‎pygeos-tools/src/geos/pygeos_tools/utilities/model/SepModel.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
import sys
33
import numpy as np
44
import argparse
5+
import mpi4py
6+
mpi4py.rc.initialize = False
57
from mpi4py import MPI
68

79

‎pygeos-tools/src/geos/pygeos_tools/utilities/output/SEGYTraceOutput.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
import os
22
import numpy as np
33
import segyio
4+
import mpi4py
5+
mpi4py.rc.initialize = False
46
from mpi4py import MPI
57

68

‎pygeos-tools/src/geos/pygeos_tools/utilities/output/SEPTraceOutput.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
import os
22
import numpy as np
33
from geos.pygeos_tools.utilities.model.SepModel import SEPModel
4+
import mpi4py
5+
mpi4py.rc.initialize = False
46
from mpi4py import MPI
57

68

‎pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,12 +13,16 @@
1313
# ------------------------------------------------------------------------------------------------------------
1414

1515
import numpy as np
16+
17+
import mpi4py
18+
mpi4py.rc.initialize = False
19+
from mpi4py import MPI
20+
1621
import os
1722
import pygeosx
1823
import sys
1924
from geos.pygeos_tools.utilities.input.Xml import XML
2025
from geos.pygeos_tools.utilities.input.GeosxArgs import GeosxArgs
21-
from mpi4py import MPI
2226

2327

2428
class Solver:

‎pygeos-tools/src/geos/pygeos_tools/wrapper.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
import sys
22
import numpy as np
3+
import mpi4py
4+
mpi4py.rc.initialize = False
35
from mpi4py import MPI
46
import matplotlib.pyplot as plt
57
import pylvarray

‎pygeos-tools/tests/acoustic_modeling.py

Lines changed: 101 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -1,145 +1,166 @@
11
import argparse
22
import os
3+
4+
#GEOS
35
from geos.pygeos_tools.utilities.input import XML
46
from geos.pygeos_tools.utilities.acquisition_library.EquispacedAcquisition import EQUISPACEDAcquisition
57
from geos.pygeos_tools.utilities.solvers import AcousticSolver
68
from geos.pygeos_tools.utilities.output import SeismicTraceOutput
7-
from mpi4py import MPI
89

10+
import mpi4py
11+
mpi4py.rc.initialize = True
12+
from mpi4py import MPI
913

1014
def parse_args():
1115
"""Get arguments
1216
1317
Returns:
1418
argument '--xml': Input xml file for GEOSX
1519
"""
16-
parser = argparse.ArgumentParser( description="Modeling acquisition example - Acoustic" )
17-
parser.add_argument( '--xml', type=str, required=True, help="Input xml file for GEOS" )
18-
parser.add_argument( '--soutdir', required=False, type=str, default="./", help="Path to seismogram output dir" )
19-
parser.add_argument( '--soutn', required=False, type=str, default="seismo", help="Name of output seismograms" )
20-
parser.add_argument( '--param_file',
21-
type=str,
22-
required=False,
23-
default="identity",
24-
dest="pfile",
25-
help="Optional file containing modelling parameters" )
26-
27-
args, _ = parser.parse_known_args()
20+
parser = argparse.ArgumentParser(description="Modeling acquisition example - Acoustic")
21+
parser.add_argument('--xml', type=str, required=True,
22+
help="Input xml file for GEOS")
23+
parser.add_argument('--soutdir', required=False, type=str, default="./",
24+
help="Path to seismogram output dir")
25+
parser.add_argument('--soutn', required=False, type=str, default="seismo",
26+
help="Name of output seismograms")
27+
parser.add_argument('--param_file', type=str, required=False, default="identity", dest="pfile", help="Optional file containing modelling parameters")
28+
29+
args ,_ = parser.parse_known_args()
2830
return args
2931

30-
31-
def parse_workflow_parameters( pfile ):
32-
with open( pfile, "r" ) as f:
32+
def parse_workflow_parameters(pfile):
33+
with open(pfile, "r") as f:
3334
hdrStr = f.read()
3435

3536
hdrList = []
36-
for fl in hdrStr.split( '\n' ):
37-
l = fl.split( "#" )[ 0 ]
37+
for fl in hdrStr.split('\n'):
38+
l = fl.split("#")[0]
3839
if l:
3940
# add "--" to facilitate parsing that follows
4041
l = "--" + l
41-
hdrList += l.split( "=" )
42-
43-
parser = argparse.ArgumentParser( "Modelling workflow parser" )
44-
parser.add_argument( "--mintime", dest="mintime", default=None, type=float, help="Min time for the simulation" )
45-
parser.add_argument( "--maxtime", dest="maxtime", default=None, type=float, help="Max time for the simulation" )
46-
parser.add_argument( "--dt", dest="dt", default=None, type=float, help="Time step of simulation" )
47-
parser.add_argument( "--dtSeismo", dest="dtSeismo", default=None, type=float, help="Time step for " )
48-
parser.add_argument( "--sourceType", dest="sourceType", type=str, help="Source type" )
49-
parser.add_argument( "--sourceFreq", dest="sourceFreq", type=float, help="Ricker source central frequency" )
50-
51-
args, _ = parser.parse_known_args( hdrList )
42+
hdrList += l.split("=")
43+
44+
parser = argparse.ArgumentParser("Modelling workflow parser")
45+
parser.add_argument("--mintime", dest="mintime",
46+
default=None, type=float,
47+
help="Min time for the simulation")
48+
parser.add_argument("--maxtime", dest="maxtime",
49+
default=None, type=float,
50+
help="Max time for the simulation")
51+
parser.add_argument("--dt", dest="dt",
52+
default=None, type=float,
53+
help="Time step of simulation")
54+
parser.add_argument("--dtSeismo", dest="dtSeismo",
55+
default=None, type=float,
56+
help="Time step for ")
57+
parser.add_argument("--sourceType", dest="sourceType",
58+
type=str,
59+
help="Source type")
60+
parser.add_argument("--sourceFreq", dest="sourceFreq",
61+
type=float,
62+
help="Ricker source central frequency")
63+
64+
args, _ = parser.parse_known_args(hdrList)
5265
return args
5366

5467

5568
def main():
69+
70+
if MPI.Is_initialized():
71+
print("MPI initialized")
72+
else:
73+
print("MPI not initialized. Initializing...")
74+
MPI.Init()
75+
5676
comm = MPI.COMM_WORLD
5777
rank = comm.Get_rank()
58-
78+
5979
args = parse_args()
6080
xmlfile = args.xml
6181

62-
xml = XML( xmlfile )
82+
xml = XML(xmlfile)
6383

64-
wf_args = parse_workflow_parameters( args.pfile )
84+
wf_args = parse_workflow_parameters(args.pfile)
6585

66-
# Time Parameters
86+
#Time Parameters
6787
minTime = wf_args.mintime
6888
maxTime = wf_args.maxtime
6989
dt = wf_args.dt
7090
dtSeismo = wf_args.dtSeismo
7191

72-
# Output parameters
92+
#Output parameters
7393
outdirseis = args.soutdir
74-
os.makedirs( outdirseis, exist_ok=True )
94+
os.makedirs(outdirseis, exist_ok=True)
7595
outseisname = args.soutn
76-
77-
# Source parameters
96+
97+
#Source parameters
7898
sourceType = wf_args.sourceType
7999
sourceFreq = wf_args.sourceFreq
80-
100+
101+
81102
# Read acquisition
82103
if rank == 0:
83-
# acquisition = Acquisition(xml)
84-
acquisition = EQUISPACEDAcquisition( xml=xml,
85-
startFirstSourceLine=[ 305.01, 305.01, 5.01 ],
86-
endFirstSourceLine=[ 325.01, 305.01, 5.01 ],
87-
startLastSourceLine=[ 305.01, 325.01, 5.01 ],
88-
endLastSourceLine=[ 325.01, 325.01, 5.01 ],
89-
numberOfSourceLines=2,
90-
sourcesPerLine=2,
91-
startFirstReceiversLine=[ 121.02, 255.02, 58.01 ],
92-
endFirstReceiversLine=[ 471.02, 255.02, 58.01 ],
93-
startLastReceiversLine=[ 121.02, 255.02, 58.01 ],
94-
endLastReceiversLine=[ 471.02, 255.02, 58.01 ],
95-
numberOfReceiverLines=1,
96-
receiversPerLine=8 )
104+
#acquisition = Acquisition(xml)
105+
acquisition = EQUISPACEDAcquisition(xml=xml,
106+
startFirstSourceLine=[305.01,305.01,5.01],
107+
endFirstSourceLine=[325.01,305.01,5.01],
108+
startLastSourceLine=[305.01,325.01,5.01],
109+
endLastSourceLine=[325.01,325.01,5.01],
110+
numberOfSourceLines=2,
111+
sourcesPerLine=2,
112+
startFirstReceiversLine=[121.02,255.02,58.01],
113+
endFirstReceiversLine=[471.02,255.02,58.01],
114+
startLastReceiversLine=[121.02,255.02,58.01],
115+
endLastReceiversLine=[471.02,255.02,58.01],
116+
numberOfReceiverLines=1,
117+
receiversPerLine=8)
97118

98119
else:
99120
acquisition = None
100-
acquisition = comm.bcast( acquisition, root=0 )
121+
acquisition = comm.bcast(acquisition, root=0)
101122

102-
solver = AcousticSolver( dt, minTime, maxTime, dtSeismo, sourceType, sourceFreq )
123+
solver = AcousticSolver(dt,
124+
minTime,
125+
maxTime,
126+
dtSeismo,
127+
sourceType,
128+
sourceFreq)
103129

104-
for ishot, shot in enumerate( acquisition.shots ):
130+
for ishot, shot in enumerate(acquisition.shots):
105131
xmlshot = shot.xml
106132
rank = comm.Get_rank()
107-
108-
solver.initialize( rank, xmlshot )
133+
134+
solver.initialize(rank, xmlshot)
109135
solver.applyInitialConditions()
110136

111-
solver.updateSourceAndReceivers( shot.getSourceCoords(), shot.getReceiverCoords() )
112-
solver.updateVtkOutputsName( directory=f"Shot{shot.id}" )
113-
137+
solver.updateSourceAndReceivers(shot.getSourceCoords(), shot.getReceiverCoords())
138+
solver.updateVtkOutputsName(directory=f"Shot{shot.id}")
139+
114140
t = 0
115-
cycle = 0
116-
while t < solver.maxTime:
117-
if rank == 0 and cycle % 100 == 0:
118-
print( f"time = {t:.3f}s, dt = {solver.dt:.4f}, iter = {cycle+1}" )
119-
solver.execute( t )
120-
if cycle % 50 == 0:
121-
solver.outputVtk( t )
141+
cycle = 0
142+
while t < solver.maxTime :
143+
if rank == 0 and cycle%100 == 0:
144+
print(f"time = {t:.3f}s, dt = {solver.dt:.4f}, iter = {cycle+1}")
145+
solver.execute(t)
146+
if cycle%50 == 0:
147+
solver.outputVtk(t)
122148
t += solver.dt
123149
cycle += 1
124-
150+
125151
shot.flag = "Done"
126-
if rank == 0:
127-
print( f"Shot {shot.id} done" )
128-
print( "Gathering and exporting seismos" )
152+
if rank == 0 :
153+
print(f"Shot {shot.id} done")
154+
print("Gathering and exporting seismos")
129155

130156
seismos = solver.getPressureAtReceivers()
131-
157+
132158
directory = outdirseis
133159
filename = f"{outseisname}_{shot.id}"
134-
135-
SeismicTraceOutput( seismos, format="SEP" ).export( directory=directory,
136-
rootname=filename,
137-
receiverCoords=shot.getReceiverCoords(),
138-
sourceCoords=shot.getSourceCoords()[ 0 ],
139-
dt=solver.dtSeismo )
140-
160+
161+
SeismicTraceOutput(seismos, format="SEP").export(directory=directory, rootname=filename, receiverCoords=shot.getReceiverCoords(), sourceCoords=shot.getSourceCoords()[0], dt=solver.dtSeismo)
162+
141163
solver.resetWaveField()
142164

143-
144165
if __name__ == "__main__":
145166
main()

‎pygeos-tools/tests/elastic_modeling.py

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,10 @@
33
from geos.pygeos_tools.utilities.acquisition_library.Acquisition import Acquisition
44
from geos.pygeos_tools.utilities.solvers import ElasticSolver
55
from geos.pygeos_tools.utilities.output import SeismicTraceOutput
6-
from mpi4py import MPI
76

7+
import mpi4py
8+
mpi4py.rc.initialize = True
9+
from mpi4py import MPI
810

911
def parse_args():
1012
"""Get arguments
@@ -20,6 +22,13 @@ def parse_args():
2022

2123

2224
def main():
25+
26+
if MPI.Is_initialized():
27+
print("MPI initialized")
28+
else:
29+
print("MPI not initialized. Initializing...")
30+
MPI.Init()
31+
2332
comm = MPI.COMM_WORLD
2433
rank = comm.Get_rank()
2534

0 commit comments

Comments
 (0)