Created
May 16, 2022 18:16
-
-
Save PageotD/7d211e80c18e58a8e1ba618553db65f0 to your computer and use it in GitHub Desktop.
Wrap for Geopsy-GPDC
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python | |
| # -*- coding: utf-8 -*- | |
| # ------------------------------------------------------------------- | |
| # Filename: gpdcwrap.py | |
| # Author: Damien Pageot | |
| # ------------------------------------------------------------------ | |
| """ | |
| Functions to use the Geopsy-gpdc engine. | |
| """ | |
| import numpy as np | |
| from ctypes import CDLL, c_int, c_float, byref, POINTER, c_double | |
| from numpy.ctypeslib import ndpointer, load_library | |
| QGPCOREWAVE_PATH = "/home/geopsy/Codes/geopsy-install/lib/libQGpCoreWave.so" | |
| import matplotlib.pyplot as plt | |
| libCoreWave = load_library('libQGpCoreWave', QGPCOREWAVE_PATH) | |
| def dispersion_curve_init(verbose): | |
| """ | |
| Initialize the dispersion curve calculation. | |
| :param verbose: integer, 0 minimal ouput, 1 verbose output | |
| """ | |
| libCoreWave.dispersion_curve_init_.argtypes = [POINTER(c_int)] | |
| libCoreWave.dispersion_curve_init_(byref(c_int(verbose))) | |
| return | |
| def dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group): | |
| """ | |
| Calculate the Rayleigh theoretical dispersion curve. | |
| :param nLayers: integer, number of layers | |
| :param h: double, thickness of layers (m) | |
| :param vp: double, P-wave velocity in each layer (m/s) | |
| :param vs: double, S-wave velocity in each layer (m/s) | |
| :param rho: double, density in each layer (kg/m3) | |
| :param nSamples: integer, number of frequency samples | |
| :param omega: double, angular frequencies (rad/s) | |
| :param nModes: integer, number of modes including fundamental | |
| :param slowness: double, output of slowness values | |
| :param group: integer, 0 for phase, 1 for group | |
| """ | |
| libCoreWave.dispersion_curve_rayleigh_.argtypes = [ POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int)] | |
| libCoreWave.dispersion_curve_rayleigh_(byref(c_int(nLayers)), | |
| h, | |
| vp, | |
| vs, | |
| rho, | |
| byref(c_int(nSamples)), | |
| omega, | |
| byref(c_int(nModes)), | |
| slowness, | |
| byref(c_int(group))) | |
| return | |
| def dispersion_curve_love(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group): | |
| """ | |
| Calculate the Love theoretical dispersion curve. | |
| :param nLayers: integer, number of layers | |
| :param h: double, thickness of layers (m) | |
| :param vp: double, P-wave velocity in each layer (m/s) | |
| :param vs: double, S-wave velocity in each layer (m/s) | |
| :param rho: double, density in each layer (kg/m3) | |
| :param nSamples: integer, number of frequency samples | |
| :param omega: double, angular frequencies (rad/s) | |
| :param nModes: integer, number of modes including fundamental | |
| :param slowness: double, output of slowness values | |
| :param group: integer, 0 for phase, 1 for group | |
| """ | |
| libCoreWave.dispersion_curve_love_.argtypes = [ POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int), | |
| ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'), | |
| POINTER(c_int)] | |
| libCoreWave.dispersion_curve_love_(byref(c_int(nLayers)), | |
| h, | |
| vp, | |
| vs, | |
| rho, | |
| byref(c_int(nSamples)), | |
| omega, | |
| byref(c_int(nModes)), | |
| slowness, | |
| byref(c_int(group))) | |
| return | |
| if __name__ == "__main__": | |
| # Import numpy | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| # Initialize dispersion curve calculation | |
| dispersion_curve_init(0) | |
| # Poisson ratio | |
| nu = 0.25 | |
| # Number of modes (fundamental only = 1) | |
| nModes = 1 | |
| # Phase velocity (group velocity =1) | |
| group = 0 | |
| # Number of layers | |
| nLayers = 4 | |
| # Layer thickness array | |
| h = np.array([2., 4., 8., 0.], dtype=np.float64) | |
| # P and S velocity arrays | |
| # With Vp = func(Vs, nu) | |
| vs = np.array([250., 500., 750., 1500.], dtype=np.float64) | |
| vp = vs * np.sqrt(2.*(1.-nu)/(1.-2.*nu)) | |
| # Density array | |
| rho = np.array([1900., 1900., 1900., 1900.], dtype=np.float64) | |
| # Frequency and angular frequency arrays | |
| freq = np.arange(0.5, 80.5, 0.5, dtype=np.float64) | |
| omega = 2.*np.pi*freq | |
| nSamples = len(freq) | |
| # Initialize the slowness array (output) | |
| slowness = np.zeros(len(freq), dtype=np.float64) | |
| # Calculate the dispersion curve | |
| dispersion_curve_rayleigh(nLayers, h, vp, vs, rho, nSamples, omega, nModes, slowness, group) | |
| # Plot | |
| plt.xlabel("Frequency (Hz)") | |
| plt.ylabel("Phase velocity (m/s)") | |
| plt.plot(freq, 1./slowness) | |
| plt.show() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment