diff --git a/README.md b/README.md index 0081ffb0..3a4f1288 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,8 @@ To use Spyro, you'll need to have some knowledge of Python and some basic concep Discussions about development take place on our Slack channel. Everyone is invited to join using the link: https://join.slack.com/t/spyroworkspace/shared_invite/zt-u87ih28m-2h9JobfkdArs4ku3a1wLLQ +If you want to know more or cite our code please see our open access publication: https://gmd.copernicus.org/articles/15/8639/2022/gmd-15-8639-2022.html + Functionality ============= diff --git a/benchmarks/benchmark_2d.py b/benchmarks/benchmark_2d.py index 4b037521..b900368e 100644 --- a/benchmarks/benchmark_2d.py +++ b/benchmarks/benchmark_2d.py @@ -1,10 +1,13 @@ -import SeismicMesh import meshio import firedrake as fire import spyro import time +try: + import SeismicMesh +except ImportError: + pass -## Adding model parameters: +# Adding model parameters: model = {} model["opts"] = { @@ -16,7 +19,8 @@ model["parallelism"] = { "type": "spatial", # options: automatic, spatial, or custom. - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + "custom_cores_per_shot": [], # only if the user wants a different number + # of cores for every shot. # input is a list of integers with the length of the number of shots. } @@ -31,14 +35,15 @@ model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s "R": 1e-6, # theoretical reflection coefficient - "lz": 0.286, # thickness of the PML in the z-direction (km) - always positive - "lx": 0.286, # thickness of the PML in the x-direction (km) - always positive - "ly": 0.0, # thickness of the PML in the y-direction (km) - always positive + "lz": 0.286, # thickness of the PML in z-direction (km)-always positive + "lx": 0.286, # thickness of the PML in x-direction (km)-always positive + "ly": 0.0, # thickness of the PML in y-direction (km)-always positive } model["acquisition"] = { @@ -49,14 +54,14 @@ "delay": 1.0, "num_receivers": 15, "receiver_locations": spyro.create_transect( - (-5.72, 6.29), + (-5.72, 6.29), (-5.72, 14.29), 15), } # Simulate for 1.0 seconds. model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 4.0, # Final time for event "dt": 0.0005, # timestep size "amplitude": 1, # the Ricker has an amplitude of 1. @@ -66,11 +71,12 @@ comm = spyro.utils.mpi_init(model) -## Starting meshing procedure with seismic mesh. This can be done seperately in order to not have to +# Starting meshing procedure with seismic mesh. This can be done seperately +# in order to not have to # when testing multiple cores. -print('Entering mesh generation', flush = True) -if model['opts']['degree'] == 2: +print('Entering mesh generation', flush=True) +if model['opts']['degree'] == 2: M = 5.85 elif model['opts']['degree'] == 3: M = 3.08 @@ -84,25 +90,30 @@ Lx = model["mesh"]["Lx"] pad = model["BCs"]["lz"] -bbox = (-Real_Lz, 0.0,-pad, Lx+pad) +bbox = (-Real_Lz, 0.0, -pad, Lx+pad) rec = SeismicMesh.Rectangle(bbox) if comm.comm.rank == 0: points, cells = SeismicMesh.generate_mesh( - domain=rec, - edge_length=edge_length, - comm = comm.ensemble_comm, - verbose = 0 - ) - - points, cells = SeismicMesh.geometry.delete_boundary_entities(points, cells, min_qual= 0.6) - - meshio.write_points_cells("meshes/benchmark_2d.msh", - points,[("triangle", cells)], - file_format="gmsh22", - binary = False - ) - -## Mesh generation finishes here. + domain=rec, + edge_length=edge_length, + comm=comm.ensemble_comm, + verbose=0 + ) + + points, cells = SeismicMesh.geometry.delete_boundary_entities( + points, + cells, + min_qual=0.6, + ) + + meshio.write_points_cells( + "meshes/benchmark_2d.msh", + points, [("triangle", cells)], + file_format="gmsh22", + binary=False + ) + +# Mesh generation finishes here. mesh = fire.Mesh( "meshes/benchmark_2d.msh", @@ -117,7 +128,12 @@ if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0: print(f"Setting up {method} a {degree}tetra element", flush=True) -element = fire.FiniteElement(method, mesh.ufl_cell(), degree=degree, variant="KMV") +element = fire.FiniteElement( + method, + mesh.ufl_cell(), + degree=degree, + variant="KMV" +) V = fire.FunctionSpace(mesh, element) @@ -139,5 +155,13 @@ ) t1 = time.time() -p, p_r = spyro.solvers.forward(model, mesh, comm, vp, sources, wavelet, receivers) +p, p_r = spyro.solvers.forward( + model, + mesh, + comm, + vp, + sources, + wavelet, + receivers +) print(time.time() - t1, flush=True) diff --git a/benchmarks/benchmark_3d.py b/benchmarks/benchmark_3d.py index b53e6001..22b1f0d6 100644 --- a/benchmarks/benchmark_3d.py +++ b/benchmarks/benchmark_3d.py @@ -1,10 +1,13 @@ -import SeismicMesh import meshio import firedrake as fire import spyro import time +try: + import SeismicMesh +except ImportError: + pass -## Adding model parameters: +# Adding model parameters: model = {} model["opts"] = { @@ -16,7 +19,8 @@ model["parallelism"] = { "type": "spatial", # options: automatic, spatial, or custom. - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # only if the user wants a different number of cores for every shot. + "custom_cores_per_shot": [], # input is a list of integers with the length of the number of shots. } @@ -31,14 +35,15 @@ model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 1, # damping layer has a exponent variation "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s "R": 1e-3, # theoretical reflection coefficient - "lz": 0.286, # thickness of the PML in the z-direction (km) - always positive - "lx": 0.286, # thickness of the PML in the x-direction (km) - always positive - "ly": 0.286, # thickness of the PML in the y-direction (km) - always positive + "lz": 0.286, # thickness of the PML in z-direction (km) - always positive + "lx": 0.286, # thickness of the PML in x-direction (km) - always positive + "ly": 0.286, # thickness of the PML in y-direction (km) - always positive } model["acquisition"] = { @@ -56,7 +61,7 @@ # Simulate for 1.0 seconds. model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 1.0, # Final time for event "dt": 0.005, # timestep size "amplitude": 1, # the Ricker has an amplitude of 1. @@ -66,11 +71,12 @@ comm = spyro.utils.mpi_init(model) -## Starting meshing procedure with seismic mesh. This can be done seperately in order to not have to +# Starting meshing procedure with seismic mesh. This can be done seperately +# in order to not have to # when testing multiple cores. -print('Entering mesh generation', flush = True) -if model['opts']['degree'] == 2: +print('Entering mesh generation', flush=True) +if model['opts']['degree'] == 2: M = 5.1 elif model['opts']['degree'] == 3: M = 3.1 @@ -81,25 +87,34 @@ Ly = model["mesh"]["Ly"] pad = model["BCs"]["lz"] -bbox = (-Real_Lz, 0.0,-pad, Lx+pad, -pad, Ly+pad) +bbox = (-Real_Lz, 0.0, -pad, Lx+pad, -pad, Ly+pad) cube = SeismicMesh.Cube(bbox) points, cells = SeismicMesh.generate_mesh( - domain=cube, - edge_length=edge_length, - max_iter = 80, - comm = comm.ensemble_comm, - verbose = 2 - ) + domain=cube, + edge_length=edge_length, + max_iter=80, + comm=comm.ensemble_comm, + verbose=2 +) -points, cells = SeismicMesh.sliver_removal(points=points, bbox=bbox, max_iter=100, domain=cube, edge_length=edge_length, preserve=True) +points, cells = SeismicMesh.sliver_removal( + points=points, + bbox=bbox, + max_iter=100, + domain=cube, + edge_length=edge_length, + preserve=True +) -meshio.write_points_cells("meshes/benchmark_3d.msh", - points,[("tetra", cells)], - file_format="gmsh22", - binary = False - ) +meshio.write_points_cells( + "meshes/benchmark_3d.msh", + points, + [("tetra", cells)], + file_format="gmsh22", + binary=False +) -## Mesh generation finishes here. +# Mesh generation finishes here. mesh = fire.Mesh( "meshes/benchmark_3d.msh", @@ -114,7 +129,12 @@ if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0: print(f"Setting up {method} a {degree}tetra element", flush=True) -element = fire.FiniteElement(method, mesh.ufl_cell(), degree=degree, variant="KMV") +element = fire.FiniteElement( + method, + mesh.ufl_cell(), + degree=degree, + variant="KMV" +) V = fire.FunctionSpace(mesh, element) @@ -136,5 +156,13 @@ ) t1 = time.time() -p, p_r = spyro.solvers.forward(model, mesh, comm, vp, sources, wavelet, receivers) +p, p_r = spyro.solvers.forward( + model, + mesh, + comm, + vp, + sources, + wavelet, + receivers +) print(time.time() - t1, flush=True) diff --git a/benchmarks/benchmark_cuboid.py b/benchmarks/benchmark_cuboid.py index b564bc0b..4754e2db 100644 --- a/benchmarks/benchmark_cuboid.py +++ b/benchmarks/benchmark_cuboid.py @@ -32,7 +32,7 @@ model["PML"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 1.0, # maximum acoustic wave velocity in PML - km/s @@ -42,6 +42,7 @@ "ly": 1.0, # thickness of the PML in the y-direction (km) - always positive } + model["acquisition"] = { "source_type": "Ricker", "num_sources": 1, @@ -56,7 +57,7 @@ # Simulate for 1.0 seconds. model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 1.0, # Final time for event "dt": 0.0005, # timestep size "amplitude": 1, # the Ricker has an amplitude of 1. diff --git a/demos/run_forward.py b/demos/run_forward.py index 61559ecf..a2066dab 100644 --- a/demos/run_forward.py +++ b/demos/run_forward.py @@ -21,7 +21,7 @@ } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s @@ -40,7 +40,7 @@ "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500), } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 6.00, # Final time for event "dt": 0.001, "amplitude": 1, # the Ricker has an amplitude of 1. diff --git a/demos/run_fwi.py b/demos/run_fwi.py index 0be941c1..3a3c2277 100644 --- a/demos/run_fwi.py +++ b/demos/run_fwi.py @@ -20,7 +20,7 @@ "degree": 5, # p order "dimension": 2, # dimension "regularization": True, # regularization is on? - "gamma": 1.e-6, # regularization parameter + "gamma": 1.0e-6, # regularization parameter } model["parallelism"] = { "type": "automatic", @@ -35,7 +35,7 @@ } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s @@ -54,7 +54,7 @@ "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500), } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 6.00, # Final time for event "dt": 0.001, "amplitude": 1, # the Ricker has an amplitude of 1. @@ -179,8 +179,8 @@ def gradient(self, g, x, tol): if comm.comm.size > 1: dJ /= comm.comm.size # regularize the gradient if asked. - if model['opts']['regularization']: - gamma = model['opts']['gamma'] + if model["opts"]["regularization"]: + gamma = model["opts"]["gamma"] dJ = regularize_gradient(vp, dJ, gamma) # mask the water layer dJ.dat.data[water] = 0.0 diff --git a/paper/run_forward_2d.py b/paper/run_forward_2d.py index 0f05bd31..06fdeb91 100644 --- a/paper/run_forward_2d.py +++ b/paper/run_forward_2d.py @@ -22,7 +22,7 @@ } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s @@ -41,7 +41,7 @@ "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500), } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 5.00, # Final time for event "dt": 0.00025, "amplitude": 1, # the Ricker has an amplitude of 1. diff --git a/paper/run_forward_3d.py b/paper/run_forward_3d.py index 2f66dfce..18a0047a 100644 --- a/paper/run_forward_3d.py +++ b/paper/run_forward_3d.py @@ -35,7 +35,7 @@ } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 6.0, # maximum acoustic wave velocity in PML - km/s @@ -54,7 +54,7 @@ "receiver_locations": receivers, } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 4.00, # Final time for event "dt": 0.00075, "amplitude": 1, # the Ricker has an amplitude of 1. diff --git a/paper/run_fwi_2d.py b/paper/run_fwi_2d.py index 8b3d70e2..34d5c7c9 100644 --- a/paper/run_fwi_2d.py +++ b/paper/run_fwi_2d.py @@ -15,7 +15,7 @@ def get_memory_usage(): """Return the memory usage in Mo.""" process = psutil.Process(os.getpid()) - mem = process.memory_info()[0] / float(2 ** 20) + mem = process.memory_info()[0] / float(2**20) return mem @@ -46,7 +46,7 @@ def get_memory_usage(): } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s @@ -65,7 +65,7 @@ def get_memory_usage(): "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500), } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 5.00, # Final time for event "dt": 0.001, "amplitude": 1, # the Ricker has an amplitude of 1. @@ -117,8 +117,8 @@ def regularize_gradient(vp, dJ): """Tikhonov regularization""" m_u = TrialFunction(V) m_v = TestFunction(V) - mgrad = m_u * m_v * dx(rule=qr_x) - ffG = dot(grad(vp), grad(m_v)) * dx(rule=qr_x) + mgrad = m_u * m_v * dx(rule=quad_rule) + ffG = dot(grad(vp), grad(m_v)) * dx(rule=quad_rule) G = mgrad - ffG lhsG, rhsG = lhs(G), rhs(G) gradreg = Function(V) diff --git a/paper/run_fwi_3d.py b/paper/run_fwi_3d.py index fa2d0913..483fc111 100644 --- a/paper/run_fwi_3d.py +++ b/paper/run_fwi_3d.py @@ -30,7 +30,7 @@ def get_memory_usage(): """Return the memory usage in Mo.""" process = psutil.Process(os.getpid()) - mem = process.memory_info()[0] / float(2 ** 20) + mem = process.memory_info()[0] / float(2**20) return mem @@ -61,7 +61,7 @@ def get_memory_usage(): } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 6.0, # maximum acoustic wave velocity in PML - km/s @@ -80,7 +80,7 @@ def get_memory_usage(): "receiver_locations": receivers, } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 4.00, # Final time for event "dt": 0.00075, "amplitude": 1, # the Ricker has an amplitude of 1. @@ -139,8 +139,8 @@ def regularize_gradient(vp, dJ): """Tikhonov regularization""" m_u = TrialFunction(V) m_v = TestFunction(V) - mgrad = m_u * m_v * dx(rule=qr_x) - ffG = dot(grad(vp), grad(m_v)) * dx(rule=qr_x) + mgrad = m_u * m_v * dx(rule=quad_rule) + ffG = dot(grad(vp), grad(m_v)) * dx(rule=quad_rule) G = mgrad - ffG lhsG, rhsG = lhs(G), rhs(G) gradreg = Function(V) diff --git a/setup.py b/setup.py index b04d3b58..1fa5f94a 100644 --- a/setup.py +++ b/setup.py @@ -8,5 +8,10 @@ author="Keith J. Roberts, Alexandre F. G. Olender, Lucas Franceschini", url="https://github.com/krober10nd/spyro", packages=find_packages(), - install_requires=["numpy", "scipy", "matplotlib", "exdown==0.7.0", "segyio"], + install_requires=[ + "numpy", + "scipy", + "matplotlib", + "exdown==0.7.0", + "segyio"], ) diff --git a/spyro/__init__.py b/spyro/__init__.py index 68430d60..4c45979f 100644 --- a/spyro/__init__.py +++ b/spyro/__init__.py @@ -4,11 +4,12 @@ from .receivers.Receivers import Receivers from .sources.Sources import Sources, ricker_wavelet, full_ricker_wavelet from .utils import utils -from .utils.geometry_creation import create_transect, create_2d_grid, insert_fixed_value, create_3d_grid +from .utils.geometry_creation import create_transect, create_2d_grid +from .utils.geometry_creation import insert_fixed_value, create_3d_grid from .utils.estimate_timestep import estimate_timestep from .io import io from . import solvers -from .import tools +from . import tools __all__ = [ "io", diff --git a/spyro/domains/quadrature.py b/spyro/domains/quadrature.py index 47744cf1..0beaf886 100644 --- a/spyro/domains/quadrature.py +++ b/spyro/domains/quadrature.py @@ -4,7 +4,21 @@ def quadrature_rules(V): - """ Quadrature rule - Gauss-Lobatto-Legendre, Gauss-Legendre and Equi-spaced, KMV""" + """Returns quadrature rule - Gauss-Lobatto-Legendre, Gauss-Legendre and Equi-spaced, KMV + + Returns quadradure rule to use with UFL's dx object when integrating + + Parameters + ---------- + V : obj + UFL Function Space + + Returns + ------- + qr_x, qr_s, qr_k : obj + quadrature rules for Firedrake to use + """ + degree = V.ufl_element().degree() dimension = V.mesh().geometric_dimension() cell_geometry = V.mesh().ufl_cell() @@ -59,6 +73,7 @@ def quadrature_rules(V): # Spectral method - Gauss-Lobatto-Legendre rule # 1D def gauss_lobatto_legendre_line_rule(degree): + """Returns GLL quad rule for a given degree in a line""" fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) finat_ps = finat.point_set.GaussLobattoLegendrePointSet @@ -68,6 +83,7 @@ def gauss_lobatto_legendre_line_rule(degree): # 3D def gauss_lobatto_legendre_cube_rule(dimension, degree): + """Returns GLL quad rule for a given degree in a multidimensional space""" make_tensor_rule = finat.quadrature.TensorProductQuadratureRule result = gauss_lobatto_legendre_line_rule(degree) for _ in range(1, dimension): diff --git a/spyro/domains/space.py b/spyro/domains/space.py index 5e4130cd..b2331609 100644 --- a/spyro/domains/space.py +++ b/spyro/domains/space.py @@ -6,7 +6,7 @@ def FE_method(mesh, method, degree): Space discretization - Continuous or Discontinuous Galerkin methods""" cell_geometry = mesh.ufl_cell() - if method == "CG" or method == 'spectral': + if method == "CG" or method == "spectral": # CG - Continuous Galerkin if cell_geometry == quadrilateral or cell_geometry == hexahedron: element = FiniteElement( diff --git a/spyro/io/__init__.py b/spyro/io/__init__.py index 030752d7..add90fb3 100644 --- a/spyro/io/__init__.py +++ b/spyro/io/__init__.py @@ -1,15 +1,34 @@ -from .io import write_function_to_grid, create_segy, is_owner, save_shots, load_shots, read_mesh, interpolate, ensemble_forward, ensemble_forward_ad, ensemble_forward_elastic_waves, ensemble_gradient, ensemble_gradient_elastic_waves, ensemble_plot +from .io import ( + write_function_to_grid, + create_segy, + is_owner, + save_shots, + load_shots, + read_mesh, +) +from .io import ( + interpolate, + ensemble_forward, + ensemble_forward_ad, + ensemble_forward_elastic_waves, + ensemble_gradient, + ensemble_gradient_elastic_waves, + ensemble_plot, +) -__all__ = ["write_function_to_grid", - "create_segy", - "is_owner", - "save_shots", - "load_shots", - "read_mesh", - "interpolate", - "ensemble_forward", - "ensemble_forward_ad", - "ensemble_forward_elastic_waves", - "ensemble_gradient", - "ensemble_gradient_elastic_waves", - "ensemble_plot"] + +__all__ = [ + "write_function_to_grid", + "create_segy", + "is_owner", + "save_shots", + "load_shots", + "read_mesh", + "interpolate", + "ensemble_forward", + "ensemble_forward_ad", + "ensemble_forward_elastic_waves", + "ensemble_gradient", + "ensemble_gradient_elastic_waves", + "ensemble_plot", +] diff --git a/spyro/io/io.py b/spyro/io/io.py index 1f1bf485..5ba47cd8 100644 --- a/spyro/io/io.py +++ b/spyro/io/io.py @@ -1,6 +1,4 @@ from __future__ import with_statement - -import os import pickle import firedrake as fire @@ -15,52 +13,70 @@ def ensemble_save(func): """Decorator for read and write shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) _comm = args[1] - custom_file_name = kwargs.get('file_name') + custom_file_name = kwargs.get("file_name") for snum in range(num): if is_owner(_comm, snum) and _comm.comm.rank == 0: if custom_file_name is None: - func(*args, **dict(kwargs, file_name = "shots/shot_record_"+str(snum+1)+".dat")) + func( + *args, + **dict( + kwargs, + file_name="shots/shot_record_" + str(snum + 1) + ".dat", + ) + ) else: - func(*args, **dict(kwargs, file_name = custom_file_name)) + func(*args, **dict(kwargs, file_name=custom_file_name)) + return wrapper def ensemble_load(func): - """Decorator for read and write shots for ensemble parallelism""" + """Decorator for loading shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) _comm = args[1] - custom_file_name = kwargs.get('file_name') + custom_file_name = kwargs.get("file_name") for snum in range(num): if is_owner(_comm, snum): if custom_file_name is None: - values = func(*args, **dict(kwargs, file_name = "shots/shot_record_"+str(snum+1)+".dat")) + values = func( + *args, + **dict( + kwargs, + file_name="shots/shot_record_" + str(snum + 1) + ".dat", + ) + ) else: - values = func(*args, **dict(kwargs, file_name = custom_file_name)) - return values + values = func(*args, **dict(kwargs, file_name=custom_file_name)) + return values + return wrapper def ensemble_plot(func): """Decorator for `plot_shots` to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) _comm = args[1] for snum in range(num): if is_owner(_comm, snum) and _comm.comm.rank == 0: - func(*args, **dict(kwargs, file_name = str(snum+1))) + func(*args, **dict(kwargs, file_name=str(snum + 1))) return wrapper def ensemble_forward(func): """Decorator for forward to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) @@ -72,8 +88,10 @@ def wrapper(*args, **kwargs): return wrapper + def ensemble_forward_ad(func): - """Decorator for forward to distribute shots for ensemble parallelism""" + """Decorator for forward_ad to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) @@ -89,8 +107,10 @@ def wrapper(*args, **kwargs): return wrapper + def ensemble_forward_elastic_waves(func): """Decorator for forward elastic waves to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") num = len(acq["source_pos"]) @@ -105,6 +125,7 @@ def wrapper(*args, **kwargs): def ensemble_gradient(func): """Decorator for gradient to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") save_adjoint = kwargs.get("save_adjoint") @@ -124,6 +145,7 @@ def wrapper(*args, **kwargs): def ensemble_gradient_elastic_waves(func): """Decorator for gradient (elastic waves) to distribute shots for ensemble parallelism""" + def wrapper(*args, **kwargs): acq = args[0].get("acquisition") save_adjoint = kwargs.get("save_adjoint") @@ -174,8 +196,8 @@ def create_segy(velocity, filename): velocity = np.flipud(velocity.T) - spec.sorting = 2 # not sure what this means - spec.format = 1 # not sure what this means + spec.sorting = 2 # not sure what this means + spec.format = 1 # not sure what this means spec.samples = range(velocity.shape[0]) spec.ilines = range(velocity.shape[1]) spec.xlines = range(velocity.shape[0]) @@ -248,9 +270,8 @@ def is_owner(ens_comm, rank): return ens_comm.ensemble_comm.rank == (rank % ens_comm.ensemble_comm.size) - - def _check_units(c): + """Checks if velocity is in m/s or km/s""" if min(c.dat.data[:]) > 100.0: # data is in m/s but must be in km/s if fire.COMM_WORLD.rank == 0: diff --git a/spyro/plots/plots.py b/spyro/plots/plots.py index aecb98be..dade1732 100644 --- a/spyro/plots/plots.py +++ b/spyro/plots/plots.py @@ -1,11 +1,9 @@ # from scipy.io import savemat import matplotlib.pyplot as plt - import numpy as np - from ..io import ensemble_plot -__all__ = ["plot_shots"] +__all__ = ["plot_shots"] @ensemble_plot diff --git a/spyro/pml/damping.py b/spyro/pml/damping.py index da425224..9a812df2 100644 --- a/spyro/pml/damping.py +++ b/spyro/pml/damping.py @@ -21,7 +21,7 @@ def functions( y2=None, b_pml=None, ): - """ Damping functions for the perfect matched layer for 2D and 3D""" + """Damping functions for the perfect matched layer for 2D and 3D""" damping_type = model["BCs"]["damping_type"] if damping_type == "polynomial": @@ -94,13 +94,12 @@ def functions( ) ) sigma_y = Function(V, name="sigma_y").interpolate(aux1 + aux2) - - # sgm_y = File("pmlField/sigma_y.pvd") # sgm_y.write(sigma_y) return (sigma_x, sigma_y, sigma_z) + def matrices_2D(sigma_x, sigma_y): """Damping matrices for a two-dimensional problem""" Gamma_1 = as_tensor([[sigma_x, 0.0], [0.0, sigma_y]]) diff --git a/spyro/receivers/Receivers.py b/spyro/receivers/Receivers.py index e6b17599..49d6f516 100644 --- a/spyro/receivers/Receivers.py +++ b/spyro/receivers/Receivers.py @@ -11,6 +11,43 @@ class Receivers: """Interpolate data defined on a triangular mesh to a set of 2D/3D coordinates for variable spatial order using Lagrange interpolation. + + Can interpolate receiveir values that do not coincide with + mesh or DOF points + + ... + + Attributes + ---------- + mesh : Firedrake.mesh + mesh where receivers are located + V: Firedrake.FunctionSpace object + The space of the finite elements + my_ensemble: Firedrake.ensemble_communicator + An ensemble communicator + dimension: int + The dimension of the space + degree: int + Degree of the function space + receiver_locations: list + List of tuples containing all receiver locations + num_receivers: int + Number of receivers + quadrilateral: boolean + Boolean that specifies if cells are quadrilateral + is_local: list of booleans + List that checks if receivers are present in cores + spatial paralelism + + Methods + ------- + build_maps() + Calculates and stores tabulations for interpolation + interpolate(field) + Interpolates field value at receiver locations + apply_receivers_as_source(rhs_forcing, residual, IT) + Applies receivers as source with values from residual + in timestep IT, for usage with adjoint propagation """ def __init__(self, model, mesh, V, my_ensemble): @@ -42,13 +79,13 @@ def __init__(self, model, mesh, V, my_ensemble): self.dimension = model["opts"]["dimension"] self.degree = model["opts"]["degree"] self.receiver_locations = model["acquisition"]["receiver_locations"] - - if self.dimension==3 and model["aut_dif"]['status']: + + if self.dimension == 3 and model["aut_dif"]["status"]: self.column_x = model["acquisition"]["num_rec_x_columns"] self.column_y = model["acquisition"]["num_rec_y_columns"] self.column_z = model["acquisition"]["num_rec_z_columns"] - self.num_receivers = self.column_x*self.column_y - + self.num_receivers = self.column_x * self.column_y + else: self.num_receivers = len(self.receiver_locations) @@ -57,7 +94,7 @@ def __init__(self, model, mesh, V, my_ensemble): self.cell_tabulations = None self.cellNodeMaps = None self.nodes_per_cell = None - self.quadrilateral = (model["opts"]['quadrature']=='GLL') + self.quadrilateral = model["opts"]["quadrature"] == "GLL" self.is_local = [0] * self.num_receivers if not self.automatic_adjoint: self.build_maps() @@ -73,14 +110,25 @@ def num_receivers(self, value): self.__num_receivers = value def build_maps(self): + """Calculates and stores tabulations for interpolation + + Is always automatticaly called when initializing the class, + therefore should only be called again if a mesh related attribute + changes. + """ + for rid in range(self.num_receivers): tolerance = 1e-6 if self.dimension == 2: receiver_z, receiver_x = self.receiver_locations[rid] - cell_id = self.mesh.locate_cell([receiver_z, receiver_x], tolerance=tolerance ) + cell_id = self.mesh.locate_cell( + [receiver_z, receiver_x], tolerance=tolerance + ) elif self.dimension == 3: receiver_z, receiver_x, receiver_y = self.receiver_locations[rid] - cell_id = self.mesh.locate_cell([receiver_z, receiver_x, receiver_y], tolerance=tolerance ) + cell_id = self.mesh.locate_cell( + [receiver_z, receiver_x, receiver_y], tolerance=tolerance + ) self.is_local[rid] = cell_id ( @@ -108,20 +156,37 @@ def interpolate(self, field): return [self.__new_at(field, rn) for rn in range(self.num_receivers)] def apply_receivers_as_source(self, rhs_forcing, residual, IT): - """ - The adjoint operation of interpolation (injection) + """The adjoint operation of interpolation (injection) + + Injects residual, and timestep IT, at receiver locations + as source and stores their value in the right hand side + operator rhs_forcing. + + Parameters + ---------- + rhs_forcing: object + Firedrake assembled right hand side operator to store values + residual: list + List of residual values at different receiver locations + and timesteps + IT: int + Desired time step number to get residual value from + Returns + ------- + rhs_forcing: object + Firedrake assembled right hand side operator with injected values """ for rid in range(self.num_receivers): value = residual[IT][rid] if self.is_local[rid]: idx = np.int_(self.cellNodeMaps[rid]) phis = self.cell_tabulations[rid] - + tmp = np.dot(phis, value) rhs_forcing.dat.data_with_halos[idx] += tmp else: tmp = rhs_forcing.dat.data_with_halos[0] - + return rhs_forcing def __func_receiver_locator(self): @@ -180,10 +245,15 @@ def __func_receiver_locator_2D(self): cellNodeMaps = np.zeros((num_recv, nodes_per_cell)) cellVertices = [] - if self.quadrilateral == True: + if self.quadrilateral is True: end_vertex_id = 4 degree = self.degree - cell_ends = [0, (degree+1)*(degree+1)-degree-1, (degree+1)*(degree+1)-1, degree] + cell_ends = [ + 0, + (degree + 1) * (degree + 1) - degree - 1, + (degree + 1) * (degree + 1) - 1, + degree, + ] else: end_vertex_id = 3 cell_ends = [0, 1, 2] @@ -198,8 +268,12 @@ def __func_receiver_locator_2D(self): cellNodeMaps[receiver_id, :] = cell_node_map[cell_id, :] for vertex_number in range(0, end_vertex_id): cellVertices[receiver_id].append([]) - z = node_locations[cell_node_map[cell_id, cell_ends[vertex_number]], 0] - x = node_locations[cell_node_map[cell_id, cell_ends[vertex_number]], 1] + z = node_locations[ + cell_node_map[cell_id, cell_ends[vertex_number]], 0 + ] + x = node_locations[ + cell_node_map[cell_id, cell_ends[vertex_number]], 1 + ] cellVertices[receiver_id][vertex_number] = (z, x) return cellId_maps, cellVertices, cellNodeMaps @@ -301,14 +375,14 @@ def __func_node_locations_3D(self): return node_locations def __func_build_cell_tabulations(self): - if self.dimension == 2 and self.quadrilateral == False: + if self.dimension == 2 and self.quadrilateral is False: return self.__func_build_cell_tabulations_2D() - elif self.dimension == 3 and self.quadrilateral == False: + elif self.dimension == 3 and self.quadrilateral is False: return self.__func_build_cell_tabulations_3D() - elif self.dimension == 2 and self.quadrilateral == True: + elif self.dimension == 2 and self.quadrilateral is True: return self.__func_build_cell_tabulations_2D_quad() - elif self.dimension == 3 and self.quadrilateral == True: - raise ValueError('3D GLL hexas not yet supported.') + elif self.dimension == 3 and self.quadrilateral is True: + raise ValueError("3D GLL hexas not yet supported.") else: raise ValueError @@ -359,10 +433,12 @@ def __func_build_cell_tabulations_3D(self): return cell_tabulations def __func_build_cell_tabulations_2D_quad(self): - finatelement = FiniteElement('CG', self.mesh.ufl_cell(), degree=self.degree, variant='spectral') + finatelement = FiniteElement( + "CG", self.mesh.ufl_cell(), degree=self.degree, variant="spectral" + ) V = FunctionSpace(self.mesh, finatelement) u = TrialFunction(V) - Q=u.function_space() + Q = u.function_space() element = Q.finat_element.fiat_equivalent cell_tabulations = np.zeros((self.num_receivers, self.nodes_per_cell)) @@ -371,7 +447,7 @@ def __func_build_cell_tabulations_2D_quad(self): cell_id = self.is_local[receiver_id] if cell_id is not None: # getting coordinates to change to reference element - p = self.receiver_locations[receiver_id] + p = self.receiver_locations[receiver_id] v0 = self.cellVertices[receiver_id][0] v1 = self.cellVertices[receiver_id][1] v2 = self.cellVertices[receiver_id][2] @@ -383,43 +459,44 @@ def __func_build_cell_tabulations_2D_quad(self): cell_tabulations[receiver_id, :] = phi_tab.transpose() - return cell_tabulations - + def set_point_cloud(self, comm): # Receivers always parallel to z-axis rec_pos = self.receiver_locations - - # 2D -- - if self.dimension==2: + + # 2D -- + if self.dimension == 2: num_rec = self.num_receivers - δz = np.linspace(rec_pos[0,0], rec_pos[num_rec-1,0], 1) - δx = np.linspace(rec_pos[0,1], rec_pos[num_rec-1,1], num_rec) - - Z, X = np.meshgrid(δz,δx) - xs = np.vstack((Z.flatten(), X.flatten())).T - - #3D - elif self.dimension==3: - δz = np.linspace(rec_pos[0][0], rec_pos[1][0], self.column_z) - δx = np.linspace(rec_pos[0][1], rec_pos[1][1], self.column_x) - δy = np.linspace(rec_pos[0][2], rec_pos[1][2], self.column_y) - - Z, X, Y = np.meshgrid(δz,δx,δy) - xs = np.vstack((Z.flatten(),X.flatten(), Y.flatten())).T + δz = np.linspace(rec_pos[0, 0], rec_pos[num_rec - 1, 0], 1) + δx = np.linspace(rec_pos[0, 1], rec_pos[num_rec - 1, 1], num_rec) + + Z, X = np.meshgrid(δz, δx) + xs = np.vstack((Z.flatten(), X.flatten())).T + + # 3D + elif self.dimension == 3: + δz = np.linspace(rec_pos[0][0], rec_pos[1][0], self.column_z) + δx = np.linspace(rec_pos[0][1], rec_pos[1][1], self.column_x) + δy = np.linspace(rec_pos[0][2], rec_pos[1][2], self.column_y) + + Z, X, Y = np.meshgrid(δz, δx, δy) + xs = np.vstack((Z.flatten(), X.flatten(), Y.flatten())).T else: - print("This dimension is not accepted.") - quit() - - point_cloud = VertexOnlyMesh(self.mesh, xs) - + print("This dimension is not accepted.") + quit() + + point_cloud = VertexOnlyMesh(self.mesh, xs) + return point_cloud -## Some helper functions + +# Some helper functions def choosing_element(V, degree): + """Chooses UFL element based on desired function space""" cell_geometry = V.mesh().ufl_cell() if cell_geometry == quadrilateral: T = UFCQuadrilateral() @@ -447,6 +524,7 @@ def choosing_element(V, degree): def change_to_reference_triangle(p, a, b, c): + """Changes variables to reference triangle""" (xa, ya) = a (xb, yb) = b (xc, yc) = c @@ -488,6 +566,7 @@ def change_to_reference_triangle(p, a, b, c): def change_to_reference_tetrahedron(p, a, b, c, d): + """Changes variables to reference tetrahedron""" (xa, ya, za) = a (xb, yb, zb) = b (xc, yc, zc) = c @@ -747,7 +826,9 @@ def change_to_reference_tetrahedron(p, a, b, c, d): return (pnx, pny, pnz) + def change_to_reference_quad(p, v0, v1, v2, v3): + """Changes varibales to reference quadrilateral""" (px, py) = p # Irregular quad (x0, y0) = v0 @@ -772,41 +853,38 @@ def change_to_reference_quad(p, v0, v1, v2, v3): sumx = x0 - x1 + x2 - x3 sumy = y0 - y1 + y2 - y3 - gover = np.array([[sumx, dx2], - [sumy, dy2]]) + gover = np.array([[sumx, dx2], [sumy, dy2]]) - g_under= np.array([[dx1, dx2 ], - [dy1, dy2 ]]) + g_under = np.array([[dx1, dx2], [dy1, dy2]]) gunder = np.linalg.det(g_under) - - hover = np.array([[dx1, sumx], - [dy1, sumy]]) - hunder= gunder + hover = np.array([[dx1, sumx], [dy1, sumy]]) + + hunder = gunder - g = np.linalg.det(gover)/gunder - h = np.linalg.det(hover)/hunder + g = np.linalg.det(gover) / gunder + h = np.linalg.det(hover) / hunder i = 1.0 - a = x1 - x0 + g*x1 - b = x3 - x0 + h*x3 + a = x1 - x0 + g * x1 + b = x3 - x0 + h * x3 c = x0 - d = y1 - y0 + g*y1 - e = y3 - y0 + h*y3 + d = y1 - y0 + g * y1 + e = y3 - y0 + h * y3 f = y0 - A = e*i - f*h - B = c*h - b*i - C = b*f - c*e - D = f*g - d*i - E = a*i - c*g - F = c*d - a*f - G = d*h - e*g - H = b*g - a*h - I = a*e - b*d + A = e * i - f * h + B = c * h - b * i + C = b * f - c * e + D = f * g - d * i + E = a * i - c * g + F = c * d - a * f + G = d * h - e * g + H = b * g - a * h + I = a * e - b * d - pnx = (A*px + B*py + C)/(G*px + H*py + I) - pny = (D*px + E*py + F)/(G*px + H*py + I) + pnx = (A * px + B * py + C) / (G * px + H * py + I) + pny = (D * px + E * py + F) / (G * px + H * py + I) - return (pnx, pny) \ No newline at end of file + return (pnx, pny) diff --git a/spyro/solvers/forward.py b/spyro/solvers/forward.py index 627f48a2..6ca3f025 100644 --- a/spyro/solvers/forward.py +++ b/spyro/solvers/forward.py @@ -178,7 +178,7 @@ def forward( t = 0.0 # ------------------------------------------------------- - m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt ** 2)) * v * dx(rule=qr_x) + m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x) a = c * c * dot(grad(u_n), grad(v)) * dx(rule=qr_x) # explicit nf = 0 diff --git a/spyro/solvers/forward_AD.py b/spyro/solvers/forward_AD.py index 8ab55954..1ee1a995 100644 --- a/spyro/solvers/forward_AD.py +++ b/spyro/solvers/forward_AD.py @@ -1,10 +1,10 @@ from firedrake import * -from firedrake.assemble import create_assembly_callable -from .. import utils +# from .. import utils from ..domains import quadrature, space -from ..pml import damping -from ..io import ensemble_forward + +# from ..pml import damping +# from ..io import ensemble_forward from . import helpers # Note this turns off non-fatal warnings @@ -59,17 +59,13 @@ def forward( method = model["opts"]["method"] degree = model["opts"]["degree"] - dim = model["opts"]["dimension"] - dt = model["timeaxis"]["dt"] - tf = model["timeaxis"]["tf"] + dim = model["opts"]["dimension"] + dt = model["timeaxis"]["dt"] + tf = model["timeaxis"]["tf"] nspool = model["timeaxis"]["nspool"] - fspool = model["timeaxis"]["fspool"] - delay = model["acquisition"]["delay"] - dstep = int(delay / dt) # number of timesteps with source - PML = model["BCs"]["status"] - nt = int(tf / dt) # number of timesteps + nt = int(tf / dt) # number of timesteps excitations.current_source = source_num - params = set_params(method) + params = set_params(method) element = space.FE_method(mesh, method, degree) V = FunctionSpace(mesh, element) @@ -91,56 +87,49 @@ def forward( if output: outfile = helpers.create_output_file("forward.pvd", comm, source_num) - t = 0.0 - m = 1/(c*c) - m1 = m*((u - 2.0 * u_n + u_nm1) / Constant(dt ** 2)) * v * dx(rule=qr_x) - a = dot(grad(u_n), grad(v)) * dx(rule=qr_x) # explicit - f = Function(V) + t = 0.0 + m = 1 / (c * c) + m1 = m * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x) + a = dot(grad(u_n), grad(v)) * dx(rule=qr_x) # explicit + f = Function(V) nf = 0 if model["BCs"]["outer_bc"] == "non-reflective": nf = c * ((u_n - u_nm1) / dt) * v * ds(rule=qr_s) - - h = CellSize(mesh) - FF = m1 + a + nf - (1/(h/degree*h/degree))*f * v * dx(rule=qr_x) - X = Function(V) - B = Function(V) + + h = CellSize(mesh) + FF = m1 + a + nf - (1 / (h / degree * h / degree)) * f * v * dx(rule=qr_x) + X = Function(V) lhs_ = lhs(FF) rhs_ = rhs(FF) problem = LinearVariationalProblem(lhs_, rhs_, X) - solver = LinearVariationalSolver(problem, solver_parameters=params) + solver = LinearVariationalSolver(problem, solver_parameters=params) usol_recv = [] - save_step = 0 - - P = FunctionSpace(receivers, "DG", 0) + + P = FunctionSpace(receivers, "DG", 0) interpolator = Interpolator(u_np1, P) - J0 = 0.0 - + J0 = 0.0 + for step in range(nt): excitations.apply_source(f, wavelet[step]) - + solver.solve() u_np1.assign(X) rec = Function(P) interpolator.interpolate(output=rec) - - fwi = kwargs.get("fwi") + + fwi = kwargs.get("fwi") p_true_rec = kwargs.get("true_rec") - - usol_recv.append(rec.dat.data) + + usol_recv.append(rec.dat.data) if fwi: - J0 += calc_objective_func( - rec, - p_true_rec[step], - step, - dt, - P) + J0 += calc_objective_func(rec, p_true_rec[step], step, dt, P) if step % nspool == 0: assert ( @@ -155,19 +144,20 @@ def forward( u_n.assign(u_np1) t = step * float(dt) - + if fwi: return usol_recv, J0 else: return usol_recv -def calc_objective_func(p_rec,p_true_rec, IT, dt,P): - true_rec = Function(P) +def calc_objective_func(p_rec, p_true_rec, IT, dt, P): + true_rec = Function(P) true_rec.dat.data[:] = p_true_rec - J = 0.5 * assemble(inner(true_rec-p_rec, true_rec-p_rec) * dx) + J = 0.5 * assemble(inner(true_rec - p_rec, true_rec - p_rec) * dx) return J + def set_params(method): if method == "KMV": params = {"ksp_type": "preonly", "pc_type": "jacobi"} @@ -183,5 +173,5 @@ def set_params(method): params = {"ksp_type": "preonly", "pc_type": "jacobi"} else: raise ValueError("method is not yet supported") - + return params diff --git a/spyro/solvers/gradient.py b/spyro/solvers/gradient.py index f7f4d4ec..987888b9 100644 --- a/spyro/solvers/gradient.py +++ b/spyro/solvers/gradient.py @@ -91,7 +91,7 @@ def gradient( nt = int(tf / dt) # number of timesteps - receiver_locations = model["acquisition"]["receiver_locations"] + # receiver_locations = model["acquisition"]["receiver_locations"] dJ = Function(V, name="gradient") @@ -161,7 +161,7 @@ def gradient( t = 0.0 # ------------------------------------------------------- - m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt ** 2)) * v * dx(rule=qr_x) + m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x) a = c * c * dot(grad(u_n), grad(v)) * dx(rule=qr_x) # explicit nf = 0 @@ -235,7 +235,7 @@ def gradient( gradi = Function(V) grad_prob = LinearVariationalProblem(lhsG, rhsG, gradi) - + if method == "KMV": grad_solver = LinearVariationalSolver( grad_prob, @@ -279,10 +279,10 @@ def gradient( u_np1, psi_np1, pp_np1 = X.split() psi_nm1.assign(psi_n) - psi_n .assign(psi_np1) + psi_n.assign(psi_np1) pp_nm1.assign(pp_n) - pp_n .assign(pp_np1) + pp_n.assign(pp_np1) else: u_np1.assign(X) diff --git a/spyro/solvers/helpers.py b/spyro/solvers/helpers.py index 3c4ffa66..948a0d12 100644 --- a/spyro/solvers/helpers.py +++ b/spyro/solvers/helpers.py @@ -14,8 +14,9 @@ ] - def fill(usol_recv, is_local, nt, nr): + """Fills usol_recv with -99999 value + when it isn't local to any core""" usol_recv = np.asarray(usol_recv) for ti in range(nt): for rn in range(nr): @@ -25,6 +26,7 @@ def fill(usol_recv, is_local, nt, nr): def create_output_file(name, comm, source_num): + """Saves shots in output file""" if io.is_owner(comm, source_num): outfile = File( os.getcwd() @@ -39,6 +41,7 @@ def create_output_file(name, comm, source_num): def display(comm, source_num): + """Displays current shot and ensemble in terminal""" if comm.comm.rank == 0: print( "Timestepping for shot #", @@ -51,15 +54,23 @@ def display(comm, source_num): def display_progress(comm, t): + """Displays progress time""" if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0: print(f"Simulation time is: {t:{10}.{4}} seconds", flush=True) + def parallel_print(string, comm): + """Prints in parallel""" if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0: - print(string, flush = True) + print(string, flush=True) + def receivers_local(mesh, dimension, receiver_locations): + """Locates receiveirs in cells""" if dimension == 2: - return [mesh.locate_cell([z, x],tolerance=0.01) for z, x in receiver_locations] + return [mesh.locate_cell([z, x], tolerance=0.01) for z, x in receiver_locations] elif dimension == 3: - return [mesh.locate_cell([z, x, y],tolerance=0.01) for z, x, y in receiver_locations] + return [ + mesh.locate_cell([z, x, y], tolerance=0.01) + for z, x, y in receiver_locations + ] diff --git a/spyro/sources/Sources.py b/spyro/sources/Sources.py index 6b740c85..308e3f67 100644 --- a/spyro/sources/Sources.py +++ b/spyro/sources/Sources.py @@ -1,12 +1,45 @@ import math import numpy as np -from firedrake import * from scipy.signal import butter, filtfilt import spyro class Sources(spyro.receivers.Receivers.Receivers): - """Methods that inject a wavelet into a mesh""" + """Methods that inject a wavelet into a mesh + + ... + + Attributes + ---------- + mesh : Firedrake.mesh + mesh where receivers are located + V: Firedrake.FunctionSpace object + The space of the finite elements + my_ensemble: Firedrake.ensemble_communicator + An ensemble communicator + dimension: int + The dimension of the space + degree: int + Degree of the function space + source_locations: list + List of tuples containing all source locations + num_sources: int + Number of sources + quadrilateral: boolean + Boolean that specifies if cells are quadrilateral + is_local: list of booleans + List that checks if sources are present in cores + spatial paralelism + + Methods + ------- + build_maps() + Calculates and stores tabulations for interpolation + interpolate(field) + Interpolates field value at receiver locations + apply_source(rhs_forcing, value) + Applies value at source locations in rhs_forcing operator + """ def __init__(self, model, mesh, V, my_ensemble): """Initializes class and gets all receiver parameters from @@ -25,7 +58,7 @@ def __init__(self, model, mesh, V, my_ensemble): Returns ------- - Receivers: :class: 'Receiver' object + Sources: :class: 'Source' object """ @@ -43,24 +76,23 @@ def __init__(self, model, mesh, V, my_ensemble): self.cell_tabulations = None self.cellNodeMaps = None self.nodes_per_cell = None - self.is_local = [0]*self.num_receivers + self.is_local = [0] * self.num_receivers self.current_source = None - self.quadrilateral = (model["opts"]['quadrature']=='GLL') + self.quadrilateral = model["opts"]["quadrature"] == "GLL" super().build_maps() - def apply_source(self, rhs_forcing, value): """Applies source in a assembled right hand side.""" for source_id in range(self.num_receivers): - if self.is_local[source_id] and source_id==self.current_source: + if self.is_local[source_id] and source_id == self.current_source: for i in range(len(self.cellNodeMaps[source_id])): - rhs_forcing.dat.data_with_halos[int(self.cellNodeMaps[source_id][i])] = ( - value * self.cell_tabulations[source_id][i] - ) - else: + rhs_forcing.dat.data_with_halos[ + int(self.cellNodeMaps[source_id][i]) + ] = (value * self.cell_tabulations[source_id][i]) + else: for i in range(len(self.cellNodeMaps[source_id])): - tmp = rhs_forcing.dat.data_with_halos[0] + tmp = rhs_forcing.dat.data_with_halos[0] # noqa: F841 return rhs_forcing diff --git a/spyro/tools/__init__.py b/spyro/tools/__init__.py index e0e241f9..84f49c55 100644 --- a/spyro/tools/__init__.py +++ b/spyro/tools/__init__.py @@ -5,14 +5,15 @@ from .grid_point_calculator import grid_point_to_mesh_point_converter_for_seismicmesh from .grid_point_calculator import error_calc_line from .gradient_test_ad import gradient_test_acoustic as gradient_test_acoustic_ad + __all__ = [ "wave_solver", "generate_mesh", - 'error_calc', - 'create_model_for_grid_point_calculation', - 'minimum_grid_point_calculator', - 'time_interpolation', - 'grid_point_to_mesh_point_converter_for_seismicmesh', - 'error_calc_line', - 'gradient_test_acoustic_ad' + "error_calc", + "create_model_for_grid_point_calculation", + "minimum_grid_point_calculator", + "time_interpolation", + "grid_point_to_mesh_point_converter_for_seismicmesh", + "error_calc_line", + "gradient_test_acoustic_ad", ] diff --git a/spyro/tools/demo.py b/spyro/tools/demo.py index ec05d045..b037ad28 100644 --- a/spyro/tools/demo.py +++ b/spyro/tools/demo.py @@ -1,27 +1,26 @@ -### Demo to illustrate how a grid point density calcultor runs the experiments +# Demo to illustrate how a grid point density calcultor runs the experiments import spyro + # First we need to define experiment parameters: grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'homogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'velocity_model_file_name': 'vel_z6.25m_x12.5m_exact.segy', - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'near', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "homogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "velocity_model_file_name": "vel_z6.25m_x12.5m_exact.segy", + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "near", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 5, # Degree to use in the reference case (int) - 'G_reference': 15.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 6.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.05, - 'g_accuracy': 1e-1 - } + # Line search parameters + "reference_degree": 5, # Degree to use in the reference case (int) + "G_reference": 15.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 6.0, # Initial G for line search (float) + "accepted_error_threshold": 0.05, + "g_accuracy": 1e-1, +} G = spyro.tools.minimum_grid_point_calculator(grid_point_calculator_parameters) - diff --git a/spyro/tools/gradient_test_ad.py b/spyro/tools/gradient_test_ad.py index 191daa5c..a1c3ce51 100644 --- a/spyro/tools/gradient_test_ad.py +++ b/spyro/tools/gradient_test_ad.py @@ -2,17 +2,16 @@ from firedrake import * from pyadjoint import enlisting import spyro -from spyro.domains import quadrature -import matplotlib.pyplot as plt -import sys forward = spyro.solvers.forward_AD -def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): #{{{ + +def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): import firedrake_adjoint as fire_adj + with fire_adj.stop_annotating(): if comm.comm.rank == 0: - print('######## Starting gradient test ########', flush = True) + print("######## Starting gradient test ########", flush=True) sources = spyro.Sources(model, mesh, V, comm) receivers = spyro.Receivers(model, mesh, V, comm) @@ -25,50 +24,61 @@ def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): point_cloud = receivers.set_point_cloud(comm) # simulate the exact model if comm.comm.rank == 0: - print('######## Running the exact model ########', flush = True) + print("######## Running the exact model ########", flush=True) p_exact_recv = forward( model, mesh, comm, vp_exact, sources, wavelet, point_cloud ) - # simulate the guess model if comm.comm.rank == 0: - print('######## Running the guess model ########', flush = True) + print("######## Running the guess model ########", flush=True) p_guess_recv, Jm = forward( - model, mesh, comm, vp_guess, sources, wavelet, - point_cloud, fwi=True, true_rec=p_exact_recv + model, + mesh, + comm, + vp_guess, + sources, + wavelet, + point_cloud, + fwi=True, + true_rec=p_exact_recv, ) if comm.comm.rank == 0: - print("\n Cost functional at fixed point : " + str(Jm) + " \n ", flush = True) + print("\n Cost functional at fixed point : " + str(Jm) + " \n ", flush=True) # compute the gradient of the control (to be verified) if comm.comm.rank == 0: - print('######## Computing the gradient by automatic differentiation ########', flush = True) + print( + "######## Computing the gradient by automatic differentiation ########", + flush=True, + ) control = fire_adj.Control(vp_guess) - dJ = fire_adj.compute_gradient(Jm, control) + dJ = fire_adj.compute_gradient(Jm, control) if mask: dJ *= mask - + # File("gradient.pvd").write(dJ) - #steps = [1e-3, 1e-4, 1e-5, 1e-6, 1e-7] # step length - #steps = [1e-4, 1e-5, 1e-6, 1e-7] # step length + # steps = [1e-3, 1e-4, 1e-5, 1e-6, 1e-7] # step length + # steps = [1e-4, 1e-5, 1e-6, 1e-7] # step length steps = [1e-5, 1e-6, 1e-7, 1e-8] # step length with fire_adj.stop_annotating(): delta_m = Function(V) # model direction (random) delta_m.assign(dJ) - Jhat = fire_adj.ReducedFunctional(Jm, control) + Jhat = fire_adj.ReducedFunctional(Jm, control) derivative = enlisting.Enlist(Jhat.derivative()) hs = enlisting.Enlist(delta_m) - + projnorm = sum(hi._ad_dot(di) for hi, di in zip(hs, derivative)) - - # this deepcopy is important otherwise pertubations accumulate + # this deepcopy is important otherwise pertubations accumulate vp_original = vp_guess.copy(deepcopy=True) if comm.comm.rank == 0: - print('######## Computing the gradient by finite diferences ########', flush = True) + print( + "######## Computing the gradient by finite diferences ########", + flush=True, + ) errors = [] for step in steps: # range(3): # steps.append(step) @@ -76,10 +86,17 @@ def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): # J(m + delta_m*h) vp_guess = vp_original + step * delta_m p_guess_recv, Jp = forward( - model, mesh, comm, vp_guess, sources, wavelet, - point_cloud, fwi=True, true_rec=p_exact_recv + model, + mesh, + comm, + vp_guess, + sources, + wavelet, + point_cloud, + fwi=True, + true_rec=p_exact_recv, ) - + fd_grad = (Jp - Jm) / step if comm.comm.rank == 0: print( @@ -91,15 +108,14 @@ def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): + str(fd_grad) + ", grad'*dir : " + str(projnorm) - + " \n ", flush = True + + " \n ", + flush=True, ) errors.append(100 * ((fd_grad - projnorm) / projnorm)) - - + fire_adj.get_working_tape().clear_tape() # all errors less than 1 % errors = np.array(errors) assert (np.abs(errors) < 5.0).all() -#}}} diff --git a/spyro/tools/grid_point_calculator.py b/spyro/tools/grid_point_calculator.py index e3581940..7991f4d7 100644 --- a/spyro/tools/grid_point_calculator.py +++ b/spyro/tools/grid_point_calculator.py @@ -2,14 +2,14 @@ import numpy as np from scipy import interpolate import meshio -import SeismicMesh import firedrake as fire import time import copy import spyro + def minimum_grid_point_calculator(grid_point_calculator_parameters): - """ Function to calculate necessary grid point density. + """Function to calculate necessary grid point density. Parameters ---------- @@ -19,61 +19,67 @@ def minimum_grid_point_calculator(grid_point_calculator_parameters): Returns ------- G: `float` - Minimum grid point density necessary for a `experiment_type` mesh with a FEM whith + Minimum grid point density necessary for a `experiment_type` mesh with a FEM whith the degree and method specified within the specified error tolerance - """ - G_reference = grid_point_calculator_parameters['G_reference'] - degree_reference = grid_point_calculator_parameters['reference_degree'] - G_initial = grid_point_calculator_parameters['G_initial'] - desired_degree = grid_point_calculator_parameters['desired_degree'] - TOL = grid_point_calculator_parameters['accepted_error_threshold'] - - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, degree_reference) - #print("Model built at time "+str(time.time()-start_time), flush = True) + G_reference = grid_point_calculator_parameters["G_reference"] + degree_reference = grid_point_calculator_parameters["reference_degree"] + G_initial = grid_point_calculator_parameters["G_initial"] + desired_degree = grid_point_calculator_parameters["desired_degree"] + TOL = grid_point_calculator_parameters["accepted_error_threshold"] + + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, degree_reference + ) + # print("Model built at time "+str(time.time()-start_time), flush = True) comm = spyro.utils.mpi_init(model) - #print("Comm built at time "+str(time.time()-start_time), flush = True) + # print("Comm built at time "+str(time.time()-start_time), flush = True) if comm.comm.rank == 0: - print("Entering search", flush = True) - p_exact = wave_solver(model, G =G_reference, comm = comm) + print("Entering search", flush=True) + p_exact = wave_solver(model, G=G_reference, comm=comm) if comm.comm.rank == 0: - print("p_exact calculation finished", flush = True) + print("p_exact calculation finished", flush=True) comm.comm.barrier() - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, desired_degree) - G = searching_for_minimum(model, p_exact, TOL, starting_G=G_initial, comm = comm) + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, desired_degree + ) + G = searching_for_minimum(model, p_exact, TOL, starting_G=G_initial, comm=comm) return G -def wave_solver(model, G, comm = False): - minimum_mesh_velocity = model['testing_parameters']['minimum_mesh_velocity'] - model["mesh"]["meshfile"] = "meshes/2Dhomogeneous"+str(G)+".msh" + +def wave_solver(model, G, comm=False): + minimum_mesh_velocity = model["testing_parameters"]["minimum_mesh_velocity"] + model["mesh"]["meshfile"] = "meshes/2Dhomogeneous" + str(G) + ".msh" try: mesh, V = spyro.io.read_mesh(model, comm) - except: + except: # noqa E722 model = generate_mesh(model, G, comm) mesh, V = spyro.io.read_mesh(model, comm) - if model['testing_parameters']['experiment_type'] == 'homogeneous': + if model["testing_parameters"]["experiment_type"] == "homogeneous": vp_exact = fire.Constant(minimum_mesh_velocity) - elif model['testing_parameters']['experiment_type'] == 'heterogeneous': + elif model["testing_parameters"]["experiment_type"] == "heterogeneous": vp_exact = spyro.io.interpolate(model, mesh, V, guess=False) - if model["opts"]["method"] == 'KMV': - estimate_max_eigenvalue=True + if model["opts"]["method"] == "KMV": + estimate_max_eigenvalue = True else: - estimate_max_eigenvalue=False - new_dt = 0.2*spyro.estimate_timestep(mesh, V, vp_exact,estimate_max_eigenvalue=estimate_max_eigenvalue) + estimate_max_eigenvalue = False + new_dt = 0.2 * spyro.estimate_timestep( + mesh, V, vp_exact, estimate_max_eigenvalue=estimate_max_eigenvalue + ) - model['timeaxis']['dt'] = comm.comm.allreduce(new_dt, op=MPI.MIN) + model["timeaxis"]["dt"] = comm.comm.allreduce(new_dt, op=MPI.MIN) if comm.comm.rank == 0: print( f"Maximum stable timestep is: {model['timeaxis']['dt']} seconds", flush=True, ) - if model['timeaxis']['dt'] > 0.001: - model['timeaxis']['dt'] = 0.001 + if model["timeaxis"]["dt"] > 0.001: + model["timeaxis"]["dt"] = 0.001 if comm.comm.rank == 0: print( f"Timestep used is: {model['timeaxis']['dt']} seconds", @@ -83,76 +89,90 @@ def wave_solver(model, G, comm = False): sources = spyro.Sources(model, mesh, V, comm) receivers = spyro.Receivers(model, mesh, V, comm) wavelet = spyro.full_ricker_wavelet( - dt=model["timeaxis"]["dt"], - tf=model["timeaxis"]["tf"], - freq=model["acquisition"]["frequency"], - ) + dt=model["timeaxis"]["dt"], + tf=model["timeaxis"]["tf"], + freq=model["acquisition"]["frequency"], + ) for sn in range(model["acquisition"]["num_sources"]): if spyro.io.is_owner(comm, sn): t1 = time.time() p_field, p_recv = spyro.solvers.forward( - model, mesh, comm, vp_exact, sources, wavelet, receivers, source_num=sn, output= False) + model, + mesh, + comm, + vp_exact, + sources, + wavelet, + receivers, + source_num=sn, + output=False, + ) print(time.time() - t1) return p_recv -def generate_mesh(model,G, comm): + +def generate_mesh(model, G, comm): if model["opts"]["dimension"] == 2: - mesh = generate_mesh2D(model,G, comm) + mesh = generate_mesh2D(model, G, comm) elif model["opts"]["dimension"] == 3: - mesh = generate_mesh3D(model,G, comm) + mesh = generate_mesh3D(model, G, comm) else: raise ValueError("Wrong dimension in input model.") return mesh -def searching_for_minimum(model, p_exact, TOL, accuracy = 0.1, starting_G = 7.0, comm=False): + +def searching_for_minimum( + model, p_exact, TOL, accuracy=0.1, starting_G=7.0, comm=False +): error = 100.0 G = starting_G # fast loop - print("Entering fast loop", flush = True) + print("Entering fast loop", flush=True) while error > TOL: - dif = max(G*0.1, accuracy) + dif = max(G * 0.1, accuracy) G = G + dif - print('With G equal to '+str(G) ) - print("Entering wave solver", flush = True) - p0 = wave_solver(model,G, comm) - error = error_calc(p_exact, p0, model, comm = comm) - print('Error of '+str(error)) + print("With G equal to " + str(G)) + print("Entering wave solver", flush=True) + p0 = wave_solver(model, G, comm) + error = error_calc(p_exact, p0, model, comm=comm) + print("Error of " + str(error)) G -= dif - G = np.round(G,1)-accuracy + G = np.round(G, 1) - accuracy # slow loop - if dif > accuracy : - print("Entering slow loop", flush = True) + if dif > accuracy: + print("Entering slow loop", flush=True) error = 100.0 while error > TOL: dif = accuracy G = G + dif - print('With G equal to '+str(G) ) - print("Entering wave solver", flush = True) - p0 = wave_solver(model,G, comm ) - error = error_calc(p_exact, p0, model, comm = comm) - print('Error of '+str(error)) + print("With G equal to " + str(G)) + print("Entering wave solver", flush=True) + p0 = wave_solver(model, G, comm) + error = error_calc(p_exact, p0, model, comm=comm) + print("Error of " + str(error)) return G + def grid_point_to_mesh_point_converter_for_seismicmesh(model, G): - degree = model["opts"]['degree'] - if model["opts"]["method"] == 'KMV': + degree = model["opts"]["degree"] + if model["opts"]["method"] == "KMV": if degree == 1: - M = G/0.707813887967734 + M = G / 0.707813887967734 if degree == 2: - M = 0.5*G/0.8663141029672784 + M = 0.5 * G / 0.8663141029672784 if degree == 3: - M = 0.2934695559090401*G/0.7483761673104953 + M = 0.2934695559090401 * G / 0.7483761673104953 if degree == 4: - M = 0.21132486540518713*G/0.7010127254535244 + M = 0.21132486540518713 * G / 0.7010127254535244 if degree == 5: - M = 0.20231237605867816*G/0.9381929803311276 + M = 0.20231237605867816 * G / 0.9381929803311276 - if model["opts"]["method"] == 'CG': + if model["opts"]["method"] == "CG": raise ValueError("Correct M to G conversion to be inputed for CG") # if degree == 1: # M = G @@ -165,7 +185,7 @@ def grid_point_to_mesh_point_converter_for_seismicmesh(model, G): # if degree == 5: # M = 0.2*G - if model["opts"]["method"] == 'spectral': + if model["opts"]["method"] == "spectral": raise ValueError("Correct M to G conversion to be inputed for spectral") # if degree == 1: # M = G @@ -180,66 +200,65 @@ def grid_point_to_mesh_point_converter_for_seismicmesh(model, G): return M -def error_calc(p_exact, p, model, comm = False): + +def error_calc(p_exact, p, model, comm=False): # p0 doesn't necessarily have the same dt as p_exact # therefore we have to interpolate the missing points # to have them at the same length # testing shape - times_p_exact, r_p_exact = p_exact.shape - times_p, r_p = p.shape - if times_p_exact > times_p: #then we interpolate p_exact + times_p_exact, _ = p_exact.shape + times_p, _ = p.shape + if times_p_exact > times_p: # then we interpolate p_exact times, receivers = p.shape - dt = model["timeaxis"]['tf']/times + dt = model["timeaxis"]["tf"] / times p_exact = time_interpolation(p_exact, p, model) - elif times_p_exact < times_p: #then we interpolate p + elif times_p_exact < times_p: # then we interpolate p times, receivers = p_exact.shape - dt = model["timeaxis"]['tf']/times + dt = model["timeaxis"]["tf"] / times p = time_interpolation(p, p_exact, model) - else: #then we dont need to interpolate + else: # then we dont need to interpolate times, receivers = p.shape - dt = model["timeaxis"]['tf']/times - #p = time_interpolation(p, p_exact, model) + dt = model["timeaxis"]["tf"] / times + # p = time_interpolation(p, p_exact, model) - p_diff = p_exact-p max_absolute_diff = 0.0 max_percentage_diff = 0.0 - if comm.ensemble_comm.rank ==0: + if comm.ensemble_comm.rank == 0: numerator = 0.0 denominator = 0.0 for receiver in range(receivers): numerator_time_int = 0.0 denominator_time_int = 0.0 - for t in range(times-1): - top_integration = (p_exact[t,receiver]-p[t,receiver])**2*dt - bot_integration = (p_exact[t,receiver])**2*dt + for t in range(times - 1): + top_integration = (p_exact[t, receiver] - p[t, receiver]) ** 2 * dt + bot_integration = (p_exact[t, receiver]) ** 2 * dt # Adding 1e-25 filter to receivers to eliminate noise - numerator_time_int += top_integration + numerator_time_int += top_integration denominator_time_int += bot_integration - - diff = p_exact[t,receiver]-p[t,receiver] + diff = p_exact[t, receiver] - p[t, receiver] if abs(diff) > 1e-15 and abs(diff) > max_absolute_diff: max_absolute_diff = copy.deepcopy(diff) - - if abs(diff) > 1e-15 and abs(p_exact[t,receiver]) > 1e-15: - percentage_diff = abs( diff/p_exact[t,receiver] )*100 + + if abs(diff) > 1e-15 and abs(p_exact[t, receiver]) > 1e-15: + percentage_diff = abs(diff / p_exact[t, receiver]) * 100 if percentage_diff > max_percentage_diff: max_percentage_diff = copy.deepcopy(percentage_diff) - numerator += numerator_time_int + numerator += numerator_time_int denominator += denominator_time_int - + if denominator > 1e-15: - error = np.sqrt(numerator/denominator) + error = np.sqrt(numerator / denominator) # if numerator < 1e-15: # print('Warning: error too small to measure correctly.', flush = True) # error = 0.0 if denominator < 1e-15: - print("Warning: receivers don't appear to register a shot.", flush = True) + print("Warning: receivers don't appear to register a shot.", flush=True) error = 0.0 # print("ERROR IS ", flush = True) @@ -250,158 +269,165 @@ def error_calc(p_exact, p, model, comm = False): # print(max_percentage_diff, flush = True) return error -def error_calc_line(p_exact, p, model, comm = False): + +def error_calc_line(p_exact, p, model, comm=False): # p0 doesn't necessarily have the same dt as p_exact # therefore we have to interpolate the missing points # to have them at the same length # testing shape - times_p_exact, = p_exact.shape - times_p, = p.shape - if times_p_exact > times_p: #then we interpolate p_exact - times,= p.shape - dt = model["timeaxis"]['tf']/times + (times_p_exact,) = p_exact.shape + (times_p,) = p.shape + if times_p_exact > times_p: # then we interpolate p_exact + (times,) = p.shape + dt = model["timeaxis"]["tf"] / times p_exact = time_interpolation_line(p_exact, p, model) - elif times_p_exact < times_p: #then we interpolate p - times,= p_exact.shape - dt = model["timeaxis"]['tf']/times + elif times_p_exact < times_p: # then we interpolate p + (times,) = p_exact.shape + dt = model["timeaxis"]["tf"] / times p = time_interpolation_line(p, p_exact, model) - else: #then we dont need to interpolate - times, = p.shape - dt = model["timeaxis"]['tf']/times - + else: # then we dont need to interpolate + (times,) = p.shape + dt = model["timeaxis"]["tf"] / times - if comm.ensemble_comm.rank ==0: + if comm.ensemble_comm.rank == 0: numerator_time_int = 0.0 denominator_time_int = 0.0 # Integrating with trapezoidal rule - for t in range(times-1): - numerator_time_int += (p_exact[t]-p[t])**2 - denominator_time_int += (p_exact[t])**2 - numerator_time_int -= ((p_exact[0]-p[0])**2 + (p_exact[times-1]-p[times-1])**2)/2 + for t in range(times - 1): + numerator_time_int += (p_exact[t] - p[t]) ** 2 + denominator_time_int += (p_exact[t]) ** 2 + numerator_time_int -= ( + (p_exact[0] - p[0]) ** 2 + (p_exact[times - 1] - p[times - 1]) ** 2 + ) / 2 numerator_time_int *= dt - denominator_time_int -= (p_exact[0]**2+p_exact[times-1]**2)/2 + denominator_time_int -= (p_exact[0] ** 2 + p_exact[times - 1] ** 2) / 2 denominator_time_int *= dt - - #if denominator_time_int > 1e-15: - error = np.sqrt(numerator_time_int/denominator_time_int) - #if numerator_time_int < 1e-15: - # print('Warning: error too small to measure correctly.', flush = True) - #error = 0.0 + # if denominator_time_int > 1e-15: + error = np.sqrt(numerator_time_int / denominator_time_int) + if denominator_time_int < 1e-15: - print("Warning: receivers don't appear to register a shot.", flush = True) + print("Warning: receivers don't appear to register a shot.", flush=True) error = 0.0 return error + def time_interpolation(p_old, p_exact, model): times, receivers = p_exact.shape - dt = model["timeaxis"]['tf']/times + dt = model["timeaxis"]["tf"] / times times_old, rec = p_old.shape - dt_old = model["timeaxis"]['tf']/times_old - time_vector_old = np.zeros((1,times_old)) + dt_old = model["timeaxis"]["tf"] / times_old + time_vector_old = np.zeros((1, times_old)) for ite in range(times_old): - time_vector_old[0,ite] = dt_old*ite + time_vector_old[0, ite] = dt_old * ite - time_vector_new = np.zeros((1,times)) + time_vector_new = np.zeros((1, times)) for ite in range(times): - time_vector_new[0,ite] = dt*ite + time_vector_new[0, ite] = dt * ite p = np.zeros((times, receivers)) for receiver in range(receivers): - f = interpolate.interp1d(time_vector_old[0,:], p_old[:,receiver] ) - p[:,receiver] = f(time_vector_new[0,:]) + f = interpolate.interp1d(time_vector_old[0, :], p_old[:, receiver]) + p[:, receiver] = f(time_vector_new[0, :]) return p + def time_interpolation_line(p_old, p_exact, model): - times, = p_exact.shape - dt = model["timeaxis"]['tf']/times + (times,) = p_exact.shape + dt = model["timeaxis"]["tf"] / times - times_old, = p_old.shape - dt_old = model["timeaxis"]['tf']/times_old - time_vector_old = np.zeros((1,times_old)) + (times_old,) = p_old.shape + dt_old = model["timeaxis"]["tf"] / times_old + time_vector_old = np.zeros((1, times_old)) for ite in range(times_old): - time_vector_old[0,ite] = dt_old*ite + time_vector_old[0, ite] = dt_old * ite - time_vector_new = np.zeros((1,times)) + time_vector_new = np.zeros((1, times)) for ite in range(times): - time_vector_new[0,ite] = dt*ite + time_vector_new[0, ite] = dt * ite p = np.zeros((times,)) - f = interpolate.interp1d(time_vector_old[0,:], p_old[:] ) - p[:] = f(time_vector_new[0,:]) + f = interpolate.interp1d(time_vector_old[0, :], p_old[:]) + p[:] = f(time_vector_new[0, :]) return p -def generate_mesh2D(model,G, comm): + +def generate_mesh2D(model, G, comm): + import SeismicMesh if comm.comm.rank == 0: - print('Entering mesh generation', flush = True) + print("Entering mesh generation", flush=True) M = grid_point_to_mesh_point_converter_for_seismicmesh(model, G) - method = model["opts"]["method"] - Lz = model["mesh"]['Lz'] - lz = model['BCs']['lz'] - Lx = model["mesh"]['Lx'] - lx = model['BCs']['lx'] + Lz = model["mesh"]["Lz"] + lz = model["BCs"]["lz"] + Lx = model["mesh"]["Lx"] + lx = model["BCs"]["lx"] Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx + Real_Lx = Lx + 2 * lx + + if model["testing_parameters"]["experiment_type"] == "homogeneous": - if model['testing_parameters']['experiment_type']== 'homogeneous': - - minimum_mesh_velocity = model['testing_parameters']['minimum_mesh_velocity'] - frequency = model["acquisition"]['frequency'] - lbda = minimum_mesh_velocity/frequency + minimum_mesh_velocity = model["testing_parameters"]["minimum_mesh_velocity"] + frequency = model["acquisition"]["frequency"] + lbda = minimum_mesh_velocity / frequency Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx - edge_length = lbda/M + Real_Lx = Lx + 2 * lx + edge_length = lbda / M - bbox = (-Real_Lz, 0.0, -lx, Real_Lx-lx) + bbox = (-Real_Lz, 0.0, -lx, Real_Lx - lx) rec = SeismicMesh.Rectangle(bbox) if comm.comm.rank == 0: # Creating rectangular mesh points, cells = SeismicMesh.generate_mesh( - domain=rec, - edge_length=edge_length, - mesh_improvement = False, - comm = comm.ensemble_comm, - verbose = 0 + domain=rec, + edge_length=edge_length, + mesh_improvement=False, + comm=comm.ensemble_comm, + verbose=0, + ) + print("entering spatial rank 0 after mesh generation") + + points, cells = SeismicMesh.geometry.delete_boundary_entities( + points, cells, min_qual=0.6 + ) + + meshio.write_points_cells( + "meshes/2Dhomogeneous" + str(G) + ".msh", + points, + [("triangle", cells)], + file_format="gmsh22", + binary=False, + ) + meshio.write_points_cells( + "meshes/2Dhomogeneous" + str(G) + ".vtk", + points, + [("triangle", cells)], + file_format="vtk", ) - print('entering spatial rank 0 after mesh generation') - - points, cells = SeismicMesh.geometry.delete_boundary_entities(points, cells, min_qual= 0.6) - a=np.amin(SeismicMesh.geometry.simp_qual(points, cells)) - - meshio.write_points_cells("meshes/2Dhomogeneous"+str(G)+".msh", - points,[("triangle", cells)], - file_format="gmsh22", - binary = False - ) - meshio.write_points_cells("meshes/2Dhomogeneous"+str(G)+".vtk", - points,[("triangle", cells)], - file_format="vtk" - ) if comm.comm.rank == 0: - print('Finishing mesh generation', flush = True) + print("Finishing mesh generation", flush=True) - elif model['testing_parameters']['experiment_type']== 'heterogeneous': + elif model["testing_parameters"]["experiment_type"] == "heterogeneous": # Name of SEG-Y file containg velocity model. fname = "vel_z6.25m_x12.5m_exact.segy" # Bounding box describing domain extents (corner coordinates) bbox = (-12000.0, 0.0, 0.0, 67000.0) - rectangle =SeismicMesh.Rectangle(bbox) + rectangle = SeismicMesh.Rectangle(bbox) # Desired minimum mesh size in domain - frequency = model["acquisition"]['frequency'] - hmin = 1429.0/(M*frequency) + frequency = model["acquisition"]["frequency"] + hmin = 1429.0 / (M * frequency) # Construct mesh sizing object from velocity model ef = SeismicMesh.get_sizing_function_from_segy( @@ -415,96 +441,108 @@ def generate_mesh2D(model,G, comm): pad_style="edge", ) - points, cells = SeismicMesh.generate_mesh(domain=rectangle, edge_length=ef, verbose = 0, mesh_improvement=False ) + points, cells = SeismicMesh.generate_mesh( + domain=rectangle, edge_length=ef, verbose=0, mesh_improvement=False + ) - meshio.write_points_cells("meshes/2Dheterogeneous"+str(G)+".msh", - points/1000,[("triangle", cells)], - file_format="gmsh22", - binary = False - ) - meshio.write_points_cells("meshes/2Dheterogeneous"+str(G)+".vtk", - points/1000,[("triangle", cells)], - file_format="vtk" - ) + meshio.write_points_cells( + "meshes/2Dheterogeneous" + str(G) + ".msh", + points / 1000, + [("triangle", cells)], + file_format="gmsh22", + binary=False, + ) + meshio.write_points_cells( + "meshes/2Dheterogeneous" + str(G) + ".vtk", + points / 1000, + [("triangle", cells)], + file_format="vtk", + ) comm.comm.barrier() - if method == "CG" or method == "KMV": - mesh = fire.Mesh( - "meshes/2Dheterogeneous"+str(G)+".msh", - distribution_parameters={ - "overlap_type": (fire.DistributedMeshOverlapType.NONE, 0) - }, - ) - + return model + def generate_mesh3D(model, G, comm): + import SeismicMesh - print('Entering mesh generation', flush = True) + print("Entering mesh generation", flush=True) M = grid_point_to_mesh_point_converter_for_seismicmesh(model, G) method = model["opts"]["method"] - Lz = model["mesh"]['Lz'] - lz = model['BCs']['lz'] - Lx = model["mesh"]['Lx'] - lx = model['BCs']['lx'] - Ly = model["mesh"]['Ly'] - ly= model['BCs']['ly'] + Lz = model["mesh"]["Lz"] + lz = model["BCs"]["lz"] + Lx = model["mesh"]["Lx"] + lx = model["BCs"]["lx"] + Ly = model["mesh"]["Ly"] + ly = model["BCs"]["ly"] Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx - Real_Ly = Ly + 2*ly + Real_Lx = Lx + 2 * lx + Real_Ly = Ly + 2 * ly - minimum_mesh_velocity = model['testing_parameters']['minimum_mesh_velocity'] - frequency = model["acquisition"]['frequency'] - lbda = minimum_mesh_velocity/frequency + minimum_mesh_velocity = model["testing_parameters"]["minimum_mesh_velocity"] + frequency = model["acquisition"]["frequency"] + lbda = minimum_mesh_velocity / frequency - edge_length = lbda/M - #print(Real_Lz) + edge_length = lbda / M + # print(Real_Lz) - bbox = (-Real_Lz, 0.0, -lx, Real_Lx-lx, -ly, Real_Ly-ly) + bbox = (-Real_Lz, 0.0, -lx, Real_Lx - lx, -ly, Real_Ly - ly) cube = SeismicMesh.Cube(bbox) if comm.comm.rank == 0: # Creating rectangular mesh points, cells = SeismicMesh.generate_mesh( - domain=cube, - edge_length=edge_length, - mesh_improvement = False, - max_iter = 75, - comm = comm.ensemble_comm, - verbose = 0 + domain=cube, + edge_length=edge_length, + mesh_improvement=False, + max_iter=75, + comm=comm.ensemble_comm, + verbose=0, ) - points, cells = SeismicMesh.sliver_removal(points=points, bbox=bbox, domain=cube, edge_length=edge_length, preserve=True, max_iter=200) - - print('entering spatial rank 0 after mesh generation') - - meshio.write_points_cells("meshes/3Dhomogeneous"+str(G)+".msh", - points,[("tetra", cells)], - file_format="gmsh22", - binary = False - ) - meshio.write_points_cells("meshes/3Dhomogeneous"+str(G)+".vtk", - points,[("tetra", cells)], - file_format="vtk" - ) + points, cells = SeismicMesh.sliver_removal( + points=points, + bbox=bbox, + domain=cube, + edge_length=edge_length, + preserve=True, + max_iter=200, + ) + + print("entering spatial rank 0 after mesh generation") + + meshio.write_points_cells( + "meshes/3Dhomogeneous" + str(G) + ".msh", + points, + [("tetra", cells)], + file_format="gmsh22", + binary=False, + ) + meshio.write_points_cells( + "meshes/3Dhomogeneous" + str(G) + ".vtk", + points, + [("tetra", cells)], + file_format="vtk", + ) comm.comm.barrier() if method == "CG" or method == "KMV": mesh = fire.Mesh( - "meshes/3Dhomogeneous"+str(G)+".msh", + "meshes/3Dhomogeneous" + str(G) + ".msh", distribution_parameters={ "overlap_type": (fire.DistributedMeshOverlapType.NONE, 0) }, ) - - print('Finishing mesh generation', flush = True) + print("Finishing mesh generation", flush=True) return mesh + def mesh_generation(model, Gs, comm): for G in Gs: - mesh = generate_mesh(model, G, comm) - + _ = generate_mesh(model, G, comm) + return True diff --git a/spyro/tools/input_models.py b/spyro/tools/input_models.py index 3bbaf05e..fa435958 100644 --- a/spyro/tools/input_models.py +++ b/spyro/tools/input_models.py @@ -1,7 +1,7 @@ import numpy as np import spyro import shutil -import SeismicMesh + def create_3d_grid(start, end, num): """Create a 3d grid of `num**3` points between `start1` @@ -22,7 +22,7 @@ def create_3d_grid(start, end, num): """ (start1, start2, start3) = start - (end1, end2, end3) = end + (end1, end2, end3) = end x = np.linspace(start1, end1, num) y = np.linspace(start2, end2, num) z = np.linspace(start3, end3, num) @@ -30,10 +30,12 @@ def create_3d_grid(start, end, num): points = np.vstack((X.flatten(), Y.flatten(), Z.flatten())).T return [tuple(point) for point in points] + def create_model_2D_homogeneous(grid_point_calculator_parameters, degree): - ''' Creates models with the correct parameters for for grid point calculation experiments + """Creates models with the correct parameters for for grid point + calculation experiments on the 2D homogeneous case with a grid of receivers near the source. - + Parameters ---------- grid_point_calculator_parameters: Python 'dictionary' @@ -42,58 +44,66 @@ def create_model_2D_homogeneous(grid_point_calculator_parameters, degree): ------- model: Python `dictionary` Contains model options and parameters for use in Spyro - - ''' - minimum_mesh_velocity = grid_point_calculator_parameters['minimum_velocity_in_the_domain'] - frequency = grid_point_calculator_parameters['source_frequency'] - dimension = grid_point_calculator_parameters['dimension'] - receiver_type = grid_point_calculator_parameters['receiver_setup'] - method = grid_point_calculator_parameters['FEM_method_to_evaluate'] + """ + minimum_mesh_velocity = grid_point_calculator_parameters[ + "minimum_velocity_in_the_domain" + ] + frequency = grid_point_calculator_parameters["source_frequency"] + dimension = grid_point_calculator_parameters["dimension"] + receiver_type = grid_point_calculator_parameters["receiver_setup"] + + method = grid_point_calculator_parameters["FEM_method_to_evaluate"] model = {} if minimum_mesh_velocity > 500: - print("Warning: minimum mesh velocity seems to be in m/s, input should be in km/s", flush = True) + print( + "Warning: minimum mesh velocity seems to be in m/s, input should be in km/s", + flush=True, + ) # domain calculations pady = 0.0 Ly = 0.0 - lbda = minimum_mesh_velocity/frequency - pml_fraction = lbda + lbda = minimum_mesh_velocity / frequency pad = lbda - Lz = 40*lbda#100*lbda - Real_Lz = Lz+ pad - #print(Real_Lz) - Lx = 30*lbda#90*lbda - Real_Lx = Lx+ 2*pad + Lz = 40 * lbda # 100*lbda + Real_Lz = Lz + pad + # print(Real_Lz) + Lx = 30 * lbda # 90*lbda + Real_Lx = Lx + 2 * pad # source location - source_z = -Real_Lz/2.#1.0 - #print(source_z) - source_x = Real_Lx/2. - source_coordinates = [(source_z, source_x)] #Source at the center. If this is changes receiver's bin has to also be changed. + source_z = -Real_Lz / 2.0 # 1.0 + # print(source_z) + source_x = Real_Lx / 2.0 + # Source at the center. If this is changes receiver's bin has to also be + # changed. + source_coordinates = [(source_z, source_x)] padz = pad padx = pad # time calculations - tmin = 1./frequency - final_time = 20*tmin #should be 35 + tmin = 1.0 / frequency + final_time = 20 * tmin # Should be 35 # receiver calculations - receiver_bin_center1 = 10*lbda#20*lbda - receiver_bin_width = 5*lbda#15*lbda - receiver_quantity = 36#2500 # 50 squared + receiver_bin_center1 = 10 * lbda # 20*lbda + receiver_bin_width = 5 * lbda # 15*lbda + receiver_quantity = 36 # 2500 # 50 squared + + bin1_startZ = source_z + receiver_bin_center1 - receiver_bin_width / 2.0 + bin1_endZ = source_z + receiver_bin_center1 + receiver_bin_width / 2.0 + bin1_startX = source_x - receiver_bin_width / 2.0 + bin1_endX = source_x + receiver_bin_width / 2.0 - bin1_startZ = source_z + receiver_bin_center1 - receiver_bin_width/2. - bin1_endZ = source_z + receiver_bin_center1 + receiver_bin_width/2. - bin1_startX = source_x - receiver_bin_width/2. - bin1_endX = source_x + receiver_bin_width/2. + receiver_coordinates = spyro.create_2d_grid( + bin1_startZ, bin1_endZ, bin1_startX, bin1_endX, int(np.sqrt(receiver_quantity)) + ) - receiver_coordinates = spyro.create_2d_grid(bin1_startZ, bin1_endZ, bin1_startX, bin1_endX, int(np.sqrt(receiver_quantity))) - # Choose method and parameters model["opts"] = { "method": method, @@ -106,7 +116,7 @@ def create_model_2D_homogeneous(grid_point_calculator_parameters, degree): model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -138,29 +148,30 @@ def create_model_2D_homogeneous(grid_point_calculator_parameters, degree): } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": final_time, # Final time for event "dt": 0.001, # timestep size "nspool": 200, # how frequently to output solution to pvds "fspool": 100, # how frequently to save solution to RAM - } + } model["parallelism"] = { - "type": "spatial", # options: automatic (same number of cores for evey processor), custom, off. - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. + "type": "spatial", # options: automatic (same number of cores for evey processor), custom, off. + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. } - model['testing_parameters'] = { - 'minimum_mesh_velocity': minimum_mesh_velocity, - 'pml_fraction': padz/Lz, - 'receiver_type': receiver_type, - 'experiment_type': 'homogeneous' + model["testing_parameters"] = { + "minimum_mesh_velocity": minimum_mesh_velocity, + "pml_fraction": padz / Lz, + "receiver_type": receiver_type, + "experiment_type": "homogeneous", } return model + def create_model_2D_heterogeneous(grid_point_calculator_parameters, degree): - ''' Creates models with the correct parameters for for grid point calculation experiments. - + """Creates models with the correct parameters for for grid point calculation experiments. + Parameters ---------- frequency: `float` @@ -180,99 +191,121 @@ def create_model_2D_heterogeneous(grid_point_calculator_parameters, degree): ------- model: Python `dictionary` Contains model options and parameters for use in Spyro - - ''' - minimum_mesh_velocity = grid_point_calculator_parameters['minimum_velocity_in_the_domain'] - frequency = grid_point_calculator_parameters['source_frequency'] - dimension = grid_point_calculator_parameters['dimension'] - receiver_type = grid_point_calculator_parameters['receiver_setup'] - method = grid_point_calculator_parameters['FEM_method_to_evaluate'] - velocity_model = grid_point_calculator_parameters['velocity_model_file_name'] + """ + import SeismicMesh + + minimum_mesh_velocity = grid_point_calculator_parameters[ + "minimum_velocity_in_the_domain" + ] + frequency = grid_point_calculator_parameters["source_frequency"] + dimension = grid_point_calculator_parameters["dimension"] + receiver_type = grid_point_calculator_parameters["receiver_setup"] + + method = grid_point_calculator_parameters["FEM_method_to_evaluate"] + velocity_model = grid_point_calculator_parameters["velocity_model_file_name"] model = {} if minimum_mesh_velocity > 500: - print("Warning: minimum mesh velocity seems to be in m/s, input should be in km/s", flush = True) + print( + "Warning: minimum mesh velocity seems to be in m/s, input should be in km/s", + flush=True, + ) # domain calculations pady = 0.0 Ly = 0.0 - #using the BP2004 velocity model - - Lz = 12000.0/1000. - Lx = 67000.0/1000. - pad = 1000./1000. - Real_Lz = Lz+ pad - Real_Lx = Lx+ 2*pad + # using the BP2004 velocity model + + Lz = 12000.0 / 1000.0 + Lx = 67000.0 / 1000.0 + pad = 1000.0 / 1000.0 + Real_Lx = Lx + 2 * pad source_z = -1.0 - source_x = Real_Lx/2. - source_coordinates = [(source_z,source_x)] - if velocity_model != None: + source_x = Real_Lx / 2.0 + source_coordinates = [(source_z, source_x)] + if velocity_model is not None: if velocity_model[-4:] == "segy": - SeismicMesh.write_velocity_model(velocity_model, ofname = 'velocity_models/gridsweepcalc') + SeismicMesh.write_velocity_model( + velocity_model, ofname="velocity_models/gridsweepcalc" + ) elif velocity_model[-4:] == "hdf5": - shutil.copy(velocity_model,'velocity_models/gridsweepcalc.hdf5' ) + shutil.copy(velocity_model, "velocity_models/gridsweepcalc.hdf5") else: raise ValueError("Velocity model filetype not recognized.") else: - print("Warning: running without a velocity model is suitable for testing purposes only.", flush = True) + print( + "Warning: running without a velocity model is suitable for testing purposes only.", + flush=True, + ) padz = pad padx = pad - - if receiver_type == 'bins': + if receiver_type == "bins": # time calculations - tmin = 1./frequency - final_time = 25*tmin #should be 35 + tmin = 1.0 / frequency + final_time = 25 * tmin # should be 35 # receiver calculations - receiver_bin_center1 = 2.5*750.0/1000 - receiver_bin_width = 500.0/1000 - receiver_quantity_in_bin = 100#2500 # 50 squared + receiver_bin_center1 = 2.5 * 750.0 / 1000 + receiver_bin_width = 500.0 / 1000 + receiver_quantity_in_bin = 100 # 2500 # 50 squared - bin1_startZ = source_z - receiver_bin_width/2. - bin1_endZ = source_z + receiver_bin_width/2. - bin1_startX = source_x + receiver_bin_center1 - receiver_bin_width/2. - bin1_endX = source_x + receiver_bin_center1 + receiver_bin_width/2. + bin1_startZ = source_z - receiver_bin_width / 2.0 + bin1_endZ = source_z + receiver_bin_width / 2.0 + bin1_startX = source_x + receiver_bin_center1 - receiver_bin_width / 2.0 + bin1_endX = source_x + receiver_bin_center1 + receiver_bin_width / 2.0 - receiver_coordinates = spyro.create_2d_grid(bin1_startZ, bin1_endZ, bin1_startX, bin1_endX, int(np.sqrt(receiver_quantity_in_bin))) + receiver_coordinates = spyro.create_2d_grid( + bin1_startZ, + bin1_endZ, + bin1_startX, + bin1_endX, + int(np.sqrt(receiver_quantity_in_bin)), + ) - receiver_bin_center2 = 6500.0/1000 - receiver_bin_width = 500.0/1000 + receiver_bin_center2 = 6500.0 / 1000 + receiver_bin_width = 500.0 / 1000 - bin2_startZ = source_z - receiver_bin_width/2. - bin2_endZ = source_z + receiver_bin_width/2. - bin2_startX = source_x + receiver_bin_center2 - receiver_bin_width/2. - bin2_endX = source_x + receiver_bin_center2 + receiver_bin_width/2. + bin2_startZ = source_z - receiver_bin_width / 2.0 + bin2_endZ = source_z + receiver_bin_width / 2.0 + bin2_startX = source_x + receiver_bin_center2 - receiver_bin_width / 2.0 + bin2_endX = source_x + receiver_bin_center2 + receiver_bin_width / 2.0 - receiver_coordinates= receiver_coordinates + spyro.create_2d_grid(bin2_startZ, bin2_endZ, bin2_startX, bin2_endX, int(np.sqrt(receiver_quantity_in_bin))) + receiver_coordinates = receiver_coordinates + spyro.create_2d_grid( + bin2_startZ, + bin2_endZ, + bin2_startX, + bin2_endX, + int(np.sqrt(receiver_quantity_in_bin)), + ) - receiver_quantity = 2*receiver_quantity_in_bin + receiver_quantity = 2 * receiver_quantity_in_bin - - if receiver_type == 'line': + if receiver_type == "line": # time calculations - tmin = 1./frequency - final_time = 2*10*tmin + 5.0 #should be 35 + tmin = 1.0 / frequency + final_time = 2 * 10 * tmin + 5.0 # should be 35 # receiver calculations - receiver_bin_center1 = 2000.0/1000 - receiver_bin_center2 = 10000.0/1000 + receiver_bin_center1 = 2000.0 / 1000 + receiver_bin_center2 = 10000.0 / 1000 receiver_quantity = 500 - bin1_startZ = source_z - bin1_endZ = source_z + bin1_startZ = source_z + bin1_endZ = source_z bin1_startX = source_x + receiver_bin_center1 - bin1_endX = source_x + receiver_bin_center2 + bin1_endX = source_x + receiver_bin_center2 - receiver_coordinates = spyro.create_transect( (bin1_startZ, bin1_startX), (bin1_endZ, bin1_endX), receiver_quantity) + receiver_coordinates = spyro.create_transect( + (bin1_startZ, bin1_startX), (bin1_endZ, bin1_endX), receiver_quantity + ) - # Choose method and parameters model["opts"] = { "method": method, @@ -284,7 +317,7 @@ def create_model_2D_heterogeneous(grid_point_calculator_parameters, degree): model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -316,88 +349,93 @@ def create_model_2D_heterogeneous(grid_point_calculator_parameters, degree): } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": final_time, # Final time for event "dt": 0.001, # timestep size "nspool": 200, # how frequently to output solution to pvds "fspool": 100, # how frequently to save solution to RAM - } + } model["parallelism"] = { - "type": "off", # options: automatic (same number of cores for evey processor), custom, off. - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. + "type": "off", # options: automatic (same number of cores for evey processor), custom, off. + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. } - model['testing_parameters'] = { - 'minimum_mesh_velocity': minimum_mesh_velocity, - 'pml_fraction': padz/Lz, - 'receiver_type': receiver_type + model["testing_parameters"] = { + "minimum_mesh_velocity": minimum_mesh_velocity, + "pml_fraction": padz / Lz, + "receiver_type": receiver_type, } # print(source_coordinates) # print(receiver_coordinates) return model + def create_model_3D_homogeneous(grid_point_calculator_parameters, degree): - minimum_mesh_velocity = grid_point_calculator_parameters['minimum_velocity_in_the_domain'] - frequency = grid_point_calculator_parameters['source_frequency'] - dimension = grid_point_calculator_parameters['dimension'] - receiver_type = grid_point_calculator_parameters['receiver_setup'] + minimum_mesh_velocity = grid_point_calculator_parameters[ + "minimum_velocity_in_the_domain" + ] + frequency = grid_point_calculator_parameters["source_frequency"] + dimension = grid_point_calculator_parameters["dimension"] + + method = grid_point_calculator_parameters["FEM_method_to_evaluate"] - method = grid_point_calculator_parameters['FEM_method_to_evaluate'] - model = {} - lbda = minimum_mesh_velocity/frequency + lbda = minimum_mesh_velocity / frequency pad = lbda - Lz = 15*lbda#100*lbda - Real_Lz = Lz+ pad - #print(Real_Lz) - Lx = 30*lbda#90*lbda - Real_Lx = Lx+ 2*pad + Lz = 15 * lbda # 100*lbda + Real_Lz = Lz + pad + # print(Real_Lz) + Lx = 30 * lbda # 90*lbda Ly = Lx - Real_Ly = Ly + 2*pad + Real_Ly = Ly + 2 * pad # source location - source_z = -Real_Lz/2.#1.0 - #print(source_z) - source_x = lbda*1.5 - source_y = Real_Ly/2.0 - source_coordinates = [(source_z, source_x, source_y)] #Source at the center. If this is changes receiver's bin has to also be changed. + source_z = -Real_Lz / 2.0 # 1.0 + # print(source_z) + source_x = lbda * 1.5 + source_y = Real_Ly / 2.0 + source_coordinates = [ + (source_z, source_x, source_y) + ] # Source at the center. If this is changes receiver's bin has to also be changed. padz = pad padx = pad pady = pad # time calculations - tmin = 1./frequency - final_time = 20*tmin #should be 35 + tmin = 1.0 / frequency + final_time = 20 * tmin # should be 35 # receiver calculations - receiver_bin_center1 = 10*lbda#20*lbda - receiver_bin_width = 5*lbda#15*lbda - receiver_quantity = 36#2500 # 50 squared + receiver_bin_center1 = 10 * lbda # 20*lbda + receiver_bin_width = 5 * lbda # 15*lbda + receiver_quantity = 36 # 2500 # 50 squared - bin1_startZ = source_z - receiver_bin_width/2. - bin1_endZ = source_z + receiver_bin_width/2. - bin1_startX = source_x + receiver_bin_center1 - receiver_bin_width/2. - bin1_endX = source_x + receiver_bin_center1 + receiver_bin_width/2. - bin1_startY = source_y - receiver_bin_width/2. - bin1_endY = source_y + receiver_bin_width/2. + bin1_startZ = source_z - receiver_bin_width / 2.0 + bin1_endZ = source_z + receiver_bin_width / 2.0 + bin1_startX = source_x + receiver_bin_center1 - receiver_bin_width / 2.0 + bin1_endX = source_x + receiver_bin_center1 + receiver_bin_width / 2.0 + bin1_startY = source_y - receiver_bin_width / 2.0 + bin1_endY = source_y + receiver_bin_width / 2.0 - receiver_coordinates = create_3d_grid( (bin1_startZ,bin1_startX,bin1_startY) , (bin1_endZ,bin1_endX,bin1_endY) , 6) + receiver_coordinates = create_3d_grid( + (bin1_startZ, bin1_startX, bin1_startY), (bin1_endZ, bin1_endX, bin1_endY), 6 + ) # Choose method and parameters model["opts"] = { "method": method, "variant": None, "element": "tetra", # tria or tetra - 'quadrature': 'KMV', + "quadrature": "KMV", "degree": degree, # p order "dimension": dimension, # dimension } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -424,24 +462,25 @@ def create_model_3D_homogeneous(grid_point_calculator_parameters, degree): } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": final_time, # Final time for event "dt": 0.0002, # timestep size "nspool": 200, # how frequently to output solution to pvds "fspool": 100, # how frequently to save solution to RAM - } + } model["parallelism"] = { - "type": "spatial", + "type": "spatial", } # print(source_coordinates) # print(receiver_coordinates) return model + def create_model_for_grid_point_calculation(grid_point_calculator_parameters, degree): - ''' Creates models with the correct parameters for for grid point calculation experiments + """Creates models with the correct parameters for for grid point calculation experiments on the 2D homogeneous case with a grid of receivers near the source. - + Parameters ---------- grid_point_calculator_parameters: Python 'dictionary' @@ -450,16 +489,16 @@ def create_model_for_grid_point_calculation(grid_point_calculator_parameters, de ------- model: Python `dictionary` Contains model options and parameters for use in Spyro - - ''' - dimension = grid_point_calculator_parameters['dimension'] - experiment_type = grid_point_calculator_parameters['velocity_profile_type'] - if dimension == 2 and experiment_type == 'homogeneous': + + """ + dimension = grid_point_calculator_parameters["dimension"] + experiment_type = grid_point_calculator_parameters["velocity_profile_type"] + if dimension == 2 and experiment_type == "homogeneous": model = create_model_2D_homogeneous(grid_point_calculator_parameters, degree) - elif dimension == 2 and experiment_type == 'heterogeneous': + elif dimension == 2 and experiment_type == "heterogeneous": model = create_model_2D_heterogeneous(grid_point_calculator_parameters, degree) elif dimension == 3: model = create_model_3D_homogeneous(grid_point_calculator_parameters, degree) - - return model \ No newline at end of file + + return model diff --git a/spyro/utils/estimate_timestep.py b/spyro/utils/estimate_timestep.py index 6e3e4053..d77399cc 100644 --- a/spyro/utils/estimate_timestep.py +++ b/spyro/utils/estimate_timestep.py @@ -31,7 +31,7 @@ def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True): Asp = scipy.sparse.csr_matrix((av, aj, ai)) Asp_inv = scipy.sparse.csr_matrix((av_inv, aj, ai)) - K = fd.assemble(c*c*dot(grad(u), grad(v)) * dxlump) + K = fd.assemble(c * c * dot(grad(u), grad(v)) * dxlump) ai, aj, av = K.petscmat.getValuesCSR() Ksp = scipy.sparse.csr_matrix((av, aj, ai)) @@ -51,11 +51,11 @@ def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True): # print(max_eigval) if np.sqrt(max_eigval) > 0.0: - max_dt = np.float(2 / np.sqrt(max_eigval)) + max_dt = np.float(2 / np.sqrt(max_eigval)) else: max_dt = 100000000 - #print( + # print( # f"Maximum stable timestep should be about: {np.float(2 / np.sqrt(max_eigval))} seconds", # flush=True, - #) + # ) return max_dt diff --git a/spyro/utils/geometry_creation.py b/spyro/utils/geometry_creation.py index 1d700a02..4ec22636 100644 --- a/spyro/utils/geometry_creation.py +++ b/spyro/utils/geometry_creation.py @@ -50,6 +50,7 @@ def create_2d_grid(start1, end1, start2, end2, num): points = np.vstack((X.flatten(), Y.flatten())).T return [tuple(point) for point in points] + def create_3d_grid(start, end, num): """Create a 3d grid of `num**3` points between `start1` and `end1` and `start2` and `end2` @@ -69,7 +70,7 @@ def create_3d_grid(start, end, num): """ (start1, start2, start3) = start - (end1, end2, end3) = end + (end1, end2, end3) = end x = np.linspace(start1, end1, num) y = np.linspace(start2, end2, num) z = np.linspace(start3, end3, num) @@ -77,6 +78,7 @@ def create_3d_grid(start, end, num): points = np.vstack((X.flatten(), Y.flatten(), Z.flatten())).T return [tuple(point) for point in points] + def insert_fixed_value(points, value, insert): """Insert a fixed `value` in each tuple at index `insert`. diff --git a/spyro/utils/utils.py b/spyro/utils/utils.py index 6ab0ec7b..354fe9ba 100644 --- a/spyro/utils/utils.py +++ b/spyro/utils/utils.py @@ -45,9 +45,9 @@ def compute_functional(model, residual, velocity=None): J += residual[ti][rn] ** 2 J *= 0.5 - if regularize: - Jreg = assemble(0.5 * gamma * dot(grad(vp), grad(vp)) * dx) - J += Jreg + # if regularize: + # Jreg = assemble(0.5 * gamma * dot(grad(vp), grad(vp)) * dx) + # J += Jreg return J @@ -75,8 +75,8 @@ def mysize(COMM=COMM_SELF): def mpi_init(model): """Initialize computing environment""" - rank = myrank() - size = mysize() + # rank = myrank() + # size = mysize() available_cores = COMM_WORLD.size if model["parallelism"]["type"] == "automatic": num_cores_per_shot = available_cores / len(model["acquisition"]["source_pos"]) @@ -112,7 +112,7 @@ def communicate(array, my_ensemble): """ array_reduced = copy.copy(array) - + if my_ensemble.comm.size > 1: if my_ensemble.comm.rank == 0 and my_ensemble.ensemble_comm.rank == 0: print("Spatial parallelism, reducing to comm 0", flush=True) @@ -125,5 +125,5 @@ def analytical_solution_for_pressure_based_on_MMS(model, mesh, time): degree = model["opts"]["degree"] V = FunctionSpace(mesh, "CG", degree) z, x = SpatialCoordinate(mesh) - p = Function(V).interpolate((time ** 2) * sin(pi * z) * sin(pi * x)) + p = Function(V).interpolate((time**2) * sin(pi * z) * sin(pi * x)) return p diff --git a/test/inputfiles/Model1_2d_CG.py b/test/inputfiles/Model1_2d_CG.py index b0cb7553..7fb9e4fe 100644 --- a/test/inputfiles/Model1_2d_CG.py +++ b/test/inputfiles/Model1_2d_CG.py @@ -14,7 +14,7 @@ # Choose method and parameters opts = { "method": "CG", - "quadrature": 'KMV', + "quadrature": "KMV", "variant": None, "type": "SIP", # for DG only - SIP, NIP and IIP "degree": 1, # p order @@ -25,10 +25,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. # input is a list of integers with the length of the number of shots. - } +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -41,7 +41,7 @@ BCs = { "status": False, # True, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -60,7 +60,7 @@ } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.10, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -72,7 +72,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options model = { diff --git a/test/inputfiles/Model1_2d_DG.py b/test/inputfiles/Model1_2d_DG.py index 9ec9ab8b..e5f405a9 100644 --- a/test/inputfiles/Model1_2d_DG.py +++ b/test/inputfiles/Model1_2d_DG.py @@ -25,10 +25,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -41,7 +41,7 @@ PML = { "status": False, # True, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -60,7 +60,7 @@ } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.10, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -72,7 +72,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/Model1_3d_CG.py b/test/inputfiles/Model1_3d_CG.py index e2b44b14..1678f8f4 100644 --- a/test/inputfiles/Model1_3d_CG.py +++ b/test/inputfiles/Model1_3d_CG.py @@ -14,7 +14,7 @@ # Choose method and parameters opts = { "method": "CG", - "quadrature": 'KMV', + "quadrature": "KMV", "variant": None, "type": "SIP", # for DG only - SIP, NIP and IIP "degree": 1, # p order @@ -25,10 +25,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. # input is a list of integers with the length of the number of shots. - } +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -41,7 +41,7 @@ BCs = { "status": False, # True or false - "outer_bc": "non-reflective", # Neuumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # Neuumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -56,11 +56,11 @@ "frequency": 2.0, "delay": 1.0, "source_pos": [(-0.5, 0.5, 0.5)], # z,x,y - "receiver_locations":[()], + "receiver_locations": [()], } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.20, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -72,7 +72,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options model = { diff --git a/test/inputfiles/Model1_3d_DG.py b/test/inputfiles/Model1_3d_DG.py index 8207cc12..4438c8ba 100644 --- a/test/inputfiles/Model1_3d_DG.py +++ b/test/inputfiles/Model1_3d_DG.py @@ -25,10 +25,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -41,7 +41,7 @@ PML = { "status": False, # True or false - "outer_bc": "non-reflective", # Neuumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # Neuumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -56,11 +56,11 @@ "frequency": 2.0, "delay": 1.0, "source_pos": [(-0.5, 0.5, 0.5)], # z,x,y - "receiver_locations":[()], + "receiver_locations": [()], } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.20, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -72,7 +72,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/Model1_gradient_2d.py b/test/inputfiles/Model1_gradient_2d.py index 3dac78d1..fc3ec86a 100644 --- a/test/inputfiles/Model1_gradient_2d.py +++ b/test/inputfiles/Model1_gradient_2d.py @@ -54,7 +54,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 1.0, # Final time for event "dt": 0.001, # timestep size "nspool": 9999, # how frequently to output solution to pvds @@ -66,7 +66,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options model = { diff --git a/test/inputfiles/Model1_gradient_2d_pml.py b/test/inputfiles/Model1_gradient_2d_pml.py index b923ee80..a457aee1 100644 --- a/test/inputfiles/Model1_gradient_2d_pml.py +++ b/test/inputfiles/Model1_gradient_2d_pml.py @@ -31,7 +31,7 @@ } BCs = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 1.0, # maximum acoustic wave velocity in PML - km/s @@ -49,7 +49,7 @@ "receiver_locations": spyro.create_transect((-0.95, 0.1), (-0.95, 0.9), 100), } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 1.0, # Final time for event "dt": 0.001, # timestep size "amplitude": 1, # the Ricker has an amplitude of 1. @@ -58,7 +58,7 @@ } aut_dif = { - "status": False, + "status": False, } model_pml = { diff --git a/test/inputfiles/Model1_gradient_3d.py b/test/inputfiles/Model1_gradient_3d.py index 747d1c4b..2beecdf4 100644 --- a/test/inputfiles/Model1_gradient_3d.py +++ b/test/inputfiles/Model1_gradient_3d.py @@ -35,7 +35,7 @@ BCs = { "status": True, # True, # True or false - "outer_bc": "none", # neumann, non-reflective (outer boundary condition) + "outer_bc": "none", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 2, "cmax": 3.0, # maximum acoustic wave velocity in PML - km/s @@ -57,7 +57,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.5, # Final time for event "dt": 0.0005, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -69,7 +69,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/Model1_gradient_3d_pml.py b/test/inputfiles/Model1_gradient_3d_pml.py index c39518b8..f07cba1d 100644 --- a/test/inputfiles/Model1_gradient_3d_pml.py +++ b/test/inputfiles/Model1_gradient_3d_pml.py @@ -31,7 +31,7 @@ } BCs = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.0, # maximum acoustic wave velocity in PML - km/s @@ -41,7 +41,9 @@ "ly": 0.20, # thickness of the PML in the y-direction (km) - always positive } -receivers = spyro.insert_fixed_value(spyro.create_2d_grid(0.1,0.9,0.1,0.9,10),-0.9, 0) +receivers = spyro.insert_fixed_value( + spyro.create_2d_grid(0.1, 0.9, 0.1, 0.9, 10), -0.9, 0 +) print(len(receivers)) acquisition = { @@ -49,10 +51,10 @@ "source_pos": [(-0.1, 0.5, 0.5)], "frequency": 5.0, "delay": 1.0, - "receiver_locations": receivers, + "receiver_locations": receivers, } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 1.00, # Final time for event "dt": 0.0001, # timestep size "amplitude": 1, # the Ricker has an amplitude of 1. @@ -61,7 +63,7 @@ } aut_dif = { - "status": False, + "status": False, } model_pml = { diff --git a/test/inputfiles/Model1_parallel_2d.py b/test/inputfiles/Model1_parallel_2d.py index 296195ea..d74852ee 100644 --- a/test/inputfiles/Model1_parallel_2d.py +++ b/test/inputfiles/Model1_parallel_2d.py @@ -54,7 +54,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 2.0, # Final time for event "dt": 0.001, # timestep size "nspool": 9999, # how frequently to output solution to pvds @@ -66,7 +66,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/ModelMMS_2d_CG.py b/test/inputfiles/ModelMMS_2d_CG.py index d4a9cde1..27ebcae0 100644 --- a/test/inputfiles/ModelMMS_2d_CG.py +++ b/test/inputfiles/ModelMMS_2d_CG.py @@ -13,7 +13,7 @@ # Choose method and parameters opts = { "method": "CG", - "quadrature": 'KMV', + "quadrature": "KMV", "variant": None, "type": "SIP", # for DG only - SIP, NIP and IIP "degree": 1, # p order @@ -24,10 +24,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -40,7 +40,7 @@ PML = { "status": False, # True or false - "outer_bc": "dirichlet", # neumann, non-reflective (outer boundary condition) + "outer_bc": "dirichlet", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -59,7 +59,7 @@ } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.4, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -71,7 +71,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/ModelMMS_2d_DG.py b/test/inputfiles/ModelMMS_2d_DG.py index bf8c97bb..e4fd0c11 100644 --- a/test/inputfiles/ModelMMS_2d_DG.py +++ b/test/inputfiles/ModelMMS_2d_DG.py @@ -24,10 +24,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -40,7 +40,7 @@ PML = { "status": False, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -59,7 +59,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.4, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -71,7 +71,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/ModelMMS_3d_CG.py b/test/inputfiles/ModelMMS_3d_CG.py index 6842b62f..6790234e 100644 --- a/test/inputfiles/ModelMMS_3d_CG.py +++ b/test/inputfiles/ModelMMS_3d_CG.py @@ -13,7 +13,7 @@ # Choose method and parameters opts = { "method": "CG", - "quadrature": 'KMV', + "quadrature": "KMV", "variant": None, "type": "SIP", # for DG only - SIP, NIP and IIP "degree": 1, # p order @@ -24,10 +24,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -40,7 +40,7 @@ PML = { "status": False, # True or false - "outer_bc": "dirichlet", # neumann, non-reflective (outer boundary condition) + "outer_bc": "dirichlet", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -59,7 +59,7 @@ } # equi-spaced for now timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.4, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -71,7 +71,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/inputfiles/ModelMMS_3d_DG.py b/test/inputfiles/ModelMMS_3d_DG.py index d718e44d..ee545bee 100644 --- a/test/inputfiles/ModelMMS_3d_DG.py +++ b/test/inputfiles/ModelMMS_3d_DG.py @@ -24,10 +24,10 @@ } parallelism = { - "type": "automatic", # options: automatic, custom, off - "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. - # input is a list of integers with the length of the number of shots. - } + "type": "automatic", # options: automatic, custom, off + "custom_cores_per_shot": [], # only if the user wants a different number of cores for every shot. + # input is a list of integers with the length of the number of shots. +} mesh = { "Lz": 2.000, # depth in km - always positive @@ -40,7 +40,7 @@ PML = { "status": False, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -59,7 +59,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.4, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -71,7 +71,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } # Create your model with all the options diff --git a/test/model.py b/test/model.py index 6dcb679e..222474fc 100644 --- a/test/model.py +++ b/test/model.py @@ -9,7 +9,7 @@ # Choose method and parameters opts = { "method": "DG", - "quadrature": 'KMV', + "quadrature": "KMV", "variant": None, "type": "SIP", # for DG only - SIP, NIP and IIP "degree": 1, # p order @@ -34,7 +34,7 @@ BCs = { "status": False, # True or false - "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # neumann, non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial. hyperbolic, shifted_hyperbolic "exponent": 1, "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s @@ -55,7 +55,7 @@ } timeaxis = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.4, # Final time for event "dt": 0.001, # timestep size "nspool": 20, # how frequently to output solution to pvds @@ -67,7 +67,7 @@ } # cutoff frequencies (Hz) for Ricker source and to low-pass the observed shot record aut_dif = { - "status": False, + "status": False, } diff --git a/test/not_a_test.py b/test/not_a_test.py index b9bfbe2b..38a388a2 100644 --- a/test/not_a_test.py +++ b/test/not_a_test.py @@ -106,14 +106,14 @@ def gradient(self, g, x, tol): def update(self, x, flag, iteration): vp_guess.assign(Function(V, x.vec, name="velocity")) - paramsDict = { - "Step": { - "Line Search": {"Descent Method": {"Type": "Quasi-Newton Method"}}, - "Type": "Line Search", - }, - "Status Test": {"Gradient Tolerance": 1e-12, "Iteration Limit": 20}, - } - params = ROL.ParameterList(paramsDict, "Parameters") + # paramsDict = { + # "Step": { + # "Line Search": {"Descent Method": {"Type": "Quasi-Newton Method"}}, + # "Type": "Line Search", + # }, + # "Status Test": {"Gradient Tolerance": 1e-12, "Iteration Limit": 20}, + # } + # params = ROL.ParameterList(paramsDict, "Parameters") inner_product = L2Inner() obj = Objective(inner_product) diff --git a/test/test_MMS.py b/test/test_MMS.py index 3ce4977e..facbecec 100644 --- a/test/test_MMS.py +++ b/test/test_MMS.py @@ -6,7 +6,8 @@ from .model import model -@pytest.fixture(params=["triangle", "tetrahedral","square"]) + +@pytest.fixture(params=["triangle", "tetrahedral", "square"]) def mesh_type(request): return request.param @@ -20,7 +21,7 @@ def mesh(mesh_type): ) model["acquisition"]["source_pos"] = [(-0.05, 1.5)] elif mesh_type == "square": - model["opts"]['quadrature']=='GLL' + model["opts"]["quadrature"] == "GLL" model["opts"]["dimension"] = 2 model["acquisition"]["receiver_locations"] = spyro.create_transect( (0.0, 1.0), (0.0, 0.9), 256 @@ -34,9 +35,9 @@ def mesh(mesh_type): model["acquisition"]["source_pos"] = [(-0.05, 1.5, 1.5)] return { - "triangle": lambda n: UnitSquareMesh(2 ** n, 2 ** n), - "tetrahedral": lambda n: UnitCubeMesh(2 ** n, 2 ** n, 2 ** n), - "square": lambda n: UnitSquareMesh(2 ** n, 2 ** n, quadrilateral = True), + "triangle": lambda n: UnitSquareMesh(2**n, 2**n), + "tetrahedral": lambda n: UnitCubeMesh(2**n, 2**n, 2**n), + "square": lambda n: UnitSquareMesh(2**n, 2**n, quadrilateral=True), }[mesh_type] @@ -64,21 +65,24 @@ def timestep_method(timestep_method_type): @pytest.fixture def interpolation_expr(mesh_type): return { - "square": lambda x, y: (0.10 ** 2) * sin(pi * x) * sin(pi * y), - "triangle": lambda x, y: (0.10 ** 2) * sin(pi * x) * sin(pi * y), - "tetrahedral": lambda x, y, z: (0.10 ** 2) * sin(pi * x) * sin(pi * y) * sin(pi * z), + "square": lambda x, y: (0.10**2) * sin(pi * x) * sin(pi * y), + "triangle": lambda x, y: (0.10**2) * sin(pi * x) * sin(pi * y), + "tetrahedral": lambda x, y, z: (0.10**2) + * sin(pi * x) + * sin(pi * y) + * sin(pi * z), }[mesh_type] def run_solve(timestep_method, method, model, mesh, expr): testmodel = deepcopy(model) cell_geometry = mesh.ufl_cell() - if method =="CG" or method == 'spectral': + if method == "CG" or method == "spectral": if cell_geometry == quadrilateral or cell_geometry == hexahedron: - variant="spectral" - testmodel["opts"]["quadrature"]="GLL" + variant = "spectral" + testmodel["opts"]["quadrature"] = "GLL" else: - variant="equispaced" + variant = "equispaced" elif method == "KMV": variant = "KMV" @@ -104,6 +108,7 @@ def run_solve(timestep_method, method, model, mesh, expr): expr = expr(*SpatialCoordinate(mesh)) return errornorm(interpolate(expr, V), p[-1]) + def test_method(mesh, timestep_method, spatial_method, interpolation_expr): if mesh(3).ufl_cell() == quadrilateral and spatial_method == "KMV": pytest.skip("KMV isn't possible in quadrilaterals.") diff --git a/test/test_gradient.py b/test/test_gradient.py index 304a9a38..9422112f 100644 --- a/test/test_gradient.py +++ b/test/test_gradient.py @@ -1,5 +1,4 @@ import numpy as np -import pytest from firedrake import * import spyro from spyro.domains import quadrature @@ -43,12 +42,15 @@ def _make_vp_guess(V, mesh): File("guess_vel.pvd").write(vp_guess) return vp_guess + def test_gradient(): _test_gradient(model) + def test_gradient_pml(): _test_gradient(model_pml, pml=True) + def _test_gradient(options, pml=False): comm = spyro.utils.mpi_init(options) @@ -61,8 +63,6 @@ def _test_gradient(options, pml=False): Lx = model_pml["mesh"]["Lx"] Lz = model_pml["mesh"]["Lz"] x1 = 0.0 - x2 = Lx - z1 = 0.0 z2 = -Lz boxx1 = Function(V).interpolate(conditional(x > x1, 1.0, 0.0)) boxx2 = Function(V).interpolate(conditional(x < Lx, 1.0, 0.0)) diff --git a/test/test_gradient_3d.py b/test/test_gradient_3d.py index 80becc64..8ba0dbba 100644 --- a/test/test_gradient_3d.py +++ b/test/test_gradient_3d.py @@ -1,5 +1,5 @@ import numpy as np -import pytest +import pytest from firedrake import * @@ -31,40 +31,38 @@ def _make_vp_guess(V, mesh): File("guess_vel.pvd").write(vp_guess) return vp_guess + @pytest.mark.skip(reason="no way of currently testing this") def test_gradient_3d(): _test_gradient(model_pml, pml=True) + def _test_gradient(options, pml=False): comm = spyro.utils.mpi_init(options) mesh, V = spyro.io.read_mesh(options, comm) - vp_exact = _make_vp_exact_pml(V, mesh) - #dt = spyro.estimate_timestep(mesh, V, vp_exact) - #print(dt, flush=True) + # dt = spyro.estimate_timestep(mesh, V, vp_exact) + # print(dt, flush=True) z, x, y = SpatialCoordinate(mesh) Lx = model_pml["mesh"]["Lx"] Lz = model_pml["mesh"]["Lz"] Ly = model_pml["mesh"]["Ly"] x1 = 0.0 - x2 = Lx - z1 = 0.0 z2 = -Lz - y1 = 0.0 - y2 = Ly - boxx1 = Function(V).interpolate(conditional(x>x1 ,1.0,0.0)) - boxx2 = Function(V).interpolate( conditional(x x1, 1.0, 0.0)) + boxx2 = Function(V).interpolate(conditional(x < Lx, 1.0, 0.0)) - boxy1 = Function(V).interpolate(conditional(y>y1 ,1.0,0.0)) - boxy2 = Function(V).interpolate( conditional(y y1, 1.0, 0.0)) + boxy2 = Function(V).interpolate(conditional(y < Ly, 1.0, 0.0)) - boxz1 = Function(V).interpolate( conditional(z>z2,1.0,0.0) ) - box1 = Function(V).interpolate( boxx1 * boxx2 * boxz1 * boxy1 * boxy2 ) + boxz1 = Function(V).interpolate(conditional(z > z2, 1.0, 0.0)) + box1 = Function(V).interpolate(boxx1 * boxx2 * boxz1 * boxy1 * boxy2) File("inner-product-energy.pvd").write(box1) vp_guess = _make_vp_guess(V, mesh) @@ -112,18 +110,18 @@ def _test_gradient(options, pml=False): dJ = gradient(options, mesh, comm, vp_guess, receivers, p_guess, misfit) File("gradient.pvd").write(dJ) - steps = [1e-3, 1e-4, 1e-5]#, 1e-6] # step length + steps = [1e-3, 1e-4, 1e-5] # , 1e-6] # step length delta_m = Function(V) # model direction (random) - #delta_m.vector()[:] = 0.2 # np.random.rand(V.dim()) + # delta_m.vector()[:] = 0.2 # np.random.rand(V.dim()) delta_m.assign(dJ) # this deepcopy is important otherwise pertubations accumulate vp_original = vp_guess.copy(deepcopy=True) errors = [] - for step in steps: #range(3): - #steps.append(step) + for step in steps: # range(3): + # steps.append(step) # perturb the model and calculate the functional (again) # J(m + delta_m*h) vp_guess = vp_original + step * delta_m @@ -138,7 +136,7 @@ def _test_gradient(options, pml=False): ) Jp = functional(options, p_exact_recv - p_guess_recv) - projnorm = assemble(box1*dJ * delta_m * dx(rule=qr_x)) + projnorm = assemble(box1 * dJ * delta_m * dx(rule=qr_x)) fd_grad = (Jp - Jm) / step print( "\n Cost functional for step " @@ -153,12 +151,12 @@ def _test_gradient(options, pml=False): ) errors.append(100 * ((fd_grad - projnorm) / projnorm)) - #step /= 2 + # step /= 2 # all errors less than 1 % errors = np.array(errors) assert (np.abs(errors) < 1.0).all() -if __name__ == "__main__": - test_gradient_3d(model) +# if __name__ == "__main__": +# test_gradient_3d(model) diff --git a/test/test_io.py b/test/test_io.py index 9e7a2165..b8e473b8 100644 --- a/test/test_io.py +++ b/test/test_io.py @@ -2,29 +2,37 @@ import math import pytest import spyro -from SeismicMesh import write_velocity_model +try: + from SeismicMesh import write_velocity_model + + seismic_mesh = True +except ImportError: + seismic_mesh = False + + +@pytest.mark.skipif(not seismic_mesh, reason="No SeismicMesh") def test_read_and_write_segy(): vp_name = "velocity_models/test" - segy_file = vp_name+".segy" - hdf5_file = vp_name+".hdf5" - mesh = fire.UnitSquareMesh(10,10) - mesh.coordinates.dat.data[:,0] *= -1 + segy_file = vp_name + ".segy" + hdf5_file = vp_name + ".hdf5" + mesh = fire.UnitSquareMesh(10, 10) + mesh.coordinates.dat.data[:, 0] *= -1 - V = fire.FunctionSpace(mesh, 'CG', 3) + V = fire.FunctionSpace(mesh, "CG", 3) x, y = fire.SpatialCoordinate(mesh) r = 0.2 xc = -0.5 - yc = 0.5 + yc = 0.5 vp = fire.Function(V) - c = fire.conditional( (x-xc)**2+(y-yc)**2 < r**2 , 3.0 , 1.5) + c = fire.conditional((x - xc) ** 2 + (y - yc) ** 2 < r**2, 3.0, 1.5) vp.interpolate(c) - xi, yi, zi = spyro.io.write_function_to_grid(vp,V, 10.0/1000.0) - spyro.io.create_segy(zi,segy_file) + xi, yi, zi = spyro.io.write_function_to_grid(vp, V, 10.0 / 1000.0) + spyro.io.create_segy(zi, segy_file) write_velocity_model(segy_file, vp_name) model = {} @@ -47,18 +55,16 @@ def test_read_and_write_segy(): "status": False, } - vp_read = spyro.io.interpolate(model, mesh,V, guess = False) + vp_read = spyro.io.interpolate(model, mesh, V, guess=False) fire.File("velocity_models/test.pvd").write(vp_read) - value_at_center = vp_read.at(xc,yc) - test1 = math.isclose(value_at_center,3.0) - value_outside_circle = vp_read.at(xc+r+0.1,yc) - test2 = math.isclose(value_outside_circle,1.5) + value_at_center = vp_read.at(xc, yc) + test1 = math.isclose(value_at_center, 3.0) + value_outside_circle = vp_read.at(xc + r + 0.1, yc) + test2 = math.isclose(value_outside_circle, 1.5) assert all([test1, test2]) - -if __name__ == "__main__": - test_read_and_write_segy() - +if __name__ == "__main__": + test_read_and_write_segy() diff --git a/test/test_newat.py b/test/test_newat.py index 985f27f3..c915afb6 100644 --- a/test/test_newat.py +++ b/test/test_newat.py @@ -1,7 +1,6 @@ import math import numpy as np from copy import deepcopy -import pytest from firedrake import * import spyro @@ -17,6 +16,7 @@ def triangle_area(p1, p2, p3): return abs(x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2 + def test_correct_receiver_location_generation2D(): """Tests if receiver locations where generated correctly""" comm = spyro.utils.mpi_init(model) @@ -27,6 +27,7 @@ def test_correct_receiver_location_generation2D(): assert np.allclose(receivers, answer) + def test_correct_receiver_to_cell_location2D(): """Tests if the receivers where located in the correct cell""" comm = spyro.utils.mpi_init(model) @@ -86,6 +87,7 @@ def test_correct_receiver_to_cell_location2D(): assert all([test1, test2, test3]) + def test_correct_at_value2D(): comm = spyro.utils.mpi_init(model) model["opts"]["degree"] = 3 @@ -117,6 +119,7 @@ def test_correct_at_value2D(): assert all([test1, test2]) + def test_correct_at_value2D_quad(): model_quad = deepcopy(model) comm = spyro.utils.mpi_init(model_quad) @@ -150,6 +153,7 @@ def test_correct_at_value2D_quad(): assert all([test1, test2]) + def tetrahedral_volume(p1, p2, p3, p4): (x1, y1, z1) = p1 (x2, y2, z2) = p2 @@ -165,6 +169,7 @@ def tetrahedral_volume(p1, p2, p3, p4): return volume + def test_correct_receiver_location_generation3D(): """Tests if receiver locations where generated correctly""" @@ -179,6 +184,7 @@ def test_correct_receiver_location_generation3D(): assert np.allclose(receivers.receiver_locations, answer) + def test_correct_receiver_to_cell_location3D(): """Tests if the receivers where located in the correct cell""" @@ -252,6 +258,7 @@ def test_correct_receiver_to_cell_location3D(): assert all([test1, test2, test3]) + def test_correct_at_value3D(): test_model2 = deepcopy(model3D) test_model2["acquisition"]["num_receivers"] = 3 diff --git a/test/test_parallel_source.py b/test/test_parallel_source.py index 34e143c5..e46289d0 100644 --- a/test/test_parallel_source.py +++ b/test/test_parallel_source.py @@ -44,9 +44,9 @@ def test_parallel_source(): receivers, ) - #print(np.amax(np.abs(r))) - #spyro.io.save_shots('serial_shot.dat', r) - r_s = spyro.io.load_shots(os.getcwd()+'/test/serial_shot.dat') + # print(np.amax(np.abs(r))) + # spyro.io.save_shots('serial_shot.dat', r) + r_s = spyro.io.load_shots(os.getcwd() + "/test/serial_shot.dat") assert np.amax(np.abs(r - r_s)) < 1e-16 diff --git a/test/test_readmesh.py b/test/test_readmesh.py index 3d0ee4f8..f34b68c7 100644 --- a/test/test_readmesh.py +++ b/test/test_readmesh.py @@ -1,10 +1,10 @@ import numpy as np -import pytest import spyro """Read in an external mesh and interpolate velocity to it""" + def test_readmesh2(): from .inputfiles.Model1_2d_CG import model @@ -16,6 +16,7 @@ def test_readmesh2(): assert not np.isnan(np.min(vp.dat.data[:])) + def test_readmesh3(): from .inputfiles.Model1_3d_CG import model diff --git a/test/test_sources.py b/test/test_sources.py index c80195fc..dd707b09 100644 --- a/test/test_sources.py +++ b/test/test_sources.py @@ -1,22 +1,19 @@ import math from copy import deepcopy -import pytest -import firedrake as fire -import numpy as np -from firedrake import Constant, dot, dx, grad import spyro """Read in an external mesh and interpolate velocity to it""" from .inputfiles.Model1_2d_CG import model as model + def test_ricker_varies_in_time(): """This test ricker time variation when applied to a time- dependent PDE (acoustic wave second order in pressure) in firedrake. It tests if the right hand side varies in time and if the applied ricker function behaves correctly """ - ### initial ricker tests + # initial ricker tests modelRicker = deepcopy(model) frequency = 2 amplitude = 3 diff --git a/test/test_tools.py b/test/test_tools.py index 50e5b94c..3c79a9f8 100644 --- a/test/test_tools.py +++ b/test/test_tools.py @@ -1,9 +1,8 @@ -from numbers import Real -import pytest import numpy as np import math import spyro + def tetrahedral_volume(p1, p2, p3, p4): (x1, y1, z1) = p1 (x2, y2, z2) = p2 @@ -19,6 +18,7 @@ def tetrahedral_volume(p1, p2, p3, p4): return volume + def triangle_area(p1, p2, p3): """Simple function to calculate triangle area based on its 3 vertices.""" (x1, y1) = p1 @@ -27,75 +27,79 @@ def triangle_area(p1, p2, p3): return abs(x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2 + def test_mesh_generation_for_grid_calc(): grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'homogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'velocity_model_file_name': 'vel_z6.25m_x12.5m_exact.segy', - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'near', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "homogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "velocity_model_file_name": "vel_z6.25m_x12.5m_exact.segy", + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "near", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.07, - 'g_accuracy': 1e-1 - } + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.07, + "g_accuracy": 1e-1, + } Gs = [7.0, 7.1, 7.7, 8.0] - degree_reference = grid_point_calculator_parameters['reference_degree'] + degree_reference = grid_point_calculator_parameters["reference_degree"] - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, degree_reference) + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, degree_reference + ) comm = spyro.utils.mpi_init(model) for G in Gs: - model["mesh"]["meshfile"] = "meshes/2Dhomogeneous"+str(G)+".msh" + model["mesh"]["meshfile"] = "meshes/2Dhomogeneous" + str(G) + ".msh" model = spyro.tools.generate_mesh(model, G, comm) + def test_input_models_receivers(): - test1 = True #testing if 2D receivers are inside the domain on an homogeneous case + test1 = True # testing if 2D receivers are inside the domain on an homogeneous case grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'homogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'near', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "homogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "near", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.05, - 'g_accuracy': 1e-1 - } - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, 4) - - Lz = model["mesh"]['Lz'] - lz = model['BCs']['lz'] - Lx = model["mesh"]['Lx'] - lx = model['BCs']['lx'] + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.05, + "g_accuracy": 1e-1, + } + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, 4 + ) + + Lz = model["mesh"]["Lz"] + lz = model["BCs"]["lz"] + Lx = model["mesh"]["Lx"] + lx = model["BCs"]["lx"] Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx + Real_Lx = Lx + 2 * lx - p1 = (-Real_Lz,-lx) - p2 = (-Real_Lz,Real_Lx-lx) - p3 = (0.0,-lx) - p4 = (0.0,Real_Lx-lx) + p1 = (-Real_Lz, -lx) + p2 = (-Real_Lz, Real_Lx - lx) + p3 = (0.0, -lx) + p4 = (0.0, Real_Lx - lx) - areaSquare = Real_Lz*Real_Lx + areaSquare = Real_Lz * Real_Lx for r in model["acquisition"]["receiver_locations"]: area1 = triangle_area(p1, p2, r) @@ -103,31 +107,32 @@ def test_input_models_receivers(): area3 = triangle_area(p3, p4, r) area4 = triangle_area(p2, p4, r) test = math.isclose((area1 + area2 + area3 + area4), areaSquare, rel_tol=1e-09) - if test == False: + if test is False: test1 = False - test2 =True # For 3D case + test2 = True # For 3D case grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'homogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 3, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'near', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "homogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 3, # Domain dimension. Either 2 or 3. + "receiver_setup": "near", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.05, - 'g_accuracy': 1e-1 - } - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters,4) + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.05, + "g_accuracy": 1e-1, + } + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, 4 + ) # FInish volume test later # Lz = model["mesh"]['Lz'] @@ -163,48 +168,49 @@ def test_input_models_receivers(): # if test == False: # test1 = False - assert all([test1, test2]) + def test_input_models_receivers_heterogeneous(): - test1 = True #testing if 2D receivers bins are inside the domain on an heterogeneous case + test1 = True # testing if 2D receivers bins are inside the domain on an heterogeneous case grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'heterogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'velocity_model_file_name':None, - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'bins', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "heterogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "velocity_model_file_name": None, + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "bins", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.05, - 'g_accuracy': 1e-1 - } - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, 4) - - Lz = model["mesh"]['Lz'] - lz = model['BCs']['lz'] - Lx = model["mesh"]['Lx'] - lx = model['BCs']['lx'] + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.05, + "g_accuracy": 1e-1, + } + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, 4 + ) + + Lz = model["mesh"]["Lz"] + lz = model["BCs"]["lz"] + Lx = model["mesh"]["Lx"] + lx = model["BCs"]["lx"] Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx + Real_Lx = Lx + 2 * lx - p1 = (-Real_Lz,-lx) - p2 = (-Real_Lz,Real_Lx-lx) - p3 = (0.0,-lx) - p4 = (0.0,Real_Lx-lx) + p1 = (-Real_Lz, -lx) + p2 = (-Real_Lz, Real_Lx - lx) + p3 = (0.0, -lx) + p4 = (0.0, Real_Lx - lx) - areaSquare = Real_Lz*Real_Lx + areaSquare = Real_Lz * Real_Lx for r in model["acquisition"]["receiver_locations"]: area1 = triangle_area(p1, p2, r) @@ -212,47 +218,48 @@ def test_input_models_receivers_heterogeneous(): area3 = triangle_area(p3, p4, r) area4 = triangle_area(p2, p4, r) test = math.isclose((area1 + area2 + area3 + area4), areaSquare, rel_tol=1e-09) - if test == False: + if test is False: test1 = False - - test2 = True #testing if 2D receivers line are inside the domain on an heterogeneous case + + test2 = True # testing if 2D receivers line are inside the domain on an heterogeneous case grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'heterogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'velocity_model_file_name':None, - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'line', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "heterogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "velocity_model_file_name": None, + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "line", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.05, - 'g_accuracy': 1e-1 - } - model = spyro.tools.create_model_for_grid_point_calculation(grid_point_calculator_parameters, 4) - - Lz = model["mesh"]['Lz'] - lz = model['BCs']['lz'] - Lx = model["mesh"]['Lx'] - lx = model['BCs']['lx'] + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.05, + "g_accuracy": 1e-1, + } + model = spyro.tools.create_model_for_grid_point_calculation( + grid_point_calculator_parameters, 4 + ) + + Lz = model["mesh"]["Lz"] + lz = model["BCs"]["lz"] + Lx = model["mesh"]["Lx"] + lx = model["BCs"]["lx"] Real_Lz = Lz + lz - Real_Lx = Lx + 2*lx + Real_Lx = Lx + 2 * lx - p1 = (-Real_Lz,-lx) - p2 = (-Real_Lz,Real_Lx-lx) - p3 = (0.0,-lx) - p4 = (0.0,Real_Lx-lx) + p1 = (-Real_Lz, -lx) + p2 = (-Real_Lz, Real_Lx - lx) + p3 = (0.0, -lx) + p4 = (0.0, Real_Lx - lx) - areaSquare = Real_Lz*Real_Lx + areaSquare = Real_Lz * Real_Lx for r in model["acquisition"]["receiver_locations"]: area1 = triangle_area(p1, p2, r) @@ -260,39 +267,36 @@ def test_input_models_receivers_heterogeneous(): area3 = triangle_area(p3, p4, r) area4 = triangle_area(p2, p4, r) test = math.isclose((area1 + area2 + area3 + area4), areaSquare, rel_tol=1e-09) - if test == False: + if test is False: test2 = False - - assert all([test1, test2]) + def test_grid_calc2d(): grid_point_calculator_parameters = { - ## Experiment parameters - 'source_frequency' : 5.0, # Here we define the frequency of the Ricker wavelet source - 'minimum_velocity_in_the_domain' : 1.429, # The minimum velocity present in the domain. + # Experiment parameters + "source_frequency": 5.0, # Here we define the frequency of the Ricker wavelet source + "minimum_velocity_in_the_domain": 1.429, # The minimum velocity present in the domain. # if an homogeneous test case is used this velocity will be defined in the whole domain. - 'velocity_profile_type': 'homogeneous', # Either or heterogeneous. If heterogeneous is - #chosen be careful to have the desired velocity model below. - 'velocity_model_file_name': 'vel_z6.25m_x12.5m_exact.segy', - 'FEM_method_to_evaluate' : 'KMV', # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) - 'dimension' : 2, # Domain dimension. Either 2 or 3. - 'receiver_setup' : 'near', #Either near or line. Near defines a receiver grid near to the source, + "velocity_profile_type": "homogeneous", # Either or heterogeneous. If heterogeneous is + # chosen be careful to have the desired velocity model below. + "velocity_model_file_name": "vel_z6.25m_x12.5m_exact.segy", + "FEM_method_to_evaluate": "KMV", # FEM to evaluate such as `KMV` or `spectral` (GLL nodes on quads and hexas) + "dimension": 2, # Domain dimension. Either 2 or 3. + "receiver_setup": "near", # Either near or line. Near defines a receiver grid near to the source, # line defines a line of point receivers with pre-established near and far offsets. - - ## Line search parameters - 'reference_degree': 4, # Degree to use in the reference case (int) - 'G_reference': 8.0, # grid point density to use in the reference case (float) - 'desired_degree': 4, # degree we are calculating G for. (int) - 'G_initial': 7.0, # Initial G for line search (float) - 'accepted_error_threshold': 0.07, - 'g_accuracy': 1e-1 - } - + # Line search parameters + "reference_degree": 4, # Degree to use in the reference case (int) + "G_reference": 8.0, # grid point density to use in the reference case (float) + "desired_degree": 4, # degree we are calculating G for. (int) + "G_initial": 7.0, # Initial G for line search (float) + "accepted_error_threshold": 0.07, + "g_accuracy": 1e-1, + } G = spyro.tools.minimum_grid_point_calculator(grid_point_calculator_parameters) - inside = (6.9< G and G<8.0) + inside = 6.9 < G and G < 8.0 print(G) assert inside diff --git a/test_3d/test_forward_3d.py b/test_3d/test_forward_3d.py index 38cd63a0..f2219eff 100644 --- a/test_3d/test_forward_3d.py +++ b/test_3d/test_forward_3d.py @@ -1,30 +1,32 @@ from firedrake import File -import matplotlib.pyplot as plt import numpy as np -import math import spyro -import pytest -def compare_velocity(p_r, receiver_in_source_index, receiver_comparison_index, model,dt): - receiver_0 = p_r[:,receiver_in_source_index] - receiver_1 = p_r[:,receiver_comparison_index] + +def compare_velocity( + p_r, receiver_in_source_index, receiver_comparison_index, model, dt +): + receiver_0 = p_r[:, receiver_in_source_index] + receiver_1 = p_r[:, receiver_comparison_index] pos = model["acquisition"]["receiver_locations"] - time0 = np.argmax(receiver_0)*dt - time1 = np.argmax(receiver_1)*dt + time0 = np.argmax(receiver_0) * dt + time1 = np.argmax(receiver_1) * dt x0 = pos[receiver_in_source_index][1] x1 = pos[receiver_comparison_index][1] - measured_velocity = np.abs(x1-x0)/(time1-time0) + measured_velocity = np.abs(x1 - x0) / (time1 - time0) minimum_velocity = 1.5 - error_percent = 100*np.abs(measured_velocity-minimum_velocity)/minimum_velocity + error_percent = ( + 100 * np.abs(measured_velocity - minimum_velocity) / minimum_velocity + ) return error_percent -def test_forward_3d(tf = 0.6): +def test_forward_3d(tf=0.6): model = {} model["opts"] = { @@ -44,7 +46,7 @@ def test_forward_3d(tf = 0.6): } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 6.0, # maximum acoustic wave velocity in PML - km/s @@ -58,13 +60,20 @@ def test_forward_3d(tf = 0.6): "source_pos": [(-0.15, 0.25, 0.25)], "frequency": 5.0, "delay": 1.0, - "receiver_locations": [(-0.15, 0.25, 0.25), (-0.15, 0.3, 0.25), (-0.15, 0.35, 0.25), (-0.15, 0.4, 0.25), (-0.15, 0.45, 0.25), (-0.15, 0.5, 0.25), (-0.15, 0.55, 0.25), (-0.15, 0.6, 0.25)], - } - model["aut_dif"] ={ - "status": False + "receiver_locations": [ + (-0.15, 0.25, 0.25), + (-0.15, 0.3, 0.25), + (-0.15, 0.35, 0.25), + (-0.15, 0.4, 0.25), + (-0.15, 0.45, 0.25), + (-0.15, 0.5, 0.25), + (-0.15, 0.55, 0.25), + (-0.15, 0.6, 0.25), + ], } + model["aut_dif"] = {"status": False} model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": tf, # Final time for event "dt": 0.00075, "amplitude": 1, # the Ricker has an amplitude of 1. @@ -91,16 +100,15 @@ def test_forward_3d(tf = 0.6): model, mesh, comm, vp, sources, wavelet, receivers, output=False ) - dt=model["timeaxis"]["dt"] - final_time=model["timeaxis"]["tf"] + dt = model["timeaxis"]["dt"] pass_error_test = True if comm.comm.rank == 0: - error_percent = compare_velocity(p_r, 0, 7, model,dt) + error_percent = compare_velocity(p_r, 0, 7, model, dt) if error_percent < 5: pass_error_test = True else: pass_error_test = False - - assert pass_error_test \ No newline at end of file + + assert pass_error_test diff --git a/test_ad/test_gradient_3d_AD.py b/test_ad/test_gradient_3d_AD.py index 968590fc..05a6f58d 100644 --- a/test_ad/test_gradient_3d_AD.py +++ b/test_ad/test_gradient_3d_AD.py @@ -1,27 +1,21 @@ from firedrake import * from scipy.optimize import * import spyro -import time -import sys -import matplotlib.pyplot as plt -import numpy as np -import meshio -import SeismicMesh -import finat import pytest -#from ..domains import quadrature, space + +# from ..domains import quadrature, space @pytest.mark.skip(reason="no way of currently testing this with cicleCI resources") def test_gradient_3d_AD(): model = {} model["opts"] = { "method": "KMV", # either CG or KMV - "quadrature": "KMV", # Equi or KMV + "quadrature": "KMV", # Equi or KMV "degree": 2, # p order "dimension": 3, # dimension "regularization": False, # regularization is on? - "gamma": 1e-5, # regularization parameter + "gamma": 1e-5, # regularization parameter } model["parallelism"] = { @@ -41,7 +35,7 @@ def test_gradient_3d_AD(): # Specify a 250-m Absorbing Boundary Layer (ABL) on the three sides of the domain to damp outgoing waves. model["BCs"] = { "status": False, # True or False, used to turn on any type of BC - "outer_bc": "non-reflective", # none or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # none or non-reflective (outer boundary condition) "abl_bc": "none", # none, gaussian-taper, or alid "lz": 0.25, # thickness of the ABL in the z-direction (km) - always positive "lx": 0.25, # thickness of the ABL in the x-direction (km) - always positive @@ -61,39 +55,32 @@ def test_gradient_3d_AD(): "receiver_locations": [(-0.1, 0.1, 0.5), (-0.1, 0.9, 0.5)], } model["aut_dif"] = { - "status": True, + "status": True, } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.8, # Final time for event (for test 7) "dt": 0.001, # timestep size (divided by 2 in the test 4. dt for test 3 is 0.00050) "amplitude": 1, # the Ricker has an amplitude of 1. - "nspool": 20, # (20 for dt=0.00050) how frequently to output solution to pvds + "nspool": 20, # (20 for dt=0.00050) how frequently to output solution to pvds "fspool": 1, # how frequently to save solution to RAM } - comm = spyro.utils.mpi_init(model) + comm = spyro.utils.mpi_init(model) mesh, V = spyro.io.read_mesh(model, comm) element = spyro.domains.space.FE_method( mesh, model["opts"]["method"], model["opts"]["degree"] ) - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, element) z, x, y = SpatialCoordinate(mesh) - vp_exact = Function(V).interpolate( 1.0 + 0.0*x) - vp_guess = Function(V).interpolate( 0.8 + 0.0*x) - + vp_exact = Function(V).interpolate(1.0 + 0.0 * x) + vp_guess = Function(V).interpolate(0.8 + 0.0 * x) - spyro.tools.gradient_test_acoustic_ad( - model, - mesh, - V, - comm, - vp_exact, - vp_guess) + spyro.tools.gradient_test_acoustic_ad(model, mesh, V, comm, vp_exact, vp_guess) -# test_gradient_3d_AD() \ No newline at end of file +# test_gradient_3d_AD() diff --git a/test_ad/test_gradient_AD.py b/test_ad/test_gradient_AD.py index d1d40a4b..34f5839a 100644 --- a/test_ad/test_gradient_AD.py +++ b/test_ad/test_gradient_AD.py @@ -1,28 +1,20 @@ from firedrake import * from scipy.optimize import * -import pytest import spyro -import time -import sys -import matplotlib.pyplot as plt -import numpy as np -import meshio -import SeismicMesh -import finat -import pytest -#from ..domains import quadrature, space + +# from ..domains import quadrature, space # @pytest.mark.skip(reason="no way of currently testing this") def test_gradient_AD(): model = {} model["opts"] = { "method": "KMV", # either CG or KMV - "quadrature": "KMV", # Equi or KMV + "quadrature": "KMV", # Equi or KMV "degree": 1, # p order "dimension": 2, # dimension "regularization": False, # regularization is on? - "gamma": 1e-5, # regularization parameter + "gamma": 1e-5, # regularization parameter } model["parallelism"] = { @@ -42,7 +34,7 @@ def test_gradient_AD(): # Specify a 250-m Absorbing Boundary Layer (ABL) on the three sides of the domain to damp outgoing waves. model["BCs"] = { "status": False, # True or False, used to turn on any type of BC - "outer_bc": "non-reflective", # none or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # none or non-reflective (outer boundary condition) "abl_bc": "none", # none, gaussian-taper, or alid "lz": 0.25, # thickness of the ABL in the z-direction (km) - always positive "lx": 0.25, # thickness of the ABL in the x-direction (km) - always positive @@ -54,41 +46,32 @@ def test_gradient_AD(): "source_pos": [(0.75, 0.75)], "frequency": 10.0, "delay": 1.0, - "receiver_locations": spyro.create_transect( - (0.9, 0.2), (0.9, 0.8), 10 - ), + "receiver_locations": spyro.create_transect((0.9, 0.2), (0.9, 0.8), 10), } model["aut_dif"] = { - "status": True, + "status": True, } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 0.8, # Final time for event (for test 7) "dt": 0.001, # timestep size (divided by 2 in the test 4. dt for test 3 is 0.00050) "amplitude": 1, # the Ricker has an amplitude of 1. - "nspool": 20, # (20 for dt=0.00050) how frequently to output solution to pvds + "nspool": 20, # (20 for dt=0.00050) how frequently to output solution to pvds "fspool": 1, # how frequently to save solution to RAM } comm = spyro.utils.mpi_init(model) - mesh = RectangleMesh(100, 100, 1.5, 1.5) # to test FWI, mesh aligned with interface + mesh = RectangleMesh(100, 100, 1.5, 1.5) # to test FWI, mesh aligned with interface element = spyro.domains.space.FE_method( mesh, model["opts"]["method"], model["opts"]["degree"] ) - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, element) z, x = SpatialCoordinate(mesh) - vp_exact = Function(V).interpolate( 1.0 + 0.0*x) - vp_guess = Function(V).interpolate( 0.8 + 0.0*x) - - spyro.tools.gradient_test_acoustic_ad( - model, - mesh, - V, - comm, - vp_exact, - vp_guess) + vp_exact = Function(V).interpolate(1.0 + 0.0 * x) + vp_guess = Function(V).interpolate(0.8 + 0.0 * x) + spyro.tools.gradient_test_acoustic_ad(model, mesh, V, comm, vp_exact, vp_guess) diff --git a/test_parallel/test_forward.py b/test_parallel/test_forward.py index 92bf1a81..dc7aaae9 100644 --- a/test_parallel/test_forward.py +++ b/test_parallel/test_forward.py @@ -3,7 +3,7 @@ import numpy as np import math import spyro -import pytest + def plot_receiver( receiver, @@ -13,15 +13,15 @@ def plot_receiver( show=False, file_format="pdf", ): - """Plot a + """Plot a Returns ------- None """ - receiver_data = receiver[:,receiver_id] + receiver_data = receiver[:, receiver_id] - nt = int(final_time/ dt) # number of timesteps + nt = int(final_time / dt) # number of timesteps times = np.linspace(0.0, final_time, nt) plt.plot(times, receiver_data) @@ -32,38 +32,45 @@ def plot_receiver( plt.yticks(fontsize=18) # plt.xlim(start_index, end_index) # plt.ylim(tf, 0) - plt.savefig("receiver"+str(receiver_id) + "." + file_format, format=file_format) + plt.savefig("receiver" + str(receiver_id) + "." + file_format, format=file_format) if show: plt.show() plt.close() return None -def compare_velocity(p_r, receiver_in_source_index, receiver_comparison_index, model,dt): - receiver_0 = p_r[:,receiver_in_source_index] - receiver_1 = p_r[:,receiver_comparison_index] + +def compare_velocity( + p_r, receiver_in_source_index, receiver_comparison_index, model, dt +): + receiver_0 = p_r[:, receiver_in_source_index] + receiver_1 = p_r[:, receiver_comparison_index] pos = model["acquisition"]["receiver_locations"] - time0 = np.argmax(receiver_0)*dt - time1 = np.argmax(receiver_1)*dt - x0 = pos[receiver_in_source_index,1] - x1 = pos[receiver_comparison_index,1] - measured_velocity = np.abs(x1-x0)/(time1-time0) + time0 = np.argmax(receiver_0) * dt + time1 = np.argmax(receiver_1) * dt + x0 = pos[receiver_in_source_index, 1] + x1 = pos[receiver_comparison_index, 1] + measured_velocity = np.abs(x1 - x0) / (time1 - time0) minimum_velocity = 1.5 - error_percent = 100*np.abs(measured_velocity-minimum_velocity)/minimum_velocity - print(f"Velocity error of {error_percent}%.", flush = True) + error_percent = ( + 100 * np.abs(measured_velocity - minimum_velocity) / minimum_velocity + ) + print(f"Velocity error of {error_percent}%.", flush=True) return error_percent + def get_receiver_in_source_location(source_id, model): receiver_locations = model["acquisition"]["receiver_locations"] source_locations = model["acquisition"]["source_pos"] - source_x = source_locations[source_id,1] + source_x = source_locations[source_id, 1] cont = 0 for receiver_location in receiver_locations: if math.isclose(source_x, receiver_location[1]): return cont cont += 1 - return ValueError("Couldn't find a receiver whose location coincides with a source within the standard tolerance.") - + return ValueError( + "Couldn't find a receiver whose location coincides with a source within the standard tolerance." + ) def test_forward_5shots(): @@ -88,7 +95,7 @@ def test_forward_5shots(): } model["BCs"] = { "status": True, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) + "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic "exponent": 2, # damping layer has a exponent variation "cmax": 4.5, # maximum acoustic wave velocity in PML - km/s @@ -102,10 +109,10 @@ def test_forward_5shots(): "source_pos": spyro.create_transect((-0.1, 1.0), (-0.1, 15.0), 5), "frequency": 5.0, "delay": 1.0, - "receiver_locations": spyro.create_transect((-0.1, 1.0), (-0.1, 15.0),13), + "receiver_locations": spyro.create_transect((-0.1, 1.0), (-0.1, 15.0), 13), } model["timeaxis"] = { - "t0": 0.0, # Initial time for event + "t0": 0.0, # Initial time for event "tf": 3.00, # Final time for event "dt": 0.001, "amplitude": 1, # the Ricker has an amplitude of 1. @@ -113,8 +120,7 @@ def test_forward_5shots(): "fspool": 99999, # how frequently to save solution to RAM } - dt=model["timeaxis"]["dt"] - final_time=model["timeaxis"]["tf"] + dt = model["timeaxis"]["dt"] comm = spyro.utils.mpi_init(model) @@ -135,19 +141,25 @@ def test_forward_5shots(): pass_error_test = False for source_id in range(len(model["acquisition"]["source_pos"])): if comm.ensemble_comm.rank == (source_id % comm.ensemble_comm.size): - receiver_in_source_index= get_receiver_in_source_location(source_id, model) - if source_id != len(model["acquisition"]["source_pos"])-1 or source_id == 0: + receiver_in_source_index = get_receiver_in_source_location(source_id, model) + if ( + source_id != len(model["acquisition"]["source_pos"]) - 1 + or source_id == 0 + ): receiver_comparison_index = receiver_in_source_index + 1 else: receiver_comparison_index = receiver_in_source_index - 1 - error_percent = compare_velocity(p_r, receiver_in_source_index, receiver_comparison_index, model,dt) + error_percent = compare_velocity( + p_r, receiver_in_source_index, receiver_comparison_index, model, dt + ) if error_percent < 5: pass_error_test = True - print(f"For source = {source_id}: test = {pass_error_test}", flush = True) + print(f"For source = {source_id}: test = {pass_error_test}", flush=True) spyro.plots.plot_shots(model, comm, p_r, vmin=-1e-3, vmax=1e-3) spyro.io.save_shots(model, comm, p_r) assert pass_error_test + if __name__ == "__main__": test_forward_5shots()