diff --git a/spyro/abc/abc_layer.py b/spyro/abc/abc_layer.py index 890a66e1..f43cf40b 100644 --- a/spyro/abc/abc_layer.py +++ b/spyro/abc/abc_layer.py @@ -14,9 +14,10 @@ from ..plots.plots_habc import plot_function_layer_size from ..tools.habc_tools import clipping_coordinates_lay_field, extend_scalar_field_profile from ..utils.error_management import (enum_parameter_error, value_numerical_error, - value_parameter_error) + value_parameter_error, value_string_error) from ..utils.freq_tools import freq_response -from ..utils.typing import HyperLayerDegreeType, LayerShapeType, LayerSizeRefFrequency +from ..utils.typing import (BoundaryConditionsType, HyperLayerDegreeType, + LayerDampingType, LayerShapeType, LayerSizeRefFrequency) # Work from Ruben Andres Salas, Andre Luis Ferreira da Silva, # Luis Fernando Nogueira de Sá, Emilio Carlos Nelli Silva. @@ -36,11 +37,13 @@ class ABCLayer(NRBC): abc_boundary_layer_shape : `typing.LayerShapeType`, optional Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. - abc_boundary_layer_type : `str` - Type of the boundary layer. Options: 'hybrid' or 'PML'. - Default is 'hybrid'. Option 'hybrid' is based on paper of Salas et al. (2022). + abc_boundary_layer_type : `typing.LayerDampingType` + Type of the boundary layer. Options: `LayerDampingType.LOCAL`, + `LayerDampingType.HYBRID`, `LayerDampingType.PML` or `LayerDampingType.NOABCS`. + Default is `LayerDampingType.NOABCS` where no absorbing BCs are applied. + Option `LayerDampingType.HYBRID` is based on paper of Salas et al. (2022). doi: https://doi.org/10.1016/j.apm.2022.09.014 - TODO: Add reference + TODO: Add citation abc_pad_length : `float` Size of the absorbing layer abc_reference_freq : `typing.LayerSizeRefFrequency`, optional @@ -136,7 +139,7 @@ class ABCLayer(NRBC): def __init__(self, domain_dim, frequency, freq_Nyquist, dimension=2, quadrilateral=False, func_space_type=None, abc_boundary_layer_shape=LayerShapeType.RECTANGULAR, - abc_boundary_layer_type="hybrid", + abc_boundary_layer_type=LayerDampingType.HYBRID, abc_reference_freq=LayerSizeRefFrequency.SOURCE, abc_degree_type=HyperLayerDegreeType.REAL, abc_deg_layer=None, output_folder=None, comm=None): @@ -189,38 +192,22 @@ def __init__(self, domain_dim, frequency, freq_Nyquist, dimension=2, # Validate input arguments if not isinstance(domain_dim, tuple): - raise TypeError("domain_dim must be a tuple, " - f"got {type(domain_dim).__name__}.") - - if not isinstance(frequency, (float)): - raise TypeError("frequency must be a float number, " - f"got {type(frequency).__name__}.") - - if dimension not in [2, 3]: - value_parameter_error('dimension', dimension, [2, 3]) - - if abc_boundary_layer_type not in ["hybrid", "PML"]: - value_parameter_error( - 'abc_boundary_layer_type', abc_boundary_layer_type, ["hybrid", "PML"]) - - if abc_deg_layer is not None and abc_deg_layer < 2.: - raise ValueError(f"abc_deg_layer must be >= 2, got {abc_deg_layer}.") - - if output_folder is not None and not isinstance(output_folder, str): - raise TypeError("output_folder must be a string, " - f"got {type(output_folder).__name__}.") + raise TypeError(f"'domain_dim' must be a tuple, got {type(domain_dim).__name__}.") # Original domain dimensions self.domain_dim = domain_dim # Source frequency - self.frequency = frequency + self.frequency = value_numerical_error("frequency", frequency, float_num=True, + integer_num=True, lower_bound=0.) # Nyquist frequency - self.freq_Nyquist = freq_Nyquist + self.freq_Nyquist = value_numerical_error("freq_Nyquist", freq_Nyquist, + float_num=True, integer_num=True, + lower_bound=0.) # Model dimension - self.dimension = dimension + self.dimension = value_parameter_error("dimension", dimension, [2, 3]) # Quadrilateral/hexahedral elements self.quadrilateral = quadrilateral @@ -229,24 +216,30 @@ def __init__(self, domain_dim, frequency, freq_Nyquist, dimension=2, self.func_space_type = func_space_type # ABC layer parameters - self.abc_boundary_layer_type = abc_boundary_layer_type - self.abc_boundary_layer_shape = enum_parameter_error('abc_boundary_layer_shape', + self.abc_boundary_layer_type = enum_parameter_error("abc_boundary_layer_type", + abc_boundary_layer_type, + LayerDampingType) + if abc_boundary_layer_type == LayerDampingType.NOABCS: + value_parameter_error("abc_boundary_layer_type", abc_boundary_layer_type, + [LayerDampingType.HYBRID, LayerDampingType.PML]) + + self.abc_boundary_layer_shape = enum_parameter_error("abc_boundary_layer_shape", abc_boundary_layer_shape, LayerShapeType) - self.abc_reference_freq = enum_parameter_error('abc_reference_freq', + self.abc_reference_freq = enum_parameter_error("abc_reference_freq", abc_reference_freq, LayerSizeRefFrequency) - self.abc_degree_type = enum_parameter_error('abc_degree_type', abc_degree_type, + self.abc_degree_type = enum_parameter_error("abc_degree_type", abc_degree_type, HyperLayerDegreeType) # Layer degree if self.abc_boundary_layer_shape == LayerShapeType.RECTANGULAR: self.abc_deg_layer = None elif self.abc_boundary_layer_shape == LayerShapeType.HYPERSHAPE: - self.abc_deg_layer = abc_deg_layer - value_numerical_error( - 'abc_deg_layer', self.abc_deg_layer, float_num=True, - integer_num=True, lower_bound=2., include_lower_bound=True) + self.abc_deg_layer = value_numerical_error('abc_deg_layer', abc_deg_layer, + float_num=True, integer_num=True, + lower_bound=2., + include_lower_bound=True) # Communicator MPI self.comm = comm @@ -257,18 +250,19 @@ def __init__(self, domain_dim, frequency, freq_Nyquist, dimension=2, # Create the path to save data self.path_to_save_abc_layer_case(output_folder=output_folder) + # Initializing the NRBC class + NRBC.__init__(self, self.domain_dim, + self.abc_boundary_layer_shape, + dimension=self.dimension, + output_folder=self.path_case_abc, + comm=self.comm) + # # Initializing the error measure class # HABCError.__init__(self, self.dt, self.freq_Nyquist, # self.receiver_locations, # output_folder=self.path_save, # output_case=self.path_case_abc) - # # Initializing the NRBC class - # NRBC.__init__(self, domain_dim, - # self.abc_boundary_layer_shape, - # dimension=self.dimension, - # output_folder=self.path_case_abc) - def _define_layer_shape(self): """Define the shape of the absorbing layer. @@ -324,14 +318,10 @@ def formatting_abc_layer_type(self, str_to_format, for_prints=True): """ # Layer type - if self.abc_boundary_layer_type == "hybrid": + if self.abc_boundary_layer_type == LayerDampingType.HYBRID: abc_layer_str = "Absorbing" if for_prints else "habc" - elif self.abc_boundary_layer_type == "PML": + elif self.abc_boundary_layer_type == LayerDampingType.PML: abc_layer_str = "PML" if for_prints else "pml" - else: - value_parameter_error('abc_boundary_layer_type', - self.abc_boundary_layer_type, - ["hybrid", "PML"]) formatted_str = str_to_format.format(abc_layer_str) @@ -393,6 +383,9 @@ def path_to_save_abc_layer_case(self, output_folder=None): None """ + # Validate the output folder parameter + value_string_error("output_folder", output_folder) + # Identify the case of the ABC scheme for output labeling self.case_abc = self.identify_abc_layer_case() @@ -414,6 +407,7 @@ def critical_boundary_points(self, Wave): See Salas et al (2022): Hybrid absorbing scheme based on hyperelliptical layers with non-reflecting boundary conditions in scalar wave equations. doi: https://doi.org/10.1016/j.apm.2022.09.014 + TODO: Add citation Parameters ---------- @@ -471,7 +465,7 @@ def det_reference_freq(self, fpad=4): # Get the minimum frequency excited at each critical point freq_ref = freq_response(histPcrit, self.freq_Nyquist, - fpad=fpad, get_max_freq=True) + fpad=fpad, get_dominant_freq=True) pprint(f"Frequency at Critical Point {n_crit:>2.0f}: {freq_ref:.5f}", comm=self.comm) @@ -774,8 +768,7 @@ def velocity_abc(self, Wave, inf_model=False, method="point_cloud", save_file=Tr del layer_mask, lay_field # Interpolating in the space function of the problem - Wave.c = Function(Wave.function_space, - name="c[km/s])").interpolate(Wave.c) + Wave.c = Function(Wave.function_space, name="c[km/s])").interpolate(Wave.c) # Save new velocity model if save_file: @@ -789,36 +782,57 @@ def velocity_abc(self, Wave, inf_model=False, method="point_cloud", save_file=Tr outfile = VTKFile(self.path_save + file_name) outfile.write(Wave.c) - # def nrbc_on_boundary_layer(self, sommerfeld_bc=False): - # """ - # Apply the Higdon ABCs on the outer boundary of the absorbing layer - - # Parameters - # ---------- - # sommerfeld_bc : `bool`, optional - # If `True`, use Sommerfeld BC instead of Higdon BC. Default is `False` - - # Returns - # ------- - # None - # """ - - # pprint("\nApplying Non-Reflecting Boundary Conditions", comm=self.comm) + def nrbc_on_boundary_layer(self, Wave_object, non_reflect_bc, save_file=True): + """Apply Non-Reflective BCs on the outer boundary of the absorbing layer. - # # Getting boundary data from the layer boundaries - # bnd_nfs, bnd_nodes_nfs = self.layer_boundary_data(self.function_space) + Parameters + ---------- + Wave_object : `acoustic_wave.AcousticWave` + An instance of the :class:`~spyro.solvers.acoustic_wave.AcousticWave`. + non_reflect_bc : `typing.BoundaryConditionsType` + Type of boundary condition to apply on the outer absorbing layer boundaries. + - Options for Non-Reflecting BCs: + 'BoundaryConditionsType.HIGDON' or 'BoundaryConditionsType.SOMMERFELD'. + save_file : `bool`, optional + If `True`, save the velocity model with absorbing layer in a .pvd file. + Default is `True`. - # # Hypershape parameters - # layer_shape = self.abc_boundary_layer_shape - # if layer_shape == LayerShapeType.HYPERSHAPE: - # hyp_par = (self.n_hyp, *self.hyper_axes) - # else: - # hyp_par = None + Returns + ------- + None + """ - # # Applying Higdon ABCs - # self.cos_ang_HigdonBC(self.function_space, self.crit_source, bnd_nfs, - # bnd_nodes_nfs, hyp_par=hyp_par, - # sommerfeld_bc=sommerfeld_bc) + # Applying NRBCs on outer boundary layer + crit_source = bnd_nod_ids_nfs = bnd_nodes_nfs = None + if non_reflect_bc == BoundaryConditionsType.SOMMERFELD or \ + non_reflect_bc == BoundaryConditionsType.HIGDON: + + pprint("\nApplying Non-Reflecting Boundary Conditions", comm=self.comm) + + # Getting boundary data from the layer boundaries + if non_reflect_bc == BoundaryConditionsType.SOMMERFELD: + bnd_nod_ids_nfs = \ + Wave_object.mesh_ops.layer_boundary_data(Wave_object.mesh, + Wave_object.function_space, + Wave_object.mesh_parameters)[0] + + if non_reflect_bc == BoundaryConditionsType.HIGDON: + crit_source = self.crit_source + bnd_nod_ids_nfs, bnd_nodes_nfs = \ + Wave_object.mesh_ops.layer_boundary_data(Wave_object.mesh, + Wave_object.function_space, + Wave_object.mesh_parameters) + + # Hypershape parameters + hyp_par = (self.layer_geometry.n_hyp, *self.layer_geometry.hyper_axes) \ + if self.abc_boundary_layer_shape == LayerShapeType.HYPERSHAPE else None + + # Applying Higdon ABCs + self.cos_ang_HigdonBC(Wave_object.function_space, crit_source, + bnd_nod_ids_nfs, bnd_nodes_nfs, non_reflect_bc, + hyp_par=hyp_par, save_file=save_file) + else: + pprint("\nNot Non-Reflecting Boundary Conditions Prescribed", comm=self.comm) # def check_timestep_abc(self, max_divisor_tf=1, set_max_dt=True, # method='ANALYTICAL', mag_add=3): diff --git a/spyro/abc/hyp_lay.py b/spyro/abc/hyp_lay.py index 6721a280..3c1204f0 100644 --- a/spyro/abc/hyp_lay.py +++ b/spyro/abc/hyp_lay.py @@ -3,7 +3,7 @@ from scipy.special import beta, betainc, gamma from sys import float_info from ..io.basicio import parallel_print as pprint -from ..utils.error_management import (enum_parameter_error, value_dimension_error, +from ..utils.error_management import (enum_parameter_error, value_model_dimension_error, value_numerical_error, value_parameter_error) from ..utils.typing import HyperLayerDegreeType @@ -125,8 +125,7 @@ def __init__(self, domain_dim, n_hyp=2., n_type=HyperLayerDegreeType.REAL, if n_hyp < 2.: raise ValueError(f"n_hyp must be >= 2, got {n_hyp}.") - if dimension not in [2, 3]: - value_parameter_error('dimension', dimension, [2, 3]) + value_parameter_error('dimension', dimension, [2, 3]) # Original domain dimensions self.domain_dim = domain_dim @@ -789,11 +788,8 @@ def calc_hyp_geom_prop(self, domain_hyp, pad_len, lmin): pprint("Determining Hypershape Layer Parameters") # Domain dimensions w/o layer - chk_domd = len(self.domain_dim) - chk_habc = len(domain_hyp) - if self.dimension != chk_domd or self.dimension != chk_habc: - value_dimension_error(('domain_dim', 'domain_hyp'), - (chk_domd, chk_habc), self.dimension) + value_model_dimension_error(('domain_dim', 'domain_hyp'), + (self.domain_dim, domain_hyp), self.dimension) # Defining the hypershape semi-axes self.define_hyperaxes(domain_hyp) diff --git a/spyro/abc/nrbc.py b/spyro/abc/nrbc.py index 5ed8b029..ab088909 100644 --- a/spyro/abc/nrbc.py +++ b/spyro/abc/nrbc.py @@ -1,6 +1,12 @@ -# import firedrake as fire -import numpy as np +"""Non-reflecting boundary condition helpers for ABCs.""" + +from firedrake import Function, VTKFile +from numpy import abs, asarray, cos, maximum, pi, sign, sqrt, sum +from numpy.linalg import norm from os import getcwd +from ..io.basicio import parallel_print as pprint +from ..utils.error_management import value_parameter_error +from ..utils.typing import BoundaryConditionsType, LayerShapeType # Work from Ruben Andres Salas, Andre Luis Ferreira da Silva, # Luis Fernando Nogueira de Sá, Emilio Carlos Nelli Silva. @@ -17,6 +23,9 @@ class NRBC(): Attributes ---------- + abc_boundary_layer_shape : `typing.LayerShapeType`, optional + Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or + `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. angle_max : `float` Maximum incidence angle considered. Default is `numpy.pi/4`. cos_Hig : `firedrake function` @@ -29,8 +38,6 @@ class NRBC(): domain_dim : `tuple` Original domain dimensions: (length_z, length_x) for 2D or (length_z, length_x, length_y) for 3D. - layer_shape : `string` - Shape type of pad layer. Options: 'rectangular' or 'hypershape'. nrbc : `str` Type of NRBC used. Either "Higdon" or "Sommerfeld". path_save_nrbc : `str` @@ -43,27 +50,30 @@ class NRBC(): hypershape_normal_vector() Compute the normal vector to a hypershape at a boundary point. source_to_bnd_reference_vector() - Compute a unitary reference vector from the source to a boundary point. + Compute a unitary direction vector from the source to a boundary point. """ - def __init__(self, domain_dim, layer_shape, angle_max=np.pi/4., - dimension=2, output_folder=None): - """ - Initialize the NRBC class. + def __init__(self, domain_dim, abc_boundary_layer_shape, angle_max=pi/4., + dimension=2, output_folder=None, comm=None): + """Initialize the NRBC class. Parameters ---------- domain_dim : `tuple` Original domain dimensions: (length_z, length_x) for 2D or (length_z, length_x, length_y) for 3D. - layer_shape : `string` - Shape type of pad layer. Options: 'rectangular' or 'hypershape'. + abc_boundary_layer_shape : `typing.LayerShapeType`, optional + Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or + `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. angle_max : `float`, optional Maximum incidence angle considered. Default is `numpy.pi/4` (45°). dimension : `int`, optional Model dimension (2D or 3D). Default is 2D. - output_folder : str, optional + output_folder : `str`, optional The folder where output data will be saved. Default is `None`. + comm : `object`, optional + An object representing the communication interface for parallel processing. + Default is `None`. Returns ------- @@ -74,13 +84,13 @@ def __init__(self, domain_dim, layer_shape, angle_max=np.pi/4., self.domain_dim = domain_dim # Shape type of pad layer - self.layer_shape = layer_shape + self.abc_boundary_layer_shape = abc_boundary_layer_shape # Maximum incidence angle considered self.angle_max = angle_max # Maximum value of the cosine of the incidence angle - self.cos_min = np.cos(angle_max) + self.cos_min = cos(angle_max) # Model dimension self.dimension = dimension @@ -91,204 +101,210 @@ def __init__(self, domain_dim, layer_shape, angle_max=np.pi/4., else: self.path_save_nrbc = output_folder - # def source_to_bnd_reference_vector(self, source_coord, bnd_nodes_nfs): - # """ - # Compute a unitary reference vector from the source to a boundary point + # Communicator MPI + self.comm = comm + + def source_to_bnd_reference_vector(self, source_coord, bnd_nodes_nfs): + """Compute a unitary direction vector from the source to a boundary point. + + Parameters + ---------- + source_coord : `tuple` + Source coordinates. + bnd_nodes_nfs : `tuple` + Mesh node coordinates on non-free surfaces. + - (z_data[nfs_idx], x_data[nfs_idx]) for 2D. + - (z_data[nfs_idx], x_data[nfs_idx], y_data[nfs_idx]) for 3D. + + Returns + ------- + unit_ref_vct : `array` + Unit direction vector from the source to a boundary point. + """ + + # Boundary node data + bnd_z, bnd_x = bnd_nodes_nfs[:2] + + # Source coordinates + psouz = source_coord[0] + psoux = source_coord[1] + + # Components of the vector pointing to the boundary point + ref_x = bnd_x - psoux + ref_z = bnd_z - psouz + ref_vct = [ref_x, ref_z] + + if self.dimension == 3: # 3D + + # Third component of the vector pointing to the boundary point + bnd_y = bnd_nodes_nfs[2] + psouy = source_coord[2] + ref_y = bnd_y - psouy + ref_vct.append(ref_y) + + # Unitary vector pointing to the boundary point + unit_ref_vct = asarray(ref_vct) / norm(ref_vct, axis=0) + + return unit_ref_vct + + def hypershape_normal_vector(self, bnd_pnts, hyper_axes, n): + """Compute the normal vector to a boundary point of a hypershape. + + Compute the normal vector to a hyperellipse (|x/a|^n + |y/b|^n = 1) or + a hyperellipsoid (|x/a|^n + |y/b|^n + |z/c|^n = 1) at a boundary point. + The hypershape must have the center at the origin. + + Parameters + ---------- + bnd_pnts : `list` + Boundary hypershape points where the normal vector is computed. + Structure: [x, y] for 2D and [x, y, z] for 3D. + hyper_axes : `list` + Semi-axes of the hyperellipse [a, b] or hyperellipsoid [a, b, c]. + n : `float` + Degree of the hyperellipse. + + Returns + ------- + unit_nrm_vct : `array` + Unitary normal vector to the hypershape at the boundary point. + + Notes + ----- + Let f(x, y) = |x/a|^n - |y/b|^n -1 = 0 a level curve (level set for + two variables) for f(x, y, z) at z = 0. The gradient of the function + f given by ∇f(x,y) = [∂f/∂x, ∂f/∂y] is a normal vector to the curve. + The normal vector is given by the partial derivatives of the function. + """ + + # Point coordinates + x, y = bnd_pnts[:2] + + # Hypershape semi-axes + a, b = hyper_axes[:2] + + # Compute partial derivatives + df_dx = (n / (a**n)) * sign(x) * abs(x)**(n - 1) + df_dy = (n / (b**n)) * sign(y) * abs(y)**(n - 1) + + nrm_vct = [df_dx, df_dy] + + if self.dimension == 3: # 3D + + # Third coordinate + z = bnd_pnts[2] + + # Third hypershape semi-axis + c = hyper_axes[2] + + # Partial derivative with respect to third coordinate + df_dz = (n / (c**n)) * sign(z) * abs(z)**(n - 1) + + nrm_vct.append(df_dz) + + # Unitary hypershape normal vector + unit_nrm_vct = asarray(nrm_vct) / norm(nrm_vct, axis=0) + + return unit_nrm_vct + + def cos_ang_HigdonBC(self, V, source_coord, bnd_nfs, bnd_nodes_nfs, + non_reflect_bc, hyp_par=None, save_file=True): + """Compute the cosine of the incidence angle for first-order Higdon BC. + + Parameters + ---------- + V : `firedrake function space` + Function space where the Non-Reflective BCs are defined. + source_coord : `tuple` + Source coordinates. + bnd_nfs : 'array' + Mesh node indices on non-free surfaces. + bnd_nodes_nfs : `tuple` + Mesh node coordinates on non-free surfaces. + - (z_data[nfs_idx], x_data[nfs_idx]) for 2D. + - (z_data[nfs_idx], x_data[nfs_idx], y_data[nfs_idx]) for 3D. + non_reflect_bc : `typing.BoundaryConditionsType` + Type of boundary condition to apply on the outer absorbing layer boundaries. + - Options for Non-Reflecting BCs: + 'BoundaryConditionsType.HIGDON' or 'BoundaryConditionsType.SOMMERFELD'. + hyp_par : `tuple`, optional + Hyperellipse parameters. Structure: + (n_hyp, a_hyp, b_hyp) for 2D or (n_hyp, a_hyp, b_hyp, b_hyp) for 3D. + - n_hyp : `float` + Degree of the hyperellipse. + - a_hyp : `float` + Hyperellipse semi-axis in direction x. + - b_hyp : `float` + Hyperellipse semi-axis in direction z. + - c_hyp : `float` + Hyperellipse semi-axis in direction y (3D only). + save_file : `bool`, optional + If `True`, save the velocity model with absorbing layer in a .pvd file. + Default is `True`. + + Returns + ------- + None + """ + + # Check if the non-reflective BC type is valid + self.nrbc = value_parameter_error('non_reflect_bc', non_reflect_bc, + [BoundaryConditionsType.HIGDON, + BoundaryConditionsType.SOMMERFELD]) + + pprint(f"Creating Field for NRBC: {non_reflect_bc.value}", comm=self.comm) + + # Initialize field for the cosine of the incidence angle + self.cosHig = Function(V, name='cosHig') + + if self.nrbc == BoundaryConditionsType.SOMMERFELD: # Sommerfeld BC + cos_Hig = 1. + + else: # Higdon BC + + # Unitary reference vector pointing to the boundary point + unit_ref_vct = self.source_to_bnd_reference_vector(source_coord, + bnd_nodes_nfs) + + # Normal vector to the boundary + if self.abc_boundary_layer_shape == LayerShapeType.RECTANGULAR: + # Normal vector to the boundary is a orthonormal vector, then + # cosine on incidence angle can be estimated from a projection + # of the reference vector to boundary onto the orthonormal + # vectors ([1, 0, 0] (2D), [0, 1, 0] (2D), [0, 0, 1] (3D)) + cos_Hig = maximum.reduce(abs(unit_ref_vct)) + + if self.abc_boundary_layer_shape == LayerShapeType.HYPERSHAPE: + + # Original domain dimensions + length_z, length_x = self.domain_dim[:2] + + # Hypershape degree and semi-axes + n_hyp, hyp_axes = hyp_par[0], hyp_par[1:] + + # Boundary points of the hypershape centered at the origin + bnd_z, bnd_x = bnd_nodes_nfs[:2] # Boundary node data + bnd_pnts = [bnd_x - length_x / 2, bnd_z + length_z / 2] + if self.dimension == 3: # 3D + length_y = self.domain_dim[2] + bnd_y = bnd_nodes_nfs[2] + bnd_pnts.append(bnd_y - length_y / 2) + + # Normal vector to the boundary + unit_nrm_vct = self.hypershape_normal_vector(bnd_pnts, hyp_axes, n_hyp) + + # Cosine of the incidence angle + cos_Hig = sum(unit_ref_vct * unit_nrm_vct, axis=0) + + # Adjust values to minimum cosine of incidence angle + cos_Hig[cos_Hig < self.cos_min] = sqrt(1. - cos_Hig[cos_Hig < self.cos_min]**2) + + self.cosHig.dat.data_with_halos[bnd_nfs] = cos_Hig - # Parameters - # ---------- - # source_coord : `tuple` - # Source coordinates - # bnd_nodes_nfs : `tuple` - # Mesh node coordinates on non-free surfaces. - # - (z_data[nfs_idx], x_data[nfs_idx]) for 2D - # - (z_data[nfs_idx], x_data[nfs_idx], y_data[nfs_idx]) for 3D - - # Returns - # ------- - # unit_ref_vct : `array` - # Unit reference vector from the source to a boundary point - # """ - - # # Boundary node data - # bnd_z, bnd_x = bnd_nodes_nfs[:2] - - # # Source coordinates - # psouz = source_coord[0] - # psoux = source_coord[1] - - # # Components of the vector pointing to the boundary point - # ref_x = bnd_x - psoux - # ref_z = bnd_z - psouz - # ref_vct = [ref_x, ref_z] - - # if self.dimension == 3: # 3D - - # # Third component of the vector pointing to the boundary point - # bnd_y = bnd_nodes_nfs[2] - # psouy = source_coord[2] - # ref_y = bnd_y - psouy - # ref_vct.append(ref_y) - - # # Unitary vector pointing to the boundary point - # unit_ref_vct = np.asarray(ref_vct) / np.linalg.norm(ref_vct, axis=0) - - # return unit_ref_vct - - # def hypershape_normal_vector(self, bnd_pnts, hyper_axes, n): - # """ - # Compute the normal vector to a hyperellipse (|x/a|^n + |y/b|^n = 1) - # or a hyperellipsoid (|x/a|^n + |y/b|^n + |z/c|^n = 1) at a boundary - # point. The hypershape must have the center at the origin. - - # Observations: - # Let f(x, y) = |x/a|^n - |y/b|^n -1 = 0 a level curve (level set for - # two variables) for f(x, y, z) at z = 0. The gradient of the function - # f given by ∇f(x,y) = [∂f/∂x, ∂f/∂y] is a normal vector to the curve. - # The normal vector is given by the partial derivatives of the function - - # Parameters - # ---------- - # bnd_pnts : `list` - # Boundary hypershape points where the normal vector is computed. - # Structure: [x, y] for 2D and [x, y, z] for 3D - # hyper_axes : `list` - # Semi-axes of the hyperellipse [a, b] or hyperellipsoid [a, b, c] - # n : `float` - # Degree of the hyperellipse - - # Returns - # ------- - # unit_nrm_vct : `array` - # Unitary normal vector to the hypershape at the boundary point - # """ - - # # Point coordinates - # x, y = bnd_pnts[:2] - - # # Hypershape semi-axes - # a, b = hyper_axes[:2] - - # # Compute partial derivatives - # df_dx = (n / (a**n)) * np.sign(x) * abs(x)**(n - 1) - # df_dy = (n / (b**n)) * np.sign(y) * abs(y)**(n - 1) - - # nrm_vct = [df_dx, df_dy] - - # if self.dimension == 3: # 3D - - # # Third coordinate - # z = bnd_pnts[2] - - # # Third hypershape semi-axis - # c = hyper_axes[2] - - # # Partial derivative with respect to third coordinate - # df_dz = (n / (c**n)) * np.sign(z) * abs(z)**(n - 1) - - # nrm_vct.append(df_dz) - - # # Unitary hypershape normal vector - # unit_nrm_vct = np.asarray(nrm_vct) / np.linalg.norm(nrm_vct, axis=0) - - # return unit_nrm_vct - - # def cos_ang_HigdonBC(self, V, source_coord, bnd_nfs, bnd_nodes_nfs, - # hyp_par=None, sommerfeld_bc=False): - # """ - # Compute the cosine of the incidence angle for first-order Higdon BC. - - # Parameters - # ---------- - # V : `firedrake function space` - # Function space where the cosine of the incidence angle is defined - # source_coord : `tuple` - # Source coordinates - # bnd_nfs : 'array' - # Mesh node indices on non-free surfaces - # bnd_nodes_nfs : `tuple` - # Mesh node coordinates on non-free surfaces. - # - (z_data[nfs_idx], x_data[nfs_idx]) for 2D - # - (z_data[nfs_idx], x_data[nfs_idx], y_data[nfs_idx]) for 3D - # hyp_par : `tuple`, optional - # Hyperellipse parameters. Structure: - # (n_hyp, a_hyp, b_hyp) for 2D or (n_hyp, a_hyp, b_hyp, b_hyp) for 3D - # - n_hyp : `float` - # Degree of the hyperellipse - # - a_hyp : `float` - # Hyperellipse semi-axis in direction x - # - b_hyp : `float` - # Hyperellipse semi-axis in direction z - # - c_hyp : `float` - # Hyperellipse semi-axis in direction y (3D only) - # sommerfeld_bc : `bool`, optional - # If True, use Sommerfeld BC instead of Higdon BC. Default is False - - # Returns - # ------- - # None - # """ - - # nrbc_str = "Sommerfeld" if sommerfeld_bc else "Higdon" - # print("Creating Field for NRBC:", nrbc_str, flush=True) - - # # Initialize field for the cosine of the incidence angle - # self.cosHig = fire.Function(V, name='cosHig') - - # if sommerfeld_bc: # Sommerfeld BC - # self.nrbc = "Sommerfeld" - # cos_Hig = 1. - - # else: # Higdon BC - - # self.nrbc = "Higdon" - - # # Unitary reference vector pointing to the boundary point - # unit_ref_vct = self.source_to_bnd_reference_vector(source_coord, - # bnd_nodes_nfs) - - # # Normal vector to the boundary - # if self.layer_shape == 'rectangular': - # # Normal vector to the boundary is a orthonormal vector, then - # # cosine on incidence angle can be estimated from a projection - # # of the reference vector to boundary onto the orthonormal - # # vectors ([1, 0, 0] (2D), [0, 1, 0] (2D), [0, 0, 1] (3D)) - # cos_Hig = np.maximum.reduce(abs(unit_ref_vct)) - - # if self.layer_shape == 'hypershape': - - # # Original domain dimensions - # Lx, Lz = self.domain_dim[:2] - - # # Hypershape degree and semi-axes - # n_hyp, hyp_axes = hyp_par[0], hyp_par[1:] - - # # Boundary points of the hypershape centered at the origin - # bnd_z, bnd_x = bnd_nodes_nfs[:2] # Boundary node data - # bnd_pnts = [bnd_x - Lx / 2, bnd_z + Lz / 2] - # if self.dimension == 3: # 3D - # Ly = self.domain_dim[2] - # bnd_y = bnd_nodes_nfs[2] - # bnd_pnts.append(bnd_y - Ly / 2) - - # # Normal vector to the boundary - # unit_nrm_vct = self.hypershape_normal_vector(bnd_pnts, - # hyp_axes, - # n_hyp) - - # # Cosine of the incidence angle - # cos_Hig = np.sum(unit_ref_vct * unit_nrm_vct, axis=0) - - # # Adjust values to minimum cosine of incidence angle - # cos_Hig[cos_Hig < self.cos_min] = ( - # 1. - cos_Hig[cos_Hig < self.cos_min]**2)**0.5 - - # self.cosHig.dat.data_with_halos[bnd_nfs] = cos_Hig - - # # Save boundary profile of cosine of incidence angle - # if hasattr(self, 'path_save_nrbc'): - # outfile = fire.VTKFile(self.path_save_nrbc + "cosHig.pvd") - # outfile.write(self.cosHig) + # Save boundary profile of cosine of incidence angle + if save_file: + outfile = VTKFile(self.path_save_nrbc + "cosHig.pvd") + outfile.write(self.cosHig) # dx = 0.1 km REC # W/O = 107.72% - 0.80% diff --git a/spyro/abc/rec_lay.py b/spyro/abc/rec_lay.py index 894b7d4a..b9d0b028 100644 --- a/spyro/abc/rec_lay.py +++ b/spyro/abc/rec_lay.py @@ -1,5 +1,5 @@ from ..io.basicio import parallel_print as pprint -from ..utils.error_management import (value_dimension_error, value_numerical_error, +from ..utils.error_management import (value_model_dimension_error, value_numerical_error, value_parameter_error) # Work from Ruben Andres Salas, Andre Luis Ferreira da Silva, @@ -76,8 +76,7 @@ def __init__(self, domain_dim, dimension=2, comm=None): raise TypeError("'domain_dim' must be a tuple, " f"got {type(domain_dim).__name__}.") - if dimension not in [2, 3]: - value_parameter_error('dimension', dimension, [2, 3]) + value_parameter_error('dimension', dimension, [2, 3]) # Original domain dimensions self.domain_dim = domain_dim @@ -139,12 +138,8 @@ def calc_rec_geom_prop(self, domain_layer, pad_len): pprint("Determining Rectangular Layer Parameters", comm=self.comm) # Checking inputs - chk_domain = len(self.domain_dim) - chk_layer = len(domain_layer) - if self.dimension != chk_domain or self.dimension != chk_layer: - value_dimension_error(('domain_dim', 'domain_layer'), - (chk_domain, chk_layer), - self.dimension) + value_model_dimension_error(('domain_dim', 'domain_layer'), + (self.domain_dim, domain_layer), self.dimension) # Domain dimensions w/o layer length_z, length_x = self.domain_dim[:2] diff --git a/spyro/domains/space.py b/spyro/domains/space.py index fa30250c..300dce1b 100644 --- a/spyro/domains/space.py +++ b/spyro/domains/space.py @@ -37,23 +37,23 @@ def create_function_space(mesh, method, degree=None, dim=1): Parameters: ----------- - mesh: Firedrake Mesh + mesh: `Firedrake.Mesh` Mesh to be used in the finite element space. - method: str or FiniteElement + method: `str` or `FiniteElement` Method to be used for the finite element space, or an already constructed finite element. - degree: int or None + degree: `int` or `None` Degree of the finite element space. Required when ``method`` is a supported method name. Ignored only when ``method`` is an already constructed finite element. - dim: int + dim: `int` Number of vector components. If ``dim`` is 1, a scalar function space is created. If ``dim`` is greater than 1, the selected element is wrapped in a vector element with ``dim`` components. Returns: -------- - function_space: Firedrake FunctionSpace + function_space: `Firedrake.FunctionSpace` Function space. """ diff --git a/spyro/habc/damp_profile.py b/spyro/habc/damp_profile.py index de96f239..e8e1ac12 100644 --- a/spyro/habc/damp_profile.py +++ b/spyro/habc/damp_profile.py @@ -217,9 +217,7 @@ def min_reflection(self, psi_damp=None, CR_err=None, typ_CR='CR_PSI'): Heuristic factor for the minimum damping ratio ''' - if typ_CR not in ['CR_PSI', 'CR_FEM', 'CR_ERR']: - value_parameter_error('typ_CR', typ_CR, - ['CR_PSI', 'CR_FEM', 'CR_ERR']) + value_parameter_error('typ_CR', typ_CR, ['CR_PSI', 'CR_FEM', 'CR_ERR']) if typ_CR == 'CR_PSI': # Minimum coefficient reflection @@ -252,15 +250,14 @@ def Zi(p, alpha, ele_type): Z_fem : `float` Parameter for the spurious reflection coefficient in FEM ''' + value_parameter_error('ele_type', ele_type, ['lumped', 'consistent']) + if ele_type == 'lumped': m1 = 1 / 2 m2 = 0. elif ele_type == 'consistent': m1 = 1 / 3 m2 = 1 / 6 - else: - value_parameter_error('ele_type', ele_type, - ['lumped', 'consistent']) Z_fem = m2 * (np.cos(alpha * p) - 1) / ( m1 * (np.cos(alpha * p) + 1)) diff --git a/spyro/habc/error_measure.py b/spyro/habc/error_measure.py index 847c607d..5761206f 100644 --- a/spyro/habc/error_measure.py +++ b/spyro/habc/error_measure.py @@ -360,6 +360,9 @@ def __init__(self, dt, freq_Nyq, receiver_locations, forward_solution_receivers= # max_errIt = dat_reg_xCR[1] # max_errPK = dat_reg_xCR[2] + # value_parameter_error('crit_opt', crit_opt, + # ['err_difference', 'err_integral', 'err_sum']) + # if crit_opt == 'err_difference': # y_err = [eI - eP for eI, eP in zip(max_errIt, max_errPK)] @@ -369,11 +372,6 @@ def __init__(self, dt, freq_Nyq, receiver_locations, forward_solution_receivers= # elif crit_opt == 'err_sum': # y_err = [eI + eP for eI, eP in zip(max_errIt, max_errPK)] - # else: - # value_parameter_error( - # 'crit_opt', crit_opt, - # ['err_difference', 'err_integral', 'err_sum']) - # # Limits for the heuristic factor # xCR_inf, xCR_sup = self.xCR_lim diff --git a/spyro/habc/habc.py b/spyro/habc/habc.py index 839c50aa..f76f1b60 100644 --- a/spyro/habc/habc.py +++ b/spyro/habc/habc.py @@ -6,7 +6,8 @@ from ..abc.abc_layer import ABCLayer from .damp_profile import HABC_Damping from ..solvers.modal.modal_sol import Modal_Solver -from ..utils.typing import HyperLayerDegreeType, LayerShapeType, LayerSizeRefFrequency +from ..utils.typing import (HyperLayerDegreeType, LayerDampingType, + LayerShapeType, LayerSizeRefFrequency) from ..io.basicio import parallel_print as pprint # from sympy import divisors # from spyro.utils.error_management import value_parameter_error @@ -29,6 +30,13 @@ class HABCLayer(ABCLayer, HABC_Damping): abc_boundary_layer_shape : `typing.LayerShapeType`, optional Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. + abc_boundary_layer_type : `typing.LayerDampingType` + Type of the boundary layer. Options: `LayerDampingType.LOCAL`, + `LayerDampingType.HYBRID`, `LayerDampingType.PML` or `LayerDampingType.NOABCS`. + Default is `LayerDampingType.NOABCS` where no absorbing BCs are applied. + Option `LayerDampingType.HYBRID` is based on paper of Salas et al. (2022). + doi: https://doi.org/10.1016/j.apm.2022.09.014 + TODO: Add citation abc_deg_layer : `int` or `float` or `None`, optional Hypershape degree. For hypershape layers, the degree must be greater than or equal to 2. `None` is used only for rectangular layers. Default is `None`. @@ -200,7 +208,7 @@ def __init__(self, domain_dim, frequency, f_Nyquist, abc_deg_layer, ABCLayer.__init__(self, domain_dim, frequency, f_Nyquist, dimension=dimension, quadrilateral=quadrilateral, func_space_type=func_space_type, abc_boundary_layer_shape=abc_boundary_layer_shape, - abc_boundary_layer_type="hybrid", + abc_boundary_layer_type=LayerDampingType.HYBRID, abc_reference_freq=abc_reference_freq, abc_degree_type=abc_degree_type, abc_deg_layer=abc_deg_layer, output_folder=output_folder, comm=comm) diff --git a/spyro/io/boundary_layer_io.py b/spyro/io/boundary_layer_io.py index 5b6dc5f4..58c45752 100644 --- a/spyro/io/boundary_layer_io.py +++ b/spyro/io/boundary_layer_io.py @@ -1,6 +1,7 @@ -from ..utils.error_management import (enum_parameter_error, value_numerical_error, - value_parameter_error) -from ..utils.typing import HyperLayerDegreeType, LayerShapeType, LayerSizeRefFrequency +from ..io.basicio import parallel_print as pprint +from ..utils.error_management import enum_parameter_error, value_numerical_error +from ..utils.typing import (HyperLayerDegreeType, LayerDampingType, + LayerShapeType, LayerSizeRefFrequency) class Read_boundary_layer: @@ -9,15 +10,16 @@ class Read_boundary_layer: Attributes ---------- - abc_boundary_layer_shape : `typing.LayerShapeType`, optional + abc_boundary_layer_shape : `typing.LayerShapeType` Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. - abc_boundary_layer_type : `str` - Type of the boundary layer. Options: 'hybrid' or 'PML'. - Option 'hybrid' is based on paper of Salas et al. (2022). + abc_boundary_layer_type : `typing.LayerDampingType` + Type of the boundary layer. Options: `LayerDampingType.LOCAL`, + `LayerDampingType.HYBRID`, `LayerDampingType.PML` or `LayerDampingType.NOABCS`. + Default is `LayerDampingType.NOABCS` where no absorbing BCs are applied. + Option `LayerDampingType.HYBRID` is based on paper of Salas et al. (2022). doi: https://doi.org/10.1016/j.apm.2022.09.014 - abc_cmax : float - Maximum acoustic wave velocity in PML - km/s + TODO: Add citation abc_deg_eikonal : `int` Finite element order for the Eikonal analysis abc_deg_layer : `int` or `float` @@ -25,8 +27,6 @@ class Read_boundary_layer: abc_degree_type : `typing.HyperLayerDegreeType`, optional Type of the hypereshape degree. Options: 'HyperLayerDegreeType.REAL' or 'HyperLayerDegreeType.INTEGER'. Default is 'HyperLayerDegreeType.REAL' - abc_exponent : `float` - Exponent of the polynomial damping abc_extend_properties : `str` Mode to extend the properties into the absorbing layer. Options: 'abc_driven' (performed by a specific method) or @@ -35,8 +35,18 @@ class Read_boundary_layer: If True, the infinite model is created abc_pad_length : `float` Thickness of the PML in the z-direction (km) - always positive - abc_R : `float` - Theoretical reflection coefficient + abc_pml_cmax: float + Maximum propagation speed (km/s) in the PML layer. + Default is 4.7 km/s because of the usual value of a salt layer in Brazil. + abc_pml_exponent: int + Exponent for the polynomial damping profile of the PML layer. Default is 2. + abc_pml_R: float + Theoretical reflection coefficient of the PML layer. Default is 1e-6 assuming an + amplitude attenuation of 120 dB. If the source's amplitude is extremely small, + this attenuation might be excessive, resulting in a large PML layer. See: + https://ccrma.stanford.edu/~jos/mdft/Exponentials.html#fig:exponential + https://ccrma.stanford.edu/~jos/mdft/Audio_Decay_Time_T60.html + TODO: Add citation abc_reference_freq : `typing.LayerSizeRefFrequency`, optional Reference frequency for sizing the absorbing layer. Options: 'LayerSizeRefFrequency.SOURCE' or 'LayerSizeRefFrequency.BOUNDARY'. @@ -44,8 +54,8 @@ class Read_boundary_layer: abc_user_pad_length : `bool` If True, the pad length is provided by the user. If False, the pad length is determined with the HABC criterion. - damping_type : `str` - Type of the boundary layer + abc_user_pml_cmax : `bool` + If True, the maximum propagation speed in the PML layer is provided by the user. dictionary : `dict` Dictionary containing the boundary layer information """ @@ -73,10 +83,10 @@ def __init__(self, comm=None): "absorving_boundary_conditions"]["status"] self.input_dictionary[ "absorving_boundary_conditions"].setdefault("damping_type", None) - self.input_dictionary[ - "absorving_boundary_conditions"].setdefault("pad_length", None) self.abc_boundary_layer_type = self.input_dictionary[ "absorving_boundary_conditions"]["damping_type"] + self.input_dictionary[ + "absorving_boundary_conditions"].setdefault("pad_length", None) self.abc_pad_length = self.input_dictionary[ "absorving_boundary_conditions"]["pad_length"] @@ -96,18 +106,18 @@ def __init__(self, comm=None): @property def abc_boundary_layer_shape(self): - if not hasattr(self, '_abc_boundary_layer_shape'): + if not hasattr(self, "_abc_boundary_layer_shape"): self._abc_boundary_layer_shape = LayerShapeType.NOLAYER return self._abc_boundary_layer_shape @abc_boundary_layer_shape.setter def abc_boundary_layer_shape(self, value): """Set boundary layer shape with enum validation.""" - shape_enum = enum_parameter_error('abc_boundary_layer_shape', + shape_enum = enum_parameter_error("abc_boundary_layer_shape", value, LayerShapeType) if shape_enum == LayerShapeType.NOLAYER: - raise ValueError("NOLAYER not allowed for active ABC.") + raise ValueError("'NOLAYER' option not allowed for active ABC.") self._abc_boundary_layer_shape = shape_enum @@ -118,7 +128,7 @@ def abc_reference_freq(self): @abc_reference_freq.setter def abc_reference_freq(self, value): """Set reference frequency for sizing the absorbing layer with enum validation.""" - reference_freq_enum = enum_parameter_error('abc_reference_freq', value, + reference_freq_enum = enum_parameter_error("abc_reference_freq", value, LayerSizeRefFrequency) self._abc_reference_freq = reference_freq_enum @@ -129,10 +139,48 @@ def abc_degree_type(self): @abc_degree_type.setter def abc_degree_type(self, value): """Set hypershape degree type for hypershape layers with enum validation.""" - degree_type_enum = enum_parameter_error('abc_degree_type', value, + degree_type_enum = enum_parameter_error("abc_degree_type", value, HyperLayerDegreeType) self._abc_degree_type = degree_type_enum + @property + def abc_pml_exponent(self): + return self._abc_pml_exponent + + @abc_pml_exponent.setter + def abc_pml_exponent(self, value): + """Set the exponent for the polynomial damping profile in PML with validation.""" + pml_exponent = value_numerical_error("abc_pml_exponent", value, integer_num=True, + lower_bound=1, include_lower_bound=True) + self._abc_pml_exponent = pml_exponent + + @property + def abc_pml_R(self): + return self._abc_pml_R + + @abc_pml_R.setter + def abc_pml_R(self, value): + """Set the theoretical reflection coefficient in the PML layer with validation.""" + pml_R = value_numerical_error("abc_pml_R", value, float_num=True, lower_bound=0.) + self._abc_pml_R = pml_R + + @property + def abc_pml_cmax(self): + return self._abc_pml_cmax + + @abc_pml_cmax.setter + def abc_pml_cmax(self, value): + """Set the maximum propagation speed in the PML layer with validation.""" + self.abc_user_pml_cmax = True + if value is None: + pprint("Maximum propagation speed will get from model", comm=self.comm) + self.abc_user_pml_cmax = False + pml_cmax = value + else: + pml_cmax = value_numerical_error("abc_pml_cmax", value, float_num=True, + integer_num=True, lower_bound=0.) + self._abc_pml_cmax = pml_cmax + @property def abc_boundary_layer_type(self): return self._abc_boundary_layer_type @@ -141,30 +189,32 @@ def abc_boundary_layer_type(self): def abc_boundary_layer_type(self, value): """Set the type of absorbing boundary layer with validation.""" abc_dictionary = self.input_dictionary['absorving_boundary_conditions'] - accepted_damping_types = [ - "PML", - "local", - "hybrid", - None, - ] # Cheking damping type input - if value not in accepted_damping_types: - value_parameter_error('abc_boundary_layer_type', value, - accepted_damping_types) + if value is None: + pprint("No Absorbing Boundary Conditions (ABCs) applied.", comm=self.comm) + value = "no_abcs" + + self._abc_boundary_layer_type = enum_parameter_error( + "abc_boundary_layer_type", value, LayerDampingType) if value == "PML": # PML forces rectangular shape self.abc_boundary_layer_shape = LayerShapeType.RECTANGULAR + + # PML-specific parameters with defaults abc_dictionary.setdefault("exponent", 2) - self.abc_exponent = abc_dictionary["exponent"] + self.abc_pml_exponent = abc_dictionary["exponent"] abc_dictionary.setdefault("R", 1e-6) - self.abc_R = abc_dictionary["R"] - abc_dictionary.setdefault("cmax", 4.7) - self.abc_cmax = abc_dictionary["cmax"] + self.abc_pml_R = abc_dictionary["R"] + abc_dictionary.setdefault("cmax", None) + self.abc_pml_cmax = abc_dictionary["cmax"] + if value == "hybrid": + # Get shape from dictionary, default to rectangular - self.abc_boundary_layer_shape = abc_dictionary.get("layer_shape", "rectangular") + self.abc_boundary_layer_shape = abc_dictionary.get("layer_shape", + "rectangular") # Hypershape-specific validation self.abc_deg_layer = None @@ -177,7 +227,6 @@ def abc_boundary_layer_type(self, value): self.abc_degree_type = abc_dictionary.get("degree_type", "real") # Common parameters for both PML and hybrid - self._abc_boundary_layer_type = value self.abc_reference_freq = abc_dictionary.get("abc_reference_freq", "source") self.abc_deg_eikonal = abc_dictionary.get("degree_eikonal", 2) self.abc_get_ref_model = abc_dictionary.get("get_ref_model", False) @@ -200,16 +249,22 @@ def abc_pad_length(self, value): Returns ------- None - """ - if isinstance(value, (int, float)) and value < 0: - raise ValueError("Pad length must be positive") + Notes + ----- + For the HABC criterion see Salas et al (2022): Hybrid absorbing scheme based on + hyperelliptical layers with non-reflecting boundary conditions in scalar wave + equations. doi: https://doi.org/10.1016/j.apm.2022.09.014 + TODO: Add citation + """ self.abc_user_pad_len = True if value is None: - print("Pad length will be determined with HABC criterion", flush=True) - value = 0. + pprint("Pad length will be determined with HABC criterion", comm=self.comm) + pad_length = 0. self.abc_user_pad_len = False - - print(f"Pad length provided by user (km): {value:.4f}", flush=True) - self._abc_pad_length = value + else: + pad_length = value_numerical_error("abc_pad_length", value, float_num=True, + integer_num=True, lower_bound=0.) + pprint(f"Pad length provided by user (km): {pad_length:.4f}", comm=self.comm) + self._abc_pad_length = pad_length diff --git a/spyro/io/dictionaryio.py b/spyro/io/dictionaryio.py index 01189347..7cb382c9 100644 --- a/spyro/io/dictionaryio.py +++ b/spyro/io/dictionaryio.py @@ -1,4 +1,4 @@ -from ..utils.error_management import value_parameter_error +from ..utils.error_management import value_numerical_error, value_parameter_error class Read_options: @@ -51,8 +51,7 @@ def __init__(self, dictionary={}): self.variant = options_dictionary["variant"] self.method = options_dictionary["method"] - if options_dictionary["cell_type"] is not None: - self.cell_type = options_dictionary["cell_type"] + self.cell_type = options_dictionary["cell_type"] self.degree = options_dictionary["degree"] self.dimension = options_dictionary["dimension"] self.analysis = options_dictionary["analysis"] @@ -64,9 +63,7 @@ def variant(self): @variant.setter def variant(self, value): accepted_variants = ["lumped", "equispaced", "DG", None] - if value not in accepted_variants: - raise ValueError(f"Variant of {value} is not valid.") - self._variant = value + self._variant = value_parameter_error("variant", value, accepted_variants) @property def method(self): @@ -91,6 +88,11 @@ def method(self, value): "DGQ", "discontinuous_galerkin_quadrilateral", ] + + # Check if the provided value is in any of the accepted methods + value_parameter_error("method", value, mlt_equivalents + sem_equivalents + + dg_t_equivalents + dg_q_equivalents + ["CG", None]) + if value in mlt_equivalents: self._method = "mass_lumped_triangle" self.cell_type = "triangle" @@ -103,21 +105,13 @@ def method(self, value): elif value in dg_q_equivalents: self._method = "DG_quadrilateral" self.cell_type = "quadrilateral" - elif value == "DG": - raise ValueError( - "DG is not a valid method. Please specify \ - either DG_triangle or DG_quadrilateral." - ) elif value == "CG": - if "variant" in self.input_dictionary["options"] and "cell_type" \ - in self.input_dictionary["options"]: - self._method = "CG" - else: + options = self.input_dictionary["options"] + if not ("variant" in options and "cell_type" in options): raise ValueError("Cant use CG without specifying cell type and variant.") + self._method = "CG" elif value is None: self._method = None - else: - raise ValueError(f"Method of {value} is not valid.") @property def cell_type(self): @@ -183,9 +177,9 @@ def degree(self): @degree.setter def degree(self, value): - if not isinstance(value, int): - raise ValueError("Degree has to be integer") - self._degree = value + self._degree = value_numerical_error('degree', value, float_num=False, + integer_num=True, lower_bound=0, + include_lower_bound=False,) @property def dimension(self): @@ -193,9 +187,7 @@ def dimension(self): @dimension.setter def dimension(self, value): - if value not in {2, 3}: - raise ValueError(f"Dimension of {value} not 2 or 3.") - self._dimension = value + self._dimension = value_parameter_error('dimension', value, [2, 3]) @property def analysis(self): @@ -204,9 +196,7 @@ def analysis(self): @analysis.setter def analysis(self, value): allowed_analyses = ["transient", "modal", "eikonal"] - if value not in allowed_analyses: - raise value_parameter_error('analysis', value, allowed_analyses) - self._analysis = value + self._analysis = value_parameter_error('analysis', value, allowed_analyses) class Read_outputs: diff --git a/spyro/io/material_properties_io.py b/spyro/io/material_properties_io.py index a2bc7835..61188b1c 100644 --- a/spyro/io/material_properties_io.py +++ b/spyro/io/material_properties_io.py @@ -43,10 +43,8 @@ def define_property_function_space(wave, func_space_type, dg_property, # Checking input arguments opts_func_space_type = ["scalar", "vector", "tensor"] - if func_space_type not in opts_func_space_type: - error_management.value_parameter_error("func_space_type", - func_space_type, - opts_func_space_type) + error_management.value_parameter_error("func_space_type", func_space_type, + opts_func_space_type) if dg_property is False and func_space_type == "scalar": return point_to_scalar_wave_function_space(wave) diff --git a/spyro/io/model_parameters.py b/spyro/io/model_parameters.py index c299f3fa..debbf801 100644 --- a/spyro/io/model_parameters.py +++ b/spyro/io/model_parameters.py @@ -79,31 +79,46 @@ class Model_parameters(Read_options, Read_boundary_layer, Firedrake mesh. abc_active: bool Whether or not the absorbing boundary conditions are used. - abc_exponent: int - Exponent of the absorbing boundary conditions. - abc_cmax: float - Maximum acoustic wave velocity in the absorbing boundary conditions. - abc_R: float - Theoretical reflection coefficient of the absorbing boundary - conditions. - abc_pad_length: float - Thickness of the absorbing boundary conditions. - abc_boundary_layer_type : `str` - Type of the boundary layer. Option 'hybrid' is based on paper - of Salas et al. (2022). doi: https://doi.org/10.1016/j.apm.2022.09.014 - abc_boundary_layer_shape : str - Shape type of pad layer. Options: 'rectangular' or 'hypershape' - abc_deg_layer : `float` - Hypershape degree - abc_degree_type : `str` - Type of the hypereshape degree. Options: 'real' or 'integer' - abc_reference_freq : `str` - Reference frequency for sizing the hybrid absorbing layer. - Options: 'source' or 'boundary' + abc_boundary_layer_shape : `typing.LayerShapeType` + Shape type of the pad layer. Options: `LayerShapeType.RECTANGULAR` or + `LayerShapeType.HYPERSHAPE`. Default is `LayerShapeType.RECTANGULAR`. + abc_boundary_layer_type : `typing.LayerDampingType` + Type of the boundary layer. Options: `LayerDampingType.LOCAL`, + `LayerDampingType.HYBRID`, `LayerDampingType.PML` or `LayerDampingType.NOABCS`. + Default is `LayerDampingType.NOABCS` where no absorbing BCs are applied. + Option `LayerDampingType.HYBRID` is based on paper of Salas et al. (2022). + doi: https://doi.org/10.1016/j.apm.2022.09.014 + TODo: Add citation abc_deg_eikonal : `int` Finite element order for the Eikonal analysis + abc_deg_layer : `int` or `float` + Hypershape degree + abc_degree_type : `typing.HyperLayerDegreeType`, optional + Type of the hypereshape degree. Options: 'HyperLayerDegreeType.REAL' or + 'HyperLayerDegreeType.INTEGER'. Default is 'HyperLayerDegreeType.REAL' + abc_extend_properties : `str` + Mode to extend the properties into the absorbing layer. + Options: 'abc_driven' (performed by a specific method) or + 'builtin' (automatic at field definition) abc_get_ref_model : `bool` If True, the infinite model is created + abc_pad_length : `float` + Thickness of the PML in the z-direction (km) - always positive + abc_pml_cmax: float + Maximum propagation speed (km/s) in the PML layer. Default is 4.7 km/s. + abc_pml_exponent: int + Exponent for the polynomial damping profile of the PML layer. Default is 2. + abc_pml_R: float + Theoretical reflection coefficient of the PML layer. Default is 1e-6. + abc_reference_freq : `typing.LayerSizeRefFrequency`, optional + Reference frequency for sizing the absorbing layer. + Options: 'LayerSizeRefFrequency.SOURCE' or 'LayerSizeRefFrequency.BOUNDARY'. + Default is 'LayerSizeRefFrequency.SOURCE'. + abc_user_pad_length : `bool` + If True, the pad length is provided by the user. If False, + the pad length is determined with the HABC criterion. + abc_user_pml_cmax : `bool` + If True, the maximum propagation speed in the PML layer is provided by the user. source_type: str Type of source used in the simulation. Can be "ricker" for a Ricker wavelet or "MMS" for a manufactured solution. diff --git a/spyro/meshing/meshing_habc.py b/spyro/meshing/meshing_habc.py index 4e2ecf7a..e76ff62a 100644 --- a/spyro/meshing/meshing_habc.py +++ b/spyro/meshing/meshing_habc.py @@ -3,14 +3,14 @@ from netgen.geom2d import SplineGeometry from netgen.meshing import Element2D, Element3D, FaceDescriptor, Mesh, MeshPoint from scipy.spatial import cKDTree -from spyro.domains.space import create_function_space -from spyro.meshing.meshing_functions import AutomaticMesh -from spyro.meshing.meshing_operations import MeshOps -from spyro.tools.habc_tools import point_cloud_field -from spyro.utils.error_management import value_parameter_error +from ..domains.space import create_function_space +from ..io.basicio import parallel_print as pprint +from .meshing_functions import AutomaticMesh +from .meshing_operations import MeshOps +from ..tools.habc_tools import point_cloud_field +from ..utils.error_management import value_parameter_error from ..tools.version_control import is_firedrake_new - if is_firedrake_new() is False: from firedrake.__future__ import interpolate fire.interpolate = interpolate @@ -124,7 +124,7 @@ def original_boundary_data(self, mesh, function_space, mesh_parameters, function_space : `Firedrake.FunctionSpace` Function space for the current mesh operations. mesh_parameters : `meshing_parameters.MeshingParameters` - Contains mesh parameters + Contains mesh parameters. initial_velocity_model : `Firedrake.Function` Initial velocity model. @@ -138,18 +138,16 @@ def original_boundary_data(self, mesh, function_space, mesh_parameters, Mesh node coordinates on boundaries of the original domain. """ - print("Getting Boundary Mesh Data from Original Domain", flush=True) + pprint("Getting Boundary Mesh Data from Original Domain", comm=self.comm) # Extract node positions node_positions = self.extract_node_positions(mesh, function_space, output_type="array") # Extract boundary node positions - all_bnd_nodes = [] - for (bnd_ids, status) in mesh_parameters.boundary_nodes_ids.values(): - if status: - all_bnd_nodes.append(bnd_ids) - all_bnd_nodes = np.unique(np.concatenate(all_bnd_nodes)) + all_bnd_nodes = np.unique(np.concatenate([ + bnd_ids for bnd_ids, status + in mesh_parameters.boundary_nodes_ids.values() if status])) coord_msh = mesh.coordinates.dat.data_with_halos coord_bnd_nodes = node_positions[all_bnd_nodes, :] @@ -171,12 +169,11 @@ def original_boundary_data(self, mesh, function_space, mesh_parameters, # Print on screen cbnd_str = "Boundary Velocity Range (km/s): {:.3f} - {:.3f}" - print(cbnd_str.format(c_bnd_min, c_bnd_max), flush=True) + pprint(cbnd_str.format(c_bnd_min, c_bnd_max), comm=self.comm) return c_bnd_min, c_bnd_max, coord_bnd_nodes - def creating_velocity_profile(self, function_space, - initial_velocity_model, path_save): + def creating_velocity_profile(self, function_space, initial_velocity_model, path_save): """Create the velocity profile for the original domain. Parameters @@ -209,7 +206,7 @@ def creating_velocity_profile(self, function_space, # Print on screen cdom_str = "Domain Velocity Range (km/s): {:.3f} - {:.3f}" - print(cdom_str.format(c_min, c_max), flush=True) + pprint(cdom_str.format(c_min, c_max), comm=self.comm) # Save initial velocity model vel_c = fire.VTKFile(path_save + "preamble/c_vel.pvd") @@ -236,19 +233,18 @@ def create_function_space_eik(self, mesh, degree_eik, ele_type_eik='consistent') Function space for the Eikonal modeling. """ - print("Setting Mesh Properties for Eikonal Analysis", flush=True) + pprint("Setting Mesh Properties for Eikonal Analysis", comm=self.comm) - allowed_ele_types = ['consistent', 'underintegrated'] - if ele_type_eik not in allowed_ele_types: - value_parameter_error('ele_type_eik', ele_type_eik, allowed_ele_types) + allowed_ele_types = ["consistent", "underintegrated"] + value_parameter_error('ele_type_eik', ele_type_eik, allowed_ele_types) # Function space for the Eikonal modeling - if ele_type_eik == 'consistent': - funct_space_eik = create_function_space(mesh, 'CG', degree_eik) + if ele_type_eik == "consistent": + funct_space_eik = create_function_space(mesh, "CG", degree_eik) - if ele_type_eik == 'underintegrated': - method = 'spectral_quadrilateral' if self.quadrilateral \ - else 'mass_lumped_triangle' + if ele_type_eik == "underintegrated": + method = "spectral_quadrilateral" if self.quadrilateral \ + else "mass_lumped_triangle" degree = min(degree_eik, 4 if self.dimension == 2 else 3) funct_space_eik = create_function_space(mesh, method, degree) @@ -299,11 +295,11 @@ def preamble_mesh_operations(self, Wave, ele_type_eik='consistent', f_est=0.03): Function space for the Eikonal modeling. """ - print("\nCreating Mesh and Initial Velocity Model", flush=True) + pprint("\nCreating Mesh and Initial Velocity Model", comm=self.comm) # Mesh data - print(f"Original Mesh with {Wave.mesh.num_vertices()} Nodes " - f"and {Wave.mesh.num_cells()} Volume Elements", flush=True) + pprint(f"Original Mesh with {Wave.mesh.num_vertices()} Nodes " + f"and {Wave.mesh.num_cells()} Volume Elements", comm=self.comm) # Save a copy of the original mesh Wave.mesh_original = Wave.mesh @@ -428,14 +424,14 @@ def trunc_hyp_bndpts_2D(self, hyp_par, xdom, z0): or pnt_bef_trunc < 3 or pnt_aft_trunc < 3: num_bnd_pts += 1 - print(pnt_str, f"Complete Hyperellipse: {num_bnd_pts}", flush=True) + pprint(f"{pnt_str} Complete Hyperellipse: {num_bnd_pts}", comm=self.comm) bnd_pts = self.bnd_pnts_hyp_2D(a_hyp, b_hyp, n_hyp, num_bnd_pts) # Filter hyperellipse points based on the truncation plane z0 filt_bnd_pts = np.array([point for point in bnd_pts if point[1] <= z0]) - print(pnt_str, f"Truncated Hyperellipse: {len(filt_bnd_pts)}", - flush=True) + pprint(f"{pnt_str} Truncated Hyperellipse: {len(filt_bnd_pts)}", + comm=self.comm) # Identify truncation index ini_trunc = max(np.where(bnd_pts[:, 1] > z0)[0][0] - 1, 0) @@ -509,12 +505,12 @@ def create_bnd_mesh_2D(self, geo, bnd_pts, trunc_feat, spln): p2 = geo.PointData()[2][idp + 1] p3 = geo.PointData()[2][idp + 2] curves.append(["spline3", p1, p2, p3, ltrunc]) - # print(p1, p2, p3) + # pprint(p1, p2, p3) for idp in range(ini_trunc, end_trunc): p1 = geo.PointData()[2][idp] p2 = geo.PointData()[2][idp + 1] - # print(p1, p2) + # pprint(p1, p2) if idp == ini_trunc or idp == end_trunc - 1: curves.append(["line", p1, p2, ltrunc]) @@ -526,14 +522,14 @@ def create_bnd_mesh_2D(self, geo, bnd_pts, trunc_feat, spln): p2 = geo.PointData()[2][idp + 1] p3 = geo.PointData()[2][idp + 2] curves.append(["spline3", p1, p2, p3, ltrunc]) - # print(p1, p2, p3) + # pprint(p1, p2, p3) else: # Mesh with line segments for idp in range(0, num_bnd_pts - 1, 2): p1 = geo.PointData()[2][idp] p2 = geo.PointData()[2][idp + 1] - # print(p1, p2) + # pprint(p1, p2) if ini_trunc + 1 <= idp <= end_trunc - 2: curves.append(["line", p1, p2, self.lmin]) @@ -595,21 +591,19 @@ def create_hyp_trunc_mesh_2D(self, hyp_par, spln=True): optsteps2d=10, # Optimize mesh ) hyp_mesh.Compress() - print("Hyperelliptical Mesh Generated Successfully", - flush=True) + pprint("Hyperelliptical Mesh Generated Successfully", comm=self.comm) break except Exception as e: # Retry with lines if splines fail if spln: - print(f"Error Meshing with Splines: {e}", flush=True) - print("Now meshing with Lines", flush=True) + pprint(f"Error Meshing with Splines: {e}", comm=self.comm) + pprint("Now meshing with Lines", comm=self.comm) spln = False else: - print(f"Error Meshing with Lines: {e}. Exiting.", - flush=True) + pprint(f"Error Meshing with Lines: {e}. Exiting.", comm=self.comm) break # Mesh is transformed into a firedrake mesh @@ -723,17 +717,17 @@ def merge_mesh_2D(self, rec_mesh, hyp_mesh): try: # Mesh data final_mesh.Compress() - print(f"Mesh created with {len(final_mesh.Points())} points " - f"and {len(final_mesh.Elements2D())} elements", flush=True) + pprint(f"Mesh created with {len(final_mesh.Points())} points " + f"and {len(final_mesh.Elements2D())} elements", comm=self.comm) # Mesh is transformed into a firedrake mesh q = {"overlap_type": (fire.DistributedMeshOverlapType.NONE, 0)} final_mesh = fire.Mesh( final_mesh, distribution_parameters=q, comm=self.comm.comm) - print("Merged Mesh Generated Successfully", flush=True) + pprint("Merged Mesh Generated Successfully", comm=self.comm) except Exception as e: - print(f"Error Generating Merged Mesh: {e}. Exiting.", flush=True) + pprint(f"Error Generating Merged Mesh: {e}. Exiting.", comm=self.comm) return final_mesh @@ -832,18 +826,18 @@ def sharp_mesh_3D(self, rec_mesh, hyp_par, centroid): try: # Mesh data sharp_mesh.Compress() - print(f"Mesh created with {len(sharp_mesh.Points())} points " - f"and {len(sharp_mesh.Elements3D())} elements", flush=True) + pprint(f"Mesh created with {len(sharp_mesh.Points())} points " + f"and {len(sharp_mesh.Elements3D())} elements", comm=self.comm) # Mesh is transformed into a firedrake mesh q = {"overlap_type": (fire.DistributedMeshOverlapType.NONE, 0)} sharp_mesh = fire.Mesh( sharp_mesh, distribution_parameters=q, comm=self.comm.comm) - print("Sharp Mesh Generated Successfully", flush=True) + pprint("Sharp Mesh Generated Successfully", comm=self.comm) # fire.VTKFile("output/sharp_mesh.pvd").write(sharp_mesh) except Exception as e: - print(f"Error Generating Merged Mesh: {e}. Exiting.", flush=True) + pprint(f"Error Generating Merged Mesh: {e}. Exiting.", comm=self.comm) return sharp_mesh @@ -926,10 +920,10 @@ def snap_nodes_to_hyp(self, mesh, hyp_par, centroid, plane_tol=1e-5): coords = mesh.coordinates.dat.data_with_halos min_z, min_x, min_y = np.min(coords, axis=0) max_z, max_x, max_y = np.max(coords, axis=0) - print("Mesh Bounds Detected:", flush=True) - print(f" X: [{min_x:.4f}, {max_x:.4f}]", flush=True) - print(f" Y: [{min_y:.4f}, {max_y:.4f}]", flush=True) - print(f" Z: [{min_z:.4f}, {max_z:.4f}]", flush=True) + pprint("Mesh Bounds Detected:", comm=self.comm) + pprint(f" X: [{min_x:.4f}, {max_x:.4f}]", comm=self.comm) + pprint(f" Y: [{min_y:.4f}, {max_y:.4f}]", comm=self.comm) + pprint(f" Z: [{min_z:.4f}, {max_z:.4f}]", comm=self.comm) # Select nodes to snap mask_min_z = np.isclose(coords[:, 0], min_z, atol=plane_tol) @@ -948,8 +942,8 @@ def snap_nodes_to_hyp(self, mesh, hyp_par, centroid, plane_tol=1e-5): coords[pnt, :], centroid, b_hyp, a_hyp, c_hyp, n_hyp) coords[pnt, 0] = np.clip(coords[pnt, 0], -np.inf, max_z) - print(f"Boundary Nodes Snapped: {len(pnts_to_snap)}", flush=True) - print("Snapped Mesh Generated Successfully", flush=True) + pprint(f"Boundary Nodes Snapped: {len(pnts_to_snap)}", comm=self.comm) + pprint("Snapped Mesh Generated Successfully", comm=self.comm) return mesh @@ -1063,48 +1057,44 @@ def hypershape_mesh_habc(self, hyp_par, mesh_original, mesh_parameters, spln=Tru return mesh_habc - def layer_boundary_data(self, V): + def layer_boundary_data(self, mesh, V, mesh_parameters): """Generate the boundary data from the domain with the absorbing layer. Parameters ---------- + mesh : `Firedrake.Mesh` + Current mesh. V : `Firedrake.FunctionSpace` Function space for the boundary of the domain with absorbing layer. + mesh_parameters : `meshing_parameters.MeshingParameters` + Contains mesh parameters. Returns ------- - bnd_nfs : 'array' - Mesh node indices on non-free surfaces. + bnd_nod_ids_nfs : 'array' + Mesh node indices on non-free surfaces of the domain with absorbing layer. bnd_nodes_nfs : `tuple` - Mesh node coordinates on non-free surfaces. - - (z_data[nfs_idx], x_data[nfs_idx]) for 2D. - - (z_data[nfs_idx], x_data[nfs_idx], y_data[nfs_idx]) for 3D. + Mesh node coordinates on non-free surfaces of the domain with absorbing layer. """ # Boundary nodes indices - bnd_nod = fire.DirichletBC(V, 0., "on_boundary").nodes + bnd_nod_ids_nfs = np.unique(np.concatenate([ + bnd_ids for bnd_ids, status + in mesh_parameters.boundary_nodes_ids.values() if status])) # Extract node positions - node_positions = self.extract_node_positions(self.mesh, V) - - # Boundary node coordinates - z_f, x_f = node_positions[:2] - bnd_z = z_f[bnd_nod] - bnd_x = x_f[bnd_nod] - - # Identify non-free surfaces (remain unchanged) - no_free_surf = ~(abs(bnd_z) <= self.mesh_parameters.tol) - - bnd_nodes_nfs = (bnd_z[no_free_surf], bnd_x[no_free_surf]) - if self.dimension == 3: # 3D - y_f = node_positions[2] - bnd_y = y_f[bnd_nod] - bnd_nodes_nfs += (bnd_y[no_free_surf],) + node_positions = self.extract_node_positions(mesh, V, output_type="array") + coord_msh = mesh.coordinates.dat.data_with_halos + coord_bnd_nodes = node_positions[bnd_nod_ids_nfs, :] - # Boundary node indices on non-free surfaces - bnd_nfs = bnd_nod[no_free_surf] + # Identify the boundary nodes + tree = cKDTree(coord_msh) + indices = tree.query(coord_bnd_nodes, k=1, + distance_upper_bound=mesh_parameters.tol)[1] + mask_boundary = indices[indices < len(coord_msh)] + bnd_nodes_nfs = mesh.coordinates.dat.data_with_halos[mask_boundary, :] - return bnd_nfs, bnd_nodes_nfs + return bnd_nod_ids_nfs, bnd_nodes_nfs def get_spatial_coordinates_abc(self, mesh, domain_layer): """Get the ufl coordinates of the mesh with absorbing layer. diff --git a/spyro/meshing/meshing_operations.py b/spyro/meshing/meshing_operations.py index a3374aa6..6cf28d23 100644 --- a/spyro/meshing/meshing_operations.py +++ b/spyro/meshing/meshing_operations.py @@ -219,6 +219,7 @@ def extract_node_positions(self, mesh, function_space, output_type="tuple"): coords.append(assemble(interpolate(ufl_input, V))) # Get the node positions + value_parameter_error('output_type', output_type, ["tuple", "array"]) if output_type == "tuple": z_data = coords[0].dat.data_with_halos[:] x_data = coords[1].dat.data_with_halos[:] @@ -230,9 +231,6 @@ def extract_node_positions(self, mesh, function_space, output_type="tuple"): elif output_type == "array": node_positions = column_stack([comp.dat.data_with_halos[:] for comp in coords]) - - else: - value_parameter_error('output_type', output_type, ["tuple", "array"]) del coords # TODO: either put wave's similar method here or just remove the one from wave diff --git a/spyro/meshing/meshing_parameters.py b/spyro/meshing/meshing_parameters.py index f6998770..96558a16 100644 --- a/spyro/meshing/meshing_parameters.py +++ b/spyro/meshing/meshing_parameters.py @@ -675,12 +675,10 @@ def method(self, value): "spectral_quadrilateral", "DG_quadrilateral", "CG", + None, ] - if value is not None and value not in allowed_types: - value_parameter_error("method", value, allowed_types) - - self._method = value + self._method = value_parameter_error("method", value, allowed_types) @property def mesh_file(self): diff --git a/spyro/pml/damping.py b/spyro/pml/damping.py deleted file mode 100644 index 61248de9..00000000 --- a/spyro/pml/damping.py +++ /dev/null @@ -1,130 +0,0 @@ -import math - -from firedrake import * # noqa: F403 - - -def functions(Wave_obj): - """Damping functions for the perfect matched layer for 2D and 3D - - Parameters - ---------- - Wave_obj : obj - Wave object with the parameters of the problem - - Returns - ------- - sigma_x : obj - Firedrake function with the damping function in the x direction - sigma_z : obj - Firedrake function with the damping function in the z direction - sigma_y : obj - Firedrake function with the damping function in the y direction - - """ - - ps = Wave_obj.abc_exponent - cmax = Wave_obj.abc_cmax # maximum acoustic wave velocity - R = Wave_obj.abc_R # theoretical reclection coefficient - pad_length = Wave_obj.mesh_parameters.abc_pad_length # length of the padding - V = Wave_obj.function_space - dimension = Wave_obj.dimension - z = Wave_obj.mesh_z - x = Wave_obj.mesh_x - x1 = 0.0 - x2 = Wave_obj.mesh_parameters.length_x - z2 = -Wave_obj.mesh_parameters.length_z - - bar_sigma = ((3.0 * cmax) / (2.0 * pad_length)) * math.log10(1.0 / R) - aux1 = Function(V) - aux2 = Function(V) - - # Sigma X - sigma_max_x = bar_sigma # Max damping - aux1.interpolate( - conditional( - And((x >= x1 - pad_length), x < x1), - ((abs(x - x1) ** (ps)) / (pad_length ** (ps))) * sigma_max_x, - 0.0, - ) - ) - aux2.interpolate( - conditional( - And(x > x2, (x <= x2 + pad_length)), - ((abs(x - x2) ** (ps)) / (pad_length ** (ps))) * sigma_max_x, - 0.0, - ) - ) - sigma_x = Function(V, name="sigma_x").interpolate(aux1 + aux2) - - # Sigma Z - tol_z = 1.000001 - sigma_max_z = bar_sigma # Max damping - aux1.interpolate( - conditional( - And(z < z2, (z >= z2 - tol_z * pad_length)), - ((abs(z - z2) ** (ps)) / (pad_length ** (ps))) * sigma_max_z, - 0.0, - ) - ) - - sigma_z = Function(V, name="sigma_z").interpolate(aux1) - - if dimension == 2: - return (sigma_x, sigma_z) - - elif dimension == 3: - # Sigma Y - sigma_max_y = bar_sigma # Max damping - y = Wave_obj.mesh_y - y1 = 0.0 - y2 = Wave_obj.length_y - aux1.interpolate( - conditional( - And((y >= y1 - pad_length), y < y1), - ((abs(y - y1) ** (ps)) / (pad_length ** (ps))) * sigma_max_y, - 0.0, - ) - ) - aux2.interpolate( - conditional( - And(y > y2, (y <= y2 + pad_length)), - ((abs(y - y2) ** (ps)) / (pad_length ** (ps))) * sigma_max_y, - 0.0, - ) - ) - sigma_y = Function(V, name="sigma_y").interpolate(aux1 + aux2) - # sgm_y = VTKFile("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]]) - Gamma_2 = as_tensor([[sigma_x - sigma_y, 0.0], [0.0, sigma_y - sigma_x]]) - - return (Gamma_1, Gamma_2) - - -def matrices_3D(sigma_x, sigma_y, sigma_z): - """Damping matrices for a three-dimensional problem""" - Gamma_1 = as_tensor( - [[sigma_x, 0.0, 0.0], [0.0, sigma_y, 0.0], [0.0, 0.0, sigma_z]] - ) - Gamma_2 = as_tensor( - [ - [sigma_x - sigma_y - sigma_z, 0.0, 0.0], - [0.0, sigma_y - sigma_x - sigma_z, 0.0], - [0.0, 0.0, sigma_z - sigma_x - sigma_y], - ] - ) - Gamma_3 = as_tensor( - [ - [sigma_y * sigma_z, 0.0, 0.0], - [0.0, sigma_x * sigma_z, 0.0], - [0.0, 0.0, sigma_x * sigma_y], - ] - ) - - return (Gamma_1, Gamma_2, Gamma_3) diff --git a/spyro/pml/pml_nsnc.py b/spyro/pml/pml_nsnc.py index 02ad7c28..046040fe 100644 --- a/spyro/pml/pml_nsnc.py +++ b/spyro/pml/pml_nsnc.py @@ -1,14 +1,18 @@ -# import firedrake as fire -# import numpy as np +from firedrake import as_tensor, conditional, Function, VTKFile +from numpy import log from ..abc.abc_layer import ABCLayer -from ..utils.error_management import value_parameter_error -from ..utils.typing import LayerShapeType, LayerSizeRefFrequency +from ..io.basicio import parallel_print as pprint +from ..utils.error_management import enum_parameter_error, value_numerical_error +from ..utils.eval_functions_to_ufl import generate_ufl_functions +from ..utils.typing import (BoundaryConditionsType, LayerDampingType, + LayerShapeType, LayerSizeRefFrequency) # Work from Ruben Andres Salas and Alexandre Olender # non-split non-convolutional PML formulation # Formulation based on: # "Efficient PML for the wave equation". Grote and Sim (2010) # "A Modified PML Acoustic Wave Equation". Kim (2019) +# TODO: Add citation class PMLLayer(ABCLayer): @@ -16,12 +20,21 @@ class PMLLayer(ABCLayer): Attributes ---------- - bc_boundary_pml : `str` + abc_pad_length : `float` + Size of the absorbing layer + bc_boundary_pml : `typing.BoundaryConditionsType` Type of boundary condition to apply on the PML boundaries. - Options are "Higdon" or "Sommerfeld" for Non-Reflecting BCs, - or "Dirichlet" or "Neumann" for typical BCs. Default is "Higdon". - pml_mask : `Firedrake.Function` - Mask function to identify the PML domain. + - Options for Non-Reflecting BCs: + 'BoundaryConditionsType.HIGDON' or 'BoundaryConditionsType.SOMMERFELD'. + - Options for typical BCs: + 'BoundaryConditionsType.DIRICHLET' or 'BoundaryConditionsType.NEUMANN' + Default is 'BoundaryConditionsType.HIGDON'. + length_xabc : `float` + Length of the domain in the x-direction with absorbing layer + length_yabc : `float` + Length of the domain in the y-direction with absorbing layer (3D) + length_zabc : `float` + Length of the domain in the z-direction with absorbing layer sigma_max : `float` Maximum damping coefficient within the PML layer. sigma_x : `Firedrake.Function` @@ -30,27 +43,24 @@ class PMLLayer(ABCLayer): Damping profile in the y direction within the PML layer (3D). sigma_z : `Firedrake.Function` Damping profile in the z direction within the PML layer. - where_to_absorb : `tuple` - Boundary ids where absorption is applied. Methods ------- calc_pml_damping() Calculate the maximum damping coefficient for the PML layer. damping_pml_2d() - Build damping matrices for a two-dimensional problem using PML. + Build the damping matrices for a two-dimensional problem using PML. damping_pml_3d() - Build Damping matrices for a three-dimensional problem using PML. + Build the damping matrices for a three-dimensional problem using PML. pml_layer() Set the damping profile within the PML layer. - pml_parameters_boundary_conditions() - Set the boundary conditions for the PML layer. pml_sigma_field() Generate a damping profile for the PML. """ def __init__(self, domain_dim, frequency, f_Nyquist, dimension=2, - quadrilateral=False, func_space_type=None, bc_boundary_pml="Higdon", + quadrilateral=False, func_space_type=None, + bc_boundary_pml=BoundaryConditionsType.NEUMANN, abc_reference_freq=LayerSizeRefFrequency.SOURCE, output_folder=None, comm=None): """ @@ -73,10 +83,13 @@ def __init__(self, domain_dim, frequency, f_Nyquist, dimension=2, func_space_type, `str`, optional Type of function space for the state variable. Options: 'scalar' or 'vector'. Default is `None`. - bc_boundary_pml : `str`, optional + bc_boundary_pml : `typing.BoundaryConditionsType`, optional Type of boundary condition to apply on the PML boundaries. - Options are 'Higdon' or 'Sommerfeld for Non-Reflecting BCs, - or "Dirichlet" or "Neumann" for typical BCs. Default is "Higdon". + - Options for Non-Reflecting BCs: + 'BoundaryConditionsType.HIGDON' or 'BoundaryConditionsType.SOMMERFELD'. + - Options for typical BCs: + 'BoundaryConditionsType.DIRICHLET' or 'BoundaryConditionsType.NEUMANN' + Default is 'BoundaryConditionsType.HIGDON'. abc_reference_freq : `typing.LayerSizeRefFrequency`, optional Reference frequency for sizing the absorbing layer. Options: 'LayerSizeRefFrequency.SOURCE' or 'LayerSizeRefFrequency.BOUNDARY'. @@ -96,228 +109,236 @@ def __init__(self, domain_dim, frequency, f_Nyquist, dimension=2, ABCLayer.__init__(self, domain_dim, frequency, f_Nyquist, dimension=dimension, quadrilateral=quadrilateral, func_space_type=func_space_type, abc_boundary_layer_shape=LayerShapeType.RECTANGULAR, - abc_boundary_layer_type="PML", + abc_boundary_layer_type=LayerDampingType.PML, abc_reference_freq=abc_reference_freq, output_folder=output_folder, comm=comm) # Type of boundary condition to apply on the PML boundaries - self.bc_boundary_pml = bc_boundary_pml - - allowed_bcs = ["Higdon", "Sommerfeld", "Dirichlet", "Neumann"] - if self.bc_boundary_pml not in allowed_bcs: - value_parameter_error('bc_boundary_pml', self.bc_boundary_pml, allowed_bcs) - - # def calc_pml_damping(self, dgr_prof=2, CR_min=1e-8, CR_max=1e-3): - # """Calculate the maximum damping coefficient for the PML layer. - - # Parameters - # ---------- - # dgr_prof : `int`, optional - # Degree of the damping profile within the PML layer - # CR_min : `float`, optional - # Minimum value for the desired reflection coefficient at outer - # boundary of PML layer. Default is 1e-8 - # CR_max : `float`, optional - # Maximum value for the desired reflection coefficient at outer - # boundary of PML layer. Default is 1e-3 - - # Returns - # ------- - # None - # """ - - # # Desired reflection coefficient at outer boundary of PML layer. - # CR = np.clip(self.abc_R, CR_min, CR_max) - - # # Degree of the damping profile within the PML layer. - # dgr_prof = max(1, dgr_prof) - - # # Maximum damping coefficient within the PML layer - # self.sigma_max = 0. if self.abc_get_ref_model else \ - # self.c_max * (dgr_prof + 1.) / (2. * self.abc_pad_length) * np.log(1. / CR) - - # def pml_sigma_field(self, coords, V): - # """Generate a damping profile for the PML. - - # Parameters - # ---------- - # coords : 'ufl.geometry.SpatialCoordinate' - # Domain Coordinates including the absorbing layer - # V : `Firedrake.FunctionSpace` - # Function space for the mask field - - # Returns - # ------- - # None - # """ - - # # Domain dimensions - # Lx, Lz = self.domain_dim[:2] - - # # Domain coordinates - # z, x = coords[0], coords[1] - - # # Conditional value - # val_condz = (z + Lz)**2 - # val_condx1 = x**2 - # val_condx2 = (x - Lx)**2 - - # # Conditional expressions for the profile - # z_sqr = fire.conditional(z < -Lz, val_condz, 0.) - # x_sqr = fire.conditional(x < 0., val_condx1, 0.) + \ - # fire.conditional(x > Lx, val_condx2, 0.) - - # # Quadratic damping profiles - # ref_z = z_sqr / fire.Constant(self.abc_pad_length ** 2) - # ref_x = x_sqr / fire.Constant(self.abc_pad_length ** 2) - # self.sigma_z = fire.Function(V, name='sigma_z [1/s]') - # self.sigma_z.interpolate(self.pml_mask * self.sigma_max * ref_z) - # self.sigma_x = fire.Function(V, name='sigma_x [1/s]') - # self.sigma_x.interpolate(self.pml_mask * self.sigma_max * ref_x) - - # if self.dimension == 3: # 3D - - # # 3D dimension - # Ly = self.domain_dim[2] - # y = coords[2] - - # # Conditional value - # val_condy1 = y**2 - # val_condy2 = (y - Ly)**2 - - # # Conditional expressions for the profile - # y_sqr = fire.conditional(y < 0., val_condy1, 0.) + \ - # fire.conditional(y > Ly, val_condy2, 0.) - - # # Quadratic damping profile - # ref_y = y_sqr / fire.Constant(self.abc_pad_length ** 2) - # self.sigma_y = fire.Function(V, name='sigma_y [1/s]') - # self.sigma_y.interpolate(self.pml_mask * self.sigma_max * ref_y) - - # # Save damping profile - # if hasattr(self, 'path_case_abc'): - # outfile = fire.VTKFile(self.path_case_abc + "sigma_pml.pvd") - # if self.dimension == 2: # 2D - # outfile.write(self.sigma_z, self.sigma_x) - # if self.dimension == 3: # 3D - # outfile.write(self.sigma_z, self.sigma_x, self.sigma_y) - - # def pml_parameters_boundary_conditions(self): - # """Set the boundary conditions for the PML layer. - - # Parameters - # ---------- - # None - - # Returns - # ------- - # None - # """ - - # # ToDo: Integrate boundary mapping - # # Tuple of boundary ids for NRBC - # bnds = [self.absorb_top, self.absorb_bottom, - # self.absorb_right, self.absorb_left] - # if self.dimension == 3: - # bnds.extend([self.absorb_front, self.absorb_back]) - - # self.where_to_absorb = tuple(np.where(bnds)[0] + 1) # ds starts at 1 - - # # Apply boundary conditions to the PML boundaries. - # type_bc = self.bc_boundary_pml - # if not (type_bc == "Dirichlet" or type_bc == "Neumann"): - - # sommerfeld_bc = True if self.bc_boundary_pml == \ - # "Sommerfeld" else False - - # # Applying NRBCs on outer boundary layer - # self.nrbc_on_boundary_layer(sommerfeld_bc=sommerfeld_bc) - - # def pml_layer(self): - # """Set the damping profile within the PML layer. - - # Parameters - # ---------- - # None - - # Returns - # ------- - # None - # """ - - # print("\nCreating Damping PML Profile", flush=True) - - # # Compute the maximum damping coefficient - - # self.calc_pml_damping() - - # # Mesh coordinates - # coords = fire.SpatialCoordinate(self.mesh) - - # # Damping mask - # V_mask = fire.FunctionSpace(self.mesh, 'DG', 0) - # self.pml_mask = self.layer_mask_field(coords, V_mask, - # type_marker='mask', - # name_mask='pml_mask') - - # # Save damping mask - # if hasattr(self, 'path_case_abc'): - # outfile = fire.VTKFile(self.path_case_abc + "pml_mask.pvd") - # outfile.write(self.pml_mask) - - # # Damping fields - # self.pml_sigma_field(coords, self.function_space) - - # # Boundary conditions for the PML layer - # self.pml_parameters_boundary_conditions() - - # def damping_pml_2d(self): - # """Build damping matrices for a two-dimensional problem using PML. - - # Parameters - # ---------- - # sigma_z: Firedrake 'Function' - # Damping profile in the z direction - # sigma_x: Firedrake 'Function' - # Damping profile in the x direction - - # Returns - # ------- - # Gamma_1: Firedrake 'TensorFunction' - # First damping matrix - # Gamma_2: Firedrake 'TensorFunction' - # Second damping matrix - # """ - # Gamma_1 = fire.as_tensor([[self.sigma_z, 0.], [0., self.sigma_x]]) - # Gamma_2 = fire.as_tensor([[self.sigma_z - self.sigma_x, 0.], - # [0., self.sigma_x - self.sigma_z]]) - - # return Gamma_1, Gamma_2 - - # def damping_pml_3d(self): - # """Build Damping matrices for a three-dimensional problem using PML. - - # Parameters - # ---------- - # None - - # Returns - # ------- - # Gamma_1: Firedrake 'TensorFunction' - # First damping matrix - # Gamma_2: Firedrake 'TensorFunction' - # Second damping matrix - # Gamma_3: Firedrake 'TensorFunction' - # Third damping matrix - # """ - # Gamma_1 = fire.as_tensor([[self.sigma_z, 0., 0.], - # [0., self.sigma_x, 0.], - # [0., 0., self.sigma_y]]) - # Gamma_2 = fire.as_tensor([[self.sigma_z - self.sigma_x - self.sigma_y, 0., 0.], - # [0., self.sigma_x - self.sigma_z - self.sigma_y, 0.], - # [0., 0., self.sigma_y - self.sigma_x - self.sigma_z]]) - # Gamma_3 = fire.as_tensor([[self.sigma_x * self.sigma_y, 0., 0.], - # [0., self.sigma_z * self.sigma_y, 0.], - # [0., 0., self.sigma_z * self.sigma_x]]) - - # return Gamma_1, Gamma_2, Gamma_3 + self.bc_boundary_pml = enum_parameter_error('bc_boundary_pml', bc_boundary_pml, + BoundaryConditionsType) + + def calc_pml_damping(self, abc_pml_R, abc_pml_cmax, abc_pad_length, + degree_prof=2, CR_min=1e-8, CR_max=1e-3): + """Calculate the maximum damping coefficient for the PML layer. + + Parameters + ---------- + abc_pml_R: float + Theoretical reflection coefficient of the PML layer. The lower the reflection + coefficient value, the larger the size of the layer. + abc_pml_cmax: float + Maximum propagation speed (km/s) in the PML layer. + abc_pad_length : `float` + Size of the absorbing layer + degree_prof : `int`, optional + Degree of the damping profile within the PML layer. Default is 2. + CR_min : `float`, optional + Minimum value for the desired reflection coefficient at outer boundary of PML + layer. Default is 1e-8 assuming an amplitude attenuation of 160 dB. If the + source's amplitude is extremely small, this attenuation might be excessive, + resulting in a large PML layer. See information on amplitude attenuation in: + https://ccrma.stanford.edu/~jos/mdft/Exponentials.html#fig:exponential + https://ccrma.stanford.edu/~jos/mdft/Audio_Decay_Time_T60.html + TODO: Add citation + CR_max : `float`, optional + Maximum value for the desired reflection coefficient at outer boundary of PML + layer. Default is 1e-3 assuming an amplitude attenuation of 60 dB. + + Returns + ------- + None + """ + + # Desired reflection coefficient at outer boundary of PML layer. + CR = value_numerical_error("abc_pml_R", abc_pml_R, float_num=True, + lower_bound=CR_min, upper_bound=CR_max, + include_lower_bound=True, include_upper_bound=True) + + # Degree of the damping profile within the PML layer. + degree_prof = value_numerical_error("degree_prof", degree_prof, integer_num=True, + lower_bound=1, include_lower_bound=True) + + # Maximum damping coefficient within the PML layer + self.sigma_max = ((degree_prof + 1.) / (2. * abc_pad_length)) * log(1. / CR) + self.sigma_max *= abc_pml_cmax + + def pml_sigma_field(self, Wave_object, ufl_coordinates_pml, V, + degree_prof=2, save_file=True): + """Generate a damping profile for the PML. + + Parameters + ---------- + Wave_object : `acoustic_wave.AcousticWave` + An instance of the :class:`~spyro.solvers.acoustic_wave.AcousticWave`. + ufl_coordinates_pml : `ufl.geometry.SpatialCoordinate` + Domain Coordinates including the absorbing layer. + V : `Firedrake.FunctionSpace` + Function space for the damping field + degree_prof : `int`, optional + Degree of the damping profile within the PML layer. Default is 2. + save_file : `bool`, optional + If `True`, save the mesh with absorbing layer in a .pvd file. + Default is `True`. + + Returns + ------- + None + """ + + # Domain dimensions + length_z, length_x = self.domain_dim[:2] + + # UFL coordinates + z, x = ufl_coordinates_pml[0], ufl_coordinates_pml[1] + + # Conditional expression + condz = z < -length_z + condx1 = x < 0. + condx2 = x > length_x + + # Conditional value + exprz = f"sqrt((z + {length_z})**2)" + exprx1 = "(sqrt(x))**2" + exprx2 = f"sqrt((x - {length_x})**2)" + valz = generate_ufl_functions(Wave_object.mesh, exprz, self.dimension) + valx1 = generate_ufl_functions(Wave_object.mesh, exprx1, self.dimension) + valx2 = generate_ufl_functions(Wave_object.mesh, exprx2, self.dimension) + + # Conditional expressions for the profile + z_pd = conditional(condz, valz, 0.) + x_pd = conditional(condx1, valx1, 0.) + conditional(condx2, valx2, 0.) + + # Damping profiles + ref_z = (z_pd / self.abc_pad_length) ** degree_prof + ref_x = (x_pd / self.abc_pad_length) ** degree_prof + self.sigma_z = Function(V, name='sigma_z[1/s]') + self.sigma_z.interpolate(self.sigma_max * ref_z) + self.sigma_x = Function(V, name='sigma_x[1/s]') + self.sigma_x.interpolate(self.sigma_max * ref_x) + + if self.dimension == 3: # 3D + + # 3D dimension + length_y = self.domain_dim[2] + y = ufl_coordinates_pml[2] + + # Conditional expression + condy1 = y < 0. + condy2 = y > length_y + + # Conditional value + expry1 = "(sqrt(y))**2" + expry2 = f"sqrt((y - {length_y})**2)" + valy1 = generate_ufl_functions(Wave_object.mesh, expry1, self.dimension) + valy2 = generate_ufl_functions(Wave_object.mesh, expry2, self.dimension) + + # Conditional expressions for the mask + y_pd = conditional(condy1, valy1, 0.) + conditional(condy2, valy2, 0.) + + # Damping profile + ref_y = (y_pd / self.abc_pad_length) ** degree_prof + self.sigma_y = Function(V, name='sigma_y[1/s]') + self.sigma_y.interpolate(self.sigma_max * ref_y) + + # Save damping profile + if not Wave_object.abc_get_ref_model and save_file: + outfile = VTKFile(self.path_case_abc + "sigma_pml.pvd") + if self.dimension == 2: # 2D + outfile.write(self.sigma_z, self.sigma_x) + if self.dimension == 3: # 3D + outfile.write(self.sigma_z, self.sigma_x, self.sigma_y) + + def pml_layer(self, Wave_object, save_file=True): + """Set the damping profile within the PML layer. + + Parameters + ---------- + Wave_object : `acoustic_wave.AcousticWave` + An instance of the :class:`~spyro.solvers.acoustic_wave.AcousticWave`. + save_file : `bool`, optional + If `True`, save the mesh with absorbing layer in a .pvd file. + Default is `True`. + + Returns + ------- + None + """ + + pprint("\nCreating Damping PML Profile", comm=self.comm) + + # New geometry with layer if pad_length is provided by the user. + if Wave_object.abc_user_pad_len: + self.abc_pad_length = Wave_object.abc_pad_length + self.abc_new_geometry() + + # Compute the maximum damping coefficient + abc_pml_cmax = Wave_object.abc_pml_cmax if Wave_object.abc_user_pml_cmax \ + else Wave_object.c_bnd_max + + self.calc_pml_damping(Wave_object.abc_pml_R, abc_pml_cmax, + Wave_object.abc_pad_length, + degree_prof=Wave_object.abc_pml_exponent) + + # Mesh coordinates including the absorbing layer + domain_layer = Wave_object.layer_ops.abc_domain_dimensions(full_hyp=False) + ufl_coordinates_pml = \ + Wave_object.mesh_ops.get_spatial_coordinates_abc(Wave_object.mesh, + domain_layer) + + # Damping fields + self.pml_sigma_field(Wave_object, ufl_coordinates_pml, Wave_object.function_space, + degree_prof=Wave_object.abc_pml_exponent) + + # Non-Reflective BCs for the PML layer if prescribed + self.nrbc_on_boundary_layer(Wave_object, self.bc_boundary_pml, save_file=save_file) + + def damping_pml_2d(self): + """Build the damping matrices for a two-dimensional problem using PML. + + Parameters + ---------- + sigma_z: Firedrake 'Function' + Damping profile in the z direction + sigma_x: Firedrake 'Function' + Damping profile in the x direction + + Returns + ------- + Gamma_1: Firedrake 'TensorFunction' + First damping matrix + Gamma_2: Firedrake 'TensorFunction' + Second damping matrix + """ + Gamma_1 = as_tensor([[self.sigma_z, 0.], [0., self.sigma_x]]) + Gamma_2 = as_tensor([[self.sigma_z - self.sigma_x, 0.], + [0., self.sigma_x - self.sigma_z]]) + + return Gamma_1, Gamma_2 + + def damping_pml_3d(self): + """Build the damping matrices for a three-dimensional problem using PML. + + Parameters + ---------- + None + + Returns + ------- + Gamma_1: Firedrake 'TensorFunction' + First damping matrix + Gamma_2: Firedrake 'TensorFunction' + Second damping matrix + Gamma_3: Firedrake 'TensorFunction' + Third damping matrix + """ + Gamma_1 = as_tensor([[self.sigma_z, 0., 0.], + [0., self.sigma_x, 0.], + [0., 0., self.sigma_y]]) + Gamma_2 = as_tensor([[self.sigma_z - self.sigma_x - self.sigma_y, 0., 0.], + [0., self.sigma_x - self.sigma_z - self.sigma_y, 0.], + [0., 0., self.sigma_y - self.sigma_x - self.sigma_z]]) + Gamma_3 = as_tensor([[self.sigma_x * self.sigma_y, 0., 0.], + [0., self.sigma_z * self.sigma_y, 0.], + [0., 0., self.sigma_z * self.sigma_x]]) + + return Gamma_1, Gamma_2, Gamma_3 diff --git a/spyro/solvers/acoustic_solver_construction_no_pml.py b/spyro/solvers/acoustic_solver_construction_no_pml.py index 8eee1887..9fcc6677 100644 --- a/spyro/solvers/acoustic_solver_construction_no_pml.py +++ b/spyro/solvers/acoustic_solver_construction_no_pml.py @@ -1,5 +1,6 @@ import firedrake as fire -from firedrake import ds, dx, Constant, dot, grad +from firedrake import ds, dx, dot, grad +from ..utils.typing import LayerDampingType def construct_solver_or_matrix_no_pml(Wave_object): @@ -31,7 +32,7 @@ def construct_solver_or_matrix_no_pml(Wave_object): # ------------------------------------------------------- m1 = ( (1 / (Wave_object.c * Wave_object.c)) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * ((u - 2.0 * u_n + u_nm1) / dt**2) * v * dx(**quad_rule) ) @@ -43,12 +44,12 @@ def construct_solver_or_matrix_no_pml(Wave_object): le += - q * v * dx(**quad_rule) if Wave_object.abc_active: - weak_expr_abc = dot((u_n - u_nm1) / Constant(dt), v) + weak_expr_abc = dot((u_n - u_nm1) / dt, v) f_abc = (1 / Wave_object.c) * weak_expr_abc qr_s = Wave_object.surface_quadrature_rule - if Wave_object.abc_boundary_layer_type == "hybrid": + if Wave_object.abc_boundary_layer_type == LayerDampingType.HYBRID: # NRBC le += Wave_object.cosHig * f_abc * ds(**qr_s) @@ -76,6 +77,7 @@ def construct_solver_or_matrix_no_pml(Wave_object): # form = m1 + a - le # Signal for le is + in derivation, see Salas et al (2022) # doi: https://doi.org/10.1016/j.apm.2022.09.014 + # TODO: Add citation form = m1 + a + le Wave_object.rhs = fire.rhs(form) Wave_object.lhs = fire.lhs(form) diff --git a/spyro/solvers/acoustic_solver_construction_with_pml.py b/spyro/solvers/acoustic_solver_construction_with_pml.py index a449f26c..723696d3 100644 --- a/spyro/solvers/acoustic_solver_construction_with_pml.py +++ b/spyro/solvers/acoustic_solver_construction_with_pml.py @@ -1,226 +1,205 @@ -import firedrake as fire -from firedrake import dx, ds, Constant, dot, grad, inner +"""Constructs Firedrake solver for the acosutic wave with a PML.""" -from ..pml import damping +from firedrake import (Cofunction, DirichletBC, div, dot, ds as fire_ds, dx as fire_dx, + Function, grad, inner, lhs, LinearVariationalProblem, + LinearVariationalSolver, rhs, split, TestFunctions, TrialFunctions) +from ..domains.space import create_function_space +from ..utils.typing import BoundaryConditionsType -def construct_solver_or_matrix_with_pml(Wave_object): - """ - Builds solver operators for wave propagator with a PML. Doesn't create mass matrices if - matrix_free option is on, which it is by default. - """ - if Wave_object.dimension == 2: - return construct_solver_or_matrix_with_pml_2d(Wave_object) - elif Wave_object.dimension == 3: - return construct_solver_or_matrix_with_pml_3d(Wave_object) - +# Work from Keith Roberts, Eduardo Moscatelli, +# Ruben Andres Salas and Alexandre Olender +# Formulation based on: +# "Efficient PML for the wave equation". Grote and Sim (2010) +# "A Modified PML Acoustic Wave Equation". Kim (2019) +# TODO: Add citations -def construct_solver_or_matrix_with_pml_2d(Wave_object): +def forms_pml(Wave_object, W, X_n, X_nm1): """ - Builds solver operators for 2D wave propagator with a PML. Doesn't create mass matrices if - matrix_free option is on, which it is by default. + Build the variational form for the wave equation with a PML. + + Parameters + ---------- + Wave_object : `acoustic_wave.AcousticWave` + An instance of the :class:`~spyro.solvers.acoustic_wave.AcousticWave`. + W : `Firedrake.MixedFunctionSpace` + Mixed function space for the wave equation with PML + X_n : `Firedrake.Function` + State variable at time step n + X_nm1 : `Firedrake.Function` + State variable at time step n - 1 + + Returns + ------- + FF : `Firedrake.Form` + Variational form for the wave equation with PML + fix_bnd : `Firedrake.DirichletBC` + Dirichlet boundary conditions applied to the PML boundaries + + Notes + ----- + In the UFL forms, however, reference ``X_n``/``X_nm1`` directly via ``fire.split``. + The pyadjoint tape tracks ``X_n.assign(X_np1)`` between time steps; the subfunction + Functions are separately-tracked tape variables that are never explicitly written to + in the tape, so on replay their tape values stay at the initial (zero) state, making + ``J`` constant w.r.t. the control. Using ``fire.split(X_n)`` keeps the form dependent + on ``X_n`` itself, which is updated correctly on replay. + This is what makes the PML Taylor test converge. """ + + # Simulation parameters for PML formulation dt = Wave_object.dt c = Wave_object.c + c_sqr_inv = 1. / (c * c) + q_rule = Wave_object.quadrature_rule + dx = fire_dx(**q_rule) if q_rule else fire_dx + Wave_object.layer_ops.pml_layer(Wave_object) - V = Wave_object.function_space - Z = Wave_object.vector_function_space - W = V * Z - Wave_object.mixed_function_space = W - dxlump = dx(**Wave_object.quadrature_rule) - - u, pp = fire.TrialFunctions(W) - v, qq = fire.TestFunctions(W) - - X_np1 = fire.Function(W) - X_n = fire.Function(W) - X_nm1 = fire.Function(W) - - # Keep Function-typed subfunction views for non-form usage (in-place - # ``.assign``, ``.dat`` access for receiver output, etc.). They share - # storage with ``X_n``/``X_nm1`` so writes to the mixed Functions are - # visible through these views. - u_n_func, pp_n_func = X_n.subfunctions - u_nm1_func, _ = X_nm1.subfunctions - - # In the UFL forms, however, reference ``X_n``/``X_nm1`` directly via - # ``fire.split``. The pyadjoint tape tracks ``X_n.assign(X_np1)`` between - # time steps; the subfunction Functions are separately-tracked tape - # variables that are never explicitly written to in the tape, so on - # replay their tape values stay at the initial (zero) state, making - # ``J`` constant w.r.t. the control. Using ``fire.split(X_n)`` keeps the - # form dependent on ``X_n`` itself, which is updated correctly on - # replay. This is what makes the PML Taylor test converge. - u_n, pp_n = fire.split(X_n) - u_nm1, _ = fire.split(X_nm1) - - Wave_object.u_n = u_n_func - Wave_object.X_np1 = X_np1 - Wave_object.X_n = X_n - Wave_object.X_nm1 = X_nm1 + # Trial and test functions, and state variables + if Wave_object.dimension == 2: + u, pp = TrialFunctions(W) + v, qq = TestFunctions(W) + u_n, pp_n = split(X_n) + u_nm1, _ = split(X_nm1) - sigma_x, sigma_z = damping.functions(Wave_object) - Gamma_1, Gamma_2 = damping.matrices_2D(sigma_z, sigma_x) - pml1 = (sigma_x + sigma_z) * ((u - u_nm1) / Constant(2.0 * dt)) * v * dxlump + elif Wave_object.dimension == 3: + u, psi, pp = TrialFunctions(W) + v, phi, qq = TestFunctions(W) + u_n, psi_n, pp_n = split(X_n) + u_nm1, psi_nm1, _ = split(X_nm1) + + # Acoustic form + m1 = c_sqr_inv * ((u - 2. * u_n + u_nm1) / dt**2) * v + a = dot(grad(u_n), grad(v)) + FF = (m1 + a) * dx + + # Common PML forms (Backward difference for first time derivative) + pml3 = -div(pp_n) * v + # ------------------------------------------------------- + mm1 = dot((pp - pp_n) / dt, qq) # 1st order in error + + # Surfaces to apply boundary conditions (NRBCs or Traditional BCs) + bc_surf = tuple([non_free_surf for non_free_surf, status in + Wave_object.mesh_parameters.boundary_ids_map.items() if status]) + abc_type = Wave_object.layer_ops.bc_boundary_pml + + if not Wave_object.abc_get_ref_model: + + # Damping profiles and matrices + sigma_x, sigma_z = Wave_object.layer_ops.sigma_x, Wave_object.layer_ops.sigma_z + if Wave_object.dimension == 2: + Gamma_1, Gamma_2 = Wave_object.layer_ops.damping_pml_2d() + + elif Wave_object.dimension == 3: + sigma_y = Wave_object.layer_ops.sigma_y + Gamma_1, Gamma_2, Gamma_3 = Wave_object.layer_ops.damping_pml_3d() + + # PML forms (Backward difference for first time derivative) + mm2 = inner(dot(Gamma_1, pp_n), qq) + dd1 = inner(dot(Gamma_2, grad(u_n)), qq) + # ------------------------------------------------------- + if Wave_object.dimension == 2: + pml1 = (sigma_z + sigma_x) * ((u_n - u_nm1) / dt) * v + pml2 = sigma_z * sigma_x * u_n * v + + elif Wave_object.dimension == 3: + pml1 = (sigma_z + sigma_x + sigma_y) * ((u_n - u_nm1) / dt) * v + pml2 = (sigma_z * sigma_x + sigma_x * sigma_y + sigma_z * sigma_y) * u_n * v + pml4 = (sigma_z * sigma_x * sigma_y) * psi_n * v + FF += c_sqr_inv * pml4 * dx + # ------------------------------------------------------- + dd2 = -inner(dot(Gamma_3, grad(psi_n)), qq) + FF += c * c * dd2 * dx + # ------------------------------------------------------- + mmm1 = ((psi - psi_n) / dt) * phi + uuu1 = -u_n * phi + FF += (mmm1 + uuu1) * dx + + # Adding common PML forms to the variational form + FF += c_sqr_inv * (pml1 + pml2 + pml3) * dx + # ------------------------------------------------------- + FF += (mm1 + mm2 + c * c * dd1) * dx + + # exterior_markers = set(Wave_object.mesh.exterior_facets.unique_markers) + # print("Available boundary markers:", exterior_markers) + + # Apply NRBCs (Higdon or Sommerfeld) at PML boundaries + if abc_type in [BoundaryConditionsType.HIGDON, BoundaryConditionsType.SOMMERFELD]: + ds = fire_ds(bc_surf, **q_rule) if q_rule else fire_ds(bc_surf) + f_abc = c * ((u_n - u_nm1) / dt) * v + le = Wave_object.layer_ops.cosHig * f_abc * ds + FF += le + + # Dirichlet BCs for PML model or Neumann BCs for the reference model + get_ref_model_negat = not Wave_object.abc_get_ref_model + fix_bnd = DirichletBC(W.sub(0), 0., bc_surf) \ + if abc_type == BoundaryConditionsType.DIRICHLET and get_ref_model_negat else None + + return FF, fix_bnd - # typical CG FEM in 2d/3d - # ------------------------------------------------------- - m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dxlump - a = c * c * 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 - # firedrake.utility_meshes.RectangleMesh/TensorRectangleMesh: - # https://www.firedrakeproject.org/_modules/firedrake/utility_meshes.html#RectangleMesh - # Spyro builds RectangleMesh(nz, nx, length_z, length_x) and negates - # 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 - nf = ( - abc_expr * ds(2, **qr_s) - + abc_expr * ds(3, **qr_s) - + abc_expr * ds(4, **qr_s) - ) - if Wave_object.absorb_top: - nf += abc_expr * ds(1, **qr_s) - - FF = m1 + a + nf - - pml2 = sigma_x * sigma_z * u_n * v * dxlump - pml3 = inner(pp_n, grad(v)) * dxlump - FF += pml1 + pml2 + pml3 - # ------------------------------------------------------- - 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 - FF += mm1 + mm2 + dd +def construct_solver_or_matrix_with_pml(Wave_object): + """Build solver operators for wave propagator with a PML. - Wave_object.lhs = fire.lhs(FF) - Wave_object.rhs = fire.rhs(FF) - Wave_object.source_function = fire.Cofunction(W.dual()) + Doesn't create mass matrices if matrix_free option is on, which it is by default. - lin_var = fire.LinearVariationalProblem( - Wave_object.lhs, Wave_object.rhs + Wave_object.source_function, - X_np1, constant_jacobian=True) - solver_parameters = dict(Wave_object.solver_parameters) - solver_parameters["mat_type"] = "matfree" - Wave_object.solver = fire.LinearVariationalSolver( - lin_var, solver_parameters=solver_parameters, - ) + Parameters + ---------- + Wave_object : `acoustic_wave.AcousticWave` + An instance of the :class:`~spyro.solvers.acoustic_wave.AcousticWave`. + Returns + ------- + None -def construct_solver_or_matrix_with_pml_3d(Wave_object): - """ - Builds solver operators for 3D wave propagator with a PML. Doesn't create mass matrices if - matrix_free option is on, which it is by default. + Notes + ----- + Keep Function-typed subfunction views for non-form usage (in-place ``.assign``, + ``.dat`` access for receiver output, etc.). They share storage with ``X_n``/``X_nm1`` + so writes to the mixed Functions are visible through these views. """ - dt = Wave_object.dt - c = Wave_object.c + # Build mixed function space V = Wave_object.function_space - Z = Wave_object.vector_function_space - W = V * V * Z + Z = create_function_space(Wave_object.mesh, V.ufl_element(), + dim=Wave_object.dimension) + Wave_object.vector_function_space = Z + if Wave_object.dimension == 2: + W = V * Z + elif Wave_object.dimension == 3: + W = V * V * Z Wave_object.mixed_function_space = W - dxlump = dx(**Wave_object.quadrature_rule) - - u, psi, pp = fire.TrialFunctions(W) - v, phi, qq = fire.TestFunctions(W) - X_np1 = fire.Function(W) - X_n = fire.Function(W) - X_nm1 = fire.Function(W) + # State variables + X_np1 = Function(W) + X_n = Function(W) + X_nm1 = Function(W) - # See 2D variant: keep subfunction Functions for non-form usage but build - # the UFL forms against ``fire.split`` so the adjoint tape tracks the - # time-stepping updates of ``X_n`` / ``X_nm1`` correctly. - u_n_func, psi_n_func, pp_n_func = X_n.subfunctions - u_nm1_func, psi_nm1_func, _ = X_nm1.subfunctions + if Wave_object.dimension == 2: + u_n_func, pp_n_func = X_n.subfunctions + u_nm1_func, _ = X_nm1.subfunctions - u_n, psi_n, pp_n = fire.split(X_n) - u_nm1, psi_nm1, _ = fire.split(X_nm1) + elif Wave_object.dimension == 3: + u_n_func, psi_n_func, pp_n_func = X_n.subfunctions + u_nm1_func, psi_nm1_func, _ = X_nm1.subfunctions Wave_object.u_n = u_n_func Wave_object.X_np1 = X_np1 Wave_object.X_n = X_n Wave_object.X_nm1 = X_nm1 - sigma_x, sigma_y, sigma_z = damping.functions(Wave_object) - Gamma_1, Gamma_2, Gamma_3 = damping.matrices_3D(sigma_x, sigma_y, sigma_z) - - pml1 = ( - (sigma_x + sigma_y + sigma_z) - * ((u - u_nm1) / Constant(2.0 * dt)) - * v - * dxlump - ) + # Build variational forms + FF, fix_bnd = forms_pml(Wave_object, W, X_n, X_nm1) + Wave_object.lhs = lhs(FF) + Wave_object.rhs = rhs(FF) + Wave_object.source_function = Cofunction(W.dual()) + Wave_object.B = Cofunction(W.dual()) - pml2 = ( - (sigma_x * sigma_y + sigma_x * sigma_z + sigma_y * sigma_z) - * u_n - * v - * dxlump - ) - - pml3 = (sigma_x * sigma_y * sigma_z) * psi_n * v * dxlump - pml4 = inner(pp_n, grad(v)) * dxlump - - # typical CG FEM in 2d/3d - - # ------------------------------------------------------- - m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dxlump - a = c * c * 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 - # firedrake.utility_meshes.BoxMesh/TensorBoxMesh: - # https://www.firedrakeproject.org/_modules/firedrake/utility_meshes.html#BoxMesh - # Spyro builds BoxMesh(nz, nx, ny, length_z, length_x, length_y) and - # 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 - nf = ( - abc_expr * ds(2, **qr_s) - + abc_expr * ds(3, **qr_s) - + abc_expr * ds(4, **qr_s) - + abc_expr * ds(5, **qr_s) - + abc_expr * ds(6, **qr_s) - ) - if Wave_object.absorb_top: - nf += abc_expr * ds(1, **qr_s) - - FF = m1 + a + nf - - B = fire.Cofunction(W.dual()) - - FF += pml1 + pml2 + pml3 + pml4 - # ------------------------------------------------------- - 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 - FF += mm1 + mm2 + dd1 + dd2 - - mmm1 = (dot((psi - psi_n), phi) / Constant(dt)) * dxlump - uuu1 = (-u_n * phi) * dxlump - FF += mmm1 + uuu1 - - lhs_ = fire.lhs(FF) - rhs_ = fire.rhs(FF) - - source_function = fire.Cofunction(W.dual()) - Wave_object.source_function = source_function - - lin_var = fire.LinearVariationalProblem(lhs_, rhs_ + source_function, X_np1, constant_jacobian=True) + # Build solver + lin_var = LinearVariationalProblem( + Wave_object.lhs, Wave_object.rhs + Wave_object.source_function, + X_np1, bcs=fix_bnd, constant_jacobian=True) solver_parameters = dict(Wave_object.solver_parameters) solver_parameters["mat_type"] = "matfree" - solver = fire.LinearVariationalSolver( - lin_var, solver_parameters=solver_parameters, - ) - Wave_object.solver = solver - Wave_object.rhs = rhs_ - Wave_object.B = B - - return + Wave_object.solver = \ + LinearVariationalSolver(lin_var, solver_parameters=solver_parameters) diff --git a/spyro/solvers/acoustic_wave.py b/spyro/solvers/acoustic_wave.py index e45826c8..dae9913b 100644 --- a/spyro/solvers/acoustic_wave.py +++ b/spyro/solvers/acoustic_wave.py @@ -15,7 +15,7 @@ ) from ..domains.space import create_function_space from ..utils.typing import ( - AdjointType, RieszMapType, override, WaveType, + AdjointType, RieszMapType, override, WaveType, LayerDampingType, ) from ..utils import write_hdf5_velocity_model from .functionals import acoustic_energy @@ -71,12 +71,11 @@ def matrix_building(self): self.solver = None self.rhs = None self.B = None - if abc_type is None or abc_type == "local" or abc_type == "hybrid": + + if abc_type in [LayerDampingType.LOCAL, LayerDampingType.HYBRID, + LayerDampingType.NOABCS]: construct_solver_or_matrix_no_pml(self) - elif abc_type == "PML": - V = self.function_space - Z = fire.VectorFunctionSpace(V.ufl_domain(), V.ufl_element()) - self.vector_function_space = Z + elif abc_type == LayerDampingType.PML: self.X_np1 = None self.X_n = None self.X_nm1 = None @@ -200,7 +199,7 @@ def _automated_adjoint_gradient(self, riesz_map=RieszMapType.L2): ) def reset_pressure(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: self.X_n.assign(0.0) self.X_nm1.assign(0.0) else: @@ -244,49 +243,49 @@ def _initialize_model_parameters(self): @override def _set_vstate(self, vstate): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: self.X_n.assign(vstate) else: self.u_n.assign(vstate) @override def _get_vstate(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.X_n else: return self.u_n @override def _set_prev_vstate(self, vstate): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: self.X_nm1.assign(vstate) else: self.u_nm1.assign(vstate) @override def _get_prev_vstate(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.X_nm1 else: return self.u_nm1 @override def _set_next_vstate(self, vstate): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: self.X_np1.assign(vstate) else: self.u_np1.assign(vstate) @override def _get_next_vstate(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.X_np1 else: return self.u_np1 @override def get_forward_solution_receivers(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: data_with_halos = self.X_n.dat.data_ro_with_halos[0][:] else: data_with_halos = self.u_n.dat.data_ro_with_halos[:] @@ -313,12 +312,12 @@ def get_function(self, state: fire.Function = None) -> fire.Function: The scalar wave field corresponding to the specified `state` or the time step ``n``. """ if state is None: - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.X_n.sub(0) else: return self.u_n else: - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return state.sub(0) else: return state @@ -340,7 +339,7 @@ def _create_function_space(self): @override def rhs_no_pml(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.B.sub(0) else: return self.B @@ -349,7 +348,7 @@ def rhs_no_pml_source(self): """Return the source cofunction added to the variational right-hand side. """ - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return self.source_function.sub(0) else: return self.source_function @@ -366,7 +365,7 @@ def pressure_for_receivers(self): For the non-mixed case ``self.u_n`` is already the form coefficient Function, which is updated directly by ``u_n.assign(...)``. """ - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: return fire.split(self.X_n)[0] return self.u_n diff --git a/spyro/solvers/backward_time_integration.py b/spyro/solvers/backward_time_integration.py index c425b34f..d8331a6e 100644 --- a/spyro/solvers/backward_time_integration.py +++ b/spyro/solvers/backward_time_integration.py @@ -3,6 +3,7 @@ from .wave import Wave from ..io.basicio import parallel_print from ..receivers.Receivers import Receivers +from ..utils.typing import LayerDampingType def backward_wave_propagator(wave_obj: Wave, dt: float = None) -> fire.Function: @@ -75,7 +76,7 @@ def backward_wave_propagator(wave_obj: Wave, dt: float = None) -> fire.Function: # Assign the adjoint solution at the step `np1` to `uadj`. uadj.assign(wave_obj.get_function(state=wave_obj.next_vstate)) - if wave_obj.abc_boundary_layer_type == "PML": + if wave_obj.abc_boundary_layer_type == LayerDampingType.PML: # Pop to keep the list in sync, but use the element one # step behind so that u_fwd and u_adj are at the same # physical time (usol[k] = u^{k+1}; we need u^k). @@ -89,7 +90,7 @@ def backward_wave_propagator(wave_obj: Wave, dt: float = None) -> fire.Function: grad_solver.solve() _trapezoidal_gradient_integration(dJ, gradi, step, nt) - if wave_obj.abc_boundary_layer_type == "PML": + if wave_obj.abc_boundary_layer_type == LayerDampingType.PML: wave_obj.X_nm1.assign(wave_obj.X_n) wave_obj.X_n.assign(wave_obj.X_np1) else: @@ -168,20 +169,30 @@ def _build_gradient_solver(wave_obj: Wave, mask_available: bool) -> tuple[ forward_field = fire.Function(V) uadj = fire.Function(V) - if wave_obj.abc_boundary_layer_type == "PML": + if wave_obj.abc_boundary_layer_type == LayerDampingType.PML: # Always exclude PML region from gradient. # This is necessary once the gradient expression is not considering # the PML auxiliary variables. In addition, we are not interested # in the gradient in the PML region. indicator = _pml_interior_indicator(wave_obj) # Compute the gradient only in the physical domain. + + """ + TODO: Refactor the gradient due to new PML formulation + TODO: Add citations + Formulation based on: + "Efficient PML for the wave equation". Grote and Sim (2010) + "A Modified PML Acoustic Wave Equation". Kim (2019) + Acoustic Eq. is modified by dividing by c^2 (see implementation). + The remaining PML Eqs. remanin unchanged. + """ + ffG = ( 2.0 * wave_obj.c * indicator * fire.dot( fire.grad(uadj), fire.grad(forward_field)) * m_v * dx ) - parallel_print( - "Excluding PML region from gradient (mixed space)", wave_obj.comm - ) + raise ValueError("PML gradient calculation temporarily unavailable") + else: ffG = ( -2 * (wave_obj.c) ** (-3) * fire.dot(forward_field, uadj) * m_v * dx diff --git a/spyro/solvers/eikonal/eikonal_eq.py b/spyro/solvers/eikonal/eikonal_eq.py index 4e189648..509d9a6a 100644 --- a/spyro/solvers/eikonal/eikonal_eq.py +++ b/spyro/solvers/eikonal/eikonal_eq.py @@ -120,10 +120,8 @@ def __init__(self, dimension, source_locations, ele_type_eik='consistent', # Finite element type. allowed_ele_types = ['consistent', 'underintegrated'] - if ele_type_eik not in allowed_ele_types: - value_parameter_error('ele_type_eik', ele_type_eik, allowed_ele_types) - else: - self.ele_type_eik = ele_type_eik + self.ele_type_eik = value_parameter_error('ele_type_eik', ele_type_eik, + allowed_ele_types) # Finite element order for the Eikonal analysis self.degree_eik = degree_eik if degree_eik is not None else \ diff --git a/spyro/solvers/elastic_wave/isotropic_wave.py b/spyro/solvers/elastic_wave/isotropic_wave.py index c320de32..762d925e 100644 --- a/spyro/solvers/elastic_wave/isotropic_wave.py +++ b/spyro/solvers/elastic_wave/isotropic_wave.py @@ -7,8 +7,8 @@ from .forms import (isotropic_elastic_without_pml, isotropic_elastic_with_pml) from .functionals import mechanical_energy_form -from ...utils.typing import (ElasticMaterialParameter, - ElasticMaterialParameterization, override) +from ...utils.typing import (ElasticMaterialParameter, ElasticMaterialParameterization, + LayerDampingType, override) from ...domains.space import create_function_space @@ -231,7 +231,7 @@ def _get_next_vstate(self): @override def get_forward_solution_receivers(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: raise NotImplementedError else: data_with_halos = self.u_n.dat.data_ro_with_halos[:] @@ -506,21 +506,20 @@ def matrix_building(self): self.parse_boundary_conditions() self.parse_volumetric_forces() - if self.abc_boundary_layer_type is None or \ - self.abc_boundary_layer_type == "local": + if self.abc_boundary_layer_type in [LayerDampingType.LOCAL, LayerDampingType.NOABCS]: isotropic_elastic_without_pml(self) - elif self.abc_boundary_layer_type == "PML": + elif self.abc_boundary_layer_type == LayerDampingType.PML: isotropic_elastic_with_pml(self) @override def rhs_no_pml(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: raise NotImplementedError else: return self.B def rhs_no_pml_source(self): - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: raise NotImplementedError else: return self.source_function diff --git a/spyro/solvers/modal/modal_sol.py b/spyro/solvers/modal/modal_sol.py index 11127a13..3db65fdd 100644 --- a/spyro/solvers/modal/modal_sol.py +++ b/spyro/solvers/modal/modal_sol.py @@ -140,9 +140,7 @@ def __init__(self, dimension=2, method=None, calc_max_dt=False, comm=None): self.valid_methods.extend(['KRYLOVSCH_CH', 'KRYLOVSCH_CG', 'KRYLOVSCH_GH', 'KRYLOVSCH_GG', 'RAYLEIGH']) - - if self.method not in self.valid_methods: - value_parameter_error('method', method, self.valid_methods) + value_parameter_error('method', self.method, self.valid_methods) pprint(f"Solver Method: {self.method}", comm=self.comm) @@ -448,12 +446,12 @@ def freq_factor_rec(hyper_axes, bc="Neumann"): Fundamental frequency factor for rectangular or prismatic geometry. """ + # Check boundary condition type + value_parameter_error('bc', bc, ["Dirichlet", "Neumann"]) if bc == "Neumann": f_rec = 1. / max(hyper_axes) elif bc == "Dirichlet": f_rec = sum(1. / np.asarray(hyper_axes)**2)**0.5 - else: - value_parameter_error('bc', bc, ["Dirichlet", "Neumann"]) f_rec *= np.pi / 2. @@ -484,8 +482,8 @@ def freq_factor_ell(self, hyper_axes, bc="Neumann", all_axes_equal=False): Fundamental frequency factor for elliptical or ellipsoidal geometry. """ - if bc not in ["Dirichlet", "Neumann"]: - value_parameter_error('bc', bc, ["Dirichlet", "Neumann"]) + # Check boundary condition type + value_parameter_error('bc', bc, ["Dirichlet", "Neumann"]) def MMF(q): """Compute the Modified Mathieiu's Function (MMF) or its derivative. diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index 1e5e29f6..675cefed 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -3,7 +3,7 @@ from . import helpers from .. import utils -from ..utils.typing import FunctionalEvaluationMode, AdjointType +from ..utils.typing import FunctionalEvaluationMode, AdjointType, LayerDampingType def _propagate_forward_central_difference(wave_obj, source_ids): @@ -49,14 +49,14 @@ def _propagate_forward_central_difference(wave_obj, source_ids): # being one at a point and zero elsewhere. source_cof = wave_obj.sources.source_cofunction() - if wave_obj.abc_boundary_layer_type == "PML": + if wave_obj.abc_boundary_layer_type == LayerDampingType.PML: pressure_expr = fire.split(wave_obj.X_n)[0] else: pressure_expr = wave_obj.u_n interpolate_receivers = wave_obj.receivers.receiver_interpolator( pressure_expr) if ( - wave_obj.abc_boundary_layer_type == "PML" + wave_obj.abc_boundary_layer_type == LayerDampingType.PML and wave_obj.source_function is not None ): master_source_W = fire.Cofunction( diff --git a/spyro/solvers/wave.py b/spyro/solvers/wave.py index b8e665c6..2caa72a8 100644 --- a/spyro/solvers/wave.py +++ b/spyro/solvers/wave.py @@ -1,4 +1,5 @@ from abc import abstractmethod, ABCMeta +from numpy import inf import warnings import firedrake as fire @@ -15,7 +16,8 @@ from ..sources.Sources import Sources from .solver_parameters import get_default_parameters_for_method from ..utils import eval_functions_to_ufl -from ..utils.typing import (AdjointType, FunctionalEvaluationMode, +from ..utils.error_management import enum_parameter_error +from ..utils.typing import (AdjointType, FunctionalEvaluationMode, LayerDampingType, LayerShapeType, WaveType) from .modal.modal_sol import Modal_Solver from .automatic_differentiation_solver import AutomatedAdjoint @@ -120,12 +122,13 @@ def __init__(self, dictionary=None, wave_type=WaveType.NONE, comm=None): model_parameters : `Python object` Contains model parameters. """ + super().__init__(dictionary=dictionary, comm=comm) self.initial_velocity_model = None self.gradient_mask_available = False # Setting wave type - self.wave_type = wave_type + self.wave_type = enum_parameter_error("wave_type", wave_type, WaveType) self.function_space = None self.dg0_scalar_function_space = None @@ -185,7 +188,7 @@ def forward_solve(self): if self.function_space is None: self.force_rebuild_function_space() - if self.abc_boundary_layer_type != "hybrid": + if self.abc_boundary_layer_type != LayerDampingType.HYBRID: self._initialize_model_parameters() self.matrix_building() self.wave_propagator() @@ -698,9 +701,9 @@ def layer_manager(self): domain_dim = self.domain_dimensions() # Nyquist frequency - freq_Nyquist = None if self.analysis != "transient" else 1. / (2. * self.dt) + freq_Nyquist = inf if self.analysis != "transient" else 1. / (2. * self.dt) - if self.abc_boundary_layer_type == "PML": + if self.abc_boundary_layer_type == LayerDampingType.PML: from ..pml.pml_nsnc import PMLLayer self.layer_ops = PMLLayer(domain_dim, self.frequency, freq_Nyquist, dimension=self.dimension, @@ -709,7 +712,7 @@ def layer_manager(self): abc_reference_freq=self.abc_reference_freq, output_folder=self.output_folder, comm=self.comm) - if self.abc_boundary_layer_type == "hybrid": + if self.abc_boundary_layer_type == LayerDampingType.HYBRID: from ..habc.habc import HABCLayer self.layer_ops = HABCLayer(domain_dim, self.frequency, freq_Nyquist, self.abc_deg_layer, dimension=self.dimension, @@ -721,7 +724,7 @@ def layer_manager(self): output_folder=self.output_folder, comm=self.comm) # Identifier for the current case study - if self.abc_boundary_layer_type in ["PML", "hybrid"]: + if self.abc_boundary_layer_type in [LayerDampingType.PML, LayerDampingType.HYBRID]: self.case_abc = self.layer_ops.case_abc self.path_save = self.layer_ops.path_save self.path_case_abc = self.layer_ops.path_case_abc diff --git a/spyro/tools/cells_per_wavelength_calculator.py b/spyro/tools/cells_per_wavelength_calculator.py index ee991843..df679566 100644 --- a/spyro/tools/cells_per_wavelength_calculator.py +++ b/spyro/tools/cells_per_wavelength_calculator.py @@ -358,7 +358,7 @@ def find_minimum(self, starting_cpw=None, TOL=None, accuracy=None, savetxt=False ) else: Wave_obj.dt = self.fixed_timestep - print("Maximum dt is ", Wave_obj.dt, flush=True) + print("Maximum dt (seconds) is", Wave_obj.dt, flush=True) t0 = timinglib.time() Wave_obj.forward_solve() diff --git a/spyro/tools/habc_tools.py b/spyro/tools/habc_tools.py index 37532963..3ea49486 100644 --- a/spyro/tools/habc_tools.py +++ b/spyro/tools/habc_tools.py @@ -4,7 +4,8 @@ from firedrake import sqrt as fire_sqrt from numpy import clip, where from ..domains.space import create_function_space -from ..utils.error_management import value_parameter_error +from ..io.basicio import parallel_print as pprint +from ..utils.error_management import value_numerical_error, value_parameter_error from ..utils.eval_functions_to_ufl import generate_ufl_functions from ..tools.version_control import is_firedrake_new @@ -46,8 +47,6 @@ def generate_conditional_value_for_layer(domain_dim, mesh, dimension, # UFL coordinates z, x = ufl_coordinates_habc[0], ufl_coordinates_habc[1] - if dimension == 3: # 3D - y = ufl_coordinates_habc[2] # Conditional expression condz = z < -length_z @@ -71,6 +70,7 @@ def generate_conditional_value_for_layer(domain_dim, mesh, dimension, # 3D dimension length_y = domain_dim[2] + y = ufl_coordinates_habc[2] # Conditional expression condy1 = y < 0. @@ -112,8 +112,8 @@ def layer_mask_field(domain_dim, mesh, dimension, ufl_coordinates_habc, V, Function space for the mask field. damp_par : `tuple`, optional Damping parameters for the absorbing layer. - Structure: (pad_len, eta_crt, aq, bq). - - pad_len : `float` + Structure: (pad_length, eta_crt, aq, bq). + - pad_length : `float` Size of the absorbing layer. - eta_crt : `float` Critical damping coefficient (1/s). @@ -143,26 +143,28 @@ def layer_mask_field(domain_dim, mesh, dimension, ufl_coordinates_habc, V, dimension, ufl_coordinates_habc, type_marker=type_marker) - if type_marker == 'damping': # Damping profile for the absorbing layer + value_parameter_error('type_marker', type_marker, ["damping", "mask"]) - if damp_par is None: + if type_marker == "damping": + + # Damping profile for the absorbing layer + if damp_par is None or len(damp_par) != 4: raise ValueError("Damping parameters must be provided " "when 'type_marker' is 'damping'.") + if not isinstance(damp_par, tuple): + raise TypeError(f"'damp_par' must be a tuple, got {type(damp_par).__name__}.") + # Damping parameters - pad_len, eta_crt, aq, bq = damp_par + pad_length, eta_crt, aq, bq = damp_par - if pad_len <= 0: - raise ValueError(f"Invalid value for 'pad_len': {pad_len}. " - "'pad_len' must be greater than zero " - "when 'type_marker' is 'damping'.") - if eta_crt <= 0: - raise ValueError(f"Invalid value for 'eta_crt': {eta_crt}. " - "'eta_crt' must be greater than zero " - "when 'type_marker' is 'damping'.") + value_numerical_error("pad_length", pad_length, float_num=True, + integer_num=True, lower_bound=0.) + + value_numerical_error("eta_crt", eta_crt, float_num=True, lower_bound=0.) # Reference distance to the original boundary - ref_funct = fire_sqrt(ref_funct) / pad_len + ref_funct = fire_sqrt(ref_funct) / pad_length # Quadratic damping profile ref_funct = aq * ref_funct**2 @@ -170,13 +172,11 @@ def layer_mask_field(domain_dim, mesh, dimension, ufl_coordinates_habc, V, ref_funct += bq * ref_funct ref_funct *= eta_crt - elif type_marker == 'mask': # Mask filter for layer boundary domain + elif type_marker == "mask": + # Mask filter for layer boundary domain ref_funct = conditional(ref_funct > 0, 1., 0.) - else: - value_parameter_error('type_marker', type_marker, ['damping', 'mask']) - layer_mask = Function(V, name=name_mask) layer_mask.assign(assemble(interpolate(ref_funct, V))) @@ -211,10 +211,10 @@ def clipping_coordinates_lay_field(domain_dim, mesh, dimension, Mask for the absorbing layer. """ - print("Clipping Coordinates Inside Layer", flush=True) + pprint("Clipping Coordinates Inside Layer") # Domain dimensions - Lz, Lx = domain_dim[:2] + length_z, length_x = domain_dim[:2] # Vectorial space for auxiliar field of clipped coordinates method_element = "DQ" if quadrilateral else "DG" @@ -224,16 +224,16 @@ def clipping_coordinates_lay_field(domain_dim, mesh, dimension, lay_field = Function(W) lay_field.assign(assemble(interpolate(ufl_coordinates_habc, W))) lay_arr = lay_field.dat.data_with_halos[:] - lay_arr[:, 0] = clip(lay_arr[:, 0], -Lz, 0.) - lay_arr[:, 1] = clip(lay_arr[:, 1], 0., Lx) + lay_arr[:, 0] = clip(lay_arr[:, 0], -length_z, 0.) + lay_arr[:, 1] = clip(lay_arr[:, 1], 0., length_x) if dimension == 3: # 3D # 3D dimension - Ly = domain_dim[2] + length_y = domain_dim[2] # Clipping coordinates - lay_arr[:, 2] = clip(lay_arr[:, 2], 0., Ly) + lay_arr[:, 2] = clip(lay_arr[:, 2], 0., length_y) # Mask function to identify the absorbing layer domain layer_mask = layer_mask_field(domain_dim, mesh, dimension, @@ -317,7 +317,7 @@ def extend_scalar_field_profile(mesh_original, field_to_extend, lay_field, layer Extended scalar field defined on the same function space as `lay_field`. """ - print("Extending Profile Inside Layer", flush=True) + pprint("Extending Profile Inside Layer") # Extracting the nodes from the layer field lay_nodes = lay_field.dat.data_with_halos[:] @@ -327,9 +327,10 @@ def extend_scalar_field_profile(mesh_original, field_to_extend, lay_field, layer pts_to_extend = lay_nodes[ind_nodes] # Set the property of the nearest point on the original boundary + value_parameter_error('method', method, ["point_cloud", "nearest_point"]) if method == "point_cloud": - print(f"Using Cloud Points Method to Extend {name_prop} Profile", flush=True) + pprint(f"Using Cloud Points Method to Extend {name_prop} Profile") vel_to_extend = \ point_cloud_field(mesh_original, pts_to_extend, @@ -337,14 +338,11 @@ def extend_scalar_field_profile(mesh_original, field_to_extend, lay_field, layer elif method == "nearest_point": - print(f"Using Nearest Point Method to Extend {name_prop} Profile", flush=True) + pprint(f"Using Nearest Point Method to Extend {name_prop} Profile") vel_to_extend = field_to_extend.at(pts_to_extend, dont_raise=True) del pts_to_extend - else: - value_parameter_error('method', method, ["point_cloud", "nearest_point"]) - # Velocity profile inside the layer lay_field.dat.data_with_halos[ind_nodes, 0] = vel_to_extend del vel_to_extend, ind_nodes diff --git a/spyro/utils/error_management.py b/spyro/utils/error_management.py index 8e24fa9d..bf496540 100644 --- a/spyro/utils/error_management.py +++ b/spyro/utils/error_management.py @@ -7,31 +7,38 @@ def value_parameter_error(par_name, par_value, valid_values): - """Raise a ValueError with a specific error message. + """Validate parameter value and raise a ValueError with a specific error message. Parameters ---------- par_name : `str` - Name of the parameter that has an invalid value. + Name of the parameter to be validated (used in error messages). par_value : `str`, `int` or `float` - Value of the parameter that is invalid. + Value of the parameter to be validated. valid_values : `list` List of valid values for the parameter. + Returns + ------- + par_value : `str`, `int` or `float` + The validated parameter value. + Raises ------ ValueError If the parameter value is not in the list of valid values. """ - # Error message about the invalid parameter - err_str = f"Invalid {par_name}: '{par_value}'. Please use: " - opt_str = ", ".join([f"'{val}'" for val in valid_values]) - last_comma = opt_str.rfind(',') - opt_str = opt_str[:last_comma] + " or" + opt_str[last_comma + 1:] \ - if len(valid_values) > 1 else opt_str + if par_value not in valid_values: + err_str = f"Invalid {par_name}: '{par_value}'. Please use: " + opt_str = ", ".join([f"'{val}'" for val in valid_values]) + last_comma = opt_str.rfind(',') + opt_str = opt_str[:last_comma] + " or" + opt_str[last_comma + 1:] \ + if len(valid_values) > 1 else opt_str - raise ValueError(err_str + opt_str) + raise ValueError(err_str + opt_str) + + return par_value def mutually_exclusive_parameter_error(par_name_lst, par_value_lst): @@ -68,15 +75,15 @@ def mutually_exclusive_parameter_error(par_name_lst, par_value_lst): raise ValueError(exc_str + err_str + opt_str) -def value_dimension_error(par_names, par_values, expected_dim): - """Raise a ValueError if parameter dimensions mismatch expected dimension. +def value_model_dimension_error(par_names, parameters, expected_dim): + """Raise a ValueError if parameter dimensions mismatch expected model dimension. Parameters ---------- par_names : `tuple` Names of the parameters to check dimensions. - par_values : `tuple` - Values of the parameters to check dimensions. + parameters : `tuple` + Parameters to check dimensions. expected_dim : `int` Expected dimension of the parameters (2 or 3). @@ -87,13 +94,16 @@ def value_dimension_error(par_names, par_values, expected_dim): """ str_reference, str_comparison = par_names - chk_reference, chk_comparison = par_values + par_reference, par_comparison = parameters - dim_err = ( - f"Mismatch in domain dimensions\n" - f"{str_reference} ({chk_reference}), {str_comparison} ({chk_comparison}) " - f"do not match expected model dimension ({expected_dim}D).") - raise ValueError(dim_err) + chk_reference = len(par_reference) + chk_comparison = len(par_comparison) + if expected_dim != chk_reference or expected_dim != chk_comparison: + dim_err = (f"Mismatch in domain dimensions\n" + f"{str_reference} ({chk_reference}), " + f"{str_comparison} ({chk_comparison}) " + f"do not match expected model dimension ({expected_dim}D).") + raise ValueError(dim_err) def clean_inst_num(data_arr): @@ -116,14 +126,14 @@ def clean_inst_num(data_arr): def value_numerical_error(par_name, par_value, float_num=True, integer_num=False, lower_bound=None, upper_bound=None, include_lower_bound=False, include_upper_bound=False): - """Raise a ValueError with a specific error message for numerical parameters. + """Validate numerical parameters and raise a ValueError if invalid. Parameters ---------- par_name : `str` - Name of the parameter that has an invalid value. + Name of the parameter to be validated (used in error messages). par_value : `int` or `float` - Value of the parameter that is invalid. + Value of the parameter to be validated. float_num : `bool`, optional If `True`, the parameter can be a float. Default is `True`. integer_num : `bool`, optional @@ -137,10 +147,15 @@ def value_numerical_error(par_name, par_value, float_num=True, integer_num=False include_upper_bound : `bool`, optional If `True`, the upper bound is included in the valid range. Default is `False`. + Returns + ------- + par_value : `int` or `float` + The validated parameter value. + Raises ------ TypeError - If the parameter value is not of the expected type (float or integer). + If the parameter value is not of the expected type (`float` or `int`). ValueError If the parameter value is outside the specified bounds or the bounds are invalid. """ @@ -185,6 +200,8 @@ def value_numerical_error(par_name, par_value, float_num=True, integer_num=False raise ValueError(f"'{par_name}' must be {bound_str}, got {par_value}.") + return par_value + def enum_parameter_error(par_name, par_value, valid_enum): """Validate and convert an enum parameter, returning the enum instance. @@ -197,9 +214,9 @@ def enum_parameter_error(par_name, par_value, valid_enum): Parameters ---------- par_name : `str` - Name of the parameter being validated (used in error messages). + Name of the parameter to be validated (used in error messages). par_value : `object` - Value of the parameter to validate. Can be an `enum.EnumMeta` or a `str`. + Value of the parameter to be validated. Can be an `enum.EnumMeta` or a `str`. valid_enum : `enum.EnumMeta` Enum class containing the valid values for the parameter. @@ -223,10 +240,37 @@ def enum_parameter_error(par_name, par_value, valid_enum): # Check if string maps to valid enum value if isinstance(par_value, str): valid_values = [enum.value for enum in valid_enum] - if par_value not in valid_values: - value_parameter_error(par_name, par_value, valid_values) + value_parameter_error(par_name, par_value, valid_values) return valid_enum(par_value) # Invalid type - neither enum instance nor string raise TypeError(f"'{par_name}' must be {valid_enum.__name__} or str" f", got {type(par_value).__name__}") + + +def value_string_error(par_name, par_value): + """Validate string parameters and raise a TypeError if invalid. + + Parameters + ---------- + par_name : `str` + Name of the parameter to be validated (used in error messages). + par_value : `str` + Value of the parameter to be validated. + + Returns + ------- + par_value : `str` + The validated parameter value. + + Raises + ------ + TypeError + If the parameter value is not of the expected type (`str`). + """ + + # Checking the parameter type + if not isinstance(par_value, str): + raise TypeError(f"'{par_name}' must be a string, got {type(par_value).__name__}.") + + return par_value diff --git a/spyro/utils/eval_functions_to_ufl.py b/spyro/utils/eval_functions_to_ufl.py index fb593656..f0a72100 100644 --- a/spyro/utils/eval_functions_to_ufl.py +++ b/spyro/utils/eval_functions_to_ufl.py @@ -71,8 +71,7 @@ def generate_ufl_functions(mesh, expression, dimension): """ # Check model dimension - if dimension not in [2, 3]: - value_parameter_error('dimension', dimension, [2, 3]) + value_parameter_error('dimension', dimension, [2, 3]) # Get available functions and variables namespace = available_functions_to_eval(mesh, dimension) diff --git a/spyro/utils/freq_tools.py b/spyro/utils/freq_tools.py index 5608fc60..23eea799 100644 --- a/spyro/utils/freq_tools.py +++ b/spyro/utils/freq_tools.py @@ -1,10 +1,10 @@ """Utilities for calculating the frequency response of a signal.""" -from numpy import abs, concatenate, linspace, zeros +from numpy import abs, hanning, linspace, mean, pad from scipy.fft import fft -def freq_response(signal, f_Nyq, fpad=4, get_dominant_freq=False): +def freq_response(signal, f_Nyq, fpad=0, get_dominant_freq=False): """Calculate the response in frequency domain of a time signal via FFT. Parameters @@ -14,7 +14,7 @@ def freq_response(signal, f_Nyq, fpad=4, get_dominant_freq=False): f_Nyq : `float` Nyquist frequency according to the time step. f_Nyq = 1 / (2 * dt). fpad : `int`, optional - Padding factor for FFT. Default is 4. + Padding factor for FFT. Default is 0, which means no padding. get_dominant_freq : `bool`, optional If `True`, return only the dominant frequency of the spectrum. Default is `False`. @@ -35,14 +35,21 @@ def freq_response(signal, f_Nyq, fpad=4, get_dominant_freq=False): raise ValueError("Nyquist frequency is invalid. " "Cannot compute frequency response.") + # Remove DC offset + signal = signal - mean(signal) + + # Apply window to taper ends to zero + window = hanning(len(signal)) + signal_windowed = signal * window + # Zero padding for increasing smoothing in FFT - signal_with_padding = concatenate([zeros(fpad * len(signal)), signal]) + signal_with_padding = pad(signal_windowed, (0, fpad * len(signal)), 'constant') # Number of sample points N_samples = len(signal_with_padding) # Determine the number of samples of the spectrum - samples_fft = N_samples // 2 + N_samples % 2 + samples_fft = N_samples // 2 + 1 # Calculate the response in frequency domain of the signal (FFT) norm_magnitude = abs(fft(signal_with_padding)[0:samples_fft]) diff --git a/spyro/utils/typing.py b/spyro/utils/typing.py index 4ffc2ec1..eb4c278b 100644 --- a/spyro/utils/typing.py +++ b/spyro/utils/typing.py @@ -134,7 +134,7 @@ class LayerShapeType(Enum): class LayerSizeRefFrequency(Enum): - """Enum for different reference frequencies for sizing the hybrid absorbing layer. + """Enum for different reference frequencies for sizing the absorbing layer. SOURCE: Size based on dominant source frequency. BOUNDARY: Size based on wave frequency at the critical boundary point (Eikonal min.) @@ -151,3 +151,32 @@ class HyperLayerDegreeType(Enum): """ REAL = "real" INTEGER = "integer" + + +class BoundaryConditionsType(Enum): + """Enum for different types of boundary conditions including Non-Reflecting BCs. + + DIRICHLET: Dirichlet boundary condition + NEUMANN: Neumann boundary condition + HIGDON: 1-st order Higdon boundary condition (NRBC type) + SOMMERFELD: Sommerfeld (Radiation) boundary condition (NRBC type) + + """ + DIRICHLET = "Dirichlet" + NEUMANN = "Neumann" + HIGDON = "Higdon" + SOMMERFELD = "Sommerfeld" + + +class LayerDampingType(Enum): + """Enum for different types of damping for the absorbing layer. + + LOCAL: An Non-Reflecting Boundary Condition (NRBC) at the boundary of the domain. + HYBRID: A combination of a sponge layer and an NRBC at the outer layer boundary + PML: Perfectly Matched Layer (PML) damping + NOABCS: No absorbing boundary conditions applied. + """ + LOCAL = "local" + HYBRID = "hybrid" + PML = "PML" + NOABCS = "no_abcs" diff --git a/tests/on_one_core/test_cpw_calc.py b/tests/on_one_core/test_cpw_calc.py index 6a391865..5c3db5f7 100644 --- a/tests/on_one_core/test_cpw_calc.py +++ b/tests/on_one_core/test_cpw_calc.py @@ -77,11 +77,12 @@ def run_test_cpw_calc(FEM_method_to_evaluate, correct_cpw): assert all([test1, test2, test3]) +# @pytest.mark.xfail( +# reason="Waiting for seismicmesh update for compatibility", +# ) +# @pytest.mark.skipif(not is_seismicmesh_installed(), reason="SeismicMesh is not installed") @pytest.mark.slow -@pytest.mark.xfail( - reason="Waiting for seismicmesh update for compatibility", -) -@pytest.mark.skipif(not is_seismicmesh_installed(), reason="SeismicMesh is not installed") +@pytest.mark.skip(reason="PML formulation subject to another PR") def test_cpw_calc_triangles(): method = "mass_lumped_triangle" correct_cpw = 2.3 diff --git a/tests/on_one_core/test_gradient_2d_pml.py b/tests/on_one_core/test_gradient_2d_pml.py index 0c18e7a5..d395c3d6 100644 --- a/tests/on_one_core/test_gradient_2d_pml.py +++ b/tests/on_one_core/test_gradient_2d_pml.py @@ -3,7 +3,7 @@ from copy import deepcopy import firedrake as fire import spyro -from spyro.utils.typing import AdjointType +from spyro.utils.typing import AdjointType, LayerDampingType import pytest @@ -17,7 +17,7 @@ def check_gradient(Wave_obj_guess, dJ, rec_out_exact, Jm, plot=False, tol=3.0): size, = np.shape(dm.dat.data[:]) dm_data = np.random.default_rng(0).random(size) dm.dat.data_wo[:] = dm_data - if Wave_obj_guess.abc_boundary_layer_type == "PML": + if Wave_obj_guess.abc_boundary_layer_type == LayerDampingType.PML: x = Wave_obj_guess.mesh_x z = Wave_obj_guess.mesh_z inside = fire.And( @@ -224,6 +224,7 @@ def test_gradient_pml_auto_adjoint(): @pytest.mark.slow +@pytest.mark.skip(reason="PML formulation subject to another PR") def test_gradient_pml_implemented_adjoint(): test_gradient_implemented_adjoint(PML=True) diff --git a/tests/on_one_core/test_model_parameters.py b/tests/on_one_core/test_model_parameters.py index 0278ff7d..cbea30a8 100644 --- a/tests/on_one_core/test_model_parameters.py +++ b/tests/on_one_core/test_model_parameters.py @@ -84,67 +84,55 @@ def test_method_reader(): "dimension": 2, # dimension } # Trying out different method entries and seeing if all of them work for MLT - test1 = False test_dictionary = deepcopy(test_dictionary0) test_dictionary["options"]["method"] = "MLT" model = Model_parameters(dictionary=test_dictionary) - if model.method == "mass_lumped_triangle": - test1 = True + assert model.method == "mass_lumped_triangle", \ + f"Expected 'mass_lumped_triangle', got '{model.method}'" test_dictionary = deepcopy(test_dictionary0) - test2 = False test_dictionary["options"]["method"] = "KMV" model = Model_parameters(dictionary=test_dictionary) - if model.method == "mass_lumped_triangle": - test2 = True + assert model.method == "mass_lumped_triangle", \ + f"Expected 'mass_lumped_triangle', got '{model.method}'" test_dictionary = deepcopy(test_dictionary0) - test3 = False test_dictionary["options"]["method"] = "mass_lumped_triangle" model = Model_parameters(dictionary=test_dictionary) - if model.method == "mass_lumped_triangle": - test3 = True + assert model.method == "mass_lumped_triangle", \ + f"Expected 'mass_lumped_triangle', got '{model.method}'" # Trying out different method entries for spectral quads test_dictionary = deepcopy(test_dictionary0) - test4 = False test_dictionary["options"]["method"] = "spectral_quadrilateral" model = Model_parameters(dictionary=test_dictionary) - if model.method == "spectral_quadrilateral": - test4 = True + assert model.method == "spectral_quadrilateral", \ + f"Expected 'spectral_quadrilateral', got '{model.method}'" - test5 = False test_dictionary = deepcopy(test_dictionary0) test_dictionary["options"]["method"] = "CG" test_dictionary["options"]["variant"] = "GLL" - try: + with pytest.raises(ValueError, match="Invalid variant: 'GLL'."): model = Model_parameters(dictionary=test_dictionary) - except ValueError: - test5 = True - test6 = False test_dictionary = deepcopy(test_dictionary0) test_dictionary["options"]["method"] = "SEM" model = Model_parameters(dictionary=test_dictionary) - if model.method == "spectral_quadrilateral": - test6 = True + assert model.method == "spectral_quadrilateral", \ + f"Expected 'spectral_quadrilateral', got '{model.method}'" # Trying out some entries for other less used methods - test7 = False test_dictionary = deepcopy(test_dictionary0) test_dictionary["options"]["method"] = "DG_triangle" model = Model_parameters(dictionary=test_dictionary) - if model.method == "DG_triangle": - test7 = True + assert model.method == "DG_triangle", \ + f"Expected 'DG_triangle', got '{model.method}'" - test8 = False test_dictionary = deepcopy(test_dictionary0) test_dictionary["options"]["method"] = "DG_quadrilateral" model = Model_parameters(dictionary=test_dictionary) - if model.method == "DG_quadrilateral": - test8 = True - - assert all([test1, test2, test3, test4, test5, test6, test7, test8]) + assert model.method == "DG_quadrilateral", \ + f"Expected 'DG_quadrilateral', got '{model.method}'" def test_cell_type_reader(): @@ -159,55 +147,44 @@ def test_cell_type_reader(): # Testing lumped cases ct_dictionary0["options"]["variant"] = "lumped" - test1 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "triangle" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "mass_lumped_triangle": - test1 = True + assert model.method == "mass_lumped_triangle", \ + f"Expected 'mass_lumped_triangle', got '{model.method}'" - test2 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "quadrilateral" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "spectral_quadrilateral": - test2 = True + assert model.method == "spectral_quadrilateral", \ + f"Expected 'spectral_quadrilateral', got '{model.method}'" # Testing equispaced cases ct_dictionary0["options"]["variant"] = "equispaced" - test3 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "triangle" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "CG": - test3 = True + assert model.method == "CG", f"Expected 'CG', got '{model.method}'" - test4 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "quadrilateral" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "CG": - test4 = True + assert model.method == "CG", f"Expected 'CG', got '{model.method}'" # Testing DG cases ct_dictionary0["options"]["variant"] = "DG" - test5 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "triangle" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "DG_triangle": - test5 = True + assert model.method == "DG_triangle", f"Expected 'DG_triangle', got '{model.method}'" - test6 = False ct_dictionary = deepcopy(ct_dictionary0) ct_dictionary["options"]["cell_type"] = "quadrilateral" model = Model_parameters(dictionary=ct_dictionary) - if model.method == "DG_quadrilateral": - test6 = True - - assert all([test1, test2, test3, test4, test5, test6]) + assert model.method == "DG_quadrilateral", \ + f"Expected 'DG_quadrilateral', got '{model.method}'" def test_dictionary_conversion(): @@ -341,23 +318,13 @@ def test_dictionary_conversion(): model_from_new = Model_parameters(dictionary=new_dictionary) # checking relevant information from models - same = True - if model_from_new.method != model_from_old.method: - same = False - if model_from_new.initial_time != model_from_old.initial_time: - same = False - if model_from_new.degree != model_from_old.degree: - same = False - if model_from_new.dimension != model_from_old.dimension: - same = False - if model_from_new.dt != model_from_old.dt: - same = False - if model_from_new.final_time != model_from_old.final_time: - same = False - if model_from_new.forward_output_filename != model_from_old.forward_output_filename: - same = False - - assert same + assert model_from_new.method == model_from_old.method + assert model_from_new.initial_time == pytest.approx(model_from_old.initial_time) + assert model_from_new.degree == model_from_old.degree + assert model_from_new.dimension == model_from_old.dimension + assert model_from_new.dt == pytest.approx(model_from_old.dt) + assert model_from_new.final_time == pytest.approx(model_from_old.final_time) + assert model_from_new.forward_output_filename == model_from_old.forward_output_filename def test_time_exception(): # TODO: improve @@ -387,6 +354,17 @@ def test_receiver_exception(): # TODO: improve model = Model_parameters(dictionary=ex_dictionary) # noqa: F841 +def test_analysis_exception(): + ex_dictionary = deepcopy(dictionary) + with pytest.raises(ValueError, match="Invalid analysis: 'None'."): + ex_dictionary["options"] = { + "degree": 4, # p order + "dimension": 2, # dimension + "analysis": None, # Options: transient, modal or eikonal + } + model = Model_parameters(dictionary=ex_dictionary) # noqa: F841 + + if __name__ == "__main__": test_method_reader() test_cell_type_reader() @@ -394,5 +372,6 @@ def test_receiver_exception(): # TODO: improve test_time_exception() test_source_exception() test_receiver_exception() + test_analysis_exception() print("END") diff --git a/tests/on_one_core/test_plots.py b/tests/on_one_core/test_plots.py index e8ae95e6..b68b8c9a 100644 --- a/tests/on_one_core/test_plots.py +++ b/tests/on_one_core/test_plots.py @@ -39,7 +39,7 @@ def test_plot(): "frequency": 8.0, } rectangle_dictionary["time_axis"] = { - "final_time": 2.0, # Final time for event + "final_time": 1.5, # Final time for event } Wave_obj = spyro.examples.Rectangle_acoustic( dictionary=rectangle_dictionary diff --git a/tests/on_one_core/test_time_convergence.py b/tests/on_one_core/test_time_convergence.py index e8d34fcc..ff1c152c 100644 --- a/tests/on_one_core/test_time_convergence.py +++ b/tests/on_one_core/test_time_convergence.py @@ -11,7 +11,7 @@ def error_calc(p_numerical, p_analytical, nt): return div_error_time -def run_forward(dt): +def run_forward(dt, with_pml=False): # dt = float(sys.argv[1]) final_time = 1.0 @@ -41,6 +41,15 @@ def run_forward(dt): "mesh_file": None, "mesh_type": "firedrake_mesh", # options: firedrake_mesh or user_mesh } + if with_pml: + dictionary["absorving_boundary_conditions"] = { + "status": True, + "damping_type": "PML", + "exponent": 2, + "cmax": 4.5, + "R": 1e-6, + "pad_length": 0.25, + } # 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. @@ -86,7 +95,8 @@ def run_forward(dt): @pytest.mark.slow -def test_second_order_time_convergence(): +@pytest.mark.parametrize("with_pml", [False, True]) +def test_second_order_time_convergence(with_pml): """Test that the second order time convergence of the central difference method is achieved""" @@ -105,7 +115,7 @@ def test_second_order_time_convergence(): for i in range(len(dts)): dt = dts[i] - rec_out = run_forward(dt) + rec_out = run_forward(dt, with_pml=with_pml) rec_anal = np.load(analytical_files[i]) time = np.linspace(0.0, 1.0, int(1.0 / dts[i]) + 1) nt = len(time) @@ -115,7 +125,7 @@ def test_second_order_time_convergence(): theory = [t**2 for t in dts] theory = [errors[0] * th / theory[0] for th in theory] - assert math.isclose(np.log(theory[-1]), np.log(errors[-1]), rel_tol=1e-2) + assert math.isclose(np.log(theory[-1]), np.log(errors[-1]), rel_tol=3e-2) if __name__ == "__main__": diff --git a/tests/on_one_core/test_utils.py b/tests/on_one_core/test_utils.py index 6fe9ea55..132310ad 100644 --- a/tests/on_one_core/test_utils.py +++ b/tests/on_one_core/test_utils.py @@ -2,12 +2,29 @@ import scipy as sp import numpy as np import math +from spyro.utils.freq_tools import freq_response def test_butter_lowpast_filter(): + """ + Test the butter_lowpass_filter function from spyro.utils.utils module. + + Notes + ----- + https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.flattop.html#scipy.signal.windows.flattop + Flat top windows are used for taking accurate measurements of signal amplitude in + the frequency domain, with minimal scalloping error from the center of a frequency + bin to its edges, compared to others. This is a 5th-order cosine window, with the + 5 terms optimized to make the main lobe maximally flat. + See D’Antona, G., Ferrero A., “Digital Signal Processing for Measurement Systems”, + Springer Media, 2006, p. 70 DOI:10.1007/0-387-28666-7. + TODO: Add citation + """ Wave_obj = spyro.examples.Rectangle_acoustic( - dictionary={"absorving_boundary_conditions": {"absorb_top": True}} + dictionary={"absorving_boundary_conditions": {"absorb_top": True}, + "mesh": {"h": 0.09}} # Increasing from 0.05 to 0.09 to save time ) + layer_values = [1.5, 2.0, 2.5, 3.0] z_switches = [-0.25, -0.5, -0.75] Wave_obj.multiple_layer_velocity_model(z_switches, layer_values) @@ -19,10 +36,11 @@ def test_butter_lowpast_filter(): fs = 1.0 / Wave_obj.dt - # Checks if frequency with greater power density is close to 5 - (f, S) = sp.signal.periodogram(rec10, fs) + (f, S) = sp.signal.periodogram(rec10, fs, window='flattop', detrend='linear') peak_frequency = f[np.argmax(S)] - test1 = math.isclose(peak_frequency, 5.0, rel_tol=1e-2) + + # Checks if frequency with greater power density is close to 5 + test1 = math.isclose(peak_frequency, Wave_obj.frequency, rel_tol=1e-2) # Checks if the new frequency is lower than the cutoff cutoff_frequency = 3.0 @@ -35,10 +53,19 @@ def test_butter_lowpast_filter(): filtered_peak_frequency = filt_f[np.argmax(filt_S)] test2 = filtered_peak_frequency < cutoff_frequency - print(f"Peak frequency is close to what it is supposed to be: {test1}") - print(f"Filtered peak frequency is lower than cutoff frequency: {test2}") + # Get the dominant frequency of filtered signal by using FFT + freq_Nyquist = fs / 2. + freq_filt_by_FFT = freq_response(filtered_rec10, freq_Nyquist, get_dominant_freq=True) + test3 = freq_filt_by_FFT < cutoff_frequency + + print(f"Peak frequency ({peak_frequency:.2f} Hz) " + f"is close to what it is supposed to be ({Wave_obj.frequency:.2f}): {test1}") + print(f"Filtered peak frequency ({filtered_peak_frequency:.2f}) Hz " + f"is lower than cutoff frequency ({cutoff_frequency:.2f}): {test2}") + print(f"Filtered peak frequency ({freq_filt_by_FFT:.2f}) Hz " + f"is lower than cutoff frequency ({cutoff_frequency:.2f}): {test3}") - assert all([test1, test2]) + assert all([test1, test2, test3]) def test_geometry_creation():