Merge pull request #88 from NDF-Poli-USP/main

Updating with main
NDF-Poli-USP · Jul 3, 2023 · 4347d91 · 4347d91
2 parents 0494c2d + fef8f54
commit 4347d91
Show file tree

Hide file tree

Showing 11 changed files with 67 additions and 67 deletions.
diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml
@@ -17,38 +17,38 @@ jobs:
     - uses: actions/checkout@v3
     - name: Running serial tests in parallel (only one test per core)
       run: |
-          source /opt/firedrake/bin/activate
+          source /home/olender/Firedrakes/main/firedrake/bin/activate
           pytest -n 10 --cov-report=xml --cov=spyro test/
     - name: Running serial tests for adjoint
       run: |
-          source /opt/firedrake/bin/activate
+          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 /opt/firedrake/bin/activate
+          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 /opt/firedrake/bin/activate
+          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: Running parallel 3D forward test
       run: |
-          source /opt/firedrake/bin/activate
+          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 15 pytest test_3d/test_forward_3d.py 
+          mpiexec -n 10 pytest test_3d/test_forward_3d.py 
     - name: Covering parallel 3D forward test
       continue-on-error: true
       run: |
-          source /opt/firedrake/bin/activate
+          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 15 pytest --cov-report=xml --cov-append --cov=spyro test_3d/test_forward_3d.py 
+          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="6cd21147-54f7-4b77-94ad-4b138053401d" && bash <(curl -s https://codecov.io/bash)
 
diff --git a/demos/run_fwi.py b/demos/run_fwi.py
@@ -83,7 +83,7 @@
 quad_rule = finat.quadrature.make_quadrature(
     V.finat_element.cell, V.ufl_element().degree(), "KMV"
 )
-dxlump = dx(rule=quad_rule)
+dxlump = dx(scheme=quad_rule)
 
 water = np.where(vp.dat.data[:] < 1.51)
 
@@ -110,8 +110,8 @@ def regularize_gradient(vp, dJ, gamma):
     """Tikhonov regularization"""
     m_u = TrialFunction(V)
     m_v = TestFunction(V)
-    mgrad = m_u * m_v * dx(rule=qr_x)
-    ffG = dot(grad(vp), grad(m_v)) * dx(rule=qr_x)
+    mgrad = m_u * m_v * dx(scheme=qr_x)
+    ffG = dot(grad(vp), grad(m_v)) * dx(scheme=qr_x)
     G = mgrad - ffG
     lhsG, rhsG = lhs(G), rhs(G)
     gradreg = Function(V)

diff --git a/paper/run_fwi_2d.py b/paper/run_fwi_2d.py
@@ -93,7 +93,7 @@ def get_memory_usage():
 quad_rule = finat.quadrature.make_quadrature(
     V.finat_element.cell, V.ufl_element().degree(), "KMV"
 )
-dxlump = dx(rule=quad_rule)
+dxlump = dx(scheme=quad_rule)
 
 water = np.where(vp.dat.data[:] < 1.51)
 
@@ -117,8 +117,8 @@ def regularize_gradient(vp, dJ):
     """Tikhonov regularization"""
     m_u = TrialFunction(V)
     m_v = TestFunction(V)
-    mgrad = m_u * m_v * dx(rule=quad_rule)
-    ffG = dot(grad(vp), grad(m_v)) * dx(rule=quad_rule)
+    mgrad = m_u * m_v * dx(scheme=quad_rule)
+    ffG = dot(grad(vp), grad(m_v)) * dx(scheme=quad_rule)
     G = mgrad - ffG
     lhsG, rhsG = lhs(G), rhs(G)
     gradreg = Function(V)

diff --git a/paper/run_fwi_3d.py b/paper/run_fwi_3d.py
@@ -115,7 +115,7 @@ def get_memory_usage():
 quad_rule = finat.quadrature.make_quadrature(
     V.finat_element.cell, V.ufl_element().degree(), "KMV"
 )
-dxlump = dx(rule=quad_rule)
+dxlump = dx(scheme=quad_rule)
 
 water = np.where(vp.dat.data[:] < 1.51)
 
@@ -139,8 +139,8 @@ def regularize_gradient(vp, dJ):
     """Tikhonov regularization"""
     m_u = TrialFunction(V)
     m_v = TestFunction(V)
-    mgrad = m_u * m_v * dx(rule=quad_rule)
-    ffG = dot(grad(vp), grad(m_v)) * dx(rule=quad_rule)
+    mgrad = m_u * m_v * dx(scheme=quad_rule)
+    ffG = dot(grad(vp), grad(m_v)) * dx(scheme=quad_rule)
     G = mgrad - ffG
     lhsG, rhsG = lhs(G), rhs(G)
     gradreg = Function(V)

diff --git a/spyro/solvers/forward.py b/spyro/solvers/forward.py
@@ -140,7 +140,7 @@ def forward(
                 (sigma_x + sigma_z)
                 * ((u - u_nm1) / Constant(2.0 * dt))
                 * v
-                * dx(rule=qr_x)
+                * dx(scheme=qr_x)
             )
         elif dim == 3:
 
@@ -178,12 +178,12 @@ def forward(
     t = 0.0
 
     # -------------------------------------------------------
-    m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x)
-    a = c * c * dot(grad(u_n), grad(v)) * dx(rule=qr_x)  # explicit
+    m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(scheme=qr_x)
+    a = c * c * dot(grad(u_n), grad(v)) * dx(scheme=qr_x)  # explicit
 
     nf = 0
     if model["BCs"]["outer_bc"] == "non-reflective":
-        nf = c * ((u_n - u_nm1) / dt) * v * ds(rule=qr_s)
+        nf = c * ((u_n - u_nm1) / dt) * v * ds(scheme=qr_s)
 
     FF = m1 + a + nf
 
@@ -192,41 +192,41 @@ def forward(
         B = Function(W)
 
         if dim == 2:
-            pml2 = sigma_x * sigma_z * u_n * v * dx(rule=qr_x)
-            pml3 = inner(pp_n, grad(v)) * dx(rule=qr_x)
+            pml2 = sigma_x * sigma_z * u_n * v * dx(scheme=qr_x)
+            pml3 = inner(pp_n, grad(v)) * dx(scheme=qr_x)
             FF += pml1 + pml2 + pml3
             # -------------------------------------------------------
-            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(rule=qr_x)
-            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(rule=qr_x)
-            dd = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dx(rule=qr_x)
+            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(scheme=qr_x)
+            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(scheme=qr_x)
+            dd = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dx(scheme=qr_x)
             FF += mm1 + mm2 + dd
         elif dim == 3:
             pml1 = (
                 (sigma_x + sigma_y + sigma_z)
                 * ((u - u_n) / Constant(dt))
                 * v
-                * dx(rule=qr_x)
+                * dx(scheme=qr_x)
             )
             pml2 = (
                 (sigma_x * sigma_y + sigma_x * sigma_z + sigma_y * sigma_z)
                 * u_n
                 * v
-                * dx(rule=qr_x)
+                * dx(scheme=qr_x)
             )
-            pml3 = (sigma_x * sigma_y * sigma_z) * psi_n * v * dx(rule=qr_x)
-            pml4 = inner(pp_n, grad(v)) * dx(rule=qr_x)
+            pml3 = (sigma_x * sigma_y * sigma_z) * psi_n * v * dx(scheme=qr_x)
+            pml4 = inner(pp_n, grad(v)) * dx(scheme=qr_x)
 
             FF += pml1 + pml2 + pml3 + pml4
             # -------------------------------------------------------
-            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(rule=qr_x)
-            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(rule=qr_x)
-            dd1 = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dx(rule=qr_x)
-            dd2 = -c * c * inner(grad(psi_n), dot(Gamma_3, qq)) * dx(rule=qr_x)
+            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(scheme=qr_x)
+            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(scheme=qr_x)
+            dd1 = c * c * inner(grad(u_n), dot(Gamma_2, qq)) * dx(scheme=qr_x)
+            dd2 = -c * c * inner(grad(psi_n), dot(Gamma_3, qq)) * dx(scheme=qr_x)
 
             FF += mm1 + mm2 + dd1 + dd2
             # -------------------------------------------------------
-            mmm1 = (dot((psi - psi_n), phi) / Constant(dt)) * dx(rule=qr_x)
-            uuu1 = (-u_n * phi) * dx(rule=qr_x)
+            mmm1 = (dot((psi - psi_n), phi) / Constant(dt)) * dx(scheme=qr_x)
+            uuu1 = (-u_n * phi) * dx(scheme=qr_x)
 
             FF += mmm1 + uuu1
     else:

diff --git a/spyro/solvers/forward_AD.py b/spyro/solvers/forward_AD.py
@@ -89,16 +89,16 @@ def forward(
 
     t = 0.0
     m = 1 / (c * c)
-    m1 = m * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x)
-    a = dot(grad(u_n), grad(v)) * dx(rule=qr_x)  # explicit
+    m1 = m * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(scheme=qr_x)
+    a = dot(grad(u_n), grad(v)) * dx(scheme=qr_x)  # explicit
     f = Function(V)
     nf = 0
 
     if model["BCs"]["outer_bc"] == "non-reflective":
-        nf = c * ((u_n - u_nm1) / dt) * v * ds(rule=qr_s)
+        nf = c * ((u_n - u_nm1) / dt) * v * ds(scheme=qr_s)
 
     h = CellSize(mesh)
-    FF = m1 + a + nf - (1 / (h / degree * h / degree)) * f * v * dx(rule=qr_x)
+    FF = m1 + a + nf - (1 / (h / degree * h / degree)) * f * v * dx(scheme=qr_x)
     X = Function(V)
 
     lhs_ = lhs(FF)

diff --git a/spyro/solvers/gradient.py b/spyro/solvers/gradient.py
@@ -161,12 +161,12 @@ def gradient(
     t = 0.0
 
     # -------------------------------------------------------
-    m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(rule=qr_x)
-    a = c * c * dot(grad(u_n), grad(v)) * dx(rule=qr_x)  # explicit
+    m1 = ((u - 2.0 * u_n + u_nm1) / Constant(dt**2)) * v * dx(scheme=qr_x)
+    a = c * c * dot(grad(u_n), grad(v)) * dx(scheme=qr_x)  # explicit
 
     nf = 0
     if model["BCs"]["outer_bc"] == "non-reflective":
-        nf = c * ((u_n - u_nm1) / dt) * v * ds(rule=qr_s)
+        nf = c * ((u_n - u_nm1) / dt) * v * ds(scheme=qr_s)
 
     FF = m1 + a + nf
 
@@ -175,39 +175,39 @@ def gradient(
         B = Function(W)
 
         if dim == 2:
-            pml1 = (sigma_x + sigma_z) * ((u - u_n) / dt) * v * dx(rule=qr_x)
-            pml2 = sigma_x * sigma_z * u_n * v * dx(rule=qr_x)
-            pml3 = c * c * inner(grad(v), dot(Gamma_2, pp_n)) * dx(rule=qr_x)
+            pml1 = (sigma_x + sigma_z) * ((u - u_n) / dt) * v * dx(scheme=qr_x)
+            pml2 = sigma_x * sigma_z * u_n * v * dx(scheme=qr_x)
+            pml3 = c * c * inner(grad(v), dot(Gamma_2, pp_n)) * dx(scheme=qr_x)
 
             FF += pml1 + pml2 + pml3
             # -------------------------------------------------------
-            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(rule=qr_x)
-            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(rule=qr_x)
-            dd = inner(qq, grad(u_n)) * dx(rule=qr_x)
+            mm1 = (dot((pp - pp_n), qq) / Constant(dt)) * dx(scheme=qr_x)
+            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(scheme=qr_x)
+            dd = inner(qq, grad(u_n)) * dx(scheme=qr_x)
 
             FF += mm1 + mm2 + dd
         elif dim == 3:
-            pml1 = (sigma_x + sigma_y + sigma_z) * ((u - u_n) / dt) * v * dx(rule=qr_x)
-            uuu1 = (-v * psi_n) * dx(rule=qr_x)
+            pml1 = (sigma_x + sigma_y + sigma_z) * ((u - u_n) / dt) * v * dx(scheme=qr_x)
+            uuu1 = (-v * psi_n) * dx(scheme=qr_x)
             pml2 = (
                 (sigma_x * sigma_y + sigma_x * sigma_z + sigma_y * sigma_z)
                 * u_n
                 * v
-                * dx(rule=qr_x)
+                * dx(scheme=qr_x)
             )
-            dd1 = c * c * inner(grad(v), dot(Gamma_2, pp_n)) * dx(rule=qr_x)
+            dd1 = c * c * inner(grad(v), dot(Gamma_2, pp_n)) * dx(scheme=qr_x)
 
             FF += pml1 + pml2 + dd1 + uuu1
             # -------------------------------------------------------
-            mm1 = (dot((pp - pp_n), qq) / dt) * dx(rule=qr_x)
-            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(rule=qr_x)
-            pml4 = inner(qq, grad(u_n)) * dx(rule=qr_x)
+            mm1 = (dot((pp - pp_n), qq) / dt) * dx(scheme=qr_x)
+            mm2 = inner(dot(Gamma_1, pp_n), qq) * dx(scheme=qr_x)
+            pml4 = inner(qq, grad(u_n)) * dx(scheme=qr_x)
 
             FF += mm1 + mm2 + pml4
             # -------------------------------------------------------
-            pml3 = (sigma_x * sigma_y * sigma_z) * phi * u_n * dx(rule=qr_x)
-            mmm1 = (dot((psi - psi_n), phi) / dt) * dx(rule=qr_x)
-            mmm2 = -c * c * inner(grad(phi), dot(Gamma_3, pp_n)) * dx(rule=qr_x)
+            pml3 = (sigma_x * sigma_y * sigma_z) * phi * u_n * dx(scheme=qr_x)
+            mmm1 = (dot((psi - psi_n), phi) / dt) * dx(scheme=qr_x)
+            mmm2 = -c * c * inner(grad(phi), dot(Gamma_3, pp_n)) * dx(scheme=qr_x)
 
             FF += mmm1 + mmm2 + pml3
     else:
@@ -223,12 +223,12 @@ def gradient(
     # Define gradient problem
     m_u = TrialFunction(V)
     m_v = TestFunction(V)
-    mgrad = m_u * m_v * dx(rule=qr_x)
+    mgrad = m_u * m_v * dx(scheme=qr_x)
 
     uuadj = Function(V)  # auxiliarly function for the gradient compt.
     uufor = Function(V)  # auxiliarly function for the gradient compt.
 
-    ffG = 2.0 * c * dot(grad(uuadj), grad(uufor)) * m_v * dx(rule=qr_x)
+    ffG = 2.0 * c * dot(grad(uuadj), grad(uufor)) * m_v * dx(scheme=qr_x)
 
     G = mgrad - ffG
     lhsG, rhsG = lhs(G), rhs(G)

diff --git a/spyro/utils/estimate_timestep.py b/spyro/utils/estimate_timestep.py
@@ -20,7 +20,7 @@ def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True):
     quad_rule = finat.quadrature.make_quadrature(
         V.finat_element.cell, V.ufl_element().degree(), "KMV"
     )
-    dxlump = fd.dx(rule=quad_rule)
+    dxlump = fd.dx(scheme=quad_rule)
     A = fd.assemble(u * v * dxlump)
     ai, aj, av = A.petscmat.getValuesCSR()
     av_inv = []
@@ -52,7 +52,7 @@ def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True):
 
     # print(max_eigval)
     if np.sqrt(max_eigval) > 0.0:
-        max_dt = np.float(2 / np.sqrt(max_eigval))
+        max_dt = float(2 / np.sqrt(max_eigval))
     else:
         max_dt = 100000000
     # print(

diff --git a/test/not_a_test.py b/test/not_a_test.py
@@ -56,7 +56,7 @@ def test_gradient_talyor_remainder_v2():
     class L2Inner(object):
         def __init__(self):
             self.A = assemble(
-                TrialFunction(V) * TestFunction(V) * dx(rule=qr_x), mat_type="matfree"
+                TrialFunction(V) * TestFunction(V) * dx(scheme=qr_x), mat_type="matfree"
             )
             self.Ap = as_backend_type(self.A).mat()
 

diff --git a/test/test_gradient.py b/test/test_gradient.py
@@ -117,7 +117,7 @@ def _test_gradient(options, pml=False):
 
     # compute the gradient of the control (to be verified)
     dJ = gradient(options, mesh, comm, vp_guess, receivers, p_guess, misfit)
-    dJ *= mask
+    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
@@ -145,7 +145,7 @@ def _test_gradient(options, pml=False):
         )
 
         Jp = functional(options, p_exact_recv - p_guess_recv)
-        projnorm = assemble(mask * dJ * delta_m * dx(rule=qr_x))
+        projnorm = assemble(mask * dJ * delta_m * dx(scheme=qr_x))
         fd_grad = (Jp - Jm) / step
         print(
             "\n Cost functional for step "

diff --git a/test/test_gradient_3d.py b/test/test_gradient_3d.py
@@ -136,7 +136,7 @@ def _test_gradient(options, pml=False):
         )
 
         Jp = functional(options, p_exact_recv - p_guess_recv)
-        projnorm = assemble(box1 * dJ * delta_m * dx(rule=qr_x))
+        projnorm = assemble(box1 * dJ * delta_m * dx(scheme=qr_x))
         fd_grad = (Jp - Jm) / step
         print(
             "\n Cost functional for step "