From 110720d3588dec8f8d462efdb35000cd93af2f38 Mon Sep 17 00:00:00 2001 From: olender Date: Tue, 11 Jun 2024 13:31:54 -0300 Subject: [PATCH 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 8/8] 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):