Refactoring PML by ruansava · Pull Request #249 · NDF-Poli-USP/spyro

ruansava · 2026-03-24T19:36:16Z

Introduces a refactored PML formulation

Refactor PML solver and add cost save option

Introduce abc package and PML support

Refactor variable names

Refactor PML/ABC, add NRBC and damping helpers

HABC > PML

Add Neumann BC support and test updates

Add PML tests and initialize isotropic params

Add PML support and robustness fixes across modules

Introducing profilers

Allow skipping model init in IsotropicWave

Set wave.c to initial_velocity_model in tests

Set model velocity, add PML wrapper and test fixes

Fix gradient scaling and IsotropicWave inits

Add automatic mesh and update tests/tools

Olender · 2026-04-01T17:14:45Z

Should these files be here?

Yes! This would be an "output" folder

Olender · 2026-04-01T17:15:11Z

Remember to run make clean before opening a PR

Olender · 2026-04-01T17:15:56Z

I am not sure this should be a class. Can't it just be that one method that this calss contains?

Olender · 2026-04-01T17:16:47Z

+        '''
+
+        # Add specific attributes for run PML solver
+        self.c = self.initial_velocity_model


This shouldn't be done here

Olender · 2026-04-01T17:17:13Z

        - final_energy : Dissipated energy in the HABC scheme
-    f_Nyq : `float`
-        Nyquist frequency according to the time step. f_Nyq = 1 / (2 * dt)
+    f_Nyquist : `float`


Olender · 2026-04-01T17:18:36Z

        Finite element order for the Eikonal analysis
    abc_get_ref_model : `bool`
        If True, the infinite model is created
+    abc_user_pad_length : `bool`


I really like this new parameter and functionality

Olender · 2026-04-01T17:20:47Z

        if value == "PML":
-            abc_dictionary.setdefault("exponent", 2)
-            self.abc_exponent = abc_dictionary["exponent"]
-            abc_dictionary.setdefault("R", 1e-6)


There's no default cmax anymore? So in my opinion, during FWI c changes. However, we know previously what the cmax is going to be and can use this value for both optimization boundaries and PML calculation (I don't think the PML should change during FWI toghether with the velocity model)

This is the same with HABC, properties should change with the iterations. If there is no update, you can increase your reflections. You can set with the initial velocity model and maintain fixed. The idea is always to extract the maximum quantity of parameters from data automatically

Olender · 2026-04-01T17:21:36Z

-        options_dictionary = dictionary["options"]
+
+        # Check if dictionary is empty or doesn't have 'options'
+        options_dictionary = dictionary.get("options", {})


I like a simple fail if the options key doesn't exist

In fact, this should be different
def init(self, dictionary=None) because kargs should not be set with mutable objects: dictionary, list, tuple, etc
I use options_dictionary = dictionary.get("options", {}) because some text didn't need additional data

Olender · 2026-04-01T17:26:02Z

This should be called inside the automaticmeshing class. If there are non existing inputs in the meshing parameter class, please put them there

Olender · 2026-04-01T17:27:19Z

-            self.cellVertices,
-            self.cellNodeMaps,
-        ) = self.__point_locator()
+        (self.cellIDs, self.cellVertices,


PLease do not do this. It makes it harder to review and there's nothing worng with how it was

Olender · 2026-04-01T17:28:30Z

    '''Isotropic elastic wave propagator'''

-    def __init__(self, dictionary, comm=None):
+    def __init__(self, dictionary, empty_set_property=False, comm=None):


I don't think this parameter should be here. We deal with this in the model_parameters class

The problem is the field synthetic data of dicitionary of IsotropicWave. In fact, it should be made in an empty ELasticWave object; however this class is abstract. Another problem is the check that currently has this class. Without this test of material properties does not run and this should be provisionally because this class will be completely refactored. Then, when synthetic data will be deleted this field will disappear.
Another options is set with Wave class directly

Olender · 2026-04-01T17:28:45Z

        self.field_logger.add_functional("mechanical_energy",
                                         lambda: assemble(self.mechanical_energy))

+        if not empty_set_property:


This shouldn't be in this class

Please see previous comment

Olender · 2026-04-01T17:29:15Z

 def construct_solver_or_matrix_no_pml(Wave_object):
-    """Builds solver operators for wave object without a PML. Doesn't create mass matrices if
-    matrix_free option is on, which it is by default.
+    '''


Please don't do these formatting changes. Not only are they not numpy docs style, they make reviews harder

Olender · 2026-04-01T17:30:23Z

    Wave_object.u_np1 = u_np1

-    Wave_object.current_time = 0.0
+    Wave_object.current_time = 0.


Please don't do these formatting changes. They make reviews harder. If you want X.0 float to be X, we can discuss this on a wednsday and then apply changes on a single code wide PR that only changes that

Olender · 2026-04-01T17:31:13Z

-    )
+    m1 = ((Constant(1.) / (Wave_object.c * Wave_object.c))
+          * ((u - 2. * u_n + u_nm1) / Constant(dt**2))
+          * v * dx(**quad_rule))


I prefer the previous format since it is easier to follow open and closed parantheses. Also, same as the other formatting comments

Olender · 2026-04-01T17:31:25Z

    a = dot(grad(u_n), grad(v)) * dx(**quad_rule)  # explicit

-    le = 0.0
+    le = 0.


Please don't do these formatting changes. They make reviews harder. If you want X.0 float to be X, we can discuss this on a wednsday and then apply changes on a single code wide PR that only changes that

Olender · 2026-04-01T17:32:12Z

        weak_expr_abc = dot((u_n - u_nm1) / Constant(dt), v)
-
-        f_abc = (1 / Wave_object.c) * weak_expr_abc
+        f_abc = (Constant(1.) / Wave_object.c) * weak_expr_abc


The 1 doesn't change during the code so it doesn't really need to be a constant. At least thats how I understand it, but I don't see any benefit to this change

Olender · 2026-04-01T17:32:21Z

            # Damping
            le += Wave_object.eta_mask * weak_expr_abc * \
-                (1 / (Wave_object.c * Wave_object.c)) * \
+                (Constant(1.) / (Wave_object.c * Wave_object.c)) * \


Same as before

Olender · 2026-04-01T17:32:42Z

-            if Wave_object.absorb_left:
-                le += f_abc*ds(4, **qr_s)
+            # Only NRBC
+            bnds = [Wave_object.absorb_top, Wave_object.absorb_bottom,


I like the clarity of your approach here

Olender · 2026-04-01T17:32:58Z

    lin_var = fire.LinearVariationalProblem(
-        Wave_object.lhs,
-        Wave_object.rhs + Wave_object.source_function,
+        Wave_object.lhs, Wave_object.rhs + Wave_object.source_function,


Same as before

Olender · 2026-04-01T17:33:34Z

    Wave_object.solver = fire.LinearVariationalSolver(
-        lin_var, solver_parameters=solver_parameters,
-    )
+        lin_var, solver_parameters=solver_parameters)


Same comment as before

Please let's follow the code formatting approach we talked about in the meeting

I am going to stop here for today

Think of all the innocent merge conflicts you’re unleashing into the world! Think of the branches! Think of the rebases. Think of the developers who will stare into the void after a 500-line “purely cosmetic” diff!

Olender · 2026-04-01T17:47:22Z

            self.initial_velocity_model = self.guess_velocity_model
        if c is not None:
            self.initial_velocity_model.dat.data[:] = c
+        self.c = self.initial_velocity_model


This should not be necessary. Something wrong is going on somewhere else

Olender · 2026-04-01T17:49:26Z

-            spyro_plot_model(Wave_obj_real_velocity, filename=model_filename, abc_points=abc_points, high_resolution=high_resolution_model)
+            spyro_plot_model(Wave_obj_real_velocity, filename=model_filename, abc_points=abc_points)

+        Wave_obj_real_velocity.c = Wave_obj_real_velocity.initial_velocity_model


This shouldn't be here. Having this happen in a lot of places in the code is a major point of failure that will be hard to debug

ruansava added 6 commits March 13, 2026 19:01

Refactor PML solver and add cost save option

2c5de55

Refactor PML solver and add cost save option

Introduce abc package and PML support

c92e237

Introduce abc package and PML support

Refactor variable names

3e6a36d

Refactor variable names

Refactor PML/ABC, add NRBC and damping helpers

5fa3cf2

Refactor PML/ABC, add NRBC and damping helpers

HABC > PML

b1f0774

HABC > PML

Add Neumann BC support and test updates

f42ce89

Add Neumann BC support and test updates

ruansava requested review from Ig-dolci and Olender March 24, 2026 19:36

ruansava added 11 commits March 24, 2026 19:30

Add PML tests and initialize isotropic params

efc7cc5

Add PML tests and initialize isotropic params

Add PML support and robustness fixes across modules

c141bff

Add PML support and robustness fixes across modules

Update rectangle.py

0f1519c

Introducing profilers

08c9c82

Introducing profilers

Update habc.py

3163083

Update boundary_layer_io.py

c3dd41d

Allow skipping model init in IsotropicWave

7831a28

Allow skipping model init in IsotropicWave

Set wave.c to initial_velocity_model in tests

abd8bce

Set wave.c to initial_velocity_model in tests

Set model velocity, add PML wrapper and test fixes

9cad849

Set model velocity, add PML wrapper and test fixes

Fix gradient scaling and IsotropicWave inits

6de849c

Fix gradient scaling and IsotropicWave inits

Add automatic mesh and update tests/tools

2095142

Add automatic mesh and update tests/tools

Olender reviewed Apr 1, 2026

View reviewed changes

Comment thread results/acoustic_potential_energy.npy

Olender Apr 1, 2026

Copy link
Copy Markdown

Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remember to run make clean before opening a PR

Olender reviewed Apr 1, 2026

View reviewed changes

Olender requested changes Apr 1, 2026

View reviewed changes

Olender reviewed Apr 1, 2026

View reviewed changes

Uh oh!

Conversation

ruansava commented Mar 24, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ruansava Apr 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Olender Apr 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Olender Apr 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Olender Apr 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

ruansava Apr 1, 2026 •

edited

Loading

Olender Apr 1, 2026 •

edited

Loading

Olender Apr 1, 2026 •

edited

Loading

Olender Apr 1, 2026 •

edited

Loading