From 2f7ca2b11d2ae627e91246df99f73986bf6ddcc7 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 6 Aug 2025 11:59:50 -0700 Subject: [PATCH 1/2] IMprove smoothing at octree change (cherry picked from commit 1c0c6900406c8de2d179d9c63a573003bacc9750) --- simpeg_drivers/driver.py | 4 +- simpeg_drivers/utils/regularization.py | 58 ++++++++++++++++---------- 2 files changed, 38 insertions(+), 24 deletions(-) diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index fc4dc8f5..4f3db2fc 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -560,7 +560,7 @@ def get_regularization(self): f"aveCC2F{fun.orientation}", ) fun.set_weights(**{comp: average_op @ weight}) - fun.norm = average_op @ norm + fun.norm = np.round(average_op @ norm, decimals=3) functions.append(fun) if is_rotated: @@ -583,7 +583,7 @@ def get_regularization(self): f"aveCC2F{fun.orientation}", ) backward_fun.set_weights(**{comp: average_op @ weight}) - backward_fun.norm = average_op @ norm + backward_fun.norm = np.round(average_op @ norm, decimals=3) functions.append(backward_fun) # Will avoid recomputing operators if the regularization mesh is the same diff --git a/simpeg_drivers/utils/regularization.py b/simpeg_drivers/utils/regularization.py index f71528cc..a90663ba 100644 --- a/simpeg_drivers/utils/regularization.py +++ b/simpeg_drivers/utils/regularization.py @@ -339,20 +339,20 @@ def gradient_operator( :param volumes: Partial volume array. :param n_cells: Number of cells in mesh. """ - Grad = ssp.csr_matrix( + grad = ssp.csr_matrix( (volumes, (neighbors[:, 0], neighbors[:, 1])), shape=(n_cells, n_cells) ) # Normalize rows - Vol = mkvc(Grad.sum(axis=1)) - Vol[Vol > 0] = 1.0 / Vol[Vol > 0] - Grad = -sdiag(Vol) * Grad + vol = mkvc(grad.sum(axis=1)) + vol[vol > 0] = 1.0 / vol[vol > 0] + grad = -sdiag(vol) * grad diag = np.ones(n_cells) - diag[Vol == 0] = 0 - Grad = sdiag(diag) + Grad + diag[vol == 0] = 0 + grad = sdiag(diag) + grad - return Grad + return grad def rotated_gradient( @@ -399,8 +399,8 @@ def rotated_gradient( ) unit_grad = gradient_operator(neighbors, volumes, n_cells) - axes = "xyz" if dim == 3 else "xz" - return sdiag(1 / mesh.h_gridded[:, axes.find(axis)]) @ unit_grad + + return unit_grad def ensure_dip_direction_convention( @@ -464,26 +464,40 @@ def set_rotated_operators( :param forward: Whether to use forward or backward difference for derivative approximations. """ - grad_op = rotated_gradient( - function.regularization_mesh.mesh, neighbors, axis, dip, direction, forward - ) - grad_op_active = function.regularization_mesh.Pac.T @ ( - grad_op @ function.regularization_mesh.Pac - ) - active_faces = np.isclose( - grad_op_active @ np.ones(function.regularization_mesh.n_cells), 0 - ) + mesh = function.regularization_mesh + axes = "xyz" if mesh.dim == 3 else "xz" + + h_cell = mesh.mesh.h_gridded[:, axes.find(axis)] + + unit_grad_op = rotated_gradient(mesh.mesh, neighbors, axis, dip, direction, forward) + + count_op = unit_grad_op.copy() + count_op.data = np.ones_like(count_op.data) + + grad_op_active = mesh.Pac.T @ (unit_grad_op @ mesh.Pac) + count_op = mesh.Pac.T @ (count_op @ mesh.Pac) + active_faces = np.isclose(grad_op_active @ np.ones(mesh.n_cells), 0) active_faces &= grad_op_active.max(axis=1).toarray().ravel() != 0 + count_op = count_op[active_faces, :] + + h_op = sdiag( + (count_op @ (mesh.Pac.T @ h_cell) / np.asarray(count_op.sum(axis=1)).ravel()) + ** -1.0 + ) + + grad_op = h_op @ grad_op_active[active_faces, :] + avg_op = abs(grad_op_active[active_faces, :]) + avg_op = sdiag(np.asarray(avg_op.sum(axis=1)).ravel() ** -1.0) @ avg_op setattr( - function.regularization_mesh, + mesh, f"_cell_gradient_{function.orientation}", - grad_op_active[active_faces, :], + grad_op, ) setattr( - function.regularization_mesh, + mesh, f"_aveCC2F{function.orientation}", - sdiag(np.ones(function.regularization_mesh.n_cells))[active_faces, :], + avg_op, ) return function From 70142e162caac54171f235cbb9c67cb8ff2a7552 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 6 Aug 2025 20:16:57 -0700 Subject: [PATCH 2/2] MOre clean ups --- simpeg_drivers/utils/regularization.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/simpeg_drivers/utils/regularization.py b/simpeg_drivers/utils/regularization.py index a90663ba..56d6e164 100644 --- a/simpeg_drivers/utils/regularization.py +++ b/simpeg_drivers/utils/regularization.py @@ -364,7 +364,7 @@ def rotated_gradient( forward: bool = True, ) -> ssp.csr_matrix: """ - Calculated rotated gradient operator using partial volumes. + Calculated rotated gradient operator using unit partial volumes. :param mesh: Input TreeMesh. :param neighbors: Cell neighbors array. @@ -471,24 +471,21 @@ def set_rotated_operators( unit_grad_op = rotated_gradient(mesh.mesh, neighbors, axis, dip, direction, forward) - count_op = unit_grad_op.copy() - count_op.data = np.ones_like(count_op.data) + vol_avg_op = abs(unit_grad_op) + vol_avg_op.data = ( + vol_avg_op.data * mesh.mesh.cell_volumes[unit_grad_op.nonzero()[1]] + ) grad_op_active = mesh.Pac.T @ (unit_grad_op @ mesh.Pac) - count_op = mesh.Pac.T @ (count_op @ mesh.Pac) + vol_avg_op = mesh.Pac.T @ (vol_avg_op @ mesh.Pac) active_faces = np.isclose(grad_op_active @ np.ones(mesh.n_cells), 0) active_faces &= grad_op_active.max(axis=1).toarray().ravel() != 0 - count_op = count_op[active_faces, :] - - h_op = sdiag( - (count_op @ (mesh.Pac.T @ h_cell) / np.asarray(count_op.sum(axis=1)).ravel()) - ** -1.0 - ) - + vol_avg_op = vol_avg_op[active_faces, :] + vol_avg_op = sdiag(np.asarray(vol_avg_op.sum(axis=1)).ravel() ** -1) @ vol_avg_op + h_op = sdiag(vol_avg_op @ (mesh.Pac.T @ h_cell**-1.0)) grad_op = h_op @ grad_op_active[active_faces, :] - avg_op = abs(grad_op_active[active_faces, :]) - avg_op = sdiag(np.asarray(avg_op.sum(axis=1)).ravel() ** -1.0) @ avg_op + setattr( mesh, f"_cell_gradient_{function.orientation}", @@ -497,7 +494,7 @@ def set_rotated_operators( setattr( mesh, f"_aveCC2F{function.orientation}", - avg_op, + vol_avg_op, ) return function