From 7451b14db8eb5e18b1b576392b54715950e74e3e Mon Sep 17 00:00:00 2001 From: Justin Bonus Date: Sun, 29 Sep 2024 07:33:19 -0700 Subject: [PATCH] Lint / format shakermaker files for the sake of the green check-mark --- .../stochasticWave/StochasticWave.py | 2 +- ...Loads.py => StochasticWaveLoadsJONSWAP.py} | 692 +++++++++--------- .../tools/ShakerMaker/ShakerMakersubmitjob.py | 215 +++--- .../TapisFiles/ShakerMakermodel.py | 271 +++---- modules/tools/ShakerMaker/faultplotter.py | 440 +++++------ 5 files changed, 839 insertions(+), 781 deletions(-) mode change 100755 => 100644 modules/createEVENT/stochasticWave/StochasticWave.py rename modules/createEVENT/stochasticWave/{Ex4_WaveLoads.py => StochasticWaveLoadsJONSWAP.py} (96%) mode change 100755 => 100644 diff --git a/modules/createEVENT/stochasticWave/StochasticWave.py b/modules/createEVENT/stochasticWave/StochasticWave.py old mode 100755 new mode 100644 index 2d85a23cf..24f94e7ad --- a/modules/createEVENT/stochasticWave/StochasticWave.py +++ b/modules/createEVENT/stochasticWave/StochasticWave.py @@ -7,7 +7,7 @@ import re import sys -import Ex4_WaveLoads +import StochasticWaveLoadsJONSWAP """ Portions of this module are implemented courtesy of the welib python package: diff --git a/modules/createEVENT/stochasticWave/Ex4_WaveLoads.py b/modules/createEVENT/stochasticWave/StochasticWaveLoadsJONSWAP.py old mode 100755 new mode 100644 similarity index 96% rename from modules/createEVENT/stochasticWave/Ex4_WaveLoads.py rename to modules/createEVENT/stochasticWave/StochasticWaveLoadsJONSWAP.py index daaccaa05..cd91956f9 --- a/modules/createEVENT/stochasticWave/Ex4_WaveLoads.py +++ b/modules/createEVENT/stochasticWave/StochasticWaveLoadsJONSWAP.py @@ -1,346 +1,346 @@ -#!/usr/bin/env python3 - -"""Compute inline/total hydrodynamic force and moments on a monopile using Morison's equation""" # noqa: D400 - -import argparse -from fractions import Fraction - -import matplotlib.pyplot as plt -import numpy as np -import pandas as pd - -# Local -from welib.tools.figure import defaultRC - -defaultRC() - -from welib.hydro.morison import inline_load # noqa: E402 -from welib.hydro.wavekin import elevation2d, kinematics2d, wavenumber # noqa: E402 -from welib.tools.colors import python_colors # noqa: E402 - -# --- Parameters -g = 9.80665 # gravity [m/s^2] -h = 30.0 # water depth [m] -rho = 1000 # water density -D = 6 # monopile diameter [m] -CD = 1 # drag coefficient -CM = 2 # added-mass coefficient -a = 3.0 # wave peak amplitude [m] -T = 12.0 # period [s] -eps = 0 # phase shift [rad] -f = 1.0 / T -k = wavenumber(f, h, g) - -nz = 30 # number of points used in the z direction to compute loads - -z_ref = -h # reference point for moment calculation - -# --------------------------------------------------------------------------------} -# --- Inline force and moments as function of time, with or without Wheeler stretching -# --------------------------------------------------------------------------------{ -time = np.linspace(0, T, 9) - -fig1, axes1 = plt.subplots( - 2, 4, sharex=True, sharey=True, figsize=(12.8, 4.8) -) # (6.4,4.8) -fig1.subplots_adjust( - left=0.05, right=0.99, top=0.95, bottom=0.09, hspace=0.26, wspace=0.11 -) - -fig2, axes2 = plt.subplots( - 2, 4, sharex=True, sharey=True, figsize=(12.8, 4.8) -) # (6.4,4.8) -fig2.subplots_adjust( - left=0.05, right=0.99, top=0.95, bottom=0.09, hspace=0.26, wspace=0.11 -) - -XLIM = [-75, 75] # For inline force -XLIMM = [-2500, 2500] # For inline moment - -for it, t in enumerate(time[:-1]): - # Wave kinematics - eta = elevation2d(a, f, k, eps, t, x=0) - z = np.linspace(-h, eta, nz) - u, du = kinematics2d(a, f, k, eps, h, t, z, Wheeler=True, eta=eta) - u0, du0 = kinematics2d(a, f, k, eps, h, t, z) - # Wave loads with wheeler - p_tot = inline_load(u, du, D, CD, CM, rho) - p_inertia = inline_load(u, du, D, CD * 0, CM, rho) - p_drag = inline_load(u, du, D, CD, CM * 0, rho) - dM = p_tot * (z - z_ref) # [Nm/m] # noqa: N816 - - # Wave loads without Wheeler - p_tot0 = inline_load(u0, du0, D, CD, CM, rho) - p_inertia0 = inline_load(u0, du0, D, CD * 0, CM, rho) - p_drag0 = inline_load(u0, du0, D, CD, CM * 0, rho) - dM0 = p_tot0 * (z - z_ref) # [Nm/m] # noqa: N816 - - # Plot inline force - ax = axes1[int(it / 4), np.mod(it, 4)] - ax.plot(p_inertia / 1000, z, '-', c=python_colors(0), label=r'$f_{inertia}$') - ax.plot(p_drag / 1000, z, '-', c=python_colors(3), label=r'$f_{drag}$') - ax.plot(p_tot / 1000, z, 'k-', label=r'$f_{tot}$') - ax.plot(p_inertia0 / 1000, z, '+', c=python_colors(0)) - ax.plot(p_drag0 / 1000, z, '+', c=python_colors(3)) - ax.plot(p_tot0 / 1000, z, 'k+') - ax.set_title(f't/T={Fraction(t / T)}') - if it == 0: - ax.legend() - ax.plot(XLIM, [0, 0], 'k') - ax.plot(XLIM, [a, a], 'k--') - ax.plot(XLIM, [-a, -a], 'k--') - - # Plot inline moment - ax = axes2[int(it / 4), np.mod(it, 4)] - ax.plot(dM / 1000, z, 'k-', label=r'$dM_{tot}$ with Wheeler') - ax.plot(dM0 / 1000, z, 'k+', label=r'$dM_{tot}$ no-correction') - ax.set_title(f't/T={Fraction(t / T)}') - if it == 0: - ax.legend() - ax.plot(XLIMM, [0, 0], 'k') - ax.plot(XLIMM, [a, a], 'k--') - ax.plot(XLIMM, [-a, -a], 'k--') - - -axes1[0, 0].set_xlim(XLIM) -axes1[0, 0].set_ylim([-h, a + 1]) -axes1[0, 0].set_ylabel('Depth z [m]') -axes1[1, 0].set_ylabel('Depth z [m]') -for i in range(4): - axes1[1, i].set_xlabel('Inline force [kN/m]') - - -fig1.savefig('forces.png') -# fig1.savefig('forces.webp') -# fig1.show() - -axes2[0, 0].set_xlim(XLIMM) -axes2[0, 0].set_ylim([-h, a + 1]) -axes2[0, 0].set_ylabel('Depth z [m]') -axes2[1, 0].set_ylabel('Depth z [m]') -for i in range(4): - axes2[1, i].set_xlabel('Inline moment [kNm/m]') - -fig2.savefig('moments.png') -# fig2.savefig('moments.webp') -# fig2.show() -# --------------------------------------------------------------------------------} -# --- Integrated force and sea bed moment over a period -# --------------------------------------------------------------------------------{ -time = np.linspace(0, 60.0, 6001) - -veta = np.zeros(time.shape) -vF = np.zeros(time.shape) # noqa: N816 -vM = np.zeros(time.shape) # noqa: N816 -vF0 = np.zeros(time.shape) # noqa: N816 -vM0 = np.zeros(time.shape) # noqa: N816 - -XLIM = [-75, 75] # For inline force -XLIMM = [-2500, 2500] # For inline moment - -elevation = np.zeros((len(time), nz)) -velocity = np.zeros((len(time), nz)) -accel = np.zeros((len(time), nz)) -force = np.zeros((len(time), nz)) - -for it, t in enumerate(time): - # Wave kinematics - veta[it] = elevation2d(a, f, k, eps, t, x=0) - z = np.linspace(-h, veta[it], nz) - u, du = kinematics2d(a, f, k, eps, h, t, z, Wheeler=True, eta=veta[it]) - u0, du0 = kinematics2d(a, f, k, eps, h, t, z) - # Wave loads with Wheeler - p_tot = inline_load(u, du, D, CD, CM, rho) - vF[it] = np.trapz(p_tot, z) # [N] # noqa: NPY201 - vM[it] = np.trapz(p_tot * (z - z_ref), z) # [Nm] # noqa: NPY201 - # Wave loads without Wheeler - p_tot0 = inline_load(u0, du0, D, CD, CM, rho) - vF0[it] = np.trapz(p_tot0, z) # [N] # noqa: NPY201 - vM0[it] = np.trapz(p_tot0 * (z - z_ref), z) # [Nm] # noqa: NPY201 - - elevation[it, :] = z.copy() - velocity[it, :] = u.copy() - accel[it, :] = du.copy() - force[it, :] = p_tot.copy() - -# Plot -fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6.4, 4.8)) # (6.4,4.8) -fig.subplots_adjust( - left=0.12, right=0.95, top=0.95, bottom=0.11, hspace=0.14, wspace=0.20 -) -# axes[0] = axes[0] -axes[0].plot(time / T, veta, 'k-') -axes[0].set_ylabel('Elevation [m]') -axes[0].grid(True) # noqa: FBT003 -# axes[1] = axes[1] -axes[1].plot(time / T, vF0 / 1e6, label='Standard') -axes[1].plot(time / T, vF / 1e6, 'k-', label='Wheeler Correction') -axes[1].set_ylabel('Streamwise Load, Cumulative [MN]') -axes[1].legend() -axes[1].grid(True) # noqa: FBT003 -# axes[2] = axes[2] -axes[2].plot(time / T, vM0 / 1e6, label='Standard') -axes[2].plot(time / T, vM / 1e6, 'k-', label='Wheeler Correction') -axes[2].set_ylabel('Sea-Bed Moment [MNm]') -axes[2].set_xlabel('Dimensionless Time, t/T [-]') -axes[2].legend() -axes[2].grid(True) # noqa: FBT003 - -# fig.savefig('IntegratedPileLoads.png') -# fig.savefig('IntegratedPileLoads.webp') -fig.savefig('IntegratedPileLoads.png') -# fig.show() - -# now save csv of the velocity, acceleration, force, and moment -veta_df = pd.DataFrame() -u_df = pd.DataFrame() -du_df = pd.DataFrame() -for i in range(nz): - dof = 1 - name = 'Disp_' + str(i + 1) + '_' + str(dof) - veta_df[name] = elevation[:, i] - name = 'Vel_' + str(i + 1) + '_' + str(dof) - u_df[name] = velocity[:, i] - name = 'RMSA_' + str(i + 1) + '_' + str(dof) - du_df[name] = accel[:, i] - -# add column per each force recorder -result_df = pd.DataFrame() -for i in range(nz): - dof = 1 - # name = 'Node_' + str(i+1) + '_' + str(dof) - name = 'Force_' + str(i + 1) + '_' + str(dof) - result_df[name] = force[:, i] - -# write columns to columns in csv files -(veta_df.T).to_csv('disp.evt', sep=' ', encoding='utf-8', index=False, header=False) -(u_df.T).to_csv('vel.evt', sep=' ', encoding='utf-8', index=False, header=False) -(du_df.T).to_csv('accel.evt', sep=' ', encoding='utf-8', index=False, header=False) -(result_df.T).to_csv( - 'forces.evt', sep=' ', encoding='utf-8', index=False, header=False -) - -# write columns to columns in csv files -(veta_df.T).to_csv('disp.out', sep=' ', encoding='utf-8', index=False, header=False) -(u_df.T).to_csv('vel.out', sep=' ', encoding='utf-8', index=False, header=False) -(du_df.T).to_csv('accel.out', sep=' ', encoding='utf-8', index=False, header=False) -(result_df.T).to_csv( - 'forces.out', sep=' ', encoding='utf-8', index=False, header=False -) -(result_df.T).to_csv( - 'node.out', sep=' ', encoding='utf-8', index=False, header=False -) - -# make results.out dataframe with 3 columns and one row, no header. Each element is separated by a space -# results_df = pd.DataFrame({'total_impulse':vF[-1], 'max_force':vM[-1], 'total_disp':vF0[-1]}, index=[0]) -# results_df.to_csv('results.out', sep=' ', encoding='utf-8', header=False, index=False) - - -def main(df=None): # noqa: D103 - return df - - -if __name__ == '__main__': - """Entry point for the script.""" - parser = argparse.ArgumentParser( - description='Compute inline/total hydrodynamic force and moments on a monopile using Morisons equation' - ) - - parser.add_argument( - '-hw', - '--water_depth', - type=float, - default=30.0, - help='Water depth [m]', - ) - parser.add_argument( - '-Tp', - '--peak_period', - type=float, - default=12.7, - help='Wave period [s]', - ) - parser.add_argument( - '-Hs', - '--significant_wave_height', - type=float, - default=5.0, - help='Significant wave height [m]', - ) - parser.add_argument( - '-Dp', - '--pile_diameter', - type=float, - default=1.0, - help='Monopile diameter [m]', - ) - parser.add_argument( - '-Cd', - '--drag_coefficient', - type=float, - default=2.1, - help='Drag coefficient', - ) - parser.add_argument( - '-Cm', - '--mass_coefficient', - type=float, - default=2.0, - help='Mass coefficient', - ) - parser.add_argument( - '-nz', - '--number_of_recorders_z', - type=int, - default=4, - help='Number of points used in the z direction to compute loads', - ) - parser.add_argument('-t', '--time', type=float, default=1.0, help='Time [s]') - - arguments, unknowns = parser.parse_known_args() - - # hw = arguments.water_depth - # Tp = arguments.peak_period - # Hs = arguments.significant_wave_height - # # D = arguments.pile_diameter - # # CD = arguments.drag_coefficient - # # CM = arguments.mass_coefficient - # nz = arguments.number_of_recorders_z - # t = arguments.time - - # # --- Derived parameters - # h = hw - # # T = Tp - # f = 1./T - # g = 9.81 - # k = wavenumber(f, h, g) - # a = Hs - # z_ref = -h # reference point for moment calculation - # eps = 0 # phase shift [rad] - # rho = 1000 # water density - - # --- Wave kinematics - - # fig.show() - - # --------------------------------------------------------------------------------} - - # plt.suptitle('Hydro - Morison loads on monopile') - - # plt.savefig('MorisonLoads.png') - - # plt.show() - print('End of __main__ in Ex4_WaveLoads.py') # noqa: T201 - main() - - -if __name__ == '__test__': - pass - -if __name__ == '__export__': - # plt.close(fig1) - # plt.close(fig2) - # plt.close(fig) - from welib.tools.repo import export_figs_callback - - export_figs_callback(__file__) +#!/usr/bin/env python3 + +"""Compute inline/total hydrodynamic force and moments on a monopile using Morison's equation""" # noqa: D400 + +import argparse +from fractions import Fraction + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +# Local +from welib.tools.figure import defaultRC + +defaultRC() + +from welib.hydro.morison import inline_load # noqa: E402 +from welib.hydro.wavekin import elevation2d, kinematics2d, wavenumber # noqa: E402 +from welib.tools.colors import python_colors # noqa: E402 + +# --- Parameters +g = 9.80665 # gravity [m/s^2] +h = 30.0 # water depth [m] +rho = 1000 # water density +D = 6 # monopile diameter [m] +CD = 1 # drag coefficient +CM = 2 # added-mass coefficient +a = 3.0 # wave peak amplitude [m] +T = 12.0 # period [s] +eps = 0 # phase shift [rad] +f = 1.0 / T +k = wavenumber(f, h, g) + +nz = 30 # number of points used in the z direction to compute loads + +z_ref = -h # reference point for moment calculation + +# --------------------------------------------------------------------------------} +# --- Inline force and moments as function of time, with or without Wheeler stretching +# --------------------------------------------------------------------------------{ +time = np.linspace(0, T, 9) + +fig1, axes1 = plt.subplots( + 2, 4, sharex=True, sharey=True, figsize=(12.8, 4.8) +) # (6.4,4.8) +fig1.subplots_adjust( + left=0.05, right=0.99, top=0.95, bottom=0.09, hspace=0.26, wspace=0.11 +) + +fig2, axes2 = plt.subplots( + 2, 4, sharex=True, sharey=True, figsize=(12.8, 4.8) +) # (6.4,4.8) +fig2.subplots_adjust( + left=0.05, right=0.99, top=0.95, bottom=0.09, hspace=0.26, wspace=0.11 +) + +XLIM = [-75, 75] # For inline force +XLIMM = [-2500, 2500] # For inline moment + +for it, t in enumerate(time[:-1]): + # Wave kinematics + eta = elevation2d(a, f, k, eps, t, x=0) + z = np.linspace(-h, eta, nz) + u, du = kinematics2d(a, f, k, eps, h, t, z, Wheeler=True, eta=eta) + u0, du0 = kinematics2d(a, f, k, eps, h, t, z) + # Wave loads with wheeler + p_tot = inline_load(u, du, D, CD, CM, rho) + p_inertia = inline_load(u, du, D, CD * 0, CM, rho) + p_drag = inline_load(u, du, D, CD, CM * 0, rho) + dM = p_tot * (z - z_ref) # [Nm/m] # noqa: N816 + + # Wave loads without Wheeler + p_tot0 = inline_load(u0, du0, D, CD, CM, rho) + p_inertia0 = inline_load(u0, du0, D, CD * 0, CM, rho) + p_drag0 = inline_load(u0, du0, D, CD, CM * 0, rho) + dM0 = p_tot0 * (z - z_ref) # [Nm/m] # noqa: N816 + + # Plot inline force + ax = axes1[int(it / 4), np.mod(it, 4)] + ax.plot(p_inertia / 1000, z, '-', c=python_colors(0), label=r'$f_{inertia}$') + ax.plot(p_drag / 1000, z, '-', c=python_colors(3), label=r'$f_{drag}$') + ax.plot(p_tot / 1000, z, 'k-', label=r'$f_{tot}$') + ax.plot(p_inertia0 / 1000, z, '+', c=python_colors(0)) + ax.plot(p_drag0 / 1000, z, '+', c=python_colors(3)) + ax.plot(p_tot0 / 1000, z, 'k+') + ax.set_title(f't/T={Fraction(t / T)}') + if it == 0: + ax.legend() + ax.plot(XLIM, [0, 0], 'k') + ax.plot(XLIM, [a, a], 'k--') + ax.plot(XLIM, [-a, -a], 'k--') + + # Plot inline moment + ax = axes2[int(it / 4), np.mod(it, 4)] + ax.plot(dM / 1000, z, 'k-', label=r'$dM_{tot}$ with Wheeler') + ax.plot(dM0 / 1000, z, 'k+', label=r'$dM_{tot}$ no-correction') + ax.set_title(f't/T={Fraction(t / T)}') + if it == 0: + ax.legend() + ax.plot(XLIMM, [0, 0], 'k') + ax.plot(XLIMM, [a, a], 'k--') + ax.plot(XLIMM, [-a, -a], 'k--') + + +axes1[0, 0].set_xlim(XLIM) +axes1[0, 0].set_ylim([-h, a + 1]) +axes1[0, 0].set_ylabel('Depth z [m]') +axes1[1, 0].set_ylabel('Depth z [m]') +for i in range(4): + axes1[1, i].set_xlabel('Inline force [kN/m]') + + +fig1.savefig('forces.png') +# fig1.savefig('forces.webp') +# fig1.show() + +axes2[0, 0].set_xlim(XLIMM) +axes2[0, 0].set_ylim([-h, a + 1]) +axes2[0, 0].set_ylabel('Depth z [m]') +axes2[1, 0].set_ylabel('Depth z [m]') +for i in range(4): + axes2[1, i].set_xlabel('Inline moment [kNm/m]') + +fig2.savefig('moments.png') +# fig2.savefig('moments.webp') +# fig2.show() +# --------------------------------------------------------------------------------} +# --- Integrated force and sea bed moment over a period +# --------------------------------------------------------------------------------{ +time = np.linspace(0, 60.0, 6001) + +veta = np.zeros(time.shape) +vF = np.zeros(time.shape) # noqa: N816 +vM = np.zeros(time.shape) # noqa: N816 +vF0 = np.zeros(time.shape) # noqa: N816 +vM0 = np.zeros(time.shape) # noqa: N816 + +XLIM = [-75, 75] # For inline force +XLIMM = [-2500, 2500] # For inline moment + +elevation = np.zeros((len(time), nz)) +velocity = np.zeros((len(time), nz)) +accel = np.zeros((len(time), nz)) +force = np.zeros((len(time), nz)) + +for it, t in enumerate(time): + # Wave kinematics + veta[it] = elevation2d(a, f, k, eps, t, x=0) + z = np.linspace(-h, veta[it], nz) + u, du = kinematics2d(a, f, k, eps, h, t, z, Wheeler=True, eta=veta[it]) + u0, du0 = kinematics2d(a, f, k, eps, h, t, z) + # Wave loads with Wheeler + p_tot = inline_load(u, du, D, CD, CM, rho) + vF[it] = np.trapz(p_tot, z) # [N] # noqa: NPY201 + vM[it] = np.trapz(p_tot * (z - z_ref), z) # [Nm] # noqa: NPY201 + # Wave loads without Wheeler + p_tot0 = inline_load(u0, du0, D, CD, CM, rho) + vF0[it] = np.trapz(p_tot0, z) # [N] # noqa: NPY201 + vM0[it] = np.trapz(p_tot0 * (z - z_ref), z) # [Nm] # noqa: NPY201 + + elevation[it, :] = z.copy() + velocity[it, :] = u.copy() + accel[it, :] = du.copy() + force[it, :] = p_tot.copy() + +# Plot +fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6.4, 4.8)) # (6.4,4.8) +fig.subplots_adjust( + left=0.12, right=0.95, top=0.95, bottom=0.11, hspace=0.14, wspace=0.20 +) +# axes[0] = axes[0] +axes[0].plot(time / T, veta, 'k-') +axes[0].set_ylabel('Elevation [m]') +axes[0].grid(True) # noqa: FBT003 +# axes[1] = axes[1] +axes[1].plot(time / T, vF0 / 1e6, label='Standard') +axes[1].plot(time / T, vF / 1e6, 'k-', label='Wheeler Correction') +axes[1].set_ylabel('Streamwise Load, Cumulative [MN]') +axes[1].legend() +axes[1].grid(True) # noqa: FBT003 +# axes[2] = axes[2] +axes[2].plot(time / T, vM0 / 1e6, label='Standard') +axes[2].plot(time / T, vM / 1e6, 'k-', label='Wheeler Correction') +axes[2].set_ylabel('Sea-Bed Moment [MNm]') +axes[2].set_xlabel('Dimensionless Time, t/T [-]') +axes[2].legend() +axes[2].grid(True) # noqa: FBT003 + +# fig.savefig('IntegratedPileLoads.png') +# fig.savefig('IntegratedPileLoads.webp') +fig.savefig('IntegratedPileLoads.png') +# fig.show() + +# now save csv of the velocity, acceleration, force, and moment +veta_df = pd.DataFrame() +u_df = pd.DataFrame() +du_df = pd.DataFrame() +for i in range(nz): + dof = 1 + name = 'Disp_' + str(i + 1) + '_' + str(dof) + veta_df[name] = elevation[:, i] + name = 'Vel_' + str(i + 1) + '_' + str(dof) + u_df[name] = velocity[:, i] + name = 'RMSA_' + str(i + 1) + '_' + str(dof) + du_df[name] = accel[:, i] + +# add column per each force recorder +result_df = pd.DataFrame() +for i in range(nz): + dof = 1 + # name = 'Node_' + str(i+1) + '_' + str(dof) + name = 'Force_' + str(i + 1) + '_' + str(dof) + result_df[name] = force[:, i] + +# write columns to columns in csv files +(veta_df.T).to_csv('disp.evt', sep=' ', encoding='utf-8', index=False, header=False) +(u_df.T).to_csv('vel.evt', sep=' ', encoding='utf-8', index=False, header=False) +(du_df.T).to_csv('accel.evt', sep=' ', encoding='utf-8', index=False, header=False) +(result_df.T).to_csv( + 'forces.evt', sep=' ', encoding='utf-8', index=False, header=False +) + +# write columns to columns in csv files +(veta_df.T).to_csv('disp.out', sep=' ', encoding='utf-8', index=False, header=False) +(u_df.T).to_csv('vel.out', sep=' ', encoding='utf-8', index=False, header=False) +(du_df.T).to_csv('accel.out', sep=' ', encoding='utf-8', index=False, header=False) +(result_df.T).to_csv( + 'forces.out', sep=' ', encoding='utf-8', index=False, header=False +) +(result_df.T).to_csv( + 'node.out', sep=' ', encoding='utf-8', index=False, header=False +) + +# make results.out dataframe with 3 columns and one row, no header. Each element is separated by a space +# results_df = pd.DataFrame({'total_impulse':vF[-1], 'max_force':vM[-1], 'total_disp':vF0[-1]}, index=[0]) +# results_df.to_csv('results.out', sep=' ', encoding='utf-8', header=False, index=False) + + +def main(df=None): # noqa: D103 + return df + + +if __name__ == '__main__': + """Entry point for the script.""" + parser = argparse.ArgumentParser( + description='Compute inline/total hydrodynamic force and moments on a monopile using Morisons equation' + ) + + parser.add_argument( + '-hw', + '--water_depth', + type=float, + default=30.0, + help='Water depth [m]', + ) + parser.add_argument( + '-Tp', + '--peak_period', + type=float, + default=12.7, + help='Wave period [s]', + ) + parser.add_argument( + '-Hs', + '--significant_wave_height', + type=float, + default=5.0, + help='Significant wave height [m]', + ) + parser.add_argument( + '-Dp', + '--pile_diameter', + type=float, + default=1.0, + help='Monopile diameter [m]', + ) + parser.add_argument( + '-Cd', + '--drag_coefficient', + type=float, + default=2.1, + help='Drag coefficient', + ) + parser.add_argument( + '-Cm', + '--mass_coefficient', + type=float, + default=2.0, + help='Mass coefficient', + ) + parser.add_argument( + '-nz', + '--number_of_recorders_z', + type=int, + default=4, + help='Number of points used in the z direction to compute loads', + ) + parser.add_argument('-t', '--time', type=float, default=1.0, help='Time [s]') + + arguments, unknowns = parser.parse_known_args() + + # hw = arguments.water_depth + # Tp = arguments.peak_period + # Hs = arguments.significant_wave_height + # # D = arguments.pile_diameter + # # CD = arguments.drag_coefficient + # # CM = arguments.mass_coefficient + # nz = arguments.number_of_recorders_z + # t = arguments.time + + # # --- Derived parameters + # h = hw + # # T = Tp + # f = 1./T + # g = 9.81 + # k = wavenumber(f, h, g) + # a = Hs + # z_ref = -h # reference point for moment calculation + # eps = 0 # phase shift [rad] + # rho = 1000 # water density + + # --- Wave kinematics + + # fig.show() + + # --------------------------------------------------------------------------------} + + # plt.suptitle('Hydro - Morison loads on monopile') + + # plt.savefig('MorisonLoads.png') + + # plt.show() + print('End of __main__ in Ex4_WaveLoads.py') # noqa: T201 + main() + + +if __name__ == '__test__': + pass + +if __name__ == '__export__': + # plt.close(fig1) + # plt.close(fig2) + # plt.close(fig) + from welib.tools.repo import export_figs_callback + + export_figs_callback(__file__) diff --git a/modules/tools/ShakerMaker/ShakerMakersubmitjob.py b/modules/tools/ShakerMaker/ShakerMakersubmitjob.py index 264a75047..198ea3d99 100644 --- a/modules/tools/ShakerMaker/ShakerMakersubmitjob.py +++ b/modules/tools/ShakerMaker/ShakerMakersubmitjob.py @@ -1,13 +1,16 @@ -# %% -from tapipy.tapis import Tapis +# %% # noqa: INP001, D100 +import argparse +import json import os import time -import json -import argparse -#get the number of cores from argparser and runtime +from tapipy.tapis import Tapis + +# get the number of cores from argparser and runtime parser = argparse.ArgumentParser(description='Submit a job to DesignSafe') -parser.add_argument('--metadata', type=str, help='metadatafile for submitting the job') +parser.add_argument( + '--metadata', type=str, help='metadatafile for submitting the job' +) parser.add_argument('--tapisfolder', type=str, help='folder to upload the files to') parser.add_argument('--username', type=str, help='username for DesignSafe') parser.add_argument('--password', type=str, help='password for DesignSafe') @@ -15,56 +18,62 @@ args = parser.parse_args() # load the metadata file -with open(args.metadata + "/metadata.json") as json_file: +with open(args.metadata + '/metadata.json') as json_file: # noqa: PTH123 metadata = json.load(json_file) -jobinfo = metadata['jobdata'] -maxminutes = jobinfo['maxruntime'] -corespernode = jobinfo['numCores'] -numnodes = jobinfo['numNodes'] -totalcores = str(corespernode*numnodes) -queue = jobinfo['queue'] +jobinfo = metadata['jobdata'] +maxminutes = jobinfo['maxruntime'] +corespernode = jobinfo['numCores'] +numnodes = jobinfo['numNodes'] +totalcores = str(corespernode * numnodes) +queue = jobinfo['queue'] # chane hh:mm:ss to minutes -maxminutes = int(maxminutes.split(":")[0]) * 60 + int(maxminutes.split(":")[1]) + int(maxminutes.split(":")[2]) / 60 +maxminutes = ( + int(maxminutes.split(':')[0]) * 60 + + int(maxminutes.split(':')[1]) + + int(maxminutes.split(':')[2]) / 60 +) maxminutes = int(maxminutes) -systemcodepath = args.tapisfolder -faultinfo = metadata['faultdata'] -faultfiles = faultinfo["Faultfilenames"] +systemcodepath = args.tapisfolder +faultinfo = metadata['faultdata'] +faultfiles = faultinfo['Faultfilenames'] filenames = [] for filename in faultfiles: # // seperate the filename from the whole path - filename = os.path.basename(filename) + filename = os.path.basename(filename) # noqa: PTH119, PLW2901 filenames.append(filename) -faultinfo["Faultfilenames"] = filenames +faultinfo['Faultfilenames'] = filenames # Write the faualt info in the faultinfo folder -with open(args.metadata + "/faultInfo.json", "w") as outfile: +with open(args.metadata + '/faultInfo.json', 'w') as outfile: # noqa: PTH123 json.dump(faultinfo, outfile, indent=4) # write the metadata file -metadata["faultdata_path"] = ""; +metadata['faultdata_path'] = '' # delete the faultdata key -del metadata["faultdata"] -del metadata["jobdata"] +del metadata['faultdata'] +del metadata['jobdata'] -with open(args.metadata + "/metadata.json", "w") as outfile: +with open(args.metadata + '/metadata.json', 'w') as outfile: # noqa: PTH123 json.dump(metadata, outfile, indent=4) # print("username: ", args.username) # print("password: ", args.password) # exit() -print("Connecting to DesignSafe") -t = Tapis(base_url= "https://designsafe.tapis.io", - username=args.username, - password=args.password) +print('Connecting to DesignSafe') # noqa: T201 +t = Tapis( + base_url='https://designsafe.tapis.io', + username=args.username, + password=args.password, +) # Call to Tokens API to get access token t.get_tokens() @@ -73,97 +82,127 @@ # ============================================================================= # uploding files # ============================================================================= -print("Uploading files to DesignSafe") +print('Uploading files to DesignSafe') # noqa: T201 # create shakermaker directory if it does not exist -codepath = f"{t.username}/Shakermaker" -t.files.mkdir(systemId="designsafe.storage.default", path=codepath) +codepath = f'{t.username}/Shakermaker' +t.files.mkdir(systemId='designsafe.storage.default', path=codepath) for filename in os.listdir(systemcodepath): - filedata = open(f"{systemcodepath}/{filename}", "r").read() - t.files.insert(systemId="designsafe.storage.default", path=f"{codepath}/{filename}", file=filedata) + filedata = open(f'{systemcodepath}/{filename}').read() # noqa: SIM115, PTH123 + t.files.insert( + systemId='designsafe.storage.default', + path=f'{codepath}/{filename}', + file=filedata, + ) for filename in faultfiles: - filedata = open(filename, "r").read() - t.files.insert(systemId="designsafe.storage.default", path=f"{codepath}/{os.path.basename(filename)}", file=filedata) + filedata = open(filename).read() # noqa: SIM115, PTH123 + t.files.insert( + systemId='designsafe.storage.default', + path=f'{codepath}/{os.path.basename(filename)}', # noqa: PTH119 + file=filedata, + ) # upload the source time function -t.files.insert(systemId="designsafe.storage.default", path=f"{codepath}/SourceTimeFunction.py", file=open(f"{args.metadata}fault/SourceTimeFunction.py", "r").read()) -t.files.insert(systemId="designsafe.storage.default", path=f"{codepath}/metadata.json", file=open(args.metadata + "metadata.json", "r").read()) -t.files.insert(systemId="designsafe.storage.default", path=f"{codepath}/faultInfo.json", file=open(args.metadata + "faultInfo.json", "r").read()) +t.files.insert( + systemId='designsafe.storage.default', + path=f'{codepath}/SourceTimeFunction.py', + file=open(f'{args.metadata}fault/SourceTimeFunction.py').read(), # noqa: SIM115, PTH123 +) +t.files.insert( + systemId='designsafe.storage.default', + path=f'{codepath}/metadata.json', + file=open(args.metadata + 'metadata.json').read(), # noqa: SIM115, PTH123 +) +t.files.insert( + systemId='designsafe.storage.default', + path=f'{codepath}/faultInfo.json', + file=open(args.metadata + 'faultInfo.json').read(), # noqa: SIM115, PTH123 +) # ============================================================================= # submit a tapi job # ============================================================================= -print("Submitting a job to DesignSafe") -day = time.strftime("%Y_%m_%d") -time_ = time.strftime("%H_%M_%S") +print('Submitting a job to DesignSafe') # noqa: T201 +day = time.strftime('%Y_%m_%d') +time_ = time.strftime('%H_%M_%S') # attach the date and time to the username -jobname = f"EEUQ_DesignSafe_{t.username}_{day}_{time_}" -archivepath = f"{t.username}/tapis-jobs-archive/{jobname}" +jobname = f'EEUQ_DesignSafe_{t.username}_{day}_{time_}' +archivepath = f'{t.username}/tapis-jobs-archive/{jobname}' urls = [ - f"tapis://designsafe.storage.default/{t.username}/Shakermaker/ShakerMakermodel.py", - f"tapis://designsafe.storage.default/{t.username}/Shakermaker/metadata.json", - f"tapis://designsafe.storage.default/{t.username}/Shakermaker/faultInfo.json", - f"tapis://designsafe.storage.default/{t.username}/Shakermaker/SourceTimeFunction.py" + f'tapis://designsafe.storage.default/{t.username}/Shakermaker/ShakerMakermodel.py', + f'tapis://designsafe.storage.default/{t.username}/Shakermaker/metadata.json', + f'tapis://designsafe.storage.default/{t.username}/Shakermaker/faultInfo.json', + f'tapis://designsafe.storage.default/{t.username}/Shakermaker/SourceTimeFunction.py', ] for filename in filenames: - urls.append(f"tapis://designsafe.storage.default/{t.username}/Shakermaker/{filename}") + urls.append( # noqa: PERF401 + f'tapis://designsafe.storage.default/{t.username}/Shakermaker/{filename}' + ) my_job = { - "name": jobname, - "appId": "Shakermaker-app-amnp95", - "appVersion": "1.0.0", - "jobType": "BATCH", - "execSystemId": "frontera", - "execSystemExecDir":"${JobWorkingDir}/jobs/", - "execSystemInputDir":"${JobWorkingDir}/inputs/", - "execSystemOutputDir":"${JobWorkingDir}/inputs/results/", - "archiveSystemId": "designsafe.storage.default", - "archiveSystemDir": archivepath, - "nodeCount": numnodes, - "coresPerNode": corespernode, - "maxMinutes": maxminutes, - "execSystemLogicalQueue": queue, - "archiveOnAppError": False, - "fileInputArrays":[{"sourceUrls":urls,"targetDir":"*"}], - "parameterSet": { - "schedulerOptions": [{ "arg": "-A DesignSafe-SimCenter" }], - "envVariables": [ - {"key":"inputFile","value":"ShakerMakermodel.py"}, - {"key":"numProcessors","value":totalcores} - ] - } + 'name': jobname, + 'appId': 'Shakermaker-app-amnp95', + 'appVersion': '1.0.0', + 'jobType': 'BATCH', + 'execSystemId': 'frontera', + 'execSystemExecDir': '${JobWorkingDir}/jobs/', + 'execSystemInputDir': '${JobWorkingDir}/inputs/', + 'execSystemOutputDir': '${JobWorkingDir}/inputs/results/', + 'archiveSystemId': 'designsafe.storage.default', + 'archiveSystemDir': archivepath, + 'nodeCount': numnodes, + 'coresPerNode': corespernode, + 'maxMinutes': maxminutes, + 'execSystemLogicalQueue': queue, + 'archiveOnAppError': False, + 'fileInputArrays': [{'sourceUrls': urls, 'targetDir': '*'}], + 'parameterSet': { + 'schedulerOptions': [{'arg': '-A DesignSafe-SimCenter'}], + 'envVariables': [ + {'key': 'inputFile', 'value': 'ShakerMakermodel.py'}, + {'key': 'numProcessors', 'value': totalcores}, + ], + }, } job_info = t.jobs.submitJob(**my_job) # ============================================================================= # get the job status # ============================================================================= -jobs = t.jobs.getJobList(limit=1, orderBy='lastUpdated(desc),name(asc)', computeTotal=True) +jobs = t.jobs.getJobList( + limit=1, orderBy='lastUpdated(desc),name(asc)', computeTotal=True +) job_info = jobs[0] uuid = job_info.uuid # ge the status every 10 seconds job_info = t.jobs.getJob(jobUuid=uuid) -previous_status = job_info.status -while job_info.status not in['FINISHED','FAILED', 'CANCELLED','CANCELLED', 'QUEUED', 'RUNNING']: +previous_status = job_info.status +while job_info.status not in [ + 'FINISHED', + 'FAILED', + 'CANCELLED', + 'CANCELLED', + 'QUEUED', + 'RUNNING', +]: job_info = t.jobs.getJob(jobUuid=uuid) if job_info.status != previous_status: - print(f"Job status: {job_info.status}") + print(f'Job status: {job_info.status}') # noqa: T201 previous_status = job_info.status time.sleep(10.0) if job_info.status == 'FINISHED': - print("Job finished successfully") + print('Job finished successfully') # noqa: T201 if job_info.status == 'FAILED': - print("Job failed") - print("please check the job logs and contact the developers") + print('Job failed') # noqa: T201 + print('please check the job logs and contact the developers') # noqa: T201 if job_info.status == 'CANCELLED': - print("Job cancelled") + print('Job cancelled') # noqa: T201 if job_info.status == 'QUEUED': - print("Job is submitted and is in the queue") - print("This can take several days according to the queue") - print("please wait for the job to finish") - print("you can check the job status through the designsafe portal") + print('Job is submitted and is in the queue') # noqa: T201 + print('This can take several days according to the queue') # noqa: T201 + print('please wait for the job to finish') # noqa: T201 + print('you can check the job status through the designsafe portal') # noqa: T201 if job_info.status == 'RUNNING': - print("Job is running") - print("This can take several hours") - print("please wait for the job to finish") - print("you can check the job status through the designsafe portal") + print('Job is running') # noqa: T201 + print('This can take several hours') # noqa: T201 + print('please wait for the job to finish') # noqa: T201 + print('you can check the job status through the designsafe portal') # noqa: T201 # %% - diff --git a/modules/tools/ShakerMaker/TapisFiles/ShakerMakermodel.py b/modules/tools/ShakerMaker/TapisFiles/ShakerMakermodel.py index 3f32f2c5e..5490f7653 100644 --- a/modules/tools/ShakerMaker/TapisFiles/ShakerMakermodel.py +++ b/modules/tools/ShakerMaker/TapisFiles/ShakerMakermodel.py @@ -1,46 +1,55 @@ -# %% +# %% # noqa: INP001, D100 ############################################################# # This is shaker Maker model created for runing simulations # # for the EE-UQ app. # # This code created by Amin Pakzad and Pedro Arduino based # # on the initial code by Jose Abell and Jorge Crempien # # date = September 2024 # -# ########################################################### +# ########################################################### -import os import json +import os + from geopy.distance import geodesic +from mpi4py import MPI from shakermaker import shakermaker from shakermaker.crustmodel import CrustModel -from shakermaker.pointsource import PointSource from shakermaker.faultsource import FaultSource -from shakermaker.slw_extensions import DRMHDF5StationListWriter +from shakermaker.pointsource import PointSource from shakermaker.sl_extensions import DRMBox +from shakermaker.slw_extensions import DRMHDF5StationListWriter from shakermaker.station import Station from shakermaker.stationlist import StationList -from mpi4py import MPI -def calculate_distances_with_direction(lat1, lon1, lat2, lon2): + +def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # noqa: D103 # Original points point1 = (lat1, lon1) - + # Calculate north-south distance (latitude changes, longitude constant) north_point = (lat2, lon1) north_south_distance = geodesic(point1, north_point).kilometers - north_south_direction = "north" if lat2 > lat1 else "south" - + north_south_direction = 'north' if lat2 > lat1 else 'south' + # Calculate east-west distance (longitude changes, latitude constant) west_point = (lat1, lon2) west_east_distance = geodesic(point1, west_point).kilometers - west_east_direction = "east" if lon2 > lon1 else "west" - + west_east_direction = 'east' if lon2 > lon1 else 'west' + # south is negative - north_south_distance = -north_south_distance if north_south_direction == "south" else north_south_distance + north_south_distance = ( + -north_south_distance + if north_south_direction == 'south' + else north_south_distance + ) # west is negative - west_east_distance = -west_east_distance if west_east_direction == "west" else west_east_distance + west_east_distance = ( + -west_east_distance if west_east_direction == 'west' else west_east_distance + ) return north_south_distance, west_east_distance + # ====================================================================================== # Code initialization # ====================================================================================== @@ -49,19 +58,19 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): nprocs = comm.Get_size() # reading the metdadatafile -metadata_file = "metadata.json" -with open(metadata_file) as f: +metadata_file = 'metadata.json' +with open(metadata_file) as f: # noqa: PTH123 metadata = json.load(f) # make results directory if it doesn't exist -results_dir = "results" +results_dir = 'results' if rank == 0: - if not os.path.exists("results"): - os.makedirs("results") + if not os.path.exists('results'): # noqa: PTH110 + os.makedirs('results') # noqa: PTH103 comm.barrier() # Define the source parameters -# For estimating wave arrival windows, -# assume the following maximum and +# For estimating wave arrival windows, +# assume the following maximum and # minimum propagation velocities Vs_min = 3.14 # Minimum shear wave velocity Vp_max = 8.00 # Maximum primary wave velocity @@ -71,16 +80,22 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # Shakermaker configuration # ====================================================================================== -_m = 0.001 # meters (ShakerMaker uses kilometers) -dt = metadata["analysisdata"]["dt"] # Time step -nfft = metadata["analysisdata"]["nfft"] # Number of samples in the record -dk = metadata["analysisdata"]["dk"] # (Wavelength space discretization) adjust using theory -tb = 0 # How much to "advance" the simulation window... no advance -tmin = metadata["analysisdata"]["tmin"] # Time when the final results start -tmax = metadata["analysisdata"]["tmax"] # Time when the final results end -delta_h = metadata["analysisdata"]["dh"]*_m # Horizontal distance increment -delta_v_rec = metadata["analysisdata"]["delta_v_rec"]*_m # Vertical distance increment for receivers -delta_v_src = metadata["analysisdata"]["delta_v_src"]*_m # Vertical distance increment for sources +_m = 0.001 # meters (ShakerMaker uses kilometers) +dt = metadata['analysisdata']['dt'] # Time step +nfft = metadata['analysisdata']['nfft'] # Number of samples in the record +dk = metadata['analysisdata'][ + 'dk' +] # (Wavelength space discretization) adjust using theory +tb = 0 # How much to "advance" the simulation window... no advance +tmin = metadata['analysisdata']['tmin'] # Time when the final results start +tmax = metadata['analysisdata']['tmax'] # Time when the final results end +delta_h = metadata['analysisdata']['dh'] * _m # Horizontal distance increment +delta_v_rec = ( + metadata['analysisdata']['delta_v_rec'] * _m +) # Vertical distance increment for receivers +delta_v_src = ( + metadata['analysisdata']['delta_v_src'] * _m +) # Vertical distance increment for sources # nfft = 2048 # Number of samples in the record # dk = 0.2 # (Wavelength space discretization) adjust using theory @@ -100,7 +115,7 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # ====================================================================================== # Loading the crust model # ====================================================================================== -layers = metadata["crustdata"] +layers = metadata['crustdata'] num_layers = len(layers) CRUST = CrustModel(num_layers) for layer in layers: @@ -119,7 +134,7 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # faultpath = metadata["faultdata_path"] # make Fault Directory -fault_dir = "fault" +fault_dir = 'fault' # if rank == 0: # # remove the directory if it exists @@ -135,7 +150,7 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): comm.barrier() # copy the fault files to the fault directory -# src = faultpath +# src = faultpath # dst = fault_dir # # windows, linux or mac # if os.name == 'nt': # Windows @@ -156,17 +171,17 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): comm.barrier() # load the faultInfo.json file into faultdata -with open(f"faultInfo.json") as f: +with open('faultInfo.json') as f: # noqa: PTH123 faultdata = json.load(f) -faultLat = faultdata["latitude"] -faultLon = faultdata["longitude"] -filenames = faultdata["Faultfilenames"] -M0 = faultdata["M0"] -faultName = faultdata["name"] -xmean = faultdata["xmean"] -ymean = faultdata["ymean"] +faultLat = faultdata['latitude'] # noqa: N816 +faultLon = faultdata['longitude'] # noqa: N816 +filenames = faultdata['Faultfilenames'] +M0 = faultdata['M0'] +faultName = faultdata['name'] # noqa: N816 +xmean = faultdata['xmean'] +ymean = faultdata['ymean'] # if rank == 0: # for filename in filenames: @@ -182,72 +197,84 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): comm.barrier() # import SourceTimeFunction from fault file -from SourceTimeFunction import source_time_function - +from SourceTimeFunction import source_time_function # noqa: E402 for filename in filenames: sources = [] # read the json fault file - f = open(f"{filename}") + f = open(f'{filename}') # noqa: SIM115, PTH123 faultsources = json.load(f) for source in faultsources: - xsource = source["x"] - xmean - ysource = source["y"] - ymean - zsource = source["z"] - strike = source["strike"] - dip = source["dip"] - rake = source["rake"] - t0 = source["t0"] - stf = source["stf"] - stf_type = stf["type"] - params = stf["parameters"] - numparams = stf["numParameters"] - stf_func = source_time_function(*params) - sources.append(PointSource([xsource, ysource, zsource], [strike,dip,rake], tt=t0, stf=stf_func)) + xsource = source['x'] - xmean + ysource = source['y'] - ymean + zsource = source['z'] + strike = source['strike'] + dip = source['dip'] + rake = source['rake'] + t0 = source['t0'] + stf = source['stf'] + stf_type = stf['type'] + params = stf['parameters'] + numparams = stf['numParameters'] + stf_func = source_time_function(*params) + sources.append( + PointSource( + [xsource, ysource, zsource], [strike, dip, rake], tt=t0, stf=stf_func + ) + ) f.close() -FAULT = FaultSource(sources, - metadata = { - "name": f"{faultName} M0={M0}" - }) +FAULT = FaultSource(sources, metadata={'name': f'{faultName} M0={M0}'}) # ====================================================================================== # Loading the stations # ====================================================================================== -stationsType = metadata["stationdata"]["stationType"] +stationsType = metadata['stationdata']['stationType'] # noqa: N816 # single station -if stationsType.lower() in ["singlestation" ,"single"]: +if stationsType.lower() in ['singlestation', 'single']: stationslist = [] - for station in metadata["stationdata"]["Singlestations"]: - stationLat = station["latitude"] - stationLon = station["longitude"] - stationDepth = station["depth"] - meta = station["metadata"] - xstation, ystation = calculate_distances_with_direction(faultLat, faultLon, stationLat, stationLon) - stationslist.append(Station([xstation, ystation, stationDepth], metadata=meta)) - - meta = {"name": metadata["stationdata"]["name"]} - STATIONS = StationList(stationslist,metadata=meta) - -elif stationsType.lower() in ["drmbox", "drm", "drm box", "drm_box","drm station"]: - DRMdata = metadata["stationdata"]["DRMbox"] - name = DRMdata["name"] - latitude = DRMdata["latitude"] - longitude = DRMdata["longitude"] - depth = DRMdata["Depth"] - Lx = DRMdata["Width X"] - Ly = DRMdata["Width Y"] - Lz = DRMdata["Depth"] - dx = DRMdata["Mesh Size X"] - dy = DRMdata["Mesh Size Y"] - dz = DRMdata["Mesh Size Z"] - nx, ny, nz = int(Lx/dx), int(Ly/dy), int(Lz/dz) - dx = dx * _m; dy = dy * _m; dz = dz * _m - Lx = Lx * _m; Ly = Ly * _m; Lz = Lz * _m - xstation, ystation = calculate_distances_with_direction(faultLat, faultLon, latitude, longitude) - STATIONS = DRMBox([xstation, ystation, 0], [nx, ny, nz], [dx, dy, dz], metadata={"name": name}) + for station in metadata['stationdata']['Singlestations']: + stationLat = station['latitude'] # noqa: N816 + stationLon = station['longitude'] # noqa: N816 + stationDepth = station['depth'] # noqa: N816 + meta = station['metadata'] + xstation, ystation = calculate_distances_with_direction( + faultLat, faultLon, stationLat, stationLon + ) + stationslist.append( + Station([xstation, ystation, stationDepth], metadata=meta) + ) + + meta = {'name': metadata['stationdata']['name']} + STATIONS = StationList(stationslist, metadata=meta) + +elif stationsType.lower() in ['drmbox', 'drm', 'drm box', 'drm_box', 'drm station']: + DRMdata = metadata['stationdata']['DRMbox'] + name = DRMdata['name'] + latitude = DRMdata['latitude'] + longitude = DRMdata['longitude'] + depth = DRMdata['Depth'] + Lx = DRMdata['Width X'] + Ly = DRMdata['Width Y'] + Lz = DRMdata['Depth'] + dx = DRMdata['Mesh Size X'] + dy = DRMdata['Mesh Size Y'] + dz = DRMdata['Mesh Size Z'] + nx, ny, nz = int(Lx / dx), int(Ly / dy), int(Lz / dz) + dx = dx * _m + dy = dy * _m + dz = dz * _m + Lx = Lx * _m + Ly = Ly * _m + Lz = Lz * _m + xstation, ystation = calculate_distances_with_direction( + faultLat, faultLon, latitude, longitude + ) + STATIONS = DRMBox( + [xstation, ystation, 0], [nx, ny, nz], [dx, dy, dz], metadata={'name': name} + ) else: - raise ValueError(f"Unknown station type: {stationsType}") + raise ValueError(f'Unknown station type: {stationsType}') # noqa: EM102, TRY003 # ====================================================================================== @@ -255,13 +282,13 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # ====================================================================================== model = shakermaker.ShakerMaker(CRUST, FAULT, STATIONS) -if stationsType.lower() in ["drmbox", "drm", "drm box", "drm_box","drm station"]: +if stationsType.lower() in ['drmbox', 'drm', 'drm box', 'drm_box', 'drm station']: # creating the pairs model.gen_greens_function_database_pairs( - dt=dt, # Output time-step - nfft=nfft, # N timesteps - dk=dk, # wavenumber discretization - tb=tb, # Initial zero-padding + dt=dt, # Output time-step + nfft=nfft, # N timesteps + dk=dk, # wavenumber discretization + tb=tb, # Initial zero-padding tmin=tmin, tmax=tmax, smth=1, @@ -269,23 +296,23 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): verbose=True, debugMPI=False, showProgress=True, - store_here="results/greensfunctions_database", + store_here='results/greensfunctions_database', delta_h=delta_h, delta_v_rec=delta_v_rec, delta_v_src=delta_v_src, npairs_max=npairs_max, using_vectorize_manner=True, - cfactor=0.5 + cfactor=0.5, ) # wait for all processes to finish comm.barrier() model.run_create_greens_function_database( - h5_database_name="results/greensfunctions_database", - dt=dt, # Output time-step - nfft=nfft, # N timesteps - dk=dk, # wavenumber discretization - tb=tb, # Initial zero-padding + h5_database_name='results/greensfunctions_database', + dt=dt, # Output time-step + nfft=nfft, # N timesteps + dk=dk, # wavenumber discretization + tb=tb, # Initial zero-padding tmin=tmin, tmax=tmax, smth=1, @@ -297,13 +324,13 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # wait for all processes to finish comm.barrier() - writer = DRMHDF5StationListWriter("results/DRMLoad.h5drm") + writer = DRMHDF5StationListWriter('results/DRMLoad.h5drm') model.run_faster( - h5_database_name="results/greensfunctions_database", - dt=dt, # Output time-step - nfft=nfft, # N timesteps - dk=dk, # wavenumber discretization - tb=tb, # Initial zero-padding + h5_database_name='results/greensfunctions_database', + dt=dt, # Output time-step + nfft=nfft, # N timesteps + dk=dk, # wavenumber discretization + tb=tb, # Initial zero-padding tmin=tmin, tmax=tmax, smth=1, @@ -320,12 +347,12 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # single station -if stationsType.lower() in ["singlestation" ,"single"]: +if stationsType.lower() in ['singlestation', 'single']: model.run( - dt=dt, # Output time-step - nfft=nfft, # N timesteps - dk=dk, # wavenumber discretization - tb=tb, # Initial zero-padding + dt=dt, # Output time-step + nfft=nfft, # N timesteps + dk=dk, # wavenumber discretization + tb=tb, # Initial zero-padding tmin=tmin, tmax=tmax, smth=1, @@ -336,15 +363,5 @@ def calculate_distances_with_direction(lat1, lon1, lat2, lon2): ) if rank == 0: for i, s in enumerate(stationslist): - output_filename = f"results/station{i+1}.npz" + output_filename = f'results/station{i+1}.npz' s.save(output_filename) - - - - - - - - - - diff --git a/modules/tools/ShakerMaker/faultplotter.py b/modules/tools/ShakerMaker/faultplotter.py index 05b81444e..d1f065474 100644 --- a/modules/tools/ShakerMaker/faultplotter.py +++ b/modules/tools/ShakerMaker/faultplotter.py @@ -1,87 +1,95 @@ -# %% -import numpy as np -import matplotlib.pyplot as plt -import json -import os -from scipy.integrate import trapezoid -import pyvista as pv -from geopy.distance import geodesic +# %% # noqa: INP001, D100 import argparse import importlib.util -from pyproj import Transformer +import json +import os + import geopandas as gpd -from shapely.geometry import Point +import matplotlib.pyplot as plt +import numpy as np import pandas as pd import plotly.express as px +import pyvista as pv +from geopy.distance import geodesic +from pyproj import Transformer +from scipy.integrate import trapezoid +from shapely.geometry import Point -def load_function_from_file(filepath, function_name): + +def load_function_from_file(filepath, function_name): # noqa: D103 # Load the module from the given file path - spec = importlib.util.spec_from_file_location("module.name", filepath) + spec = importlib.util.spec_from_file_location('module.name', filepath) module = importlib.util.module_from_spec(spec) spec.loader.exec_module(module) - + # Get the function from the module return getattr(module, function_name) -def calculate_distances_with_direction(lat1, lon1, lat2, lon2): +def calculate_distances_with_direction(lat1, lon1, lat2, lon2): # noqa: D103 # Original points point1 = (lat1, lon1) - point2 = (lat2, lon2) - + point2 = (lat2, lon2) # noqa: F841 + # Calculate north-south distance (latitude changes, longitude constant) north_point = (lat2, lon1) north_south_distance = geodesic(point1, north_point).kilometers - north_south_direction = "north" if lat2 > lat1 else "south" - + north_south_direction = 'north' if lat2 > lat1 else 'south' + # Calculate east-west distance (longitude changes, latitude constant) west_point = (lat1, lon2) west_east_distance = geodesic(point1, west_point).kilometers - west_east_direction = "east" if lon2 > lon1 else "west" - + west_east_direction = 'east' if lon2 > lon1 else 'west' + # south is negative - north_south_distance = -north_south_distance if north_south_direction == "south" else north_south_distance + north_south_distance = ( + -north_south_distance + if north_south_direction == 'south' + else north_south_distance + ) # west is negative - west_east_distance = -west_east_distance if west_east_direction == "west" else west_east_distance + west_east_distance = ( + -west_east_distance if west_east_direction == 'west' else west_east_distance + ) return north_south_distance, west_east_distance -def PlotSources(info): - - tmpLocation = info["tmpLocation"] - faultinfopath = tmpLocation + "/fault/faultInfo.json" - faultinfo = json.load(open(faultinfopath,"r")) +def PlotSources(info): # noqa: C901, N802, D103 + tmpLocation = info['tmpLocation'] # noqa: N806 + faultinfopath = tmpLocation + '/fault/faultInfo.json' + faultinfo = json.load(open(faultinfopath)) # noqa: SIM115, PTH123 # // load the source time fuction in the faultinfo path - sourcetimeinfopath = tmpLocation + "/fault/SourceTimeFunction.py" + sourcetimeinfopath = tmpLocation + '/fault/SourceTimeFunction.py' # // load the source time function - source_time_function = load_function_from_file(sourcetimeinfopath, "source_time_function") - + source_time_function = load_function_from_file( + sourcetimeinfopath, 'source_time_function' + ) - numFaults = len(faultinfo["Faultfilenames"]) - faultFiles = faultinfo["Faultfilenames"] + numFaults = len(faultinfo['Faultfilenames']) # noqa: N806 + faultFiles = faultinfo['Faultfilenames'] # noqa: N806 - xfault = faultinfo["xmean"] - yfault = faultinfo["ymean"] + xfault = faultinfo['xmean'] + yfault = faultinfo['ymean'] # filteringData = info["filteringData"] # dataType = info["dataType"] # filteringRange = info["filteringRange"] - + if numFaults < 1: - print("No faults to plot") + print('No faults to plot') # noqa: T201 return - + if numFaults != len(faultFiles): - print("Number of files does not match number of faults") + print('Number of files does not match number of faults') # noqa: T201 return - + meshlist = [] for i in range(numFaults): - faultFile = tmpLocation + "/fault/" + faultFiles[i] - sources = json.load(open(faultFile,"r")) + faultFile = tmpLocation + '/fault/' + faultFiles[i] # noqa: N806 + sources = json.load(open(faultFile)) # noqa: SIM115, PTH123 x = np.zeros(len(sources)) y = np.zeros(len(sources)) z = np.zeros(len(sources)) @@ -90,36 +98,32 @@ def PlotSources(info): c3 = np.zeros(len(sources)) c4 = np.zeros(len(sources)) c5 = np.zeros(len(sources)) - for j,source in enumerate(sources): - - sourcetype = source["stf"]["type"] - x[j] = source["x"] - y[j] = source["y"] - z[j] = source["z"] - c1[j] = source["strike"] - c2[j] = source["dip"] - c3[j] = source["rake"] - c4[j] = source["t0"] - parameters = source["stf"]["parameters"] + for j, source in enumerate(sources): + sourcetype = source['stf']['type'] # noqa: F841 + x[j] = source['x'] + y[j] = source['y'] + z[j] = source['z'] + c1[j] = source['strike'] + c2[j] = source['dip'] + c3[j] = source['rake'] + c4[j] = source['t0'] + parameters = source['stf']['parameters'] c5[j] = source_time_function(*parameters, get_slip=True) - - - - mesh = pv.PolyData(np.c_[x,y,z]) - mesh["strike"] = c1 - mesh["dip"] = c2 - mesh["rake"] = c3 - mesh["t0"] = c4 - mesh["slip"] = c5 + mesh = pv.PolyData(np.c_[x, y, z]) + mesh['strike'] = c1 + mesh['dip'] = c2 + mesh['rake'] = c3 + mesh['t0'] = c4 + mesh['slip'] = c5 # strikes_rad = np.radians(c1) # dips_rad = np.radians(c2) # rakes_rad = np.radians(c3) # # compute the vectors # s = np.array([-np.sin(strikes_rad), np.cos(strikes_rad), np.zeros_like(strikes_rad)]) # d = np.array([ - # np.cos(strikes_rad) * np.sin(dips_rad), - # np.sin(strikes_rad) * np.sin(dips_rad), + # np.cos(strikes_rad) * np.sin(dips_rad), + # np.sin(strikes_rad) * np.sin(dips_rad), # -np.cos(dips_rad) # ]) # # Compute the rake vectors @@ -128,52 +132,50 @@ def PlotSources(info): # r_magnitudes = np.linalg.norm(r, axis=0) # r_normalized = r / r_magnitudes # mesh["rake_vector"] = r_normalized.T - - meshlist.append(mesh) - - - + meshlist.append(mesh) - Mesh = meshlist[0] - for i in range(1,numFaults): - Mesh = Mesh.merge(meshlist[i]) - - xy = Mesh.points[:,:2] - Mesh.points = Mesh.points- np.array([xfault,yfault,0]) + Mesh = meshlist[0] # noqa: N806 + for i in range(1, numFaults): + Mesh = Mesh.merge(meshlist[i]) # noqa: N806 + xy = Mesh.points[:, :2] + Mesh.points = Mesh.points - np.array([xfault, yfault, 0]) # finding hypocenter hypocenterid = np.argmin(c4) hypocenter = Mesh.points[hypocenterid].copy() # off screen rendering # remote_rendering - pl = pv.Plotter(shape=(3,2),groups=[((0, np.s_[:]))]) - pl.subplot(0,0) + pl = pv.Plotter(shape=(3, 2), groups=[((0, np.s_[:]))]) + pl.subplot(0, 0) # without scalars Mesh.set_active_scalars(None) - pl.add_mesh(Mesh,color="blue") - + pl.add_mesh(Mesh, color='blue') # plotting the hypocenter zhead = 0 ztail = hypocenter[2] # create an arrow - arrow = pv.Arrow(start=(hypocenter[0],hypocenter[1],ztail), direction=[0,0,1],scale = (zhead-ztail)) - pl.add_mesh(arrow, color="black",label="Hypocenter") + arrow = pv.Arrow( + start=(hypocenter[0], hypocenter[1], ztail), + direction=[0, 0, 1], + scale=(zhead - ztail), + ) + pl.add_mesh(arrow, color='black', label='Hypocenter') # pl.add_mesh(pv.PolyData(hypocenter),color="black",point_size=20,render_points_as_spheres=True,label="Hypocenter") # pl.subplot(1) # # strike # pl.add_mesh(Mesh.copy(), scalars="strike", cmap="viridis",label=f"Fault {info['faultname']}") - pl.subplot(1,0) + pl.subplot(1, 0) # dip - pl.add_mesh(Mesh.copy(), scalars="dip", cmap="viridis") + pl.add_mesh(Mesh.copy(), scalars='dip', cmap='viridis') - pl.subplot(1,1) + pl.subplot(1, 1) # rake - pl.add_mesh(Mesh.copy(), scalars="rake", cmap="viridis") + pl.add_mesh(Mesh.copy(), scalars='rake', cmap='viridis') # Mesh2 = Mesh.copy() # arrows = Mesh2.glyph(orient="rake_vector",scale="slip") # arrows = Mesh2.glyph(orient="rake_vector",scale="slip") @@ -181,207 +183,207 @@ def PlotSources(info): # pl.add_mesh(arrows, color="black") # computinf the vectors of the rake using slip and strike - pl.subplot(2,0) + pl.subplot(2, 0) # t0 - Mesh.set_active_scalars("t0") - pl.add_mesh(Mesh.copy(), scalars="t0", cmap="viridis") + Mesh.set_active_scalars('t0') + pl.add_mesh(Mesh.copy(), scalars='t0', cmap='viridis') - pl.subplot(2,1) + pl.subplot(2, 1) # slip - Mesh.set_active_scalars("slip") - pl.add_mesh(Mesh.copy(), scalars="slip", cmap="viridis") - + Mesh.set_active_scalars('slip') + pl.add_mesh(Mesh.copy(), scalars='slip', cmap='viridis') # get the bounds - xmin, xmax, ymin, ymax, zmin, zmax = Mesh.bounds - - + xmin, xmax, ymin, ymax, zmin, zmax = Mesh.bounds - if info["plottingStation"].lower() in ["yes","true"]: - stationCoordinates = info["stationCoordinates"] + if info['plottingStation'].lower() in ['yes', 'true']: + stationCoordinates = info['stationCoordinates'] # noqa: N806 coords = [] xy2 = [] for coord in stationCoordinates: - north, east = calculate_distances_with_direction(faultinfo["latitude"],faultinfo["longitude"],coord[0],coord[1]) - coords.append([north,east,coord[2]]) - xy2.append([north + xfault ,east + yfault]) + north, east = calculate_distances_with_direction( + faultinfo['latitude'], faultinfo['longitude'], coord[0], coord[1] + ) + coords.append([north, east, coord[2]]) + xy2.append([north + xfault, east + yfault]) - stations = pv.PolyData(np.array(coords)) - xmin0, xmax0, ymin0, ymax0, zmin0, zmax0 = stations.bounds + xmin0, xmax0, ymin0, ymax0, zmin0, zmax0 = stations.bounds # update extreme values - xmin = min(xmin,xmin0) - xmax = max(xmax,xmax0) - ymin = min(ymin,ymin0) - ymax = max(ymax,ymax0) - zmin = min(zmin,zmin0) - zmax = max(zmax,zmax0) - + xmin = min(xmin, xmin0) + xmax = max(xmax, xmax0) + ymin = min(ymin, ymin0) + ymax = max(ymax, ymax0) + zmin = min(zmin, zmin0) + zmax = max(zmax, zmax0) + colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] - pl.subplot(0,0) + pl.subplot(0, 0) # if info["plottingStation"].lower() in ["yes","true"]: - # pl.add_mesh(stations, color="red", point_size=10, render_points_as_spheres=True, label="Stations") - - - if info["plotlayers"].lower() in ["yes","true"]: - thickness = info["thickness"] + # pl.add_mesh(stations, color="red", point_size=10, render_points_as_spheres=True, label="Stations") + + if info['plotlayers'].lower() in ['yes', 'true']: + thickness = info['thickness'] # insert zero at the begining - thickness.insert(0,0) + thickness.insert(0, 0) # cummulative thickness cummulative = np.cumsum(thickness) if zmax < cummulative[-1]: - zamx = cummulative[-1] + cummulative[-1]*0.8 - cummulative[-1] = zmax+ cummulative[-1]*0.8 - for j in range(len(thickness)-1): - pl.add_mesh(pv.Cube(bounds=[xmin,xmax,ymin,ymax,cummulative[j],cummulative[j+1]]),color=colors[j],opacity=0.2,label=f"Layer {j+1}") - - - - - if not os.path.exists(tmpLocation): - os.makedirs(tmpLocation) + zamx = cummulative[-1] + cummulative[-1] * 0.8 # noqa: F841 + cummulative[-1] = zmax + cummulative[-1] * 0.8 + for j in range(len(thickness) - 1): + pl.add_mesh( + pv.Cube( + bounds=[ + xmin, + xmax, + ymin, + ymax, + cummulative[j], + cummulative[j + 1], + ] + ), + color=colors[j], + opacity=0.2, + label=f'Layer {j+1}', + ) + + if not os.path.exists(tmpLocation): # noqa: PTH110 + os.makedirs(tmpLocation) # noqa: PTH103 pl.link_views() - pl.show_grid(xtitle='X (North)',ytitle='Y (East)', ztitle='Z (Depth)') + pl.show_grid(xtitle='X (North)', ytitle='Y (East)', ztitle='Z (Depth)') pl.add_axes() - pl.camera.up = (0,0,-1) + pl.camera.up = (0, 0, -1) pl.camera.elevation = 60 pl.camera.azimuth = -90 pl.reset_camera() pl.camera.view_frustum(1.0) pl.add_legend() - pl.export_html(f"{tmpLocation}/faults.html") + pl.export_html(f'{tmpLocation}/faults.html') - transformer = Transformer.from_crs(faultinfo["epsg"], "epsg:4326") - xy = xy *1000 + transformer = Transformer.from_crs(faultinfo['epsg'], 'epsg:4326') + xy = xy * 1000 # add (xfault,yfault) to the first row of xy - xy = np.vstack((xy,[xfault,yfault])) - x2,y2 = transformer.transform(xy[:,1], xy[:,0]) - geometry1 = [Point(xy) for xy in zip(x2,y2)] + xy = np.vstack((xy, [xfault, yfault])) + x2, y2 = transformer.transform(xy[:, 1], xy[:, 0]) + geometry1 = [Point(xy) for xy in zip(x2, y2)] gdf1 = gpd.GeoDataFrame(geometry=geometry1) - gdf1["type"] = "Fault" - gdf1.loc[gdf1.index[-1], "type"] = "Fault Center" + gdf1['type'] = 'Fault' + gdf1.loc[gdf1.index[-1], 'type'] = 'Fault Center' - if info["plottingStation"].lower() in ["yes","true"]: + if info['plottingStation'].lower() in ['yes', 'true']: xy2 = np.array(xy2) * 1000 - x2,y2 = transformer.transform(xy2[:,1], xy2[:,0]) - geometry2 = [Point(xy) for xy in zip(x2,y2)] + x2, y2 = transformer.transform(xy2[:, 1], xy2[:, 0]) + geometry2 = [Point(xy) for xy in zip(x2, y2)] gdf2 = gpd.GeoDataFrame(geometry=geometry2) - gdf2["type"] = "Station" + gdf2['type'] = 'Station' # make the first row of the type column of gdf2 to be "Fault Center" - - gdf = gpd.GeoDataFrame(pd.concat([gdf1,gdf2],ignore_index=True)) + gdf = gpd.GeoDataFrame(pd.concat([gdf1, gdf2], ignore_index=True)) # plot the data plotly else: gdf = gdf1 - + # use different colors for fault and stations - color_map = { - "Fault": "blue", - "Station": "red", - "Fault Center": "green" - } + color_map = {'Fault': 'blue', 'Station': 'red', 'Fault Center': 'green'} gdf['size'] = gdf['type'].apply(lambda x: 15 if x != 'Fault' else 7) # Plot with custom colors fig = px.scatter_mapbox( - gdf, - lat=gdf.geometry.x, - lon=gdf.geometry.y, - color=gdf["type"], - size=gdf["size"], - color_discrete_map=color_map # Apply the custom color map + gdf, + lat=gdf.geometry.x, + lon=gdf.geometry.y, + color=gdf['type'], + size=gdf['size'], + color_discrete_map=color_map, # Apply the custom color map ) - min_lat, max_lat = gdf.geometry.x.min(), gdf.geometry.x.max() - min_lon, max_lon = gdf.geometry.y.min(), gdf.geometry.y.max() + min_lat, max_lat = gdf.geometry.x.min(), gdf.geometry.x.max() # noqa: F841 + min_lon, max_lon = gdf.geometry.y.min(), gdf.geometry.y.max() # noqa: F841 # Update layout and display map fig.update_layout( - mapbox_style="open-street-map", - mapbox_zoom=10, # Initial zoom level (adjustable) - # mapbox_center={"lat": (min_lat + max_lat) / 2, "lon": (min_lon + max_lon) / 2}, # Center the map - # mapbox_bounds={"west": min_lon, "east": max_lon, "south": min_lat, "north": max_lat} # Set bounds + mapbox_style='open-street-map', + mapbox_zoom=10, # Initial zoom level (adjustable) + # mapbox_center={"lat": (min_lat + max_lat) / 2, "lon": (min_lon + max_lon) / 2}, # Center the map + # mapbox_bounds={"west": min_lon, "east": max_lon, "south": min_lat, "north": max_lat} # Set bounds ) # export the as html - fig.write_html(f"{tmpLocation}/faults_map.html") - + fig.write_html(f'{tmpLocation}/faults_map.html') - - -if __name__ == "__main__": +if __name__ == '__main__': info = { - "filteringData": "No", - "filteringRange": [25, 45], - "numsataions": 1, - "stationCoordinates": [[-33.482426, -70.505214,0]], - "plottingStation": "No", - "plotlayers": "No", - "thickness": [0.2,0.8,14.5,0], - "tmpLocation": "tmp", + 'filteringData': 'No', + 'filteringRange': [25, 45], + 'numsataions': 1, + 'stationCoordinates': [[-33.482426, -70.505214, 0]], + 'plottingStation': 'No', + 'plotlayers': 'No', + 'thickness': [0.2, 0.8, 14.5, 0], + 'tmpLocation': 'tmp', } - parser = argparse.ArgumentParser() - parser.add_argument("--tmplocation", help="output location for DRM",required=True) - parser.add_argument("--thickness", help="Thickness of layers",required=False) - parser.add_argument("--plotlayers", help="Plot layers",required=False) - parser.add_argument("--plotstations", help="Plot stations",required=False) - parser.add_argument("--numsataions", help="Number of stations",required=False) - parser.add_argument("--stationcoordinates", help="Station coordinates",required=False) + parser.add_argument( + '--tmplocation', help='output location for DRM', required=True + ) + parser.add_argument('--thickness', help='Thickness of layers', required=False) + parser.add_argument('--plotlayers', help='Plot layers', required=False) + parser.add_argument('--plotstations', help='Plot stations', required=False) + parser.add_argument('--numsataions', help='Number of stations', required=False) + parser.add_argument( + '--stationcoordinates', help='Station coordinates', required=False + ) # print input arguments - args = parser.parse_args() - + if args.tmplocation: - info["tmpLocation"] = args.tmplocation - + info['tmpLocation'] = args.tmplocation + if args.thickness: - info["thickness"] = [float(i) for i in args.thickness.split(";")[:-1]] + info['thickness'] = [float(i) for i in args.thickness.split(';')[:-1]] if args.plotlayers: - info["plotlayers"] = args.plotlayers - + info['plotlayers'] = args.plotlayers + if args.plotstations: - info["plottingStation"] = args.plotstations + info['plottingStation'] = args.plotstations if args.numsataions: - info["numstations"] = int(args.numsataions) + info['numstations'] = int(args.numsataions) - if info["plottingStation"].lower() in ["yes","true"] and args.stationcoordinates: + if ( + info['plottingStation'].lower() in ['yes', 'true'] + and args.stationcoordinates + ): # delete the first three characters stringcoords = args.stationcoordinates[:] - info["stationCoordinates"] = [[float(i) for i in j.split(",")] for j in stringcoords.split(";")[:-1]] - - - if info["plottingStation"].lower() in ["yes","true"] and not args.stationcoordinates: - print("Station coordinates are required") - exit() - - if info["plottingStation"].lower() in ["yes","true"] and len(info["stationCoordinates"]) != info["numstations"]: - print("Number of stations does not match number of coordinates") - exit() - - if info["plotlayers"].lower() == "yes" and not args.thickness: - print("Thickness of layers is required") - exit() - - if info["plotlayers"].lower() == "yes" and len(info["thickness"]) < 1: - print("At least one layer is required") - exit() - - + info['stationCoordinates'] = [ + [float(i) for i in j.split(',')] for j in stringcoords.split(';')[:-1] + ] + + if ( + info['plottingStation'].lower() in ['yes', 'true'] + and not args.stationcoordinates + ): + print('Station coordinates are required') # noqa: T201 + exit() # noqa: PLR1722 + + if ( + info['plottingStation'].lower() in ['yes', 'true'] + and len(info['stationCoordinates']) != info['numstations'] + ): + print('Number of stations does not match number of coordinates') # noqa: T201 + exit() # noqa: PLR1722 + + if info['plotlayers'].lower() == 'yes' and not args.thickness: + print('Thickness of layers is required') # noqa: T201 + exit() # noqa: PLR1722 + + if info['plotlayers'].lower() == 'yes' and len(info['thickness']) < 1: + print('At least one layer is required') # noqa: T201 + exit() # noqa: PLR1722 PlotSources(info) - - - - - - - - - -