Skip to content

Commit

Permalink
Updating get_rotation_to_tri_slip_rate_partials
Browse files Browse the repository at this point in the history
  • Loading branch information
@jploveless
jploveless committed Jul 10, 2024
1 parent fb1add9 commit b985c9e
Show file tree
Hide file tree
Showing 3 changed files with 184 additions and 267 deletions.
1 change: 1 addition & 0 deletions celeri/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@
get_transverse_projection,
get_segment_oblique_projection,
lsqlin_qp,
snap_segments,
)


Expand Down
70 changes: 32 additions & 38 deletions celeri/celeri.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,6 @@
N_MESH_DIM = 3



###############################################################
# #
# #
Expand Down Expand Up @@ -227,6 +226,7 @@ def get_command(command_file_name):

return command


def read_data(command: Dict):
logger.info("Reading data files")
# Read segment data
Expand Down Expand Up @@ -471,10 +471,9 @@ def read_data(command: Dict):
return segment, block, meshes, station, mogi, sar



####################################################################################################
# #
# #
# #
# 88888888ba 88888888ba ,ad8888ba, ,ad8888ba, 88888888888 ad88888ba ad88888ba #
# 88 "8b 88 "8b d8"' `"8b d8"' `"8b 88 d8" "8b d8" "8b #
# 88 ,8P 88 ,8P d8' `8b d8' 88 Y8, Y8, #
Expand All @@ -483,7 +482,7 @@ def read_data(command: Dict):
# 88 88 `8b Y8, ,8P Y8, 88 `8b `8b #
# 88 88 `8b Y8a. .a8P Y8a. .a8P 88 Y8a a8P Y8a a8P #
# 88 88 `8b `"Y8888Y"' `"Y8888Y"' 88888888888 "Y88888P" "Y88888P" #
# #
# #
# #
####################################################################################################
def triangle_area(triangles):
Expand Down Expand Up @@ -940,7 +939,6 @@ def station_row_keep(assembly):
return assembly



def get_block_centroid(segment, block_idx):
"""
Calculate centroid of a block based on boundary polygon
Expand Down Expand Up @@ -1354,6 +1352,7 @@ def get_index(assembly, station, block, meshes):
index.n_operator_cols = 3 * index.n_blocks + 2 * index.n_tde_total
return index


def get_index_no_meshes(assembly, station, block):
# NOTE: Merge with above if possible.
# Make sure empty meshes work
Expand Down Expand Up @@ -1391,7 +1390,6 @@ def get_index_no_meshes(assembly, station, block):
return index



def get_processed_data_structures(command):
# NOTE: DEPRECATE (NEVER CALLED)
data = addict.Dict()
Expand Down Expand Up @@ -2613,9 +2611,11 @@ def get_rotation_to_tri_slip_rate_partials(meshes, mesh_idx, segment, block):
for a single mesh. Called within a loop from get_tde_coupling_constraints()
"""
n_blocks = len(block)
# Calculate right-hand rule strike for each segment
rhrstrike = segment.azimuth + 180 * (segment.dip > 90)
tri_slip_rate_partials = np.zeros((3 * meshes.n_tde, 3 * n_blocks))

# Generate strikes for elements using same sign convention as segments
tristrike = meshes.strike
tristrike[meshes.strike > 180] -= 180
# Find subset of segments that are replaced by this mesh
seg_replace_idx = np.where(
(segment.patch_flag != 0) & (segment.patch_file_name == mesh_idx)
Expand Down Expand Up @@ -2683,21 +2683,27 @@ def get_rotation_to_tri_slip_rate_partials(meshes, mesh_idx, segment, block):
meshes.lat_centroid[el_idx],
)
# Project about fault strike
unit_x_parallel = np.sin(np.deg2rad(meshes.strike[el_idx]))
unit_y_parallel = np.cos(np.deg2rad(meshes.strike[el_idx]))
unit_x_perpendicular = np.cos(np.deg2rad(meshes.strike[el_idx]))
unit_y_perpendicular = -np.sin(np.deg2rad(meshes.strike[el_idx]))
unit_x_parallel = np.cos(np.deg2rad(90 - tristrike[el_idx]))
unit_y_parallel = np.sin(np.deg2rad(90 - tristrike[el_idx]))
unit_x_perpendicular = np.sin(np.deg2rad(tristrike[el_idx] - 90))
unit_y_perpendicular = np.cos(np.deg2rad(tristrike[el_idx] - 90))
# Project by fault dip
scale_factor = -1.0 / abs(np.cos(np.deg2rad(meshes.dip[el_idx])))
scale_factor = 1.0 / abs(np.cos(np.deg2rad(meshes.dip[el_idx])))
slip_rate_matrix = np.array(
[
[
unit_x_parallel * vel_east_to_omega_x
+ unit_y_parallel * vel_north_to_omega_x,
unit_x_parallel * vel_east_to_omega_y
+ unit_y_parallel * vel_north_to_omega_y,
unit_x_parallel * vel_east_to_omega_z
+ unit_y_parallel * vel_north_to_omega_z,
(
unit_x_parallel * vel_east_to_omega_x
+ unit_y_parallel * vel_north_to_omega_x
),
(
unit_x_parallel * vel_east_to_omega_y
+ unit_y_parallel * vel_north_to_omega_y
),
(
unit_x_parallel * vel_east_to_omega_z
+ unit_y_parallel * vel_north_to_omega_z
),
],
[
scale_factor
Expand All @@ -2719,15 +2725,11 @@ def get_rotation_to_tri_slip_rate_partials(meshes, mesh_idx, segment, block):
[0, 0, 0],
]
)
sign_corr = 1
# if (rhrstrike[meshes.closest_segment_idx[el_idx]] > 90) & (
# rhrstrike[meshes.closest_segment_idx[el_idx]] < 270
# ):
# sign_corr = -1
# if (rhrstrike[meshes.closest_segment_idx[el_idx]] > 0) & (
# rhrstrike[meshes.closest_segment_idx[el_idx]] < 180
# ):
# sign_corr = 1
# This correction gives -1 for strikes > 90
# Equivalent to the if statement in get_rotation_to_slip_rate_partials
sign_corr = -np.sign(meshes.strike[el_idx] - 90)

# Insert this element's partials into operator
tri_slip_rate_partials[
row_idx : row_idx + 3, column_idx_east : column_idx_east + 3
] = (sign_corr * slip_rate_matrix)
Expand All @@ -2737,8 +2739,6 @@ def get_rotation_to_tri_slip_rate_partials(meshes, mesh_idx, segment, block):
return tri_slip_rate_partials




def get_block_motion_constraints(assembly: Dict, block: pd.DataFrame, command: Dict):
"""
Applying a priori block motion constraints
Expand Down Expand Up @@ -3535,7 +3535,6 @@ def get_h_matrices_for_tde_meshes(
return H, col_norms



######################################################################
# #
# #
Expand All @@ -3547,7 +3546,7 @@ def get_h_matrices_for_tde_meshes(
# `8b Y8, ,8P 88 `8b d8' 88 #
# Y8a a8P Y8a. .a8P 88 `888' 88 #
# "Y88888P" `"Y8888Y"' 88888888888 `8' 88888888888 #
# #
# #
# #
######################################################################

Expand Down Expand Up @@ -3904,7 +3903,6 @@ def cvxopt_to_numpy_matrix(A):
return sol



def matvec_wrapper(h_matrix_solve_parameters):
def matvec_caller(x):
return matvec(x, h_matrix_solve_parameters)
Expand Down Expand Up @@ -4612,7 +4610,6 @@ def build_and_solve_dense_no_meshes(command, assembly, operators, data):
return estimation, operators, index



#########################################################################################
# #
# #
Expand Down Expand Up @@ -4765,7 +4762,6 @@ def write_output_supplemental(
pickle.dump([command, index, data, operators, estimation, assembly], f)



########################################################
# #
# #
Expand Down Expand Up @@ -5393,8 +5389,6 @@ def common_plot_elements(segment: pd.DataFrame, lon_range: Tuple, lat_range: Tup
)




################################################################################################
# #
# #
Expand Down Expand Up @@ -5425,6 +5419,7 @@ def get_logger(command):
logger.info(f"Write log file: {command.output_path}/{command.run_name}.log")
return logger


def wrap2360(lon):
lon[np.where(lon < 0.0)] += 360.0
return lon
Expand Down Expand Up @@ -5802,4 +5797,3 @@ def align_velocities(df_1, df_2, distance_threshold):
df_1_aligned["north_vel"] = df_1_aligned["north_vel"] + rotated_vels_match[1::2]

return df_1_aligned

Loading

0 comments on commit b985c9e

Please sign in to comment.