refactor wave parameters

spyro/habc/eik.py

@@ -25,9 +25,9 @@ class HABC_Eikonal(Eikonal_Modeling):
         - (absorb_top, absorb_bottom, absorb_right, absorb_left) for 2D
         - (absorb_top, absorb_bottom, absorb_right,
             absorb_left, absorb_front, absorb_back) for 3D
-    c : `firedrake function`
+    velocity_model : `firedrake function`
         Velocity model without absorbing layer
-    c_min : `float`
+    minimum_velocity : `float`
         Minimum velocity value in the model without absorbing layer
     comm : object
         An object representing the communication interface
@@ -99,10 +99,11 @@ def __init__(self, Wave):
         self.lmin = Wave.mesh_parameters.lmin
 
         # Velocity profile model
-        self.c = Wave.c
+        self.velocity_model = Wave.velocity_model
 
-        # Minimum velocity value in the model
-        self.c_min = Wave.c_min
+        # Computed by HABC_Mesh.preamble_mesh_operations() from the velocity
+        # model on the original domain.
+        self.minimum_velocity = Wave.minimum_velocity
 
         # Absorbing boundaries
         self.boundaries = Wave.get_absorbing_boundaries()
@@ -162,9 +163,12 @@ def solve_eik(self):
         '''
 
         # Eikonal solution
-        self.yp = self.eikonal_solver(self.c, self.c_min,
-                                      self.funct_space_eik,
-                                      self.diam_mesh)
+        self.yp = self.eikonal_solver(
+            self.velocity_model,
+            self.minimum_velocity,
+            self.funct_space_eik,
+            self.diam_mesh,
+        )
 
         # Save Eikonal results
         eikonal_file = fire.VTKFile(self.path_save + "Eik.pvd")
@@ -238,8 +242,17 @@ def ident_crit_eik(self):
             # Identify minimum Eikonal and critical point on the boundary
             eikmin, pnt_crit = self.ident_eik_on_bnd(bnd_ids)
 
-            # Identifying propagation speed at critical point
-            c_bnd = np.float64(self.c.at(pnt_crit).item())
+            # Interpolate the velocity to a VertexOnlyMesh containing the
+            # critical point. point_cloud_field restores input point ordering.
+            critical_velocity = point_cloud_field(
+                self.mesh,
+                np.asarray([pnt_crit]),
+                self.velocity_model,
+                self.node_tol,
+            )
+            c_bnd = np.float64(
+                critical_velocity.dat.data_ro_with_halos[0],
+            )
 
             # Print critical point coordinates
             pnt_str = "at (in km): ({2:3.3f}, {3:3.3f})"

spyro/habc/habc.py

@@ -615,8 +615,10 @@ def velocity_habc(self, inf_model=False, method='point_cloud'):
         V = create_function_space(self.mesh, method_element, 0)
 
         # Initialize velocity field and assigning the original velocity model
-        self.c = fire.Function(V).interpolate(self.initial_velocity_model,
-                                              allow_missing_dofs=True)
+        self.velocity_model = fire.Function(V).interpolate(
+            self.initial_velocity_model,
+            allow_missing_dofs=True,
+        )
 
         # Clipping coordinates to the layer domain
         ufl_coordinates_habc = self.get_spatial_coordinates_habc()
@@ -632,13 +634,15 @@ def velocity_habc(self, inf_model=False, method='point_cloud'):
                                         method=method, name_prop="Velocity")
 
         # Interpolating the velocity model in the layer
-        self.c.interpolate(extended_velocity * layer_mask + (
-            1. - layer_mask) * self.c, allow_missing_dofs=True)
+        self.velocity_model.interpolate(extended_velocity * layer_mask + (
+            1. - layer_mask) * self.velocity_model, allow_missing_dofs=True)
         del layer_mask, lay_field
 
         # Interpolating in the space function of the problem
-        self.c = fire.Function(self.function_space,
-                               name='c [km/s])').interpolate(self.c)
+        self.velocity_model = fire.Function(
+            self.function_space,
+            name='velocity [km/s]',
+        ).interpolate(self.velocity_model)
 
         # Save new velocity model
         if inf_model:
@@ -647,7 +651,7 @@ def velocity_habc(self, inf_model=False, method='point_cloud'):
             file_name = self.case_habc + "/c_habc.pvd"
 
         outfile = fire.VTKFile(self.path_save + file_name)
-        outfile.write(self.c)
+        outfile.write(self.velocity_model)
 
     def fundamental_frequency(self, method=None, monitor=False,
                               fitting_c=(0., 0., 0., 0.)):
@@ -788,7 +792,7 @@ def fundamental_frequency(self, method=None, monitor=False,
                                            quad_rule=self.quadrature_rule,
                                            static_load_for_ceq=q_ref)
 
-            Lsp = mod_sol.solve_eigenproblem(self.c, V=self.function_space,
+            Lsp = mod_sol.solve_eigenproblem(self.velocity_model, V=self.function_space,
                                              quad_rule=self.quadrature_rule,
                                              hyp_par=hyp_par, c_eqref=c_eqref,
                                              fitting_c=fitting_c,
@@ -802,12 +806,12 @@ def fundamental_frequency(self, method=None, monitor=False,
                                                       self.domain_dim,
                                                       self.hyper_axes)
 
-            Lsp = mod_sol.solve_eigenproblem(self.c, V=self.function_space,
+            Lsp = mod_sol.solve_eigenproblem(self.velocity_model, V=self.function_space,
                                              quad_rule=self.quadrature_rule,
                                              coord_norm=coord_norm)
 
         else:
-            Lsp = mod_sol.solve_eigenproblem(self.c,
+            Lsp = mod_sol.solve_eigenproblem(self.velocity_model,
                                              V=self.function_space, shift=1e-8,
                                              quad_rule=self.quadrature_rule)
 
@@ -972,7 +976,7 @@ def check_timestep_habc(self, max_divisor_tf=1, set_max_dt=True,
         dt_sol = eigsol.Modal_Solver(
             self.dimension, method=method, calc_max_dt=True)
         max_dt = dt_sol.estimate_timestep(
-            self.c, self.function_space, self.final_time, shift=1e-8,
+            self.velocity_model, self.function_space, self.final_time, shift=1e-8,
             quad_rule=self.quadrature_rule, fraction=1.)
 
         # Rounding power

spyro/meshing/meshing_habc.py

@@ -1,3 +1,5 @@
+import warnings
+
 import firedrake as fire
 import numpy as np
 import spyro.meshing.meshing_operations as mshops
@@ -6,6 +8,7 @@
     Element3D, FaceDescriptor, Mesh, MeshPoint
 from scipy.spatial import cKDTree
 from spyro.domains.space import create_function_space
+from spyro.io.basicio import parallel_print
 from spyro.meshing.meshing_functions import AutomaticMesh
 from spyro.tools.habc_tools import point_cloud_field
 from spyro.utils.error_management import value_parameter_error
@@ -37,9 +40,9 @@ class HABC_Mesh():
         Minimum velocity value on the boundary of the original domain
     c_bnd_max : `float`
         Maximum velocity value on the boundary of the original domain
-    c_min : `float`
+    minimum_velocity : `float`
         Minimum velocity value in the model without absorbing layer
-    c_max : `float`
+    maximum_velocity : `float`
         Maximum velocity value in the model without absorbing layer
     comm : object
         An object representing the communication interface
@@ -152,13 +155,53 @@ def __init__(self, domain_dim, dimension=2, quadrilateral=False,
 
         # Communicator MPI
         self.comm = comm
+        self.minimum_velocity = None
+        self.maximum_velocity = None
 
         if not hasattr(self, "mesh_ops"):
             self.mesh_ops = mshops.MeshOps(domain_dim, dimension=dimension,
                                            quadrilateral=quadrilateral,
                                            func_space_type=func_space_type,
                                            comm=comm)
 
+    @property
+    def c_min(self):
+        """Deprecated alias for :attr:`minimum_velocity`."""
+        warnings.warn(
+            "HABC_Mesh.c_min is deprecated; use minimum_velocity.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return self.minimum_velocity
+
+    @c_min.setter
+    def c_min(self, value):
+        warnings.warn(
+            "HABC_Mesh.c_min is deprecated; use minimum_velocity.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        self.minimum_velocity = value
+
+    @property
+    def c_max(self):
+        """Deprecated alias for :attr:`maximum_velocity`."""
+        warnings.warn(
+            "HABC_Mesh.c_max is deprecated; use maximum_velocity.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return self.maximum_velocity
+
+    @c_max.setter
+    def c_max(self, value):
+        warnings.warn(
+            "HABC_Mesh.c_max is deprecated; use maximum_velocity.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        self.maximum_velocity = value
+
     def original_boundary_data(self):
         """
         Generate the boundary data from the original domain mesh
@@ -294,30 +337,45 @@ def preamble_mesh_operations(self, ele_type='consistent', f_est=0.03):
         if self.dimension == 3:
             method_element = "DQ" if self.quadrilateral else "DG"
             velocity_space = create_function_space(self.mesh, method_element, 0)
-            self.c = fire.Function(velocity_space, name='c_orig [km/s])')
-            self.c.interpolate(self.initial_velocity_model,
-                               allow_missing_dofs=True)
+            self.velocity_model = fire.Function(
+                velocity_space,
+                name='velocity_orig [km/s]',
+            )
+            self.velocity_model.interpolate(
+                self.initial_velocity_model,
+                allow_missing_dofs=True,
+            )
         else:
-            self.c = fire.Function(self.function_space, name='c_orig [km/s])')
-            self.c.assign(fire.assemble(fire.interpolate(
+            self.velocity_model = fire.Function(
+                self.function_space,
+                name='velocity_orig [km/s]',
+            )
+            self.velocity_model.assign(fire.assemble(fire.interpolate(
                 self.initial_velocity_model, self.function_space)))
 
         # Get extreme values of the velocity model
-        self.c_min = self.initial_velocity_model.dat.data_with_halos.min()
-        self.c_max = self.initial_velocity_model.dat.data_with_halos.max()
+        self.minimum_velocity = (
+            self.initial_velocity_model.dat.data_with_halos.min()
+        )
+        self.maximum_velocity = (
+            self.initial_velocity_model.dat.data_with_halos.max()
+        )
 
         # Print on screen
         cdom_str = "Domain Velocity Range (km/s): {:.3f} - {:.3f}"
-        print(cdom_str.format(self.c_min, self.c_max), flush=True)
+        parallel_print(
+            cdom_str.format(self.minimum_velocity, self.maximum_velocity),
+            comm=self.comm,
+        )
 
         # Save initial velocity model
         vel_c = fire.VTKFile(self.path_save + "preamble/c_vel.pvd")
         if self.dimension == 3:
             c_vis = fire.Function(self.function_space, name='c_orig [km/s])')
-            c_vis.interpolate(self.c)
+            c_vis.interpolate(self.velocity_model)
             vel_c.write(c_vis)
         else:
-            vel_c.write(self.c)
+            vel_c.write(self.velocity_model)
 
         # Generating boundary data from the original domain mesh
         print("Getting Boundary Mesh Data from Original Domain", flush=True)

spyro/solvers/acoustic_solver_construction_no_pml.py

@@ -30,7 +30,7 @@ def construct_solver_or_matrix_no_pml(Wave_object):
 
     # -------------------------------------------------------
     m1 = (
-        (1 / (Wave_object.c * Wave_object.c))
+        (1 / (Wave_object.velocity_model * Wave_object.velocity_model))
         * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2))
         * v
         * dx(**quad_rule)
@@ -45,7 +45,7 @@ def construct_solver_or_matrix_no_pml(Wave_object):
     if Wave_object.abc_active:
         weak_expr_abc = dot((u_n - u_nm1) / Constant(dt), v)
 
-        f_abc = (1 / Wave_object.c) * weak_expr_abc
+        f_abc = (1 / Wave_object.velocity_model) * weak_expr_abc
         qr_s = Wave_object.surface_quadrature_rule
 
         if Wave_object.abc_boundary_layer_type == "hybrid":
@@ -55,7 +55,7 @@ def construct_solver_or_matrix_no_pml(Wave_object):
 
             # Damping
             le += Wave_object.eta_mask * weak_expr_abc * \
-                (1 / (Wave_object.c * Wave_object.c)) * \
+                (1 / (Wave_object.velocity_model * Wave_object.velocity_model)) * \
                 Wave_object.eta_habc * dx(**quad_rule)
 
         else:

spyro/solvers/acoustic_solver_construction_with_pml.py

@@ -21,7 +21,7 @@ def construct_solver_or_matrix_with_pml_2d(Wave_object):
     matrix_free option is on, which it is by default.
     """
     dt = Wave_object.dt
-    c = Wave_object.c
+    velocity_model = Wave_object.velocity_model
 
     V = Wave_object.function_space
     Z = Wave_object.vector_function_space
@@ -67,7 +67,7 @@ def construct_solver_or_matrix_with_pml_2d(Wave_object):
 
     # -------------------------------------------------------
     m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dxlump
-    a = c * c * dot(grad(u_n), grad(v)) * dxlump  # explicit
+    a = velocity_model * velocity_model * dot(grad(u_n), grad(v)) * dxlump  # explicit
 
     # First-order ABC on outer PML boundaries only (not the free surface).
     # Firedrake documents the base facet labels in
@@ -77,7 +77,7 @@ def construct_solver_or_matrix_with_pml_2d(Wave_object):
     # Firedrake's first coordinate, so the labels map here as:
     #   1 = top (z = 0, free surface), 2 = bottom, 3 = left, 4 = right
     qr_s = Wave_object.surface_quadrature_rule
-    abc_expr = c * ((u_n - u_nm1) / dt) * v
+    abc_expr = velocity_model * ((u_n - u_nm1) / dt) * v
     nf = (
         abc_expr * ds(2, **qr_s)
         + abc_expr * ds(3, **qr_s)
@@ -94,7 +94,7 @@ def construct_solver_or_matrix_with_pml_2d(Wave_object):
     # -------------------------------------------------------
     mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dxlump
     mm2 = inner(dot(Gamma_1, pp_n), qq) * dxlump
-    dd = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dxlump
+    dd = velocity_model * velocity_model * inner(grad(u_n), dot(Gamma_2, qq)) * dxlump
     FF += mm1 + mm2 + dd
 
     Wave_object.lhs = fire.lhs(FF)
@@ -117,7 +117,7 @@ def construct_solver_or_matrix_with_pml_3d(Wave_object):
     matrix_free option is on, which it is by default.
     """
     dt = Wave_object.dt
-    c = Wave_object.c
+    velocity_model = Wave_object.velocity_model
 
     V = Wave_object.function_space
     Z = Wave_object.vector_function_space
@@ -170,7 +170,7 @@ def construct_solver_or_matrix_with_pml_3d(Wave_object):
 
     # -------------------------------------------------------
     m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dxlump
-    a = c * c * dot(grad(u_n), grad(v)) * dxlump  # explicit
+    a = velocity_model * velocity_model * dot(grad(u_n), grad(v)) * dxlump  # explicit
 
     # First-order ABC on outer PML boundaries only (not the free surface).
     # Firedrake documents the base facet labels in
@@ -180,7 +180,7 @@ def construct_solver_or_matrix_with_pml_3d(Wave_object):
     # negates Firedrake's first coordinate, so ds(1) is the top free surface
     # (z = 0) and ds(2) is the bottom boundary.
     qr_s = Wave_object.surface_quadrature_rule
-    abc_expr = c * ((u_n - u_nm1) / dt) * v
+    abc_expr = velocity_model * ((u_n - u_nm1) / dt) * v
     nf = (
         abc_expr * ds(2, **qr_s)
         + abc_expr * ds(3, **qr_s)
@@ -199,8 +199,8 @@ def construct_solver_or_matrix_with_pml_3d(Wave_object):
     # -------------------------------------------------------
     mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dxlump
     mm2 = inner(dot(Gamma_1, pp_n), qq) * dxlump
-    dd1 = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dxlump
-    dd2 = -c * c * inner(grad(psi_n), dot(Gamma_3, qq)) * dxlump
+    dd1 = velocity_model * velocity_model * inner(grad(u_n), dot(Gamma_2, qq)) * dxlump
+    dd2 = -velocity_model * velocity_model * inner(grad(psi_n), dot(Gamma_3, qq)) * dxlump
     FF += mm1 + mm2 + dd1 + dd2
 
     mmm1 = (dot((psi - psi_n), phi) / Constant(dt)) * dxlump