From 110720d3588dec8f8d462efdb35000cd93af2f38 Mon Sep 17 00:00:00 2001 From: olender Date: Tue, 11 Jun 2024 13:31:54 -0300 Subject: [PATCH 01/11] creating first test files --- simple_mms.py | 145 ++++++++++++++++++++++++++++++++++++++++++++++++++ temp_MMS.py | 28 ++++++++++ 2 files changed, 173 insertions(+) create mode 100644 simple_mms.py create mode 100644 temp_MMS.py diff --git a/simple_mms.py b/simple_mms.py new file mode 100644 index 00000000..de70a538 --- /dev/null +++ b/simple_mms.py @@ -0,0 +1,145 @@ +import firedrake as fire +from firedrake import dot, grad, sin, pi, dx +import FIAT +import finat +import math + + +def gauss_lobatto_legendre_line_rule(degree): + fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule + fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) + finat_ps = finat.point_set.GaussLobattoLegendrePointSet + points = finat_ps(fiat_rule.get_points()) + weights = fiat_rule.get_weights() + return finat.quadrature.QuadratureRule(points, weights) + + +def gauss_lobatto_legendre_cube_rule(dimension, degree): + make_tensor_rule = finat.quadrature.TensorProductQuadratureRule + result = gauss_lobatto_legendre_line_rule(degree) + for _ in range(1, dimension): + line_rule = gauss_lobatto_legendre_line_rule(degree) + result = make_tensor_rule([result, line_rule]) + return result + + +def analytical_solution(t, V, mesh_z, mesh_x): + analytical = fire.Function(V) + x = mesh_z + y = mesh_x + analytical.interpolate(x * (x + 1) * y * (y - 1) * t) + + return analytical + + +quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=4) + +length_z = 1.0 +length_x = 1.0 + +mesh_dx = 0.02 +nz = int(length_z / mesh_dx) +nx = int(length_x / mesh_dx) +mesh = fire.RectangleMesh( + nz, + nx, + length_z, + length_x, + quadrilateral=True +) +mesh_z, mesh_x = fire.SpatialCoordinate(mesh) + +output = fire.File("debug.pvd") + +element = fire.FiniteElement( # noqa: F405 + "CG", mesh.ufl_cell(), degree=4, variant="spectral" +) +V = fire.FunctionSpace(mesh, element) + +X = fire.Function(V) +u_nm1 = fire.Function(V) +u_n = fire.Function(V) + +# Setting C +c = fire.Function(V, name="velocity") +c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) + +# setting xy part of source +q_xy = fire.Function(V) +xy = fire.project((-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x), V) +q_xy.assign(xy) + +# Starting time integration +t = 0.0 +final_time = 1.0 +dt = 0.0001 +nt = int((final_time - t) / dt) + 1 +u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) +u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) +u_np1 = fire.Function(V, name="pressure t +dt") +u = fire.TrialFunction(V) +v = fire.TestFunction(V) + +m1 = ( + (1 / (c * c)) + * ((u - 2.0 * u_n + u_nm1) / fire.Constant(dt**2)) + * v + * dx(scheme=quad_rule) +) +a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) # explicit + +B = fire.Function(V) + +form = m1 + a +lhs = fire.lhs(form) +rhs = fire.rhs(form) + +A = fire.assemble(lhs, mat_type="matfree") +solver_parameters = { + "ksp_type": "preonly", + "pc_type": "jacobi", +} +solver = fire.LinearSolver( + A, solver_parameters=solver_parameters +) + +for step in range(nt): + q = q_xy * fire.Constant(2 * t) + m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / fire.Constant(dt**2)) + * v + * dx(scheme=quad_rule) + ) + a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) + le = q * v * dx(scheme=quad_rule) + + form = m1 + a - le + rhs = fire.rhs(form) + + B = fire.assemble(rhs, tensor=B) + + solver.solve(X, B) + + u_np1.assign(X) + + if (step - 1) % 100 == 0: + assert ( + fire.norm(u_n) < 1 + ), "Numerical instability. Try reducing dt or building the \ + mesh differently" + output.write(u_n, time=t, name="Pressure") + if t > 0: + print(f"Simulation time is: {t:{10}.{4}} seconds", flush=True) + + u_nm1.assign(u_n) + u_n.assign(u_np1) + + t = step * float(dt) + +last_analytical_sol = analytical_solution(t, V, mesh_z, mesh_x) +error = fire.errornorm(u_n, last_analytical_sol) +test = math.isclose(error, 0.0, abs_tol=1e-7) + +print(F"Test is {test}") diff --git a/temp_MMS.py b/temp_MMS.py new file mode 100644 index 00000000..01a42db9 --- /dev/null +++ b/temp_MMS.py @@ -0,0 +1,28 @@ +import math +from copy import deepcopy +from firedrake import * +import spyro + +from test.model import dictionary as model + +model["acquisition"]["source_type"] = "MMS" +model["cell_type"] = "Q" +model["variant"] = "lumped" + + +def run_solve(model): + testmodel = deepcopy(model) + + Wave_obj = spyro.AcousticWaveMMS(dictionary=testmodel) + Wave_obj.set_mesh(mesh_parameters={"dx": 0.02}) + Wave_obj.set_initial_velocity_model(expression="1 + sin(pi*-z)*sin(pi*x)") + Wave_obj.forward_solve() + + u_an = Wave_obj.analytical + u_num = Wave_obj.u_n + + return errornorm(u_num, u_an) + + +error = run_solve(model) +print(error) From b62fa7af4a4992b60b694e3105db41c0c4c91e04 Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 10:17:28 -0300 Subject: [PATCH 02/11] MMS triangular lumped working --- .../acoustic_solver_construction_no_pml.py | 13 ++++-- spyro/solvers/mms_acoustic.py | 3 +- .../time_integration_central_difference.py | 14 ++++-- test/test_MMS.py | 46 ++++++++++--------- 4 files changed, 43 insertions(+), 33 deletions(-) diff --git a/spyro/solvers/acoustic_solver_construction_no_pml.py b/spyro/solvers/acoustic_solver_construction_no_pml.py index 252c0180..43bc736d 100644 --- a/spyro/solvers/acoustic_solver_construction_no_pml.py +++ b/spyro/solvers/acoustic_solver_construction_no_pml.py @@ -20,8 +20,10 @@ def construct_solver_or_matrix_no_pml(Wave_object): u_nm1 = fire.Function(V, name="pressure t-dt") u_n = fire.Function(V, name="pressure") + u_np1 = fire.Function(V, name="pressure t+dt") Wave_object.u_nm1 = u_nm1 Wave_object.u_n = u_n + Wave_object.u_np1 = u_np1 Wave_object.current_time = 0.0 dt = Wave_object.dt @@ -35,7 +37,7 @@ def construct_solver_or_matrix_no_pml(Wave_object): ) a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) # explicit - B = fire.Function(V) + B = fire.Cofunction(V.dual()) form = m1 + a lhs = fire.lhs(form) @@ -43,9 +45,12 @@ def construct_solver_or_matrix_no_pml(Wave_object): Wave_object.lhs = lhs A = fire.assemble(lhs, mat_type="matfree") - Wave_object.solver = fire.LinearSolver( - A, solver_parameters=Wave_object.solver_parameters - ) + # Wave_object.solver = fire.LinearSolver( + # A, solver_parameters=Wave_object.solver_parameters + # ) + lin_var = fire.LinearVariationalProblem(lhs, rhs + B, u_np1) + solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} + Wave_object.solver = fire.LinearVariationalSolver(lin_var,solver_parameters=solver_parameters) Wave_object.rhs = rhs Wave_object.B = B diff --git a/spyro/solvers/mms_acoustic.py b/spyro/solvers/mms_acoustic.py index ccbfb28c..ca21dc81 100644 --- a/spyro/solvers/mms_acoustic.py +++ b/spyro/solvers/mms_acoustic.py @@ -27,8 +27,7 @@ def mms_source_in_space(self): if self.dimension == 2: # xy = fire.project(sin(pi*x)*sin(pi*y), V) # self.q_xy.assign(xy) - xy = fire.project((-(x**2) - x - y**2 + y), V) - self.q_xy.assign(xy) + self.q_xy.interpolate(-(x**2) - x - y**2 + y) elif self.dimension == 3: z = self.mesh_y # xyz = fire.project(sin(pi*x)*sin(pi*y)*sin(pi*z), V) diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index 410bf5cd..378e002c 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -1,5 +1,6 @@ import firedrake as fire from firedrake import Constant, dx, dot, grad +import numpy as np from ..io.basicio import parallel_print from . import helpers @@ -47,7 +48,7 @@ def central_difference(Wave_object, source_id=0): u_n = Wave_object.u_n u_np1 = fire.Function(Wave_object.function_space) - rhs_forcing = fire.Function(Wave_object.function_space) + rhs_forcing = fire.Cofunction(Wave_object.function_space.dual()) usol = [ fire.Function(Wave_object.function_space, name="pressure") for t in range(nt) @@ -64,7 +65,7 @@ def central_difference(Wave_object, source_id=0): f = excitations.apply_source(rhs_forcing, Wave_object.wavelet[step]) B0 = B.sub(0) B0 += f - Wave_object.solver.solve(X, B) + Wave_object.solver.solve() u_np1.assign(X) @@ -241,9 +242,12 @@ def central_difference_MMS(Wave_object, source_id=0): u_nm1 = Wave_object.u_nm1 u_n = Wave_object.u_n + u_np1 = Wave_object.u_np1 u_nm1.assign(Wave_object.analytical_solution(t - 2 * dt)) u_n.assign(Wave_object.analytical_solution(t - dt)) - u_np1 = fire.Function(Wave_object.function_space, name="pressure t +dt") + # u_np1 = fire.Function(Wave_object.function_space, name="pressure t +dt") + # u_nm1.dat.data[:] = np.load("old_u_nm1.npy") + # u_n.dat.data[:] = np.load("old_u_n.npy") u = fire.TrialFunction(Wave_object.function_space) v = fire.TestFunction(Wave_object.function_space) @@ -277,9 +281,9 @@ def central_difference_MMS(Wave_object, source_id=0): B = fire.assemble(rhs, tensor=B) - Wave_object.solver.solve(X, B) + Wave_object.solver.solve(u_np1, B) - u_np1.assign(X) + # u_np1.assign(X) usol_recv.append( Wave_object.receivers.interpolate(u_np1.dat.data_ro_with_halos[:]) diff --git a/test/test_MMS.py b/test/test_MMS.py index ab2cb44c..e1cafbc1 100644 --- a/test/test_MMS.py +++ b/test/test_MMS.py @@ -4,29 +4,10 @@ from firedrake import * import spyro -from .model import dictionary as model - +from model import dictionary as model model["acquisition"]["source_type"] = "MMS" -@pytest.fixture(params=["triangle", "square"]) -def mesh_type(request): - if mesh_type == "triangle": - model["cell_type"] = "triangles" - elif mesh_type == "square": - model["cell_type"] = "quadrilaterals" - return request.param - - -@pytest.fixture(params=["lumped", "equispaced"]) -def method_type(request): - if method_type == "lumped": - model["variant"] = "lumped" - elif method_type == "equispaced": - model["variant"] = "equispaced" - return request.param - - def run_solve(model): testmodel = deepcopy(model) @@ -41,7 +22,28 @@ def run_solve(model): return errornorm(u_num, u_an) -def test_method(mesh_type, method_type): +def run_method(mesh_type, method_type): + model["options"]["cell_type"] = mesh_type + model["options"]["variant"] = method_type + print(f"For {mesh_type} and {method_type}") error = run_solve(model) + test = math.isclose(error, 0.0, abs_tol=1e-7) + print(f"Error is {error}") + print(f"Test: {test}") + + assert test + + +def test_method_triangles_lumped(): + run_method("triangles", "lumped") + + +def test_method_quads_lumped(): + run_method("quadrilaterals", "lumped") + + +if __name__ == "__main__": + test_method_triangles_lumped() + test_method_quads_lumped() - assert math.isclose(error, 0.0, abs_tol=1e-7) + print("END") From decbcd5951cad8c67858ec8171eb6454dc53a993 Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 10:32:47 -0300 Subject: [PATCH 03/11] Switching rectangle example to triangles, now passing analytical sol test --- spyro/examples/rectangle.py | 2 +- spyro/solvers/acoustic_solver_construction_no_pml.py | 12 ++++++------ spyro/solvers/time_integration_central_difference.py | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/spyro/examples/rectangle.py b/spyro/examples/rectangle.py index fc82c892..96ed883d 100644 --- a/spyro/examples/rectangle.py +++ b/spyro/examples/rectangle.py @@ -24,7 +24,7 @@ rectangle_dictionary = {} rectangle_dictionary["options"] = { # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) - "cell_type": "Q", + "cell_type": "T", "variant": "lumped", # lumped, equispaced or DG, default is lumped "degree": 4, # p order "dimension": 2, # dimension diff --git a/spyro/solvers/acoustic_solver_construction_no_pml.py b/spyro/solvers/acoustic_solver_construction_no_pml.py index 43bc736d..d5a18f3d 100644 --- a/spyro/solvers/acoustic_solver_construction_no_pml.py +++ b/spyro/solvers/acoustic_solver_construction_no_pml.py @@ -45,12 +45,12 @@ def construct_solver_or_matrix_no_pml(Wave_object): Wave_object.lhs = lhs A = fire.assemble(lhs, mat_type="matfree") - # Wave_object.solver = fire.LinearSolver( - # A, solver_parameters=Wave_object.solver_parameters - # ) - lin_var = fire.LinearVariationalProblem(lhs, rhs + B, u_np1) - solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} - Wave_object.solver = fire.LinearVariationalSolver(lin_var,solver_parameters=solver_parameters) + Wave_object.solver = fire.LinearSolver( + A, solver_parameters=Wave_object.solver_parameters + ) + # lin_var = fire.LinearVariationalProblem(lhs, rhs + B, u_np1) + # solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} + # Wave_object.solver = fire.LinearVariationalSolver(lin_var,solver_parameters=solver_parameters) Wave_object.rhs = rhs Wave_object.B = B diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index 378e002c..51ef3bb6 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -65,7 +65,7 @@ def central_difference(Wave_object, source_id=0): f = excitations.apply_source(rhs_forcing, Wave_object.wavelet[step]) B0 = B.sub(0) B0 += f - Wave_object.solver.solve() + Wave_object.solver.solve(X, B) u_np1.assign(X) From 31724f9f019423f2e030d61ea87a987c0b6fb000 Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 12:09:12 -0300 Subject: [PATCH 04/11] updated test for correct values for MLT --- test/test_cpw_calc.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_cpw_calc.py b/test/test_cpw_calc.py index bc5e99c7..35e4e825 100644 --- a/test/test_cpw_calc.py +++ b/test/test_cpw_calc.py @@ -17,7 +17,7 @@ def test_cpw_calc(): "velocity_model_file_name": None, # FEM to evaluate such as `KMV` or `spectral` # (GLL nodes on quads and hexas) - "FEM_method_to_evaluate": "spectral_quadrilateral", + "FEM_method_to_evaluate": "mass_lumped_triangle", "dimension": 2, # Domain dimension. Either 2 or 3. # Either near or line. Near defines a receiver grid near to the source, "receiver_setup": "near", @@ -32,7 +32,7 @@ def test_cpw_calc(): # grid point density to use in the reference case (float) "C_reference": None, "desired_degree": 4, # degree we are calculating G for. (int) - "C_initial": 2.4, # Initial G for line search (float) + "C_initial": 2.2, # Initial G for line search (float) "accepted_error_threshold": 0.05, "C_accuracy": 0.1, } @@ -60,7 +60,8 @@ def test_cpw_calc(): # Check if cpw is within error TOL, starting search at min min = Cpw_calc.find_minimum() - test3 = np.isclose(2.5, min) + print(f"Minimum of {min}") + test3 = np.isclose(2.3, min) print("END") assert all([test1, test2, test3]) From b44caeb07637203d70ee607a6471480e016d17cc Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 12:24:08 -0300 Subject: [PATCH 05/11] removing outdated unused test --- test/test_gradient.py | 254 ------------------------------------------ 1 file changed, 254 deletions(-) delete mode 100644 test/test_gradient.py diff --git a/test/test_gradient.py b/test/test_gradient.py deleted file mode 100644 index c64c67c4..00000000 --- a/test/test_gradient.py +++ /dev/null @@ -1,254 +0,0 @@ -import numpy as np -import os -from firedrake import * -import spyro -from spyro.domains import quadrature -import pytest - -from .inputfiles.Model1_gradient_2d import model -from .inputfiles.Model1_gradient_2d_pml import model_pml - - -dictionary = {} -dictionary["options"] = { - "cell_type": "T", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) - "variant": "lumped", # lumped, equispaced or DG, default is lumped - "method": "MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method - "degree": 4, # p order - "dimension": 2, # dimension - "automatic_adjoint": False, -} - -# Number of cores for the shot. For simplicity, we keep things serial. -# spyro however supports both spatial parallelism and "shot" parallelism. -dictionary["parallelism"] = { - "type": "automatic", # options: automatic (same number of cores for evey processor) or spatial -} - -# Define the domain size without the PML. Here we'll assume a 0.75 x 1.50 km -# domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb -# outgoing waves on three sides (eg., -z, +-x sides) of the domain. -dictionary["mesh"] = { - "Lz": 1.0, # depth in km - always positive # Como ver isso sem ler a malha? - "Lx": 1.0, # width in km - always positive - "Ly": 0.0, # thickness in km - always positive - "mesh_file": None, -} -dictionary[ - "synthetic_data" -] = { # For use only if you are using a synthetic test model or a forward only simulation -adicionar discrição para modelo direto - "real_mesh_file": None, - "real_velocity_file": None, -} -dictionary["inversion"] = { - "perform_fwi": False, # switch to true to make a FWI - "initial_guess_model_file": None, - "shot_record_file": None, - "optimization_parameters": None, -} - -# Specify a 250-m PML on the three sides of the domain to damp outgoing waves. -dictionary["absorving_boundary_conditions"] = { - "status": False, # True or false - "outer_bc": "non-reflective", # None or non-reflective (outer boundary condition) - "damping_type": "polynomial", # polynomial, hyperbolic, shifted_hyperbolic - "exponent": 2, # damping layer has a exponent variation - "cmax": 4.7, # maximum acoustic wave velocity in PML - km/s - "R": 1e-6, # theoretical reflection coefficient - "lz": 0.25, # thickness of the PML in the z-direction (km) - always positive - "lx": 0.25, # thickness of the PML in the x-direction (km) - always positive - "ly": 0.0, # thickness of the PML in the y-direction (km) - always positive -} - -# Create a source injection operator. Here we use a single source with a -# Ricker wavelet that has a peak frequency of 8 Hz injected at the center of the mesh. -# We also specify to record the solution at 101 microphones near the top of the domain. -# This transect of receivers is created with the helper function `create_transect`. -dictionary["acquisition"] = { - "source_type": "ricker", - "source_locations": [(-0.1, 0.5)], - "frequency": 5.0, - "delay": 1.0, - "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 0.9), 20), -} - -# Simulate for 2.0 seconds. -dictionary["time_axis"] = { - "initial_time": 0.0, # Initial time for event - "final_time": 2.00, # Final time for event - "dt": 0.001, # timestep size - "amplitude": 1, # the Ricker has an amplitude of 1. - "output_frequency": 100, # how frequently to output solution to pvds - Perguntar Daiane ''post_processing_frequnecy' - "gradient_sampling_frequency": 100, # how frequently to save solution to RAM - Perguntar Daiane 'gradient_sampling_frequency' -} -dictionary["visualization"] = { - "forward_output": True, - "output_filename": "results/forward_output.pvd", - "fwi_velocity_model_output": False, - "velocity_model_filename": None, - "gradient_output": False, - "gradient_filename": None, - "adjoint_output": False, - "adjoint_filename": None, -} -outfile_total_gradient = File(os.getcwd() + "/results/Gradient.pvd") - -# forward = spyro.solvers.forward -# gradient = spyro.solvers.gradient -functional = spyro.utils.compute_functional - - -def _make_vp_exact(V, mesh): - """Create a circle with higher velocity in the center""" - z, x = SpatialCoordinate(mesh) - vp_exact = Function(V).interpolate( - 4.0 + 1.0 * tanh(10.0 * (0.5 - sqrt((z - 1.5) ** 2 + (x + 1.5) ** 2))) - ) - File("exact_vel.pvd").write(vp_exact) - return vp_exact - - -def _make_vp_exact_pml(V, mesh): - """Create a half space""" - z, x = SpatialCoordinate(mesh) - velocity = conditional(z > -0.5, 1.5, 4.0) - vp_exact = Function(V, name="vp").interpolate(velocity) - File("exact_vel.pvd").write(vp_exact) - return vp_exact - - -def _make_vp_guess(V, mesh): - """The guess is a uniform velocity of 4.0 km/s""" - z, x = SpatialCoordinate(mesh) - vp_guess = Function(V).interpolate(4.0 + 0.0 * x) - File("guess_vel.pvd").write(vp_guess) - return vp_guess - - -@pytest.mark.skip(reason="not yet implemented") -def test_gradient(): - _test_gradient(model) - - -@pytest.mark.skip(reason="no way of currently testing this") -def test_gradient_pml(): - _test_gradient(model_pml, pml=True) - - -def _test_gradient(options, pml=False): - comm = spyro.utils.mpi_init(options) - - mesh, V = spyro.basicio.read_mesh(options, comm) - - if pml: - vp_exact = _make_vp_exact_pml(V, mesh) - z, x = SpatialCoordinate(mesh) - Lx = model_pml["mesh"]["Lx"] - Lz = model_pml["mesh"]["Lz"] - x1 = 0.0 - z2 = -Lz - boxx1 = Function(V).interpolate(conditional(x > x1, 1.0, 0.0)) - boxx2 = Function(V).interpolate(conditional(x < Lx, 1.0, 0.0)) - boxz1 = Function(V).interpolate(conditional(z > z2, 1.0, 0.0)) - mask = Function(V).interpolate(boxx1 * boxx2 * boxz1) - File("mask.pvd").write(mask) - else: - vp_exact = _make_vp_exact(V, mesh) - - mask = Function(V).assign(1.0) - - vp_guess = _make_vp_guess(V, mesh) - - sources = spyro.Sources(options, mesh, V, comm) - - receivers = spyro.Receivers(options, mesh, V, comm) - - wavelet = spyro.full_ricker_wavelet( - options["timeaxis"]["dt"], - options["timeaxis"]["tf"], - options["acquisition"]["frequency"], - ) - - # simulate the exact model - p_exact, p_exact_recv = forward( - options, - mesh, - comm, - vp_exact, - sources, - wavelet, - receivers, - ) - - # simulate the guess model - p_guess, p_guess_recv = forward( - options, - mesh, - comm, - vp_guess, - sources, - wavelet, - receivers, - ) - - misfit = p_exact_recv - p_guess_recv - - qr_x, _, _ = quadrature.quadrature_rules(V) - - Jm = functional(options, misfit) - print("\n Cost functional at fixed point : " + str(Jm) + " \n ") - - # compute the gradient of the control (to be verified) - dJ = gradient(options, mesh, comm, vp_guess, receivers, p_guess, misfit) - dJ.dat.data[:] = dJ.dat.data[:] * mask.dat.data[:] - File("gradient.pvd").write(dJ) - - steps = [1e-3, 1e-4, 1e-5] # , 1e-6] # step length - - delta_m = Function(V) # model direction (random) - delta_m.assign(dJ) - - # this deepcopy is important otherwise pertubations accumulate - vp_original = vp_guess.copy(deepcopy=True) - - errors = [] - for step in steps: # range(3): - # steps.append(step) - # perturb the model and calculate the functional (again) - # J(m + delta_m*h) - vp_guess = vp_original + step * delta_m - _, p_guess_recv = forward( - options, - mesh, - comm, - vp_guess, - sources, - wavelet, - receivers, - ) - - Jp = functional(options, p_exact_recv - p_guess_recv) - projnorm = assemble(mask * dJ * delta_m * dx(scheme=qr_x)) - fd_grad = (Jp - Jm) / step - print( - "\n Cost functional for step " - + str(step) - + " : " - + str(Jp) - + ", fd approx.: " - + str(fd_grad) - + ", grad'*dir : " - + str(projnorm) - + " \n ", - ) - - errors.append(100 * ((fd_grad - projnorm) / projnorm)) - # step /= 2 - - # all errors less than 1 % - errors = np.array(errors) - assert (np.abs(errors) < 5.0).all() - - -if __name__ == "__main__": - test_gradient(model) From c895faf12a9c9495c50914695cc56dbd1c94f037 Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 13:17:21 -0300 Subject: [PATCH 06/11] updating gradient --- spyro/solvers/backward_time_integration.py | 2 +- test/test_gradient_2d.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/spyro/solvers/backward_time_integration.py b/spyro/solvers/backward_time_integration.py index c22b718d..527e4f12 100644 --- a/spyro/solvers/backward_time_integration.py +++ b/spyro/solvers/backward_time_integration.py @@ -50,7 +50,7 @@ def backward_wave_propagator(Wave_obj, dt=None): u_n = Wave_obj.u_n u_np1 = fire.Function(Wave_obj.function_space) - rhs_forcing = fire.Function(Wave_obj.function_space) + rhs_forcing = fire.Cofunction(Wave_obj.function_space.dual()) B = Wave_obj.B rhs = Wave_obj.rhs diff --git a/test/test_gradient_2d.py b/test/test_gradient_2d.py index d0c3d63c..6c2ea647 100644 --- a/test/test_gradient_2d.py +++ b/test/test_gradient_2d.py @@ -67,7 +67,7 @@ def check_gradient(Wave_obj_guess, dJ, rec_out_exact, Jm, plot=False): dictionary = {} dictionary["options"] = { - "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "cell_type": "T", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) "variant": "lumped", # lumped, equispaced or DG, default is lumped "degree": 4, # p order "dimension": 2, # dimension From 22a1c8482a613f5f60e8fffcb220781f2a7c7d25 Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 20 Jun 2024 14:08:15 -0300 Subject: [PATCH 07/11] changing to mlt --- test/test_time_convergence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_time_convergence.py b/test/test_time_convergence.py index a92f66ea..8a014dfd 100644 --- a/test/test_time_convergence.py +++ b/test/test_time_convergence.py @@ -18,7 +18,7 @@ def run_forward(dt): dictionary = {} dictionary["options"] = { - "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "cell_type": "T", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) "variant": "lumped", # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method "degree": 4, # p order "dimension": 2, # dimension From 304e01bbf2dd3eec1c21fad8ea86a1b32d493ceb Mon Sep 17 00:00:00 2001 From: olender Date: Thu, 27 Jun 2024 11:22:31 -0300 Subject: [PATCH 08/11] tracking testing files --- check_mlt.py | 28 +++++++++++ check_sem.py | 55 ++++++++++++++++++++ check_u.py | 70 ++++++++++++++++++++++++++ check_u2.py | 113 +++++++++++++++++++++++++++++++++++++++++ check_u3.py | 138 +++++++++++++++++++++++++++++++++++++++++++++++++++ temp_MMS.py | 6 +-- 6 files changed, 407 insertions(+), 3 deletions(-) create mode 100644 check_mlt.py create mode 100644 check_sem.py create mode 100644 check_u.py create mode 100644 check_u2.py create mode 100644 check_u3.py diff --git a/check_mlt.py b/check_mlt.py new file mode 100644 index 00000000..623dc913 --- /dev/null +++ b/check_mlt.py @@ -0,0 +1,28 @@ +import finat +from firedrake import * +import numpy as np + + +def isDiag(M): + i, j = np.nonzero(M) + return np.all(i == j) + +degree = 4 +mesh = RectangleMesh(3, 2, 2.0, 1.0) +element = FiniteElement( # noqa: F405 + "KMV", mesh.ufl_cell(), degree=degree, variant="KMV" +) +V = FunctionSpace(mesh, element) +quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") + +u = TrialFunction(V) +v = TestFunction(V) + +form = u*v*dx(scheme=quad_rule) +A = assemble(form) +M = A.M.values +Mdiag = M.diagonal() + +print(f"Matrix is diagonal:{isDiag(M)}") +np.save("new_diag", Mdiag) +print("END") \ No newline at end of file diff --git a/check_sem.py b/check_sem.py new file mode 100644 index 00000000..d2645690 --- /dev/null +++ b/check_sem.py @@ -0,0 +1,55 @@ +import finat +import FIAT +from firedrake import * +import numpy as np + + +def gauss_lobatto_legendre_line_rule(degree): + fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule + fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) + finat_ps = finat.point_set.GaussLobattoLegendrePointSet + points = finat_ps(fiat_rule.get_points()) + weights = fiat_rule.get_weights() + return finat.quadrature.QuadratureRule(points, weights) + +def gauss_lobatto_legendre_cube_rule(dimension, degree): + make_tensor_rule = finat.quadrature.TensorProductQuadratureRule + result = gauss_lobatto_legendre_line_rule(degree) + for _ in range(1, dimension): + line_rule = gauss_lobatto_legendre_line_rule(degree) + result = make_tensor_rule([result, line_rule]) + return result + + + +def isDiag(M): + i, j = np.nonzero(M) + return np.all(i == j) + +degree = 2 +# mesh = RectangleMesh(3, 2, 2.0, 1.0, quadrilateral=True) +mesh = RectangleMesh(1, 1, 1.0, 1.0, quadrilateral=True) +element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral') +V = FunctionSpace(mesh, element) +quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree) + +u = TrialFunction(V) +v = TestFunction(V) + +form = u*v*dx(scheme=quad_rule) +A = assemble(form) +M = A.M.values +Mdiag = M.diagonal() + +x_mesh, y_mesh = SpatialCoordinate(mesh) +x_func = Function(V) +y_func = Function(V) +x_func.interpolate(x_mesh) +y_func.interpolate(y_mesh) +x = x_func.dat.data[:] +y = y_func.dat.data[:] + +print(f"Matrix is diagonal:{isDiag(M)}") +old_diag = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_sem_diag.npy") +dif = Mdiag-old_diag +print("END") \ No newline at end of file diff --git a/check_u.py b/check_u.py new file mode 100644 index 00000000..e3f3eda2 --- /dev/null +++ b/check_u.py @@ -0,0 +1,70 @@ +import finat +from firedrake import * +import numpy as np + + +def analytical_solution(t, V, mesh_z, mesh_x): + analytical = Function(V) + x = mesh_z + y = mesh_x + analytical.interpolate(x * (x + 1) * y * (y - 1) * t) + + return analytical + + +def isDiag(M): + i, j = np.nonzero(M) + return np.all(i == j) + +degree = 4 +mesh = RectangleMesh(3, 2, 2.0, 1.0) +mesh.coordinates.dat.data[:, 0] *= -1.0 +mesh_z, mesh_x = SpatialCoordinate(mesh) +V = FunctionSpace(mesh, "KMV", degree) +quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") + +u = TrialFunction(V) +v = TestFunction(V) + +c = Function(V, name="velocity") +c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) +u_n = Function(V) +u_nm1 = Function(V) +dt = 0.0005 +t = 0.0 +u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) +u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) + +q = Function(V) +q.assign(1.0) +dt = 0.001 + +m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * v + * dx(scheme=quad_rule) +) +a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) +le = q * v * dx(scheme=quad_rule) + +form = m1 + a - le + +boundary_ids = (1, 2, 3, 4) +bcs = DirichletBC(V, 0.0, boundary_ids) +A = assemble(lhs(form), bcs=bcs) +solver_parameters = { + "ksp_type": "preonly", + "pc_type": "jacobi", +} +solver = LinearSolver( + A, solver_parameters=solver_parameters +) +M = A.M.values +Mdiag = M.diagonal() + +print(f"Matrix is diagonal:{isDiag(M)}") +np.save("/home/olender/Development/issue_50_compatibility/spyro-1/new_diag", Mdiag) + +print("END") \ No newline at end of file diff --git a/check_u2.py b/check_u2.py new file mode 100644 index 00000000..30ff1ad1 --- /dev/null +++ b/check_u2.py @@ -0,0 +1,113 @@ +import finat +from firedrake import * +import numpy as np + + +def analytical_solution(t, V, mesh_z, mesh_x): + analytical = Function(V) + x = mesh_z + y = mesh_x + analytical.interpolate(x * (x + 1) * y * (y - 1) * t) + + return analytical + + +def isDiag(M): + i, j = np.nonzero(M) + return np.all(i == j) + +degree = 4 +mesh = RectangleMesh(50, 50, 1.0, 1.0) +mesh.coordinates.dat.data[:, 0] *= -1.0 +mesh_z, mesh_x = SpatialCoordinate(mesh) +V = FunctionSpace(mesh, "KMV", degree) +quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") + +u = TrialFunction(V) +v = TestFunction(V) + +c = Function(V, name="velocity") +c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) +u_n = Function(V) +u_nm1 = Function(V) +dt = 0.0005 +t = 0.0 +final_time = 1.0 +u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) +u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) + +q_xy = Function(V) +q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x) +q = q_xy * Constant(2 * t) + +nt = int((final_time - t) / dt) + 1 + +m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * v + * dx(scheme=quad_rule) +) +a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) +le = q * v * dx(scheme=quad_rule) + +form = m1 + a - le + +B = Cofunction(V.dual()) + +boundary_ids = (1, 2, 3, 4) +bcs = DirichletBC(V, 0.0, boundary_ids) +A = assemble(lhs(form), bcs=bcs) +solver_parameters = { + "ksp_type": "preonly", + "pc_type": "jacobi", +} +solver = LinearSolver( + A, solver_parameters=solver_parameters +) +As = solver.A +petsc_matrix = As.petscmat +diagonal = petsc_matrix.getDiagonal() +Mdiag = diagonal.array + +np.save("/home/olender/Development/issue_50_compatibility/spyro-1/new_diag", Mdiag) +out = File("new_firedrake_u2.pvd") + +u_np1 = Function(V) +for step in range(nt): + q = q_xy * Constant(2 * t) + m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * v + * dx(scheme=quad_rule) + ) + a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) + le = q * v * dx(scheme=quad_rule) + + form = m1 + a - le + + B = assemble(rhs(form), tensor=B) + + solver.solve(u_np1, B) + + if (step - 1) % 100 == 0: + print(f"Time : {t}") + out.write(u_n) + assert ( + norm(u_n) < 1 + ), "Numerical instability. Try reducing dt or building the \ + mesh differently" + + u_nm1.assign(u_n) + u_n.assign(u_np1) + + t = step * float(dt) + +u_an = analytical_solution(t, V, mesh_z, mesh_x) +error = errornorm(u_n, u_an) +print(f"Error: {error}") + +print("END") \ No newline at end of file diff --git a/check_u3.py b/check_u3.py new file mode 100644 index 00000000..ed72098e --- /dev/null +++ b/check_u3.py @@ -0,0 +1,138 @@ +import finat +import FIAT +from firedrake import * +import numpy as np + + +def gauss_lobatto_legendre_line_rule(degree): + fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule + fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) + finat_ps = finat.point_set.GaussLobattoLegendrePointSet + points = finat_ps(fiat_rule.get_points()) + weights = fiat_rule.get_weights() + return finat.quadrature.QuadratureRule(points, weights) + +def gauss_lobatto_legendre_cube_rule(dimension, degree): + make_tensor_rule = finat.quadrature.TensorProductQuadratureRule + result = gauss_lobatto_legendre_line_rule(degree) + for _ in range(1, dimension): + line_rule = gauss_lobatto_legendre_line_rule(degree) + result = make_tensor_rule([result, line_rule]) + return result + + + +def analytical_solution(t, V, mesh_z, mesh_x): + analytical = Function(V) + x = mesh_z + y = mesh_x + analytical.interpolate(x * (x + 1) * y * (y - 1) * t) + + return analytical + + +def isDiag(M): + i, j = np.nonzero(M) + return np.all(i == j) + +degree = 4 +mesh = RectangleMesh(50, 50, 1.0, 1.0, quadrilateral=True) +mesh.coordinates.dat.data[:, 0] *= -1.0 +mesh_z, mesh_x = SpatialCoordinate(mesh) +element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral') +V = FunctionSpace(mesh, element) +quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree) + +u = TrialFunction(V) +v = TestFunction(V) + +c = Function(V, name="velocity") +c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) +u_n = Function(V) +u_nm1 = Function(V) +dt = 0.0005 +t = 0.0 +final_time = 1.0 +u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) +u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) + +q_xy = Function(V) +q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x) +q = q_xy * Constant(2 * t) + +nt = int((final_time - t) / dt) + 1 + +m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * v + * dx(scheme=quad_rule) +) +a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) +le = q * v * dx(scheme=quad_rule) + +form = m1 + a - le + +B = Cofunction(V.dual()) + +boundary_ids = (1, 2, 3, 4) +bcs = DirichletBC(V, 0.0, boundary_ids) +A = assemble(lhs(form), bcs=bcs) +solver_parameters = { + "ksp_type": "preonly", + "pc_type": "jacobi", +} +solver = LinearSolver( + A, solver_parameters=solver_parameters +) +As = solver.A +petsc_matrix = As.petscmat +diagonal = petsc_matrix.getDiagonal() +Mdiag = diagonal.array + +# old_Mdiag = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_diag.npy") +# old_qxy = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_qxy.npy") # OK +# old_un = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_un.npy") # OK +# old_unm1 = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_unm1.npy") # OK +# old_c = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_c.npy") # OK + +# out = File("new_firedrake_u2.pvd") + +u_np1 = Function(V) +for step in range(nt): + q = q_xy * Constant(2 * t) + m1 = ( + 1 + / (c * c) + * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) + * v + * dx(scheme=quad_rule) + ) + a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) + le = q * v * dx(scheme=quad_rule) + + form = m1 + a - le + + B = assemble(rhs(form), tensor=B) + + solver.solve(u_np1, B) + + if (step - 1) % 100 == 0: + print(f"Time : {t}") + # out.write(u_n) + assert ( + norm(u_n) < 1 + ), "Numerical instability. Try reducing dt or building the \ + mesh differently" + + u_nm1.assign(u_n) + u_n.assign(u_np1) + + t = step * float(dt) + +u_an = analytical_solution(t, V, mesh_z, mesh_x) +error = errornorm(u_n, u_an) +print(f"Error: {error}") + +print("END") \ No newline at end of file diff --git a/temp_MMS.py b/temp_MMS.py index 01a42db9..241b78c3 100644 --- a/temp_MMS.py +++ b/temp_MMS.py @@ -4,10 +4,10 @@ import spyro from test.model import dictionary as model - model["acquisition"]["source_type"] = "MMS" -model["cell_type"] = "Q" -model["variant"] = "lumped" +model["options"]["cell_type"] = "Q" +model["options"]["variant"] = "lumped" + def run_solve(model): From fd7d317d482cdb2998cfc211177f40f366e6a80b Mon Sep 17 00:00:00 2001 From: olender Date: Fri, 19 Jul 2024 18:36:47 -0300 Subject: [PATCH 09/11] updating tests --- .github/workflows/python-tests.yml | 34 +++++------------------------- test/test_MMS.py | 2 +- 2 files changed, 6 insertions(+), 30 deletions(-) diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml index 2665ad7e..f7978bab 100644 --- a/.github/workflows/python-tests.yml +++ b/.github/workflows/python-tests.yml @@ -17,53 +17,29 @@ jobs: - uses: actions/checkout@v3 - name: Running serial tests run: | - source /home/olender/Firedrakes/newest3/firedrake/bin/activate + source /home/olender/firedrakes/2024_07_19/firedrake/bin/activate pytest --cov-report=xml --cov=spyro test/ - name: Running parallel 3D forward test run: | - source /home/olender/Firedrakes/newest3/firedrake/bin/activate + source /home/olender/firedrakes/2024_07_19/firedrake/bin/activate mpiexec -n 6 pytest test_3d/test_hexahedral_convergence.py mpiexec -n 6 pytest test_parallel/test_forward.py mpiexec -n 6 pytest test_parallel/test_fwi.py - name: Covering parallel 3D forward test continue-on-error: true run: | - source /home/olender/Firedrakes/newest3/firedrake/bin/activate + source /home/olender/firedrakes/2024_07_19/firedrake/bin/activate mpiexec -n 6 pytest --cov-report=xml --cov-append --cov=spyro test_3d/test_hexahedral_convergence.py - name: Covering parallel forward test continue-on-error: true run: | - source /home/olender/Firedrakes/newest3/firedrake/bin/activate + source /home/olender/firedrakes/2024_07_19/firedrake/bin/activate mpiexec -n 6 pytest --cov-report=xml --cov-append --cov=spyro test_parallel/test_forward.py - name: Covering parallel fwi test continue-on-error: true run: | - source /home/olender/Firedrakes/newest3/firedrake/bin/activate + source /home/olender/firedrakes/2024_07_19/firedrake/bin/activate mpiexec -n 6 pytest --cov-report=xml --cov-append --cov=spyro test_parallel/test_fwi.py - # - name: Running serial tests for adjoint - # run: | - # source /home/olender/Firedrakes/main/firedrake/bin/activate - # pytest -n 10 --cov-report=xml --cov-append --cov=spyro test_ad/ - # - name: Running parallel tests - # run: | - # source /home/olender/Firedrakes/main/firedrake/bin/activate - # cp /home/olender/Testing_files/velocity_models/* velocity_models/ - # cp /home/olender/Testing_files/meshes/* meshes/ - # mpiexec -n 10 pytest test_parallel/test_forward.py - # - name: Covering parallel tests - # continue-on-error: true - # run: | - # source /home/olender/Firedrakes/main/firedrake/bin/activate - # cp /home/olender/Testing_files/velocity_models/* velocity_models/ - # cp /home/olender/Testing_files/meshes/* meshes/ - # mpiexec -n 10 pytest --cov-report=xml --cov-append --cov=spyro test_parallel/test_forward.py - # - name: Covering parallel 3D forward test - # continue-on-error: true - # run: | - # source /home/olender/Firedrakes/main/firedrake/bin/activate - # cp /home/olender/Testing_files/velocity_models/* velocity_models/ - # cp /home/olender/Testing_files/meshes/* meshes/ - # mpiexec -n 10 pytest --cov-report=xml --cov-append --cov=spyro test_3d/test_forward_3d.py - name: Uploading coverage to Codecov run: export CODECOV_TOKEN="057ec853-d7ea-4277-819b-0c5ea2f9ff57" && bash <(curl -s https://codecov.io/bash) diff --git a/test/test_MMS.py b/test/test_MMS.py index e1cafbc1..cb1da365 100644 --- a/test/test_MMS.py +++ b/test/test_MMS.py @@ -4,7 +4,7 @@ from firedrake import * import spyro -from model import dictionary as model +from .model import dictionary as model model["acquisition"]["source_type"] = "MMS" From f6c4b2bcd3579a06d184df97821908b2d9ad1902 Mon Sep 17 00:00:00 2001 From: olender Date: Fri, 19 Jul 2024 18:50:44 -0300 Subject: [PATCH 10/11] cleanup --- check_mlt.py | 28 - check_sem.py | 55 - check_u.py | 70 - check_u2.py | 113 - check_u3.py | 138 - notebook_tutorials/meshing.html | 8976 ----------------- .../premade_useful_examples.html | 8177 --------------- notebook_tutorials/simple_forward.html | 8120 --------------- .../simple_forward_with_overthrust.html | 8125 --------------- simple_mms.py | 145 - temp_MMS.py | 28 - velocity_models/tutorial | Bin 518568 -> 0 bytes 12 files changed, 33975 deletions(-) delete mode 100644 check_mlt.py delete mode 100644 check_sem.py delete mode 100644 check_u.py delete mode 100644 check_u2.py delete mode 100644 check_u3.py delete mode 100644 notebook_tutorials/meshing.html delete mode 100644 notebook_tutorials/premade_useful_examples.html delete mode 100644 notebook_tutorials/simple_forward.html delete mode 100644 notebook_tutorials/simple_forward_with_overthrust.html delete mode 100644 simple_mms.py delete mode 100644 temp_MMS.py delete mode 100644 velocity_models/tutorial diff --git a/check_mlt.py b/check_mlt.py deleted file mode 100644 index 623dc913..00000000 --- a/check_mlt.py +++ /dev/null @@ -1,28 +0,0 @@ -import finat -from firedrake import * -import numpy as np - - -def isDiag(M): - i, j = np.nonzero(M) - return np.all(i == j) - -degree = 4 -mesh = RectangleMesh(3, 2, 2.0, 1.0) -element = FiniteElement( # noqa: F405 - "KMV", mesh.ufl_cell(), degree=degree, variant="KMV" -) -V = FunctionSpace(mesh, element) -quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") - -u = TrialFunction(V) -v = TestFunction(V) - -form = u*v*dx(scheme=quad_rule) -A = assemble(form) -M = A.M.values -Mdiag = M.diagonal() - -print(f"Matrix is diagonal:{isDiag(M)}") -np.save("new_diag", Mdiag) -print("END") \ No newline at end of file diff --git a/check_sem.py b/check_sem.py deleted file mode 100644 index d2645690..00000000 --- a/check_sem.py +++ /dev/null @@ -1,55 +0,0 @@ -import finat -import FIAT -from firedrake import * -import numpy as np - - -def gauss_lobatto_legendre_line_rule(degree): - fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule - fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) - finat_ps = finat.point_set.GaussLobattoLegendrePointSet - points = finat_ps(fiat_rule.get_points()) - weights = fiat_rule.get_weights() - return finat.quadrature.QuadratureRule(points, weights) - -def gauss_lobatto_legendre_cube_rule(dimension, degree): - make_tensor_rule = finat.quadrature.TensorProductQuadratureRule - result = gauss_lobatto_legendre_line_rule(degree) - for _ in range(1, dimension): - line_rule = gauss_lobatto_legendre_line_rule(degree) - result = make_tensor_rule([result, line_rule]) - return result - - - -def isDiag(M): - i, j = np.nonzero(M) - return np.all(i == j) - -degree = 2 -# mesh = RectangleMesh(3, 2, 2.0, 1.0, quadrilateral=True) -mesh = RectangleMesh(1, 1, 1.0, 1.0, quadrilateral=True) -element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral') -V = FunctionSpace(mesh, element) -quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree) - -u = TrialFunction(V) -v = TestFunction(V) - -form = u*v*dx(scheme=quad_rule) -A = assemble(form) -M = A.M.values -Mdiag = M.diagonal() - -x_mesh, y_mesh = SpatialCoordinate(mesh) -x_func = Function(V) -y_func = Function(V) -x_func.interpolate(x_mesh) -y_func.interpolate(y_mesh) -x = x_func.dat.data[:] -y = y_func.dat.data[:] - -print(f"Matrix is diagonal:{isDiag(M)}") -old_diag = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_sem_diag.npy") -dif = Mdiag-old_diag -print("END") \ No newline at end of file diff --git a/check_u.py b/check_u.py deleted file mode 100644 index e3f3eda2..00000000 --- a/check_u.py +++ /dev/null @@ -1,70 +0,0 @@ -import finat -from firedrake import * -import numpy as np - - -def analytical_solution(t, V, mesh_z, mesh_x): - analytical = Function(V) - x = mesh_z - y = mesh_x - analytical.interpolate(x * (x + 1) * y * (y - 1) * t) - - return analytical - - -def isDiag(M): - i, j = np.nonzero(M) - return np.all(i == j) - -degree = 4 -mesh = RectangleMesh(3, 2, 2.0, 1.0) -mesh.coordinates.dat.data[:, 0] *= -1.0 -mesh_z, mesh_x = SpatialCoordinate(mesh) -V = FunctionSpace(mesh, "KMV", degree) -quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") - -u = TrialFunction(V) -v = TestFunction(V) - -c = Function(V, name="velocity") -c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) -u_n = Function(V) -u_nm1 = Function(V) -dt = 0.0005 -t = 0.0 -u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) -u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) - -q = Function(V) -q.assign(1.0) -dt = 0.001 - -m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) - * v - * dx(scheme=quad_rule) -) -a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) -le = q * v * dx(scheme=quad_rule) - -form = m1 + a - le - -boundary_ids = (1, 2, 3, 4) -bcs = DirichletBC(V, 0.0, boundary_ids) -A = assemble(lhs(form), bcs=bcs) -solver_parameters = { - "ksp_type": "preonly", - "pc_type": "jacobi", -} -solver = LinearSolver( - A, solver_parameters=solver_parameters -) -M = A.M.values -Mdiag = M.diagonal() - -print(f"Matrix is diagonal:{isDiag(M)}") -np.save("/home/olender/Development/issue_50_compatibility/spyro-1/new_diag", Mdiag) - -print("END") \ No newline at end of file diff --git a/check_u2.py b/check_u2.py deleted file mode 100644 index 30ff1ad1..00000000 --- a/check_u2.py +++ /dev/null @@ -1,113 +0,0 @@ -import finat -from firedrake import * -import numpy as np - - -def analytical_solution(t, V, mesh_z, mesh_x): - analytical = Function(V) - x = mesh_z - y = mesh_x - analytical.interpolate(x * (x + 1) * y * (y - 1) * t) - - return analytical - - -def isDiag(M): - i, j = np.nonzero(M) - return np.all(i == j) - -degree = 4 -mesh = RectangleMesh(50, 50, 1.0, 1.0) -mesh.coordinates.dat.data[:, 0] *= -1.0 -mesh_z, mesh_x = SpatialCoordinate(mesh) -V = FunctionSpace(mesh, "KMV", degree) -quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") - -u = TrialFunction(V) -v = TestFunction(V) - -c = Function(V, name="velocity") -c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) -u_n = Function(V) -u_nm1 = Function(V) -dt = 0.0005 -t = 0.0 -final_time = 1.0 -u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) -u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) - -q_xy = Function(V) -q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x) -q = q_xy * Constant(2 * t) - -nt = int((final_time - t) / dt) + 1 - -m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) - * v - * dx(scheme=quad_rule) -) -a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) -le = q * v * dx(scheme=quad_rule) - -form = m1 + a - le - -B = Cofunction(V.dual()) - -boundary_ids = (1, 2, 3, 4) -bcs = DirichletBC(V, 0.0, boundary_ids) -A = assemble(lhs(form), bcs=bcs) -solver_parameters = { - "ksp_type": "preonly", - "pc_type": "jacobi", -} -solver = LinearSolver( - A, solver_parameters=solver_parameters -) -As = solver.A -petsc_matrix = As.petscmat -diagonal = petsc_matrix.getDiagonal() -Mdiag = diagonal.array - -np.save("/home/olender/Development/issue_50_compatibility/spyro-1/new_diag", Mdiag) -out = File("new_firedrake_u2.pvd") - -u_np1 = Function(V) -for step in range(nt): - q = q_xy * Constant(2 * t) - m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) - * v - * dx(scheme=quad_rule) - ) - a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) - le = q * v * dx(scheme=quad_rule) - - form = m1 + a - le - - B = assemble(rhs(form), tensor=B) - - solver.solve(u_np1, B) - - if (step - 1) % 100 == 0: - print(f"Time : {t}") - out.write(u_n) - assert ( - norm(u_n) < 1 - ), "Numerical instability. Try reducing dt or building the \ - mesh differently" - - u_nm1.assign(u_n) - u_n.assign(u_np1) - - t = step * float(dt) - -u_an = analytical_solution(t, V, mesh_z, mesh_x) -error = errornorm(u_n, u_an) -print(f"Error: {error}") - -print("END") \ No newline at end of file diff --git a/check_u3.py b/check_u3.py deleted file mode 100644 index ed72098e..00000000 --- a/check_u3.py +++ /dev/null @@ -1,138 +0,0 @@ -import finat -import FIAT -from firedrake import * -import numpy as np - - -def gauss_lobatto_legendre_line_rule(degree): - fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule - fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) - finat_ps = finat.point_set.GaussLobattoLegendrePointSet - points = finat_ps(fiat_rule.get_points()) - weights = fiat_rule.get_weights() - return finat.quadrature.QuadratureRule(points, weights) - -def gauss_lobatto_legendre_cube_rule(dimension, degree): - make_tensor_rule = finat.quadrature.TensorProductQuadratureRule - result = gauss_lobatto_legendre_line_rule(degree) - for _ in range(1, dimension): - line_rule = gauss_lobatto_legendre_line_rule(degree) - result = make_tensor_rule([result, line_rule]) - return result - - - -def analytical_solution(t, V, mesh_z, mesh_x): - analytical = Function(V) - x = mesh_z - y = mesh_x - analytical.interpolate(x * (x + 1) * y * (y - 1) * t) - - return analytical - - -def isDiag(M): - i, j = np.nonzero(M) - return np.all(i == j) - -degree = 4 -mesh = RectangleMesh(50, 50, 1.0, 1.0, quadrilateral=True) -mesh.coordinates.dat.data[:, 0] *= -1.0 -mesh_z, mesh_x = SpatialCoordinate(mesh) -element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral') -V = FunctionSpace(mesh, element) -quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree) - -u = TrialFunction(V) -v = TestFunction(V) - -c = Function(V, name="velocity") -c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) -u_n = Function(V) -u_nm1 = Function(V) -dt = 0.0005 -t = 0.0 -final_time = 1.0 -u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) -u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) - -q_xy = Function(V) -q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x) -q = q_xy * Constant(2 * t) - -nt = int((final_time - t) / dt) + 1 - -m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) - * v - * dx(scheme=quad_rule) -) -a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) -le = q * v * dx(scheme=quad_rule) - -form = m1 + a - le - -B = Cofunction(V.dual()) - -boundary_ids = (1, 2, 3, 4) -bcs = DirichletBC(V, 0.0, boundary_ids) -A = assemble(lhs(form), bcs=bcs) -solver_parameters = { - "ksp_type": "preonly", - "pc_type": "jacobi", -} -solver = LinearSolver( - A, solver_parameters=solver_parameters -) -As = solver.A -petsc_matrix = As.petscmat -diagonal = petsc_matrix.getDiagonal() -Mdiag = diagonal.array - -# old_Mdiag = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_diag.npy") -# old_qxy = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_qxy.npy") # OK -# old_un = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_un.npy") # OK -# old_unm1 = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_unm1.npy") # OK -# old_c = np.load("/home/olender/Development/issue_50_compatibility/spyro-1/old_c.npy") # OK - -# out = File("new_firedrake_u2.pvd") - -u_np1 = Function(V) -for step in range(nt): - q = q_xy * Constant(2 * t) - m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) - * v - * dx(scheme=quad_rule) - ) - a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) - le = q * v * dx(scheme=quad_rule) - - form = m1 + a - le - - B = assemble(rhs(form), tensor=B) - - solver.solve(u_np1, B) - - if (step - 1) % 100 == 0: - print(f"Time : {t}") - # out.write(u_n) - assert ( - norm(u_n) < 1 - ), "Numerical instability. Try reducing dt or building the \ - mesh differently" - - u_nm1.assign(u_n) - u_n.assign(u_np1) - - t = step * float(dt) - -u_an = analytical_solution(t, V, mesh_z, mesh_x) -error = errornorm(u_n, u_an) -print(f"Error: {error}") - -print("END") \ No newline at end of file diff --git a/notebook_tutorials/meshing.html b/notebook_tutorials/meshing.html deleted file mode 100644 index fe948904..00000000 --- a/notebook_tutorials/meshing.html +++ /dev/null @@ -1,8976 +0,0 @@ - - - - - -meshing - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - diff --git a/notebook_tutorials/premade_useful_examples.html b/notebook_tutorials/premade_useful_examples.html deleted file mode 100644 index fe6b0895..00000000 --- a/notebook_tutorials/premade_useful_examples.html +++ /dev/null @@ -1,8177 +0,0 @@ - - - - - -premade_useful_examples - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - -
- - diff --git a/notebook_tutorials/simple_forward.html b/notebook_tutorials/simple_forward.html deleted file mode 100644 index a71fb9aa..00000000 --- a/notebook_tutorials/simple_forward.html +++ /dev/null @@ -1,8120 +0,0 @@ - - - - - -simple_forward - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - -
- - diff --git a/notebook_tutorials/simple_forward_with_overthrust.html b/notebook_tutorials/simple_forward_with_overthrust.html deleted file mode 100644 index d042fc35..00000000 --- a/notebook_tutorials/simple_forward_with_overthrust.html +++ /dev/null @@ -1,8125 +0,0 @@ - - - - - -simple_forward_with_overthrust - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - -
- - diff --git a/simple_mms.py b/simple_mms.py deleted file mode 100644 index de70a538..00000000 --- a/simple_mms.py +++ /dev/null @@ -1,145 +0,0 @@ -import firedrake as fire -from firedrake import dot, grad, sin, pi, dx -import FIAT -import finat -import math - - -def gauss_lobatto_legendre_line_rule(degree): - fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule - fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1) - finat_ps = finat.point_set.GaussLobattoLegendrePointSet - points = finat_ps(fiat_rule.get_points()) - weights = fiat_rule.get_weights() - return finat.quadrature.QuadratureRule(points, weights) - - -def gauss_lobatto_legendre_cube_rule(dimension, degree): - make_tensor_rule = finat.quadrature.TensorProductQuadratureRule - result = gauss_lobatto_legendre_line_rule(degree) - for _ in range(1, dimension): - line_rule = gauss_lobatto_legendre_line_rule(degree) - result = make_tensor_rule([result, line_rule]) - return result - - -def analytical_solution(t, V, mesh_z, mesh_x): - analytical = fire.Function(V) - x = mesh_z - y = mesh_x - analytical.interpolate(x * (x + 1) * y * (y - 1) * t) - - return analytical - - -quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=4) - -length_z = 1.0 -length_x = 1.0 - -mesh_dx = 0.02 -nz = int(length_z / mesh_dx) -nx = int(length_x / mesh_dx) -mesh = fire.RectangleMesh( - nz, - nx, - length_z, - length_x, - quadrilateral=True -) -mesh_z, mesh_x = fire.SpatialCoordinate(mesh) - -output = fire.File("debug.pvd") - -element = fire.FiniteElement( # noqa: F405 - "CG", mesh.ufl_cell(), degree=4, variant="spectral" -) -V = fire.FunctionSpace(mesh, element) - -X = fire.Function(V) -u_nm1 = fire.Function(V) -u_n = fire.Function(V) - -# Setting C -c = fire.Function(V, name="velocity") -c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x)) - -# setting xy part of source -q_xy = fire.Function(V) -xy = fire.project((-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x), V) -q_xy.assign(xy) - -# Starting time integration -t = 0.0 -final_time = 1.0 -dt = 0.0001 -nt = int((final_time - t) / dt) + 1 -u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x)) -u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x)) -u_np1 = fire.Function(V, name="pressure t +dt") -u = fire.TrialFunction(V) -v = fire.TestFunction(V) - -m1 = ( - (1 / (c * c)) - * ((u - 2.0 * u_n + u_nm1) / fire.Constant(dt**2)) - * v - * dx(scheme=quad_rule) -) -a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) # explicit - -B = fire.Function(V) - -form = m1 + a -lhs = fire.lhs(form) -rhs = fire.rhs(form) - -A = fire.assemble(lhs, mat_type="matfree") -solver_parameters = { - "ksp_type": "preonly", - "pc_type": "jacobi", -} -solver = fire.LinearSolver( - A, solver_parameters=solver_parameters -) - -for step in range(nt): - q = q_xy * fire.Constant(2 * t) - m1 = ( - 1 - / (c * c) - * ((u - 2.0 * u_n + u_nm1) / fire.Constant(dt**2)) - * v - * dx(scheme=quad_rule) - ) - a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) - le = q * v * dx(scheme=quad_rule) - - form = m1 + a - le - rhs = fire.rhs(form) - - B = fire.assemble(rhs, tensor=B) - - solver.solve(X, B) - - u_np1.assign(X) - - if (step - 1) % 100 == 0: - assert ( - fire.norm(u_n) < 1 - ), "Numerical instability. Try reducing dt or building the \ - mesh differently" - output.write(u_n, time=t, name="Pressure") - if t > 0: - print(f"Simulation time is: {t:{10}.{4}} seconds", flush=True) - - u_nm1.assign(u_n) - u_n.assign(u_np1) - - t = step * float(dt) - -last_analytical_sol = analytical_solution(t, V, mesh_z, mesh_x) -error = fire.errornorm(u_n, last_analytical_sol) -test = math.isclose(error, 0.0, abs_tol=1e-7) - -print(F"Test is {test}") diff --git a/temp_MMS.py b/temp_MMS.py deleted file mode 100644 index 241b78c3..00000000 --- a/temp_MMS.py +++ /dev/null @@ -1,28 +0,0 @@ -import math -from copy import deepcopy -from firedrake import * -import spyro - -from test.model import dictionary as model -model["acquisition"]["source_type"] = "MMS" -model["options"]["cell_type"] = "Q" -model["options"]["variant"] = "lumped" - - - -def run_solve(model): - testmodel = deepcopy(model) - - Wave_obj = spyro.AcousticWaveMMS(dictionary=testmodel) - Wave_obj.set_mesh(mesh_parameters={"dx": 0.02}) - Wave_obj.set_initial_velocity_model(expression="1 + sin(pi*-z)*sin(pi*x)") - Wave_obj.forward_solve() - - u_an = Wave_obj.analytical - u_num = Wave_obj.u_n - - return errornorm(u_num, u_an) - - -error = run_solve(model) -print(error) diff --git a/velocity_models/tutorial b/velocity_models/tutorial deleted file mode 100644 index 81a389f542af76ce23026353c707775482c0aace..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 518568 zcmeI*e~g@UodEgBb^E~g@^ZhZ??`UJXv1!BPSYu|_%yVb&dhYD$|7+yAj5cN(8@4va zw~lTb+b}WKy3yEh|ASk`C#N=zHO3yA+OTEfHGdy{U}|D=e6%rs|AS+jHf$XoYpcL$ zW45v5C(Bo_{`QrPOV%|0d-)}+zWtpw%YWGT*B`A}ec2DMJgaTnPHv7ip6g_@j<*?Y zqYoJCM#=e@M9N z7blmds~U&6dJp<_43reNuumJ#fH{{qz`;Hq17QPOQpGsnoZGK+Zsc-Hid)#HnsN3M%88)ybRV<_f<_nAw_Hm3b zbAH$^Asxu(4ElT=D1Y$J`3Kv_9B5^5Ztkw1KRLj{J{|)c;K0Go6aF?F?F0^RK-PyN z>j&G<90=gRXFs{(VEf45FkfJ7ppRn%nK8q52{9%(e?0Qv|K$T^d++(R zDGqiFa3DZ?z`UXMfc8M86oqULaDW3G-~b0WzyS_$fCC)hfcXLraDW3G-~b0WpkKfN z4sd`29N+*4IKTl8aDW3G-~b1V6LNq99N+*4IKTn@0uFG1103K02ROh14sd`29N+*4 zIA9$U4sd`2<_kE$0S@RFaDW3G-~b0WzyS_$fCC)h00%h00qdA>fCC&bU%&wlaKJhy z9N+*4j1O>t103K02ROh14sd`29N+*4IKTn(1svc22ROh14sbxffCC)h00%h00S<70 z103K02ROh14j3oo00%h00S<701NsFV-~b0WzyS_$fCC)h00%h00S<7$Iwl<800+z$ zaDW3G&@bQs2ROh14sd`29N+*4IKTl8aDW5WG2s9QIAFej103Lhbxb(G0S*`+-~b0W zzyS_$fCC)h00%h00S<701Lg}jzyS_$fCC)hfPMi7IKTl8aDW3G-~b0WzyS_$fCC&b zPRIccaDW3G-~b2o3pl_54sd`29N+*4IKTl8aDW3G;DB{ZIKTl8m@nV}2RNW#zyS_$ zfCC)h00%h00S<70103K02drbl0S<7$d;teIzya%+aDW3GFh0Nm4sd`29N+*4IKTl8 zaDW3G-~b2A7jS?B9N+*4IKTn@0uFG1103K02ROh14sd`29N+*4IAENR103K02ROh1 z4(JzffCC)h00%h00S<70103K02ROh1>zHtW0~|14zyS_$K)-+k9N+*4IKTl8aDW3G z-~b0WzyS_e$AkkM;DGr84sd`2)-mA#2RLASfCC)h00%h00S<70103K02ROh14wx_C z00%h00S<701NsFV-~b0WzyS_$fCC)h00%h00S<7$I3Wi(zyS_$fCC)RFW>+NIKTl8 zaDW3G-~b0WzyS_$fCJVs;Q$9XV7`C@9N>U{0S7q10S<70103K02ROh14sd`29I%cF z2ROh1^93B>00*pN!T}C&!1w?MIKTl8aDW3G-~b0WzyS_$fCC&bU%&wlaDW3G-~b2o z3pl_54sd`29N+*4IKTl8aDW3G;DB*L4sd`29N+*4IG|s^0S<70103K02ROh14sd`2 z9N+*4tYg9f4sgJH0S7q10sR6FaDW3G-~b0WzyS_$fCC)h00%f=9TN_4fCJ_WIKTl8 zSjU6|9N>WQ0S<70103K02ROh14sd`29N+*4IAFej103K02ROh14(JzffCC)h00%h0 z0S<70103K02ROh1+NIKTl8aDW3G-~b0WzyS_$ zfCI(}Iluu9aDW3G;DCMs2ROh14sd`29N+*4IKTl8aDW3Gu#O1_IKTn(1svc22lNX# zzyS_$fCC)h00%h00S<70103Lhbxb(G0S=fi-~b0WU>y?00%h00S<6LzkmZA-~b0WzyS_$fCC)h00%h00S*`^U;OgO*+4wx_C00%grU%&wlaDW3G-~b0WzyS_$fCC)h00*pN!T}C&z0S<70103K02ROh14sd`29N+*4%olKg103K02ROh1{Q?effCC)h z00%h00S<70103K02RLAykOLgxz#uzt?H4ad51#wR=8_Y(HJ4m*NOnyxOdXwF=k7O0 zW{*r;&Kplx@4P+rIS<2Hex)G{P*}>^;1REbc>HgC>4o(8=I=F^eEF(ZGQZz~SKZGX zc=!+RNE@#BP`Y*Higd$|w-0L``}Mj5M6=Q=gfw7K#|NG{?Xg$oeIWzp+6XT8QTeYCk`>DKxA?&f}a{?K#Nttb6?nDW@K z*IH=^0~8k0fk9ixj%H2OeK@WA4l^GOw6+f8e(C>-&!BM?Qw7%mC4>v)?Im436s$XZNAq-GhNC)b^j>#=sPHirE-+g5{F!#La?!%TO5a?*AEU)kKD0$E~{n!Oi%sOzYR++tA3r8hA=>3AswjmIwl);t|~M8^Pf9S zEj>Nm@%7c|=gUW`&dMN{7a*FIRw1MVb)7f-{afChwqJJV{H&IJKkb~lfdj>k^Hdtb z0EKa&pyOxt1G>0$O1HPZyS(!Jj}P)Wy2{rUAexm{fdd6OkhRO7 z{>fzv<+SwQjU3lZe|giLDf=En!0+zWu}+nSFhF4(D9C{uTjPT}PkmSEnP2QP1&=QM zM!NQ_Q|p+yh;UKcXMQBML!%!FQ46SKO;1>B&#*-1qsyb>~i`wPPnmJipehP647>X%#{`P`CMl%z;0B7n&g<_ikgm;1+9Dh*+P!Z?ug*!ktocN>OGyL~83<`OCs~)LUpE|a~qt&_65C$lW z13Atit5;rRh-Rf#2F69p8?tXH*A-X4%B7d@Vd`_H9fles6}(xWt+aJ=Sh<%9l?Tp7I0vFjCIiK(67Dr zgGWZ1FPL*AZG6oCQI_-g;Q$9ZK9{aqHlC)w_SU)lMfsVgPqxMh2fCsxWWS&;pC9hVRz{VjRsL6ambICFH9LQdQ-xta`&yZR^4Qf8C zdR+q$?XF%zNC#@Nj!8rV{$8gyTwGON2DO~7LNovyhjgGO9Q zp!Z{DE_cy>a-jEpI@gZqDntXYaYzTu7qr$x$TtQ#cRc&qn+XFD?XF$|2YNeBvV7&u zYMGCh|K9g&_TCNKg#*1GFJzL%tThzyMAb?{` z&~XN)Po8WaqmOVqO?~aHY1Ojvj$eb%6Z$-13qA9OjRlW`j&ii|v5^SbZ?8jo>~oi; zYd?5oNA+jGqWfNZKzpF$3{9SNM0*A!+)j_KJ}O=J*{^o|8d$KP1sn+AKuzWg_GRs^ z(>FIqHoZ530a+c=#K$gTK#n%veiu2wfsS*u>Z0GIiSM2g#ejVd?D^fUj_a>2uwYvt zd;f;*VtgRSx#Q{`7iRtOC=R4u?>si$yZST^vFEzcIXnypAm4tid#<~!e@F-FG+(gsXX~3IC;lje1Fd?bR(%+d&w%Xx z8@3Av@*gwl`qOVn6L0%ybIC=`etM9#OZK%->N2(xxSkwXA7i{Q@Gz~N=qMr@wX1{WP1Kk`iY3DTuilb-y<&i**CbZa&~!GT60 z24wqx*e)T*3G1?sN&C6?_Nk-O%jZ1O7Y9aWR-}8D`TS7({>XWo0|D9tb(%NKj<;*i zIyF7I^c#IKAR8atdFs2;jXyZNt9DGzoD6zd90=e*P5wWkSF(`J7yRW-cQ!|k>&v`h zdip1qrJpY!d8Ixwp1Xg*0Swf&c{2aZm8OUwEeO>)9>H%?w;z*0S@d3FXT93-TDPto-zmS`1bA*AEVFE_vU5WqXhwxnpl?l~r^GJYQqj7>My1;_Uq2 zey!Vl!F&(yK5SXp24v2pTeh4kv&FXsviEe@F0aa8dlvh7yC&wW_V=AE1S@FAniVH zTl(#?mpZFMdFRVPFP{S(DC`)&dFzJs^u^yV%7HX<{f>0=vo{r1zuH?{o^u=s;6UBy z3)=Jb(~G~KCbxelzXc~RN%yb#1p~^eu54msOdeWcM^M{_B&jH^n&b=Sn zZSLSU2Ld=S2;+ov$J#BQx;Je<>L!0twhb23(yJ5F0rLeM z=yUt3^xy?6fCI({qMU;x)gxF;ORo+GI1nj+QR-59@PZY<0S-jTUZi>ii)rcA3F$!H z*D>jH-!@X7>ROl5gBPrTkPetHFizNKOM}I<^y-9kpziCK^!XUB>p3`5y-N>XumVCl zVBS!BpwE^Di)rcA;ec_%Nax`wbtye~!3y922cl#zQaysjwDjt5fCG{87o{$x2QOFw zAssMZV0^&%K&4V}Kzks{IXF^1Dh*+P!a_P=zCe3Gd!SMsIG{Zc={y{zE|rEbKw%u< zK$Pr7sz-omR$2uPa3E6tqSU3*5C$kLqyy#)j1L$es8k0IXb(g=2S=($r6CMZSV#xV z7ibS?4^*lH2ebzworj~;rP2@vD2xLfh?2cX^#~BnN~^#D4n)ddl)6+J!T^PZbig_$ z#s`cKRH{Qr2h1CCfCC&bU%&wlaDW3G-~b0WzyS_$fCC)h00*pN!T}C&z0S<70103K02ROh14sd`29N+*4%olKg103K02ROh1{Q?effCC)h00%h0 z0S<70103K02RLAykOLgx00%h00S@RFaDW3G-~b0WzyS_$fCC)h00%h00qdA>fCC&b zU%&wla6rF+103K02ROh14sd`29N+*4IKTl8SjU6|9N>WY0uFG11J*I&00%f=e1HQS z-~b0WzyS_$fCC)h00%h00S=fi-~b0WzyS_$fCKsk9N+*4IKTl8aDW3G-~b0WzyS_$ zz&IfXIKTl8aDW3G&@bQs2ROh14sd`29N+*4IKTl8aDW5WG2s9QIAFej103LhegOwK zzyS_$fCC)h00%h00S<70101l92?sd90rLeM-~b1#W5NLraKQKg2ROh14sd`29N+*4 wIKTl8aDW3GFkip{4sd`29N+*4^b0t^0S<70103K02ROh14sd`2eRkmg0mAntL;wH) From f01cfa2fb0550d897a8756dfd858d7ead885e9d6 Mon Sep 17 00:00:00 2001 From: olender Date: Wed, 24 Jul 2024 10:27:31 -0300 Subject: [PATCH 11/11] removing commented code --- spyro/solvers/acoustic_solver_construction_no_pml.py | 3 --- spyro/solvers/time_integration_central_difference.py | 5 ----- 2 files changed, 8 deletions(-) diff --git a/spyro/solvers/acoustic_solver_construction_no_pml.py b/spyro/solvers/acoustic_solver_construction_no_pml.py index d5a18f3d..4c4882eb 100644 --- a/spyro/solvers/acoustic_solver_construction_no_pml.py +++ b/spyro/solvers/acoustic_solver_construction_no_pml.py @@ -48,9 +48,6 @@ def construct_solver_or_matrix_no_pml(Wave_object): Wave_object.solver = fire.LinearSolver( A, solver_parameters=Wave_object.solver_parameters ) - # lin_var = fire.LinearVariationalProblem(lhs, rhs + B, u_np1) - # solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} - # Wave_object.solver = fire.LinearVariationalSolver(lin_var,solver_parameters=solver_parameters) Wave_object.rhs = rhs Wave_object.B = B diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index 51ef3bb6..2c6bf696 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -245,9 +245,6 @@ def central_difference_MMS(Wave_object, source_id=0): u_np1 = Wave_object.u_np1 u_nm1.assign(Wave_object.analytical_solution(t - 2 * dt)) u_n.assign(Wave_object.analytical_solution(t - dt)) - # u_np1 = fire.Function(Wave_object.function_space, name="pressure t +dt") - # u_nm1.dat.data[:] = np.load("old_u_nm1.npy") - # u_n.dat.data[:] = np.load("old_u_n.npy") u = fire.TrialFunction(Wave_object.function_space) v = fire.TestFunction(Wave_object.function_space) @@ -283,8 +280,6 @@ def central_difference_MMS(Wave_object, source_id=0): Wave_object.solver.solve(u_np1, B) - # u_np1.assign(X) - usol_recv.append( Wave_object.receivers.interpolate(u_np1.dat.data_ro_with_halos[:]) )