Merge pull request #113 from Olender/newHABC

New habc

create_hdf5_velocity_model.py

@@ -0,0 +1,7 @@
+import firedrake as fire
+import spyro
+
+
+# build temp space
+
+mesh = spyro.RectangleMesh(nx, ny, Lx, Ly, quadrilateral=True)

fig18_habcclass.py

@@ -0,0 +1,79 @@
+import spyro
+import math
+from generate_velocity_model_from_paper import get_paper_velocity
+
+dictionary = {}
+dictionary["options"] = {
+    "cell_type": "T",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+    "variant": 'lumped',  # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method
+    "degree": 4,  # p order
+    "dimension": 2,  # dimension
+}
+
+# Number of cores for the shot. For simplicity, we keep things serial.
+# spyro however supports both spatial parallelism and "shot" parallelism.
+dictionary["parallelism"] = {
+    "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+}
+
+# Define the domain size without the PML. Here we'll assume a 1.00 x 1.00 km
+# domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb
+# outgoing waves on three sides (eg., -z, +-x sides) of the domain.
+dictionary["mesh"] = {
+    "Lz": 2.4,  # depth in km - always positive
+    "Lx": 4.8,  # width in km - always positive
+    "Ly": 0.0,  # thickness in km - always positive
+    "mesh_file": None,
+    "user_mesh": None,
+    "mesh_type": "SeismicMesh",
+}
+
+# Create a source injection operator. Here we use a single source with a
+# Ricker wavelet that has a peak frequency of 5 Hz injected at the center of the mesh.
+# We also specify to record the solution at a microphone near the top of the domain.
+# This transect of receivers is created with the helper function `create_transect`.
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": [(-0.6, 4.8-1.68)],
+    "frequency": 5.0,
+    "delay": 1.5,
+    "receiver_locations": [(-0.6, 4.8-1.68)],
+}
+
+# Simulate for 1.0 seconds.
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": 1.00,  # Final time for event
+    "dt": 0.0005,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 100,  # how frequently to output solution to pvds
+    "gradient_sampling_frequency": 100,  # how frequently to save solution to RAM
+}
+
+dictionary["visualization"] = {
+    "forward_output": True,
+    "forward_output_filename": "results/figeigteen_forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "velocity_model_filename": None,
+    "gradient_output": False,
+    "gradient_filename": None,
+    "debug_output": True,
+}
+
+Wave_obj = spyro.HABC(dictionary=dictionary)
+
+cpw = 6.0
+lba = 1.5 / 5.0
+edge_length = lba / cpw
+Wave_obj.set_mesh(mesh_parameters={"edge_length": edge_length})
+V = Wave_obj.function_space
+mesh = Wave_obj.mesh
+c = get_paper_velocity(mesh, V)
+
+Wave_obj.set_initial_velocity_model(velocity_model_function=c)
+Wave_obj._get_initial_velocity_model()
+
+Wave_obj.c = Wave_obj.initial_velocity_model
+# Wave_obj.get_and_set_maximum_dt(fraction=0.5)
+Wave_obj.no_boundary_forward_solve()
+Wave_obj.set_damping_field()

fig8test.py

@@ -0,0 +1,106 @@
+import spyro
+from spyro.habc import HABC
+import firedrake as fire
+import math
+
+
+def test_eikonal_values_fig8():
+    dictionary = {}
+    dictionary["options"] = {
+        "cell_type": "T",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+        "variant": 'lumped',  # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method
+        "degree": 1,  # p order
+        "dimension": 2,  # dimension
+    }
+
+    # Number of cores for the shot. For simplicity, we keep things serial.
+    # spyro however supports both spatial parallelism and "shot" parallelism.
+    dictionary["parallelism"] = {
+        "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+    }
+
+    # Define the domain size without the PML. Here we'll assume a 1.00 x 1.00 km
+    # domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb
+    # outgoing waves on three sides (eg., -z, +-x sides) of the domain.
+    dictionary["mesh"] = {
+        "Lz": 1.0,  # depth in km - always positive
+        "Lx": 1.0,  # width in km - always positive
+        "Ly": 0.0,  # thickness in km - always positive
+        "mesh_file": None,
+        "user_mesh": None,
+        "mesh_type": "firedrake_mesh",
+    }
+
+    # Create a source injection operator. Here we use a single source with a
+    # Ricker wavelet that has a peak frequency of 5 Hz injected at the center of the mesh.
+    # We also specify to record the solution at a microphone near the top of the domain.
+    # This transect of receivers is created with the helper function `create_transect`.
+    dictionary["acquisition"] = {
+        "source_type": "ricker",
+        "source_locations": [(-0.5, 0.25)],
+        "frequency": 5.0,
+        "delay": 1.5,
+        "receiver_locations": spyro.create_transect(
+            (-0.10, 0.1), (-0.10, 0.9), 20
+        ),
+    }
+
+    # Simulate for 2.0 seconds.
+    dictionary["time_axis"] = {
+        "initial_time": 0.0,  # Initial time for event
+        "final_time": 2.00,  # Final time for event
+        "dt": 0.001,  # timestep size
+        "amplitude": 1,  # the Ricker has an amplitude of 1.
+        "output_frequency": 100,  # how frequently to output solution to pvds
+        "gradient_sampling_frequency": 100,  # how frequently to save solution to RAM
+    }
+
+    dictionary["visualization"] = {
+        "forward_output" : True,
+        "output_filename": "results/forward_output.pvd",
+        "fwi_velocity_model_output": False,
+        "velocity_model_filename": None,
+        "gradient_output": False,
+        "gradient_filename": None,
+    }
+
+    Wave_no_habc = spyro.AcousticWave(dictionary=dictionary)
+
+    Lx = 1
+    Lz = 1
+    user_mesh = fire.RectangleMesh(80, 80, Lz, Lx, diagonal="crossed")
+    user_mesh.coordinates.dat.data[:, 0] *= -1.0
+    z, x = fire.SpatialCoordinate(user_mesh)
+
+    cond = fire.conditional(x < 0.5, 3.0, 1.5)
+
+    Wave_no_habc.set_mesh(user_mesh=user_mesh)
+
+    Wave_no_habc.set_initial_velocity_model(conditional=cond)
+    Wave_no_habc._get_initial_velocity_model()
+
+    Wave_no_habc.c = Wave_no_habc.initial_velocity_model
+
+    habc = HABC(Wave_no_habc, h_min=0.0125)
+
+    eikonal = habc.eikonal
+
+    min_value = eikonal.min_value
+    max_value = eikonal.max_value
+
+    paper_min = 0.085
+    paper_max = 0.56
+
+    test_min = math.isclose(min_value, paper_min, rel_tol=0.1)
+    test_max = math.isclose(max_value, paper_max, rel_tol=0.2)
+    print("min_value: ", min_value)
+    print("paper_min: ", paper_min)
+    print("max_value: ", max_value)
+    print("paper_max: ", paper_max)
+
+    assert all([test_min, test_max])
+
+
+# Verificar valores das distancias como lref e velocidades
+if __name__ == "__main__":
+    test_eikonal_values_fig8()

generate_velocity_model_from_paper.py

@@ -0,0 +1,86 @@
+import firedrake as fire
+from firedrake import And, File
+
+
+def apply_box(mesh, c, x1, y1, x2, y2, value):
+    y1 = 2.4 - y1
+    y2 = 2.4 - y2
+    x1 = 4.8 - x1
+    x2 = 4.8 - x2
+    y, x = fire.SpatialCoordinate(mesh)
+    box = fire.conditional(
+        And(And(x1 >= x, x >= x2), And(y1 >= -y, -y >= y2)), value, c
+    )
+
+    c.interpolate(box)
+    return c
+
+
+def apply_slope(mesh, c, x1, y1, x3, y3, value):
+    y, x = fire.SpatialCoordinate(mesh)
+    y1 = 2.4 - y1
+    y3 = 2.4 - y3
+    x1 = 4.8 - x1
+    x3 = 4.8 - x3
+    slope = (y3 - y1) / (x3 - x1)
+
+    slope = fire.conditional(
+        And((-y - y1) / (x - x1) <= slope, x < x1), value, c
+    )
+    c.interpolate(slope)
+    return c
+
+
+def apply_vs_from_list(velmat, mesh, Lx, Ly, c):
+    # (x1, y1, x2, y2, cm)
+    for box in velmat:
+        x1 = box[0] * Lx
+        y1 = box[1] * Ly
+        x2 = box[2] * Lx
+        y2 = box[3] * Ly
+        cm = box[4]
+        c = apply_box(mesh, c, x1, y1, x2, y2, cm)
+
+    return c
+
+
+def get_paper_velocity(mesh, V, output=True, units='km/s'):
+    if units == 'km/s':
+        multiplier = 1.0
+    elif units == 'm/s':
+        multiplier = 1000.0
+    velmat = []
+    velmat.append([0.00, 0.00, 0.35, 0.10, 2.9 * multiplier])
+    velmat.append([0.00, 0.10, 0.25, 0.30, 2.9 * multiplier])
+    velmat.append([0.00, 0.30, 0.25, 0.70, 2.0 * multiplier])
+    velmat.append([0.00, 0.70, 0.10, 1.00, 3.7 * multiplier])
+    velmat.append([0.10, 0.70, 0.30, 0.90, 3.7 * multiplier])
+    velmat.append([0.25, 0.10, 0.75, 0.30, 2.5 * multiplier])
+    velmat.append([0.25, 0.30, 0.40, 0.70, 2.5 * multiplier])
+    velmat.append([0.35, 0.00, 0.70, 0.10, 2.1 * multiplier])
+    velmat.append([0.70, 0.00, 0.90, 0.10, 3.4 * multiplier])
+    velmat.append([0.80, 0.10, 0.90, 0.35, 3.4 * multiplier])
+    velmat.append([0.90, 0.00, 1.00, 0.20, 3.4 * multiplier])
+    velmat.append([0.90, 0.20, 1.00, 0.65, 2.6 * multiplier])
+    velmat.append([0.75, 0.10, 0.80, 0.50, 4.0 * multiplier])
+    velmat.append([0.80, 0.35, 0.90, 0.80, 4.0 * multiplier])
+    velmat.append([0.85, 0.80, 0.90, 0.95, 3.6 * multiplier])
+    velmat.append([0.90, 0.65, 1.00, 1.00, 3.6 * multiplier])
+    velmat.append([0.00, 0.00, 0.00, 0.00, 1.5 * multiplier])
+
+    Lx = 4.8
+    Ly = 2.4
+
+    c = fire.Function(V)
+
+    c.dat.data[:] = 1.5 * multiplier
+
+    c = apply_slope(
+        mesh, c, 0.4 * Lx, 0.3 * Ly, 0.75 * Lx, 0.65 * Ly, 3.3 * multiplier
+    )
+    c = apply_vs_from_list(velmat, mesh, Lx, Ly, c)
+
+    if output is True:
+        File("testing.pvd").write(c)
+
+    return c

habc_propagation.py

@@ -0,0 +1,63 @@
+from spyro.solvers import HABC_wave
+
+
+dictionary = {}
+dictionary["options"] = {
+    "cell_type": "T",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+    "variant": 'lumped',  # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method
+    "degree": 1,  # p order
+    "dimension": 2,  # dimension
+}
+
+# Number of cores for the shot. For simplicity, we keep things serial.
+# spyro however supports both spatial parallelism and "shot" parallelism.
+dictionary["parallelism"] = {
+    "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+}
+
+# Define the domain size without the PML. Here we'll assume a 1.00 x 1.00 km
+# domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb
+# outgoing waves on three sides (eg., -z, +-x sides) of the domain.
+dictionary["mesh"] = {
+    "Lz": 1.0,  # depth in km - always positive
+    "Lx": 1.0,  # width in km - always positive
+    "Ly": 0.0,  # thickness in km - always positive
+    "mesh_type": "firedrake_mesh",  # options: firedrake_mesh or user_mesh
+    "mesh_file": None,  # specify the mesh file
+}
+
+# Create a source injection operator. Here we use a single source with a
+# Ricker wavelet that has a peak frequency of 5 Hz injected at the center of the mesh.
+# We also specify to record the solution at a microphone near the top of the domain.
+# This transect of receivers is created with the helper function `create_transect`.
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": [(-1.0, 1.0)],
+    "frequency": 5.0,
+    "delay": 1.5,
+    "receiver_locations": [(-0.0, 0.5)],
+}
+
+# Simulate for 2.0 seconds.
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": 1.0,  # Final time for event
+    "dt": 0.0005,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 100,  # how frequently to output solution to pvds
+    "gradient_sampling_frequency": 1,  # how frequently to save solution to RAM
+}
+
+dictionary["visualization"] = {
+    "forward_output": False,
+    "output_filename": "results/forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "velocity_model_filename": None,
+    "gradient_output": False,
+    "gradient_filename": None,
+}
+
+
+Wave_obj = HABC_wave(dictionary=dictionary)
+Wave_obj.set_mesh(dx=0.02)
+Wave_obj.forward_solve()