From 0494c2d553bc37922ccd3e9a3dbf659853921808 Mon Sep 17 00:00:00 2001 From: Alexandre Olender Date: Mon, 19 Jun 2023 18:51:15 -0300 Subject: [PATCH] Adding missing documentation --- spyro/domains/quadrature.py | 28 +++++++- spyro/domains/space.py | 17 ++++- spyro/io/io.py | 50 ++++++++++++++- spyro/pml/damping.py | 46 +++++++++++++- spyro/receivers/Receivers.py | 33 +++++++++- spyro/solvers/helpers.py | 75 ++++++++++++++++++++-- spyro/sources/Sources.py | 50 ++++++++++++++- spyro/tools/gradient_test_ad.py | 23 +++++++ spyro/tools/grid_point_calculator.py | 95 ++++++++++++++++++++++++++++ spyro/utils/estimate_timestep.py | 1 + spyro/utils/utils.py | 16 +++++ 11 files changed, 420 insertions(+), 14 deletions(-) diff --git a/spyro/domains/quadrature.py b/spyro/domains/quadrature.py index 0beaf886..06b006f0 100644 --- a/spyro/domains/quadrature.py +++ b/spyro/domains/quadrature.py @@ -73,7 +73,18 @@ def quadrature_rules(V): # Spectral method - Gauss-Lobatto-Legendre rule # 1D def gauss_lobatto_legendre_line_rule(degree): - """Returns GLL quad rule for a given degree in a line""" + """Returns GLL quad rule for a given degree in a line + + Parameters + ---------- + degree : int + degree of the polynomial + + Returns + ------- + result : obj + quadrature rule + """ fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) finat_ps = finat.point_set.GaussLobattoLegendrePointSet @@ -83,7 +94,20 @@ def gauss_lobatto_legendre_line_rule(degree): # 3D def gauss_lobatto_legendre_cube_rule(dimension, degree): - """Returns GLL quad rule for a given degree in a multidimensional space""" + """Returns GLL quad rule for a given degree in a multidimensional space + + Parameters + ---------- + dimension : int + dimension of the space + degree : int + degree of the polynomial + + Returns + ------- + result : obj + quadrature rule + """ make_tensor_rule = finat.quadrature.TensorProductQuadratureRule result = gauss_lobatto_legendre_line_rule(degree) for _ in range(1, dimension): diff --git a/spyro/domains/space.py b/spyro/domains/space.py index b2331609..6ba2640f 100644 --- a/spyro/domains/space.py +++ b/spyro/domains/space.py @@ -4,7 +4,22 @@ def FE_method(mesh, method, degree): """Define the finite element method: Space discretization - Continuous - or Discontinuous Galerkin methods""" + or Discontinuous Galerkin methods + + Parameters + ---------- + mesh : obj + Firedrake mesh + method : str + Finite element method + degree : int + Degree of the finite element method + + Returns + ------- + element : obj + Firedrake finite element + """ cell_geometry = mesh.ufl_cell() if method == "CG" or method == "spectral": # CG - Continuous Galerkin diff --git a/spyro/io/io.py b/spyro/io/io.py index 218ffd24..166e6f17 100644 --- a/spyro/io/io.py +++ b/spyro/io/io.py @@ -12,7 +12,18 @@ def ensemble_save(func): - """Decorator for read and write shots for ensemble parallelism""" + """Decorator for read and write shots for ensemble parallelism + + Parameters + ---------- + func : function + Function to be decorated + + Returns + ------- + wrapper : function + Decorated function + """ def wrapper(*args, **kwargs): acq = args[0].get("acquisition") @@ -179,7 +190,26 @@ def wrapper(*args, **kwargs): def write_function_to_grid(function, V, grid_spacing): - """Interpolate a Firedrake function to a structured grid""" + """Interpolate a Firedrake function to a structured grid + + Parameters + ---------- + function : firedrake.Function + Function to interpolate + V : firedrake.FunctionSpace + Function space of function + grid_spacing : float + Spacing of grid points in metres + + Returns + ------- + x : numpy.ndarray + x coordinates of grid points + y : numpy.ndarray + y coordinates of grid points + z : numpy.ndarray + z coordinates of grid points + """ # get DoF coordinates m = V.ufl_domain() W = fire.VectorFunctionSpace(m, V.ufl_element()) @@ -206,7 +236,20 @@ def write_function_to_grid(function, V, grid_spacing): def create_segy(velocity, filename): - """Write the velocity data into a segy file named filename""" + """Write the velocity data into a segy file named filename + + Parameters + ---------- + velocity : firedrake.Function + Velocity in a firedrake function + filename : str + Name of the segy file to write + + Returns + ------- + None + + """ spec = segyio.spec() velocity = np.flipud(velocity.T) @@ -223,6 +266,7 @@ def create_segy(velocity, filename): for tr, il in enumerate(spec.ilines): f.trace[tr] = velocity[:, tr] + return None @ensemble_save def save_shots(model, comm, array, file_name=None): diff --git a/spyro/pml/damping.py b/spyro/pml/damping.py index 9a812df2..630c7594 100644 --- a/spyro/pml/damping.py +++ b/spyro/pml/damping.py @@ -21,7 +21,51 @@ def functions( y2=None, b_pml=None, ): - """Damping functions for the perfect matched layer for 2D and 3D""" + """Damping functions for the perfect matched layer for 2D and 3D + + Parameters + ---------- + model : dict + Dictionary with the model parameters + V : obj + Firedrake function space + dimension : int + Dimension of the problem + x : obj + Firedrake spatial coordinate + x1 : float + x coordinate of the left boundary of the PML + x2 : float + x coordinate of the right boundary of the PML + a_pml : float + Width of the PML in the x direction + z : obj + Firedrake spatial coordinate + z1 : float + z coordinate of the bottom boundary of the PML + z2 : float + z coordinate of the top boundary of the PML + c_pml : float + Width of the PML in the z direction + y : obj, optional + Firedrake spatial coordinate, by default None + y1 : float, optional + y coordinate of the back boundary of the PML, by default None + y2 : float, optional + y coordinate of the front boundary of the PML, by default None + b_pml : float, optional + Width of the PML in the y direction, by default None + + 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 + + """ damping_type = model["BCs"]["damping_type"] if damping_type == "polynomial": diff --git a/spyro/receivers/Receivers.py b/spyro/receivers/Receivers.py index 49d6f516..56362b27 100644 --- a/spyro/receivers/Receivers.py +++ b/spyro/receivers/Receivers.py @@ -115,6 +115,17 @@ def build_maps(self): Is always automatticaly called when initializing the class, therefore should only be called again if a mesh related attribute changes. + + Returns + ------- + cellIDs: list + List of cell IDs for each receiver + cellVertices: list + List of vertices for each receiver + cellNodeMaps: list + List of node maps for each receiver + cell_tabulations: list + List of tabulations for each receiver """ for rid in range(self.num_receivers): @@ -143,10 +154,12 @@ def build_maps(self): def interpolate(self, field): """Interpolate the solution to the receiver coordinates for one simulation timestep. + Parameters ---------- field: array-like An array of the solution at a given timestep at all nodes + Returns ------- solution_at_receivers: list @@ -171,6 +184,7 @@ def apply_receivers_as_source(self, rhs_forcing, residual, IT): and timesteps IT: int Desired time step number to get residual value from + Returns ------- rhs_forcing: object @@ -288,6 +302,7 @@ def __new_at(self, udat, receiver_id): receiver_id: a list of integers A list of receiver ids, ranging from 0 to total receivers minus one. + Returns ------- at: Function value at given receiver @@ -328,6 +343,7 @@ def __func_receiver_locator_3D(self): the receiver. The matrix has the deegres of freedom of the nodes inside same element as the receiver. + """ num_recv = self.num_receivers @@ -360,6 +376,7 @@ def __func_node_locations_3D(self): """Function that returns a list which includes a numpy matrix where line n has the x and y values of the nth degree of freedom, and a numpy matrix of the vertex coordinates. + """ x, y, z = SpatialCoordinate(self.mesh) ux = Function(self.space).interpolate(x) @@ -496,7 +513,21 @@ def set_point_cloud(self, comm): def choosing_element(V, degree): - """Chooses UFL element based on desired function space""" + """Chooses UFL element based on desired function space + and degree of interpolation. + + Parameters + ---------- + V : firedrake.FunctionSpace + Function space to be used. + degree : int + Degree of interpolation. + + Returns + ------- + element : UFL element + UFL element to be used in the interpolation. + """ cell_geometry = V.mesh().ufl_cell() if cell_geometry == quadrilateral: T = UFCQuadrilateral() diff --git a/spyro/solvers/helpers.py b/spyro/solvers/helpers.py index 948a0d12..3dc71430 100644 --- a/spyro/solvers/helpers.py +++ b/spyro/solvers/helpers.py @@ -16,7 +16,25 @@ def fill(usol_recv, is_local, nt, nr): """Fills usol_recv with -99999 value - when it isn't local to any core""" + when it isn't local to any core + + Parameters + ---------- + usol_recv : list + List of numpy arrays + is_local : list + List of booleans indicating if the receiver is local to the core + nt : int + Number of timesteps + nr : int + Number of receivers + + Returns + ------- + usol_recv : list + List of numpy arrays + + """ usol_recv = np.asarray(usol_recv) for ti in range(nt): for rn in range(nr): @@ -26,7 +44,22 @@ def fill(usol_recv, is_local, nt, nr): def create_output_file(name, comm, source_num): - """Saves shots in output file""" + """Saves shots in output file + + Parameters + ---------- + name : str + Name of the output file + comm : object + MPI communicator + source_num : int + Source number + + Returns + ------- + outfile : object + Firedrake.File object + """ if io.is_owner(comm, source_num): outfile = File( os.getcwd() @@ -41,7 +74,16 @@ def create_output_file(name, comm, source_num): def display(comm, source_num): - """Displays current shot and ensemble in terminal""" + """Displays current shot and ensemble in terminal + + Parameters + ---------- + comm : object + MPI communicator + source_num : int + Source number + + """ if comm.comm.rank == 0: print( "Timestepping for shot #", @@ -54,7 +96,15 @@ def display(comm, source_num): def display_progress(comm, t): - """Displays progress time""" + """Displays progress time + + Parameters + ---------- + comm : object + MPI communicator + t : float + Current time + """ if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0: print(f"Simulation time is: {t:{10}.{4}} seconds", flush=True) @@ -66,7 +116,22 @@ def parallel_print(string, comm): def receivers_local(mesh, dimension, receiver_locations): - """Locates receiveirs in cells""" + """Locates receivers in cells + + Parameters + ---------- + mesh : object + Firedrake mesh object + dimension : int + Dimension of the mesh + receiver_locations : list + List of receiver locations + + Returns + ------- + list + List of receiver locations in cells + """ if dimension == 2: return [mesh.locate_cell([z, x], tolerance=0.01) for z, x in receiver_locations] elif dimension == 3: diff --git a/spyro/sources/Sources.py b/spyro/sources/Sources.py index 308e3f67..a31475c4 100644 --- a/spyro/sources/Sources.py +++ b/spyro/sources/Sources.py @@ -83,7 +83,20 @@ def __init__(self, model, mesh, V, my_ensemble): super().build_maps() def apply_source(self, rhs_forcing, value): - """Applies source in a assembled right hand side.""" + """Applies source in a assembled right hand side. + + Parameters + ---------- + rhs_forcing: Firedrake.Function + The right hand side of the equation + value: float + The value to be applied at the source location + + Returns + ------- + rhs_forcing: Firedrake.Function + The right hand side of the equation with the source applied + """ for source_id in range(self.num_receivers): if self.is_local[source_id] and source_id == self.current_source: for i in range(len(self.cellNodeMaps[source_id])): @@ -110,6 +123,23 @@ def ricker_wavelet(t, freq, amp=1.0, delay=1.5): """Creates a Ricker source function with a delay in term of multiples of the distance between the minimums. + + Parameters + ---------- + t: float + Time + freq: float + Frequency of the wavelet + amp: float + Amplitude of the wavelet + delay: float + Delay in term of multiples of the distance + between the minimums + + Returns + ------- + float + Value of the wavelet at time t """ t = t - delay * math.sqrt(6.0) / (math.pi * freq) return ( @@ -124,6 +154,24 @@ def ricker_wavelet(t, freq, amp=1.0, delay=1.5): def full_ricker_wavelet(dt, tf, freq, amp=1.0, cutoff=None): """Compute the Ricker wavelet optionally applying low-pass filtering using cutoff frequency in Hertz. + + Parameters + ---------- + dt: float + Time step + tf: float + Final time + freq: float + Frequency of the wavelet + amp: float + Amplitude of the wavelet + cutoff: float + Cutoff frequency in Hertz + + Returns + ------- + full_wavelet: numpy array + Array containing the wavelet values at each time step """ nt = int(tf / dt) # number of timesteps time = 0.0 diff --git a/spyro/tools/gradient_test_ad.py b/spyro/tools/gradient_test_ad.py index a1c3ce51..1ff6b876 100644 --- a/spyro/tools/gradient_test_ad.py +++ b/spyro/tools/gradient_test_ad.py @@ -7,6 +7,29 @@ def gradient_test_acoustic(model, mesh, V, comm, vp_exact, vp_guess, mask=None): + """Gradient test for the acoustic FWI problem + + Parameters + ---------- + model : `dictionary` + Contains simulation parameters and options. + mesh : a Firedrake.mesh + 2D/3D simplicial mesh read in by Firedrake.Mesh + V : Firedrake.FunctionSpace object + The space of the finite elements + comm : Firedrake.ensemble_communicator + An ensemble communicator + vp_exact : Firedrake.Function + The exact velocity model + vp_guess : Firedrake.Function + The guess velocity model + mask : Firedrake.Function, optional + A mask for the gradient test + + Returns + ------- + None + """ import firedrake_adjoint as fire_adj with fire_adj.stop_annotating(): diff --git a/spyro/tools/grid_point_calculator.py b/spyro/tools/grid_point_calculator.py index 7991f4d7..ba14ea9f 100644 --- a/spyro/tools/grid_point_calculator.py +++ b/spyro/tools/grid_point_calculator.py @@ -51,6 +51,23 @@ def minimum_grid_point_calculator(grid_point_calculator_parameters): def wave_solver(model, G, comm=False): + """Forward solver for the acoustic wave equation + + Parameters + ---------- + model : `dictionary` + Contains simulation parameters and options. + G : `float` + Grid point density + comm : Firedrake.ensemble_communicator, optional + An ensemble communicator + + Returns + ------- + p_recv : + The pressure field at the receivers + + """ minimum_mesh_velocity = model["testing_parameters"]["minimum_mesh_velocity"] model["mesh"]["meshfile"] = "meshes/2Dhomogeneous" + str(G) + ".msh" try: @@ -114,6 +131,22 @@ def wave_solver(model, G, comm=False): def generate_mesh(model, G, comm): + """Function to generate a mesh + + Parameters + ---------- + model : `dictionary` + Contains simulation parameters and options. + G : `float` + Grid point density + comm : Firedrake.ensemble_communicator + An ensemble communicator + + Returns + ------- + mesh : `firedrake.mesh` + The mesh + """ if model["opts"]["dimension"] == 2: mesh = generate_mesh2D(model, G, comm) elif model["opts"]["dimension"] == 3: @@ -126,6 +159,29 @@ def generate_mesh(model, G, comm): def searching_for_minimum( model, p_exact, TOL, accuracy=0.1, starting_G=7.0, comm=False ): + """Function to find the minimum grid point density for a given error + + Parameters + ---------- + model : `dictionary` + Contains simulation parameters and options. + p_exact : `firedrake.Function` + The exact pressure field + TOL : `float` + The accepted error threshold + accuracy : `float`, optional + The accuracy of the search + starting_G : `float`, optional + The starting grid point density + comm : Firedrake.ensemble_communicator, optional + An ensemble communicator + + Returns + ------- + G : `float` + The minimum grid point density + + """ error = 100.0 G = starting_G @@ -202,6 +258,25 @@ def grid_point_to_mesh_point_converter_for_seismicmesh(model, G): def error_calc(p_exact, p, model, comm=False): + """ Calculates the error between the exact and the numerical solution + + Parameters + ---------- + p_exact : `firedrake.Function` + The exact pressure field + p : `firedrake.Function` + The numerical pressure field + model : `dictionary` + Contains simulation parameters and options. + comm : Firedrake.ensemble_communicator, optional + An ensemble communicator + + Returns + ------- + error : `float` + The error between the exact and the numerical solution + + """ # p0 doesn't necessarily have the same dt as p_exact # therefore we have to interpolate the missing points # to have them at the same length @@ -357,6 +432,16 @@ def time_interpolation_line(p_old, p_exact, model): def generate_mesh2D(model, G, comm): + """Generates 2D mesh using seismicmesh + Parameters + ---------- + model : dict + Dictionary containing the model parameters + G : float + Grid points per wavelength + comm : object + MPI communicator + """ import SeismicMesh if comm.comm.rank == 0: @@ -465,6 +550,16 @@ def generate_mesh2D(model, G, comm): def generate_mesh3D(model, G, comm): + """Generates 3D mesh using seismicmesh + Parameters + ---------- + model : dict + Dictionary containing the model parameters + G : float + Grid points per wavelength + comm : object + MPI communicator + """ import SeismicMesh print("Entering mesh generation", flush=True) diff --git a/spyro/utils/estimate_timestep.py b/spyro/utils/estimate_timestep.py index d77399cc..f7706eee 100644 --- a/spyro/utils/estimate_timestep.py +++ b/spyro/utils/estimate_timestep.py @@ -13,6 +13,7 @@ def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True): generalized eigenvalue exactly ONLY WORKS WITH KMV ELEMENTS + """ u, v = fd.TrialFunction(V), fd.TestFunction(V) diff --git a/spyro/utils/utils.py b/spyro/utils/utils.py index 354fe9ba..2fcc9e5b 100644 --- a/spyro/utils/utils.py +++ b/spyro/utils/utils.py @@ -8,6 +8,22 @@ def butter_lowpass_filter(shot, cutoff, fs, order=2): """Low-pass filter the shot record with sampling-rate fs Hz and cutoff freq. Hz + + Parameters + ---------- + shot : numpy array + Shot record + cutoff : float + Cutoff frequency in Hertz + fs : float + Sampling rate in Hertz + order : int + Order of the filter + + Returns + ------- + filtered_shot : numpy array + Filtered shot record """ nyq = 0.5 * fs # Nyquist Frequency normal_cutoff = cutoff / nyq