From 72c314732d55dd2b6a3474e2067c89c54941a7b5 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 24 Oct 2024 09:35:30 -0300 Subject: [PATCH 01/11] Simplify fwi test --- docs/source/advanced_tut.rst | 2 +- tests/demos/test_demos_run.py | 18 +---- tests/regression/test_fwi_demos.py | 101 +++++++++++++++++++++++++++++ 3 files changed, 105 insertions(+), 16 deletions(-) create mode 100644 tests/regression/test_fwi_demos.py diff --git a/docs/source/advanced_tut.rst b/docs/source/advanced_tut.rst index a6f4642aa1..d9fbf252ba 100644 --- a/docs/source/advanced_tut.rst +++ b/docs/source/advanced_tut.rst @@ -23,4 +23,4 @@ element systems. A pressure-convection-diffusion preconditioner for the Navier-Stokes equations. Rayleigh-Benard convection. Netgen support. - Full-waveform inversion: Full-waveform inversion: spatial and wave sources parallelism. + Full-waveform inversion: spatial and wave sources parallelism. diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index 9094bc99f7..465083ee35 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -30,11 +30,6 @@ "test_extrusion_lsw.py", ] -parallel_demos = [ - "full_waveform_inversion.py", -] - - # Discover the demo files by globbing the demo directory @pytest.fixture(params=glob.glob("%s/*/*.py.rst" % demo_dir), ids=lambda x: basename(x)) @@ -77,6 +72,8 @@ def py_file(rst_file, tmpdir, monkeypatch): @pytest.mark.skipcomplex # Will need to add a seperate case for a complex demo. +@pytest.mark.markif_fixture(pytest.mark.slow, rst_file="linear_wave_equation.py.rst") +@pytest.mark.markif_fixture(pytest.mark.slow, rst_file="full_waveform_inversion.py.rst") def test_demo_runs(py_file, env): # Add pytest skips for missing imports or packages if basename(py_file) in ("stokes.py", "rayleigh-benard.py", "saddle_point_systems.py", "navier_stokes.py", "netgen_mesh.py"): @@ -123,14 +120,5 @@ def test_demo_runs(py_file, env): import vtkmodules.vtkCommonDataModel # noqa: F401 except ImportError: pytest.skip(reason=f"VTK unavailable, skipping {basename(py_file)}") - if basename(py_file) in parallel_demos: - if basename(py_file) == "full_waveform_inversion.py": - processes = 2 - else: - raise NotImplementedError("You need to specify the number of processes for this test") - - executable = ["mpiexec", "-n", str(processes), sys.executable, py_file] - else: - executable = [sys.executable, py_file] - subprocess.check_call(executable, env=env) + subprocess.check_call([sys.executable, py_file], env=env) diff --git a/tests/regression/test_fwi_demos.py b/tests/regression/test_fwi_demos.py new file mode 100644 index 0000000000..ba01a4964a --- /dev/null +++ b/tests/regression/test_fwi_demos.py @@ -0,0 +1,101 @@ + +import finat +import pytest +import numpy as np +from firedrake import * +from firedrake.adjoint import * +from firedrake.__future__ import interpolate + + +def ricker_wavelet(t, fs, amp=1.0): + ts = 1.5 + t0 = t - ts * np.sqrt(6.0) / (np.pi * fs) + return (amp * (1.0 - (1.0 / 2.0) * (2.0 * np.pi * fs) * (2.0 * np.pi * fs) * t0 * t0) + * np.exp((-1.0 / 4.0) * (2.0 * np.pi * fs) * (2.0 * np.pi * fs) * t0 * t0)) + + +def wave_equation_solver(c, source_function, dt, V): + u = TrialFunction(V) + v = TestFunction(V) + u_np1 = Function(V) # timestep n+1 + u_n = Function(V) # timestep n + u_nm1 = Function(V) # timestep n-1 + # Quadrature rule for lumped mass matrix. + quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") + m = (1 / (c * c)) + time_term = m * (u - 2.0 * u_n + u_nm1) / Constant(dt**2) * v * dx(scheme=quad_rule) + nf = (1 / c) * ((u_n - u_nm1) / dt) * v * ds + a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) + F = time_term + a + nf + lin_var = LinearVariationalProblem(lhs(F), rhs(F) + source_function, u_np1) + # Since the linear system matrix is diagonal, the solver parameters are set to construct a solver, + # which applies a single step of Jacobi preconditioning. + solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} + solver = LinearVariationalSolver(lin_var, solver_parameters=solver_parameters) + return solver, u_np1, u_n, u_nm1 + + +@pytest.mark.parallel(nprocs=2) +def test_fwi_demos(): + M = 2 + my_ensemble = Ensemble(COMM_WORLD, M) + num_sources = my_ensemble.ensemble_comm.size + source_number = my_ensemble.ensemble_comm.rank + mesh = UnitSquareMesh(20, 20, comm=my_ensemble.comm) + + source_locations = np.linspace((0.3, 0.1), (0.7, 0.1), num_sources) + receiver_locations = np.linspace((0.2, 0.9), (0.8, 0.9), 10) + dt = 0.01 # time step in seconds + final_time = 0.8 # final time in seconds + frequency_peak = 7.0 # The dominant frequency of the Ricker wavelet in Hz. + + V = FunctionSpace(mesh, "KMV", 1) + x, z = SpatialCoordinate(mesh) + c_true = Function(V).interpolate(1.75 + 0.25 * tanh(200 * (0.125 - sqrt((x - 0.5) ** 2 + (z - 0.5) ** 2)))) + + source_mesh = VertexOnlyMesh(mesh, [source_locations[source_number]]) + + V_s = FunctionSpace(source_mesh, "DG", 0) + + d_s = Function(V_s) + d_s.interpolate(1.0) + source_cofunction = assemble(d_s * TestFunction(V_s) * dx) + q_s = Cofunction(V.dual()).interpolate(source_cofunction) + + receiver_mesh = VertexOnlyMesh(mesh, receiver_locations) + V_r = FunctionSpace(receiver_mesh, "DG", 0) + + true_data_receivers = [] + total_steps = int(final_time / dt) + 1 + f = Cofunction(V.dual()) # Wave equation forcing term. + solver, u_np1, u_n, u_nm1 = wave_equation_solver(c_true, f, dt, V) + interpolate_receivers = interpolate(u_np1, V_r) + + for step in range(total_steps): + f.assign(ricker_wavelet(step * dt, frequency_peak) * q_s) + solver.solve() + u_nm1.assign(u_n) + u_n.assign(u_np1) + true_data_receivers.append(assemble(interpolate_receivers)) + + c_guess = Function(V).interpolate(1.5) + + continue_annotation() + + f = Cofunction(V.dual()) # Wave equation forcing term. + solver, u_np1, u_n, u_nm1 = wave_equation_solver(c_guess, f, dt, V) + interpolate_receivers = interpolate(u_np1, V_r) + J_val = 0.0 + for step in range(total_steps): + f.assign(ricker_wavelet(step * dt, frequency_peak) * q_s) + solver.solve() + u_nm1.assign(u_n) + u_n.assign(u_np1) + guess_receiver = assemble(interpolate_receivers) + misfit = guess_receiver - true_data_receivers[step] + J_val += 0.5 * assemble(inner(misfit, misfit) * dx) + + J_hat = EnsembleReducedFunctional(J_val, Control(c_guess), my_ensemble) + + taylor_test(J_hat, c_guess, Function(V).interpolate(0.1)) + get_working_tape().clear_tape() From 79f4972ac4eb30040a30cecb2e0516b249d64334 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 24 Oct 2024 09:46:15 -0300 Subject: [PATCH 02/11] Remove skip slow test --- tests/demos/test_demos_run.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index 465083ee35..e2494e222a 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -72,8 +72,6 @@ def py_file(rst_file, tmpdir, monkeypatch): @pytest.mark.skipcomplex # Will need to add a seperate case for a complex demo. -@pytest.mark.markif_fixture(pytest.mark.slow, rst_file="linear_wave_equation.py.rst") -@pytest.mark.markif_fixture(pytest.mark.slow, rst_file="full_waveform_inversion.py.rst") def test_demo_runs(py_file, env): # Add pytest skips for missing imports or packages if basename(py_file) in ("stokes.py", "rayleigh-benard.py", "saddle_point_systems.py", "navier_stokes.py", "netgen_mesh.py"): From 982619d21463eb1c2e5bc7759b351d2def89e471 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 24 Oct 2024 10:02:36 -0300 Subject: [PATCH 03/11] Lint --- tests/demos/test_demos_run.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index e2494e222a..3c9ff5436a 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -30,6 +30,7 @@ "test_extrusion_lsw.py", ] + # Discover the demo files by globbing the demo directory @pytest.fixture(params=glob.glob("%s/*/*.py.rst" % demo_dir), ids=lambda x: basename(x)) From 045c47828dc8f6eb3b5b6d0f51d73ee617604a57 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 24 Oct 2024 12:03:47 -0300 Subject: [PATCH 04/11] Fix test_demos_run --- tests/demos/test_demos_run.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index 3c9ff5436a..61717a204c 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -30,6 +30,10 @@ "test_extrusion_lsw.py", ] +parallel_demos = [ + "full_waveform_inversion.py", +] + # Discover the demo files by globbing the demo directory @pytest.fixture(params=glob.glob("%s/*/*.py.rst" % demo_dir), @@ -119,5 +123,9 @@ def test_demo_runs(py_file, env): import vtkmodules.vtkCommonDataModel # noqa: F401 except ImportError: pytest.skip(reason=f"VTK unavailable, skipping {basename(py_file)}") + if basename(py_file) in parallel_demos: + # Skip this test. It is expensive and reproduced in a simpler form + # at test/regression/test_fwi_demos.py + pytest.skip("Skipping parallel full waveform inversion (FWI) test") subprocess.check_call([sys.executable, py_file], env=env) From bd5d8ce9130f2755d2b83a262b835c0c74022c32 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 24 Oct 2024 13:39:41 -0300 Subject: [PATCH 05/11] Skip complex test --- tests/regression/test_fwi_demos.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/regression/test_fwi_demos.py b/tests/regression/test_fwi_demos.py index ba01a4964a..b8f85c8954 100644 --- a/tests/regression/test_fwi_demos.py +++ b/tests/regression/test_fwi_demos.py @@ -35,6 +35,7 @@ def wave_equation_solver(c, source_function, dt, V): return solver, u_np1, u_n, u_nm1 +@pytest.mark.skipcomplex @pytest.mark.parallel(nprocs=2) def test_fwi_demos(): M = 2 From 2906ed62793adcfebaccf565155e43ac75a51279 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Mon, 4 Nov 2024 20:35:15 +0000 Subject: [PATCH 06/11] Add a RUN_MODE enviroment --- .../full_waveform_inversion.py.rst | 28 +++-- tests/demos/test_demos_run.py | 15 ++- tests/regression/test_fwi_demos.py | 102 ------------------ 3 files changed, 33 insertions(+), 112 deletions(-) delete mode 100644 tests/regression/test_fwi_demos.py diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index 8d0b415888..93cf832677 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -72,9 +72,14 @@ gradients (which we will discuss later). Instantiating an ensemble requires a communicator (usually MPI_COMM_WORLD) plus the number of MPI processes to be used in each member of the ensemble (2, in this case):: + import os from firedrake import * M = 2 my_ensemble = Ensemble(COMM_WORLD, M) + test = False + if os.getenv("RUN_MODE") == "test": + # Set the test flag to True. + test = True Each ensemble member will have the same spatial parallelism with the number of ensemble members given by dividing the size of the original communicator by the number processes in each ensemble member. @@ -100,18 +105,27 @@ The source number is defined with the ``Ensemble.ensemble_comm`` rank:: In this example, we consider a two-dimensional square domain with a side length of 1.0 km. The mesh is built over the ``my_ensemble.comm`` (spatial) communicator:: - - Lx, Lz = 1.0, 1.0 - mesh = UnitSquareMesh(80, 80, comm=my_ensemble.comm) + + if test: + # Setup for a faster test execution. + nx, ny = 15, 15 + else: + nx, ny = 80, 80 + mesh = UnitSquareMesh(nx, ny, comm=my_ensemble.comm) The basic input for the FWI problem are defined as follows:: import numpy as np + frequency_peak = 7.0 # The dominant frequency of the Ricker wavelet in Hz. source_locations = np.linspace((0.3, 0.1), (0.7, 0.1), num_sources) receiver_locations = np.linspace((0.2, 0.9), (0.8, 0.9), 20) - dt = 0.002 # time step in seconds - final_time = 1.0 # final time in seconds - frequency_peak = 7.0 # The dominant frequency of the Ricker wavelet in Hz. + if test: + # Setup for a faster test execution. + dt = 0.03 # time step in seconds + final_time = 0.6 # final time in seconds + else: + dt = 0.002 # time step in seconds + final_time = 1.0 # final time in seconds Sources and receivers locations are illustrated in the following figure: @@ -267,6 +281,8 @@ To have the step 4, we need first to tape the forward problem. That is done by c misfit = guess_receiver - true_data_receivers[step] J_val += 0.5 * assemble(inner(misfit, misfit) * dx) + print("Functional value: ", J_val) + We now instantiate :class:`~.EnsembleReducedFunctional`:: J_hat = EnsembleReducedFunctional(J_val, Control(c_guess), my_ensemble) diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index 61717a204c..3c4cb82cf5 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -46,6 +46,8 @@ def rst_file(request): def env(): env = os.environ.copy() env["MPLBACKEND"] = "pdf" + # Set environment variables for the test. + env["RUN_MODE"] = "test" return env @@ -124,8 +126,13 @@ def test_demo_runs(py_file, env): except ImportError: pytest.skip(reason=f"VTK unavailable, skipping {basename(py_file)}") if basename(py_file) in parallel_demos: - # Skip this test. It is expensive and reproduced in a simpler form - # at test/regression/test_fwi_demos.py - pytest.skip("Skipping parallel full waveform inversion (FWI) test") + if basename(py_file) == "full_waveform_inversion.py": + processes = 2 + else: + raise NotImplementedError("You need to specify the number of processes for this test") - subprocess.check_call([sys.executable, py_file], env=env) + executable = ["mpiexec", "-n", str(processes), sys.executable, py_file] + else: + executable = [sys.executable, py_file] + + subprocess.check_call(executable, env=env) diff --git a/tests/regression/test_fwi_demos.py b/tests/regression/test_fwi_demos.py deleted file mode 100644 index b8f85c8954..0000000000 --- a/tests/regression/test_fwi_demos.py +++ /dev/null @@ -1,102 +0,0 @@ - -import finat -import pytest -import numpy as np -from firedrake import * -from firedrake.adjoint import * -from firedrake.__future__ import interpolate - - -def ricker_wavelet(t, fs, amp=1.0): - ts = 1.5 - t0 = t - ts * np.sqrt(6.0) / (np.pi * fs) - return (amp * (1.0 - (1.0 / 2.0) * (2.0 * np.pi * fs) * (2.0 * np.pi * fs) * t0 * t0) - * np.exp((-1.0 / 4.0) * (2.0 * np.pi * fs) * (2.0 * np.pi * fs) * t0 * t0)) - - -def wave_equation_solver(c, source_function, dt, V): - u = TrialFunction(V) - v = TestFunction(V) - u_np1 = Function(V) # timestep n+1 - u_n = Function(V) # timestep n - u_nm1 = Function(V) # timestep n-1 - # Quadrature rule for lumped mass matrix. - quad_rule = finat.quadrature.make_quadrature(V.finat_element.cell, V.ufl_element().degree(), "KMV") - m = (1 / (c * c)) - time_term = m * (u - 2.0 * u_n + u_nm1) / Constant(dt**2) * v * dx(scheme=quad_rule) - nf = (1 / c) * ((u_n - u_nm1) / dt) * v * ds - a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule) - F = time_term + a + nf - lin_var = LinearVariationalProblem(lhs(F), rhs(F) + source_function, u_np1) - # Since the linear system matrix is diagonal, the solver parameters are set to construct a solver, - # which applies a single step of Jacobi preconditioning. - solver_parameters = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"} - solver = LinearVariationalSolver(lin_var, solver_parameters=solver_parameters) - return solver, u_np1, u_n, u_nm1 - - -@pytest.mark.skipcomplex -@pytest.mark.parallel(nprocs=2) -def test_fwi_demos(): - M = 2 - my_ensemble = Ensemble(COMM_WORLD, M) - num_sources = my_ensemble.ensemble_comm.size - source_number = my_ensemble.ensemble_comm.rank - mesh = UnitSquareMesh(20, 20, comm=my_ensemble.comm) - - source_locations = np.linspace((0.3, 0.1), (0.7, 0.1), num_sources) - receiver_locations = np.linspace((0.2, 0.9), (0.8, 0.9), 10) - dt = 0.01 # time step in seconds - final_time = 0.8 # final time in seconds - frequency_peak = 7.0 # The dominant frequency of the Ricker wavelet in Hz. - - V = FunctionSpace(mesh, "KMV", 1) - x, z = SpatialCoordinate(mesh) - c_true = Function(V).interpolate(1.75 + 0.25 * tanh(200 * (0.125 - sqrt((x - 0.5) ** 2 + (z - 0.5) ** 2)))) - - source_mesh = VertexOnlyMesh(mesh, [source_locations[source_number]]) - - V_s = FunctionSpace(source_mesh, "DG", 0) - - d_s = Function(V_s) - d_s.interpolate(1.0) - source_cofunction = assemble(d_s * TestFunction(V_s) * dx) - q_s = Cofunction(V.dual()).interpolate(source_cofunction) - - receiver_mesh = VertexOnlyMesh(mesh, receiver_locations) - V_r = FunctionSpace(receiver_mesh, "DG", 0) - - true_data_receivers = [] - total_steps = int(final_time / dt) + 1 - f = Cofunction(V.dual()) # Wave equation forcing term. - solver, u_np1, u_n, u_nm1 = wave_equation_solver(c_true, f, dt, V) - interpolate_receivers = interpolate(u_np1, V_r) - - for step in range(total_steps): - f.assign(ricker_wavelet(step * dt, frequency_peak) * q_s) - solver.solve() - u_nm1.assign(u_n) - u_n.assign(u_np1) - true_data_receivers.append(assemble(interpolate_receivers)) - - c_guess = Function(V).interpolate(1.5) - - continue_annotation() - - f = Cofunction(V.dual()) # Wave equation forcing term. - solver, u_np1, u_n, u_nm1 = wave_equation_solver(c_guess, f, dt, V) - interpolate_receivers = interpolate(u_np1, V_r) - J_val = 0.0 - for step in range(total_steps): - f.assign(ricker_wavelet(step * dt, frequency_peak) * q_s) - solver.solve() - u_nm1.assign(u_n) - u_n.assign(u_np1) - guess_receiver = assemble(interpolate_receivers) - misfit = guess_receiver - true_data_receivers[step] - J_val += 0.5 * assemble(inner(misfit, misfit) * dx) - - J_hat = EnsembleReducedFunctional(J_val, Control(c_guess), my_ensemble) - - taylor_test(J_hat, c_guess, Function(V).interpolate(0.1)) - get_working_tape().clear_tape() From c151665d6d468b50749e140848dc9976d2f4e0d2 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Wed, 6 Nov 2024 12:59:20 +0000 Subject: [PATCH 07/11] Add FIREDRAKE_CI_TESTS --- .../full_waveform_inversion.py.rst | 41 +++++++++++-------- tests/demos/test_demos_run.py | 2 - 2 files changed, 23 insertions(+), 20 deletions(-) diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index 93cf832677..54893aec50 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -69,6 +69,26 @@ To achieve this, we use ensemble parallelism, which involves solving simultaneou equation (3) with different forcing terms :math:`f_s(\mathbf{x}, t)`, different :math:`J_s` and their gradients (which we will discuss later). +We start by importing firedrake and the necessary parameters used in spatial and time execution:: + + from firedrake import * + import os + test = False + if os.getenv("FIREDRAKE_CI_TESTS") == 1: + test = True + + def spatial_time_parameters(test): + if test: + # Setup for a faster test execution. + dt = 0.03 # time step in seconds + final_time = 0.6 # final time in seconds + nx, ny = 15, 15 + else: + dt = 0.002 # time step in seconds + final_time = 1.0 # final time in seconds + nx, ny = 80, 80 + return dt, final_time, nx, ny + Instantiating an ensemble requires a communicator (usually MPI_COMM_WORLD) plus the number of MPI processes to be used in each member of the ensemble (2, in this case):: @@ -76,10 +96,7 @@ processes to be used in each member of the ensemble (2, in this case):: from firedrake import * M = 2 my_ensemble = Ensemble(COMM_WORLD, M) - test = False - if os.getenv("RUN_MODE") == "test": - # Set the test flag to True. - test = True + Each ensemble member will have the same spatial parallelism with the number of ensemble members given by dividing the size of the original communicator by the number processes in each ensemble member. @@ -106,26 +123,15 @@ The source number is defined with the ``Ensemble.ensemble_comm`` rank:: In this example, we consider a two-dimensional square domain with a side length of 1.0 km. The mesh is built over the ``my_ensemble.comm`` (spatial) communicator:: - if test: - # Setup for a faster test execution. - nx, ny = 15, 15 - else: - nx, ny = 80, 80 + dt, final_time, nx, ny = spatial_time_parameters(test) mesh = UnitSquareMesh(nx, ny, comm=my_ensemble.comm) -The basic input for the FWI problem are defined as follows:: +The frequency of the Ricker wavelet, the source and receiver locations are defined as follows:: import numpy as np frequency_peak = 7.0 # The dominant frequency of the Ricker wavelet in Hz. source_locations = np.linspace((0.3, 0.1), (0.7, 0.1), num_sources) receiver_locations = np.linspace((0.2, 0.9), (0.8, 0.9), 20) - if test: - # Setup for a faster test execution. - dt = 0.03 # time step in seconds - final_time = 0.6 # final time in seconds - else: - dt = 0.002 # time step in seconds - final_time = 1.0 # final time in seconds Sources and receivers locations are illustrated in the following figure: @@ -281,7 +287,6 @@ To have the step 4, we need first to tape the forward problem. That is done by c misfit = guess_receiver - true_data_receivers[step] J_val += 0.5 * assemble(inner(misfit, misfit) * dx) - print("Functional value: ", J_val) We now instantiate :class:`~.EnsembleReducedFunctional`:: diff --git a/tests/demos/test_demos_run.py b/tests/demos/test_demos_run.py index 3c4cb82cf5..9094bc99f7 100644 --- a/tests/demos/test_demos_run.py +++ b/tests/demos/test_demos_run.py @@ -46,8 +46,6 @@ def rst_file(request): def env(): env = os.environ.copy() env["MPLBACKEND"] = "pdf" - # Set environment variables for the test. - env["RUN_MODE"] = "test" return env From 7c964b16b603b5421c44c719b74abdb759bc2269 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Wed, 6 Nov 2024 13:08:20 +0000 Subject: [PATCH 08/11] minor changes --- .../full_waveform_inversion/full_waveform_inversion.py.rst | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index 54893aec50..82f68d0a8e 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -69,7 +69,7 @@ To achieve this, we use ensemble parallelism, which involves solving simultaneou equation (3) with different forcing terms :math:`f_s(\mathbf{x}, t)`, different :math:`J_s` and their gradients (which we will discuss later). -We start by importing firedrake and the necessary parameters used in spatial and time execution:: +We start by importing firedrake and the necessary parameters used in spatial and time executions:: from firedrake import * import os @@ -92,11 +92,9 @@ We start by importing firedrake and the necessary parameters used in spatial an Instantiating an ensemble requires a communicator (usually MPI_COMM_WORLD) plus the number of MPI processes to be used in each member of the ensemble (2, in this case):: - import os from firedrake import * M = 2 my_ensemble = Ensemble(COMM_WORLD, M) - Each ensemble member will have the same spatial parallelism with the number of ensemble members given by dividing the size of the original communicator by the number processes in each ensemble member. @@ -121,7 +119,7 @@ The source number is defined with the ``Ensemble.ensemble_comm`` rank:: source_number = my_ensemble.ensemble_comm.rank In this example, we consider a two-dimensional square domain with a side length of 1.0 km. The mesh is -built over the ``my_ensemble.comm`` (spatial) communicator:: +built over the ``my_ensemble.comm`` (spatial) communicator.:: dt, final_time, nx, ny = spatial_time_parameters(test) mesh = UnitSquareMesh(nx, ny, comm=my_ensemble.comm) @@ -287,7 +285,6 @@ To have the step 4, we need first to tape the forward problem. That is done by c misfit = guess_receiver - true_data_receivers[step] J_val += 0.5 * assemble(inner(misfit, misfit) * dx) - We now instantiate :class:`~.EnsembleReducedFunctional`:: J_hat = EnsembleReducedFunctional(J_val, Control(c_guess), my_ensemble) From 9f03e38fc0a936a583605cfb8411872a52ee881c Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Wed, 6 Nov 2024 14:13:38 +0000 Subject: [PATCH 09/11] fixing --- demos/full_waveform_inversion/full_waveform_inversion.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index 82f68d0a8e..bcd69c8b8f 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -74,7 +74,7 @@ We start by importing firedrake and the necessary parameters used in spatial an from firedrake import * import os test = False - if os.getenv("FIREDRAKE_CI_TESTS") == 1: + if os.getenv("FIREDRAKE_CI_TESTS") == "1": test = True def spatial_time_parameters(test): From 0fda2bd338b379f88b499c58c2f30c881bc01644 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 7 Nov 2024 11:31:47 +0000 Subject: [PATCH 10/11] enhacing --- .../full_waveform_inversion.py.rst | 36 ++++++++----------- 1 file changed, 14 insertions(+), 22 deletions(-) diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index bcd69c8b8f..29b89c6b8a 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -69,26 +69,6 @@ To achieve this, we use ensemble parallelism, which involves solving simultaneou equation (3) with different forcing terms :math:`f_s(\mathbf{x}, t)`, different :math:`J_s` and their gradients (which we will discuss later). -We start by importing firedrake and the necessary parameters used in spatial and time executions:: - - from firedrake import * - import os - test = False - if os.getenv("FIREDRAKE_CI_TESTS") == "1": - test = True - - def spatial_time_parameters(test): - if test: - # Setup for a faster test execution. - dt = 0.03 # time step in seconds - final_time = 0.6 # final time in seconds - nx, ny = 15, 15 - else: - dt = 0.002 # time step in seconds - final_time = 1.0 # final time in seconds - nx, ny = 80, 80 - return dt, final_time, nx, ny - Instantiating an ensemble requires a communicator (usually MPI_COMM_WORLD) plus the number of MPI processes to be used in each member of the ensemble (2, in this case):: @@ -119,9 +99,21 @@ The source number is defined with the ``Ensemble.ensemble_comm`` rank:: source_number = my_ensemble.ensemble_comm.rank In this example, we consider a two-dimensional square domain with a side length of 1.0 km. The mesh is -built over the ``my_ensemble.comm`` (spatial) communicator.:: +built over the ``my_ensemble.comm`` (spatial) communicator. + +.. code-block:: python + + import os + if os.getenv("FIREDRAKE_CI_TESTS") == "1": + # Setup for a faster test execution. + dt = 0.03 # time step in seconds + final_time = 0.6 # final time in seconds + nx, ny = 15, 15 + else: + dt = 0.002 # time step in seconds + final_time = 1.0 # final time in seconds + nx, ny = 80, 80 - dt, final_time, nx, ny = spatial_time_parameters(test) mesh = UnitSquareMesh(nx, ny, comm=my_ensemble.comm) The frequency of the Ricker wavelet, the source and receiver locations are defined as follows:: From 292d805bbfd05f14c307eeee514ea0c4b782b823 Mon Sep 17 00:00:00 2001 From: Iglesia Dolci Date: Thu, 7 Nov 2024 14:14:56 +0000 Subject: [PATCH 11/11] Fix convertion from rst to python file --- demos/full_waveform_inversion/full_waveform_inversion.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst index 29b89c6b8a..45d200a05d 100644 --- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst +++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst @@ -101,7 +101,7 @@ The source number is defined with the ``Ensemble.ensemble_comm`` rank:: In this example, we consider a two-dimensional square domain with a side length of 1.0 km. The mesh is built over the ``my_ensemble.comm`` (spatial) communicator. -.. code-block:: python +:: import os if os.getenv("FIREDRAKE_CI_TESTS") == "1":