diff --git a/.github/tests/uniform.py b/.github/tests/uniform.py new file mode 100644 index 0000000..7f11464 --- /dev/null +++ b/.github/tests/uniform.py @@ -0,0 +1,40 @@ +import pyfmm +import numpy as np + +xarr = np.arange(0, 100, 0.08) +yarr = np.arange(0, 50, 0.05) +zarr = np.array([0.0]) # 二维情况 + + +# 慢度场 +slw = np.ones((len(xarr), len(yarr), len(zarr)), dtype='f') +print(slw.shape) + +srcloc = [10, 20, 0.0] + +# FMM解 +FMMTT = pyfmm.travel_time_source( + srcloc, + xarr, yarr, zarr, slw) + +# FSM解 +FSMTT = pyfmm.travel_time_source( + srcloc, + xarr, yarr, zarr, slw, useFSM=True, FSMparallel=True) + +# 真实解 +xx, yy, zz = srcloc +real_TT = np.sqrt(((xarr-xx)**2)[:,None,None] + ((yarr-yy)**2)[None,:,None] + ((zarr-zz)**2)[None,None,:]) + +# 误差 +FMM_error = np.mean(np.abs(FMMTT - real_TT)) +FSM_error = np.mean(np.abs(FSMTT - real_TT)) +print("FMM_error = ", FMM_error) +print("FSM_error = ", FSM_error) + +tol = 0.1 + +if FMM_error > tol: + raise ValueError(f"FMM_error({FMM_error}) > tol({tol})") +if FSM_error > tol: + raise ValueError(f"FMM_error({FSM_error}) > tol({tol})") \ No newline at end of file diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..92a1273 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,329 @@ +name: Build C programs and libraries, publish to PyPI and create a release + +on: + push: + tags: + - "v*.*.*" + +jobs: + build: # 编译C库 + runs-on: ${{ matrix.os }} + strategy: + matrix: + include: + - os: windows-2022 + arch: x86_64 # windows x86_6 + - os: ubuntu-22.04 + arch: x86_64 # Ubuntu x86_64 + - os: macos-13 + arch: x86_64 # macOS Intel + - os: macos-14 + arch: arm64 # macOS Apple Silicon + + fail-fast: true + + defaults: + run: + shell: bash + + steps: + - name: Set MSYS2 (Windows) # windows平台使用MSYS2工具编译程序 + if: contains(matrix.os, 'windows') + uses: msys2/setup-msys2@v2 + with: + msystem: UCRT64 + install: >- + git + mingw-w64-ucrt-x86_64-gcc + make + + - name: Configured git attributes for line endings (Windows) # 防止接下来在windows平台checkout代码时，文本文件的换行符发生变化，导致MSYS2工作出错 + if: contains(matrix.os, 'windows') + run: git config --global core.autocrlf input + + - name: Checkout code # 下载库代码 + uses: actions/checkout@v3 + with: + fetch-depth: 1 + + # --------------------- 检查文件中的版本名是否和tag名相符 -------------------- + - name: Check version + run: | + TAG_NAME=${{ github.ref_name }} + PYVERSION_FILE="pyfmm/_version.py" + + # 提取 Python 版本号 + PYTHON_VERSION=$(sed -n 's/^__version__ = "\([0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*\)"/\1/p' ${PYVERSION_FILE}) + if [ -z "$PYTHON_VERSION" ]; then + echo "Error: Failed to extract version from ${PYVERSION_FILE}" + exit 1 + fi + + # 检查 Python 版本与 Tag 是否匹配 + if [ "$PYTHON_VERSION" != "${TAG_NAME#v}" ]; then + echo "Error: Version mismatch between tag (${TAG_NAME}) and ${PYVERSION_FILE} (${PYTHON_VERSION})" + exit 1 + fi + + echo "Version check passed:" + echo " Tag: ${TAG_NAME}" + echo " Python version: ${PYTHON_VERSION}" + + # --------------------------- 安装依赖 ------------------------------------------ + # - name: Install dependencies (Ubuntu) + # if: contains(matrix.os, 'ubuntu') + # run: | + # sudo apt update + # sudo apt install -y build-essential + + - name: Install dependencies (macOS) + # 匹配包含 "macos" 的操作系统 + if: contains(matrix.os, 'macos') + run: | + brew install libomp + + # ----------------- 编译C库和C程序 ---------------------------------- + - name: Build the project (macOS) + if: contains(matrix.os, 'macos') + working-directory: ./pyfmm/C_extension + run: | + make ARCH="-arch ${{ matrix.arch }}" \ + CC=gcc-14 + + make cleanbuild + otool -L lib/* + + - name: Build the project (Ubuntu) + if: contains(matrix.os, 'ubuntu') + working-directory: ./pyfmm/C_extension + run: | + make + make cleanbuild + ldd lib/* + + - name: Build the project (Windows) + if: contains(matrix.os, 'windows') + shell: msys2 {0} + working-directory: ./pyfmm/C_extension + run: | + make + make cleanbuild + ldd lib/* + + + # ------------------------ 定义接下来打包程序命名时的系统名后缀 --------------- + - name: Define the package OS suffix + run: | + # 符合pypi命名规范，否则上传失败 + if [[ "${{ matrix.os }}" == *"ubuntu"* ]]; then + SUFFIX_PLAT_NAME="manylinux2014_x86_64" + elif [[ "${{ matrix.os }}" == *"macos"* && "${{ matrix.arch }}" == *"x86_64"* ]]; then + SUFFIX_PLAT_NAME="macosx_10_9_x86_64" + elif [[ "${{ matrix.os }}" == *"macos"* && "${{ matrix.arch }}" == *"arm64"* ]]; then + SUFFIX_PLAT_NAME="macosx_11_0_arm64" + elif [[ "${{ matrix.os }}" == *"windows"* ]]; then + SUFFIX_PLAT_NAME="win_amd64" + else + echo " Unsupported OS: ${{ matrix.os }} (${{ matrix.arch }})" + exit 1 + fi + + echo "SUFFIX_PLAT_NAME=$SUFFIX_PLAT_NAME" >> $GITHUB_ENV + + # --------------------------- 打包整个程序 --------------------- + - name: Package the binary + run: | + PACK_NAME=pyfmm_kit-${{ github.ref_name }}-${{ env.SUFFIX_PLAT_NAME }} + echo "PACK_NAME=$PACK_NAME" >> $GITHUB_ENV + FILE_CONTENT=$(ls -A) + mkdir -p $PACK_NAME + cp -r ${FILE_CONTENT} $PACK_NAME/ + tar -czvf $PACK_NAME.tar.gz $PACK_NAME + rm -rf $PACK_NAME + + # -------------------- upload artifacts ----------------------- + - name: Upload artifact (*.tar.gz) + uses: actions/upload-artifact@v4 + with: + name: ${{ matrix.os }}-${{ matrix.arch }}_tar + path: ${{ env.PACK_NAME }}.tar.gz + + + # ======================================================================================= + test_project: # 在全新系统上测试程序，不安装其它依赖，看能否运行 + runs-on: ${{ matrix.os }} + needs: build + strategy: + matrix: + include: + - os: windows-2022 + arch: x86_64 # windows x86_6 + - os: ubuntu-22.04 + arch: x86_64 # Ubuntu x86_64 + - os: macos-13 + arch: x86_64 # macOS Intel + - os: macos-14 + arch: arm64 # macOS Apple Silicon + + fail-fast: true + + defaults: + run: + shell: bash + + steps: + - name: Download artifacts + uses: actions/download-artifact@v4 + with: + name: ${{ matrix.os }}-${{ matrix.arch }}_tar + path: artifacts + + - name: Display structure of downloaded files, and Uncompress + run: | + ls -R artifacts + echo "------------------- tar output -----------------------------" + tar -xzvf artifacts/*.tar.gz + echo "------------------------------------------------------------" + + # 获得压缩包解压出来的文件夹名 + PACK_NAME=$(ls | grep pyfmm_kit) + echo "PACK_NAME=$PACK_NAME" >> $GITHUB_ENV + + # 从解压出的文件夹命名来推断${{ env.SUFFIX_PLAT_NAME }} + SUFFIX_PLAT_NAME=$(echo $PACK_NAME | sed 's/.*-\(.*\)/\1/') + echo "SUFFIX_PLAT_NAME=$SUFFIX_PLAT_NAME" >> $GITHUB_ENV + + echo $PACK_NAME + echo $SUFFIX_PLAT_NAME + + # --------------------------- 安装依赖 ------------------------------------------ + # 实际使用时可能需要安装libomp + # - name: Install libomp (Ubuntu) + # if: contains(matrix.os, 'ubuntu') + # run: | + # sudo apt install -y libomp-dev + + # - name: Set alias (MacOS) + # if: contains(matrix.os, 'macos') + # run: | + # brew install coreutils + # echo "alias timeout=gtimeout" >> ~/.bashrc + + # --------------------搭建python环境，开始测试，并制作wheel文件 ------------------------------ + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: '3.9' + + - name: Install dependencies + working-directory: ${{ env.PACK_NAME }} + run: | + python -m pip install --upgrade pip + pip install --upgrade setuptools wheel build + pip install -v . + + - name: Clean up build and egg-info directories + working-directory: ${{ env.PACK_NAME }} + run: | + # 清理临时文件 + rm -rf build/ + rm -rf pyfmm_kit.egg-info/ + + - name: Test + working-directory: ${{ env.PACK_NAME }}/.github/tests + run: | + python uniform.py + + # --------------------------- 制作wheels --------------------- + - name: Build the Python Wheel + working-directory: ${{ env.PACK_NAME }} + run: | + python setup.py bdist_wheel --plat-name=${{ env.SUFFIX_PLAT_NAME }} # 只制作wheel，这里不打包源码 + + - name: Upload artifact (*.whl) + uses: actions/upload-artifact@v4 + with: + name: ${{ env.PACK_NAME }}_whl + path: ${{ env.PACK_NAME }}/dist/*.whl + + + # =============================================================================== + publish_pypi: # 上传pypi + runs-on: ubuntu-latest + needs: test_project + permissions: + id-token: write + contents: read + defaults: + run: + shell: bash + + + steps: + - name: Download artifacts + uses: actions/download-artifact@v4 + with: + pattern: "*_whl" + path: artifacts + + - name: Display structure of downloaded files + run: ls -R artifacts + + - name: Move wheel files to dist/ + run: | + mkdir -p dist + mv artifacts/*/*.whl dist/ + + - name: Pypi Publish + uses: pypa/gh-action-pypi-publish@release/v1 + with: + packages-dir: dist/ + + + # =============================================================================== + release: # 创建Release + runs-on: ubuntu-latest + needs: build # test_project + permissions: + contents: write + + defaults: + run: + shell: bash + + steps: + - name: Download artifacts + uses: actions/download-artifact@v4 + with: + path: artifacts + pattern: "*_tar" + + - name: Display structure of downloaded files + run: ls -R artifacts + + # 对artifacts中的每个tar.gz文件，根据其内部的文件夹名重新命名压缩包 + - name: Rename tar.gz files + run: | + mkdir -p release + for file in artifacts/*/*.tar.gz; do + TEMP_FILE=tmp + # 这里不能把tar合并到下方的管道中，在github actions中会报错 + tar -tzf $file > $TEMP_FILE + PACK_NAME=$(head -n 1 $TEMP_FILE | cut -f 1 -d '/') + mv $file release/${PACK_NAME}.tar.gz + rm $TEMP_FILE + done + + ls -R release + + - name: Create Release + id: create_release + uses: softprops/action-gh-release@v1 + with: + draft: true + files: "release/*.tar.gz" + + + + + \ No newline at end of file diff --git a/.github/workflows/publish_pypi.yml b/.github/workflows/publish_pypi.yml deleted file mode 100644 index 705022c..0000000 --- a/.github/workflows/publish_pypi.yml +++ /dev/null @@ -1,50 +0,0 @@ -name: Publish to PyPI via Trusted Publisher - -on: - push: - tags: - - "*" - -jobs: - - publish: - runs-on: ubuntu-latest - - permissions: - id-token: write - contents: read - - steps: - - name: Checkout code - uses: actions/checkout@v3 - with: - fetch-depth: 0 - - - name: Check version - run: | - TAG_NAME=$(echo ${{ github.ref_name }}) - PYVERSION_FILE="pyfmm/_version.py" - PYTHON_VERSION=$(grep -oP '__version__ = "\K[0-9]+\.[0-9]+\.[0-9]+' ${PYVERSION_FILE}) - - if [ "$PYTHON_VERSION" != "${TAG_NAME:1}" ]; then - echo "Version mismatch between tag and ${PYVERSION_FILE}" - exit 1 - fi - - - - name: Set up Python - uses: actions/setup-python@v4 - with: - python-version: "3.9" - - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install --upgrade setuptools wheel build - - - name: Build the package - run: | - python setup.py sdist # 暂不加bdist_wheel，先使用源码编译，后续会改进 - - - name: Pypi Publish - uses: pypa/gh-action-pypi-publish@release/v1 \ No newline at end of file diff --git a/.github/workflows/publish_testpypi.yml b/.github/workflows/publish_testpypi.yml deleted file mode 100644 index 9c752fb..0000000 --- a/.github/workflows/publish_testpypi.yml +++ /dev/null @@ -1,51 +0,0 @@ -name: Publish to TestPyPI - -on: - push: - tags: - - "*" - -jobs: - publish: - runs-on: ubuntu-latest - - permissions: - contents: read - - steps: - - name: Checkout code - uses: actions/checkout@v3 - with: - fetch-depth: 0 - - - name: Check version - run: | - TAG_NAME=$(echo ${{ github.ref_name }}) - PYVERSION_FILE="pyfmm/_version.py" - PYTHON_VERSION=$(grep -oP '__version__ = "\K[0-9]+\.[0-9]+\.[0-9]+' ${PYVERSION_FILE}) - - if [ "$PYTHON_VERSION" != "${TAG_NAME:1}" ]; then - echo "Version mismatch between tag and ${PYVERSION_FILE}" - exit 1 - fi - - - name: Set up Python - uses: actions/setup-python@v4 - with: - python-version: "3.9" - - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install --upgrade setuptools wheel build twine - - - name: Build the package - run: | - python setup.py sdist - - - name: Publish to TestPyPI - env: - TWINE_USERNAME: __token__ - TWINE_PASSWORD: ${{ secrets.TESTPYPI_API_TOKEN }} - run: | - python -m twine upload --repository testpypi dist/* diff --git a/.github/workflows/testbuild.yml b/.github/workflows/testbuild.yml new file mode 100644 index 0000000..1713f69 --- /dev/null +++ b/.github/workflows/testbuild.yml @@ -0,0 +1,228 @@ +name: Just test + +# 触发机制为“非主分支的提交” +on: + push: + branches-ignore: + - main # 排除主分支 + - master # 如果有 master 分支，也可以排除 + +jobs: + build: # 编译C库 + runs-on: ${{ matrix.os }} + strategy: + matrix: + include: + - os: windows-2022 + arch: x86_64 # windows x86_6 + - os: ubuntu-22.04 + arch: x86_64 # Ubuntu x86_64 + - os: macos-13 + arch: x86_64 # macOS Intel + - os: macos-14 + arch: arm64 # macOS Apple Silicon + + fail-fast: true + + defaults: + run: + shell: bash + + steps: + - name: Set MSYS2 (Windows) # windows平台使用MSYS2工具编译程序 + if: contains(matrix.os, 'windows') + uses: msys2/setup-msys2@v2 + with: + msystem: UCRT64 + install: >- + git + mingw-w64-ucrt-x86_64-gcc + make + + - name: Configured git attributes for line endings (Windows) # 防止接下来在windows平台checkout代码时，文本文件的换行符发生变化，导致MSYS2工作出错 + if: contains(matrix.os, 'windows') + run: git config --global core.autocrlf input + + - name: Checkout code # 下载库代码 + uses: actions/checkout@v3 + with: + fetch-depth: 1 + + # --------------------------- 安装依赖 ------------------------------------------ + # - name: Install dependencies (Ubuntu) + # if: contains(matrix.os, 'ubuntu') + # run: | + # sudo apt update + # sudo apt install -y build-essential + + - name: Install dependencies (macOS) + # 匹配包含 "macos" 的操作系统 + if: contains(matrix.os, 'macos') + run: | + brew install libomp + + # ----------------- 编译C库和C程序 ---------------------------------- + - name: Build the project (macOS) + if: contains(matrix.os, 'macos') + working-directory: ./pyfmm/C_extension + run: | + make ARCH="-arch ${{ matrix.arch }}" \ + CC=gcc-14 + + make cleanbuild + otool -L lib/* + + - name: Build the project (Ubuntu) + if: contains(matrix.os, 'ubuntu') + working-directory: ./pyfmm/C_extension + run: | + make + make cleanbuild + ldd lib/* + + - name: Build the project (Windows) + if: contains(matrix.os, 'windows') + shell: msys2 {0} + working-directory: ./pyfmm/C_extension + run: | + make + make cleanbuild + ldd lib/* + + + # ------------------------ 定义接下来打包程序命名时的系统名后缀 --------------- + - name: Define the package OS suffix + run: | + # 符合pypi命名规范，否则上传失败 + if [[ "${{ matrix.os }}" == *"ubuntu"* ]]; then + SUFFIX_PLAT_NAME="manylinux2014_x86_64" + elif [[ "${{ matrix.os }}" == *"macos"* && "${{ matrix.arch }}" == *"x86_64"* ]]; then + SUFFIX_PLAT_NAME="macosx_10_9_x86_64" + elif [[ "${{ matrix.os }}" == *"macos"* && "${{ matrix.arch }}" == *"arm64"* ]]; then + SUFFIX_PLAT_NAME="macosx_11_0_arm64" + elif [[ "${{ matrix.os }}" == *"windows"* ]]; then + SUFFIX_PLAT_NAME="win_amd64" + else + echo " Unsupported OS: ${{ matrix.os }} (${{ matrix.arch }})" + exit 1 + fi + + echo "SUFFIX_PLAT_NAME=$SUFFIX_PLAT_NAME" >> $GITHUB_ENV + + # --------------------------- 打包整个程序 --------------------- + - name: Package the binary + run: | + PACK_NAME=pyfmm_kit-${{ github.ref_name }}-${{ env.SUFFIX_PLAT_NAME }} + echo "PACK_NAME=$PACK_NAME" >> $GITHUB_ENV + FILE_CONTENT=$(ls -A) + mkdir -p $PACK_NAME + cp -r ${FILE_CONTENT} $PACK_NAME/ + tar -czvf $PACK_NAME.tar.gz $PACK_NAME + rm -rf $PACK_NAME + + # -------------------- upload artifacts ----------------------- + - name: Upload artifact (*.tar.gz) + uses: actions/upload-artifact@v4 + with: + name: ${{ matrix.os }}-${{ matrix.arch }}_tar + path: ${{ env.PACK_NAME }}.tar.gz + + + # ======================================================================================= + test_project: # 在全新系统上测试程序，不安装其它依赖，看能否运行 + runs-on: ${{ matrix.os }} + needs: build + strategy: + matrix: + include: + - os: windows-2022 + arch: x86_64 # windows x86_6 + - os: ubuntu-22.04 + arch: x86_64 # Ubuntu x86_64 + - os: macos-13 + arch: x86_64 # macOS Intel + - os: macos-14 + arch: arm64 # macOS Apple Silicon + + fail-fast: true + + defaults: + run: + shell: bash + + steps: + - name: Download artifacts + uses: actions/download-artifact@v4 + with: + name: ${{ matrix.os }}-${{ matrix.arch }}_tar + path: artifacts + + - name: Display structure of downloaded files, and Uncompress + run: | + ls -R artifacts + echo "------------------- tar output -----------------------------" + tar -xzvf artifacts/*.tar.gz + echo "------------------------------------------------------------" + + # 获得压缩包解压出来的文件夹名 + PACK_NAME=$(ls | grep pyfmm_kit) + echo "PACK_NAME=$PACK_NAME" >> $GITHUB_ENV + + # 从解压出的文件夹命名来推断${{ env.SUFFIX_PLAT_NAME }} + SUFFIX_PLAT_NAME=$(echo $PACK_NAME | sed 's/.*-\(.*\)/\1/') + echo "SUFFIX_PLAT_NAME=$SUFFIX_PLAT_NAME" >> $GITHUB_ENV + + echo $PACK_NAME + echo $SUFFIX_PLAT_NAME + + # --------------------------- 安装依赖 ------------------------------------------ + # 实际使用时可能需要安装libomp + # - name: Install libomp (Ubuntu) + # if: contains(matrix.os, 'ubuntu') + # run: | + # sudo apt install -y libomp-dev + + # - name: Set alias (MacOS) + # if: contains(matrix.os, 'macos') + # run: | + # brew install coreutils + # echo "alias timeout=gtimeout" >> ~/.bashrc + + # --------------------搭建python环境，开始测试，并制作wheel文件 ------------------------------ + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: '3.9' + + - name: Install dependencies + working-directory: ${{ env.PACK_NAME }} + run: | + python -m pip install --upgrade pip + pip install --upgrade setuptools wheel build + pip install -v . + + - name: Clean up build and egg-info directories + working-directory: ${{ env.PACK_NAME }} + run: | + # 清理临时文件 + rm -rf build/ + rm -rf pyfmm_kit.egg-info/ + + - name: Test + working-directory: ${{ env.PACK_NAME }}/.github/tests + run: | + python uniform.py + + + # --------------------------- 制作wheels --------------------- + # - name: Build the Python Wheel + # working-directory: ${{ env.PACK_NAME }} + # run: | + # python setup.py bdist_wheel --plat-name=${{ env.SUFFIX_PLAT_NAME }} # 只制作wheel，这里不打包源码 + + # - name: Upload artifact (*.whl) + # uses: actions/upload-artifact@v4 + # with: + # name: ${{ env.PACK_NAME }}_whl + # path: ${{ env.PACK_NAME }}/dist/*.whl + \ No newline at end of file diff --git a/pyfmm/C_extension/Makefile b/pyfmm/C_extension/Makefile index 3c37fa3..b56c9d5 100644 --- a/pyfmm/C_extension/Makefile +++ b/pyfmm/C_extension/Makefile @@ -1,87 +1,121 @@ -ifeq ($(OS),Windows_NT) -define RMDIR_FUNC - if exist $(1) rmdir /S /Q $(1) -endef -else -define RMDIR_FUNC - rm -rf $(1) -endef -endif +# ifeq ($(OS),Windows_NT) +# define RMDIR_FUNC +# if exist $(1) rmdir /S /Q $(1) +# endef +# else +# define RMDIR_FUNC +# rm -rf $(1) +# endef +# endif + +# ifeq ($(OS),Windows_NT) +# define MKDIR_FUNC +# if not exist $(1) mkdir $(1) +# endef +# else +# define MKDIR_FUNC +# mkdir -p $(1) +# endef +# endif + +# ifeq ($(OS),Windows_NT) +# FILESEP = \\ + +# else +# FILESEP = / +# endif + +SRC_DIR := src +INC_DIR := include +BUILD_DIR := build +BUILD_DIR_FLOAT = $(BUILD_DIR)/float +BUILD_DIR_DOUBLE = $(BUILD_DIR)/double +LIB_DIR = lib +LIB_NAME = $(LIB_DIR)/libfmm +LIB_EXT = .so -ifeq ($(OS),Windows_NT) -define MKDIR_FUNC - if not exist $(1) mkdir $(1) -endef -else -define MKDIR_FUNC - mkdir -p $(1) -endef -endif +CC := gcc +FOMPFLAGS := -fopenmp + +# link static library on Windows +LINK_STATIC := +LDFLAGS := $(FOMPFLAGS) -ifeq ($(OS),Windows_NT) - FILESEP = \\ +# expand stack memory for Windows +STACK_MEM := -else - FILESEP = / +ifeq ($(OS),Windows_NT) # link static oenpmp on windows + STACK_MEM := -Wl,-stack,0x1000000 + LINK_STATIC := -static + LDFLAGS := $(LINK_STATIC) $(FOMPFLAGS) endif -BIN_DIR = bin -SRC_DIR = src -INC_DIR = include -BUILD_DIR = build -BUILD_DIR_FLOAT = $(BUILD_DIR)$(FILESEP)float -BUILD_DIR_DOUBLE = $(BUILD_DIR)$(FILESEP)double -LIB_DIR = lib -LIB_NAME = $(LIB_DIR)$(FILESEP)libfmm -LIB_EXT = .so +# change architecture for macOS, from make command +ARCH = -CC = gcc -FOMPFLAGS = -fopenmp -CFLAGS = $(FOMPFLAGS) -O3 -g -ffast-math -march=native -mtune=native -fPIC -Wimplicit-function-declaration -I$(INC_DIR) -lm # -fsanitize=address -lasan -LDFLAGS = $(FOMPFLAGS) -shared -lm # -lasan +# -ffast-math -march=native -mtune=native # 如果加上这些选项，数学库-lm需要在编译动态库时再次显式指定。总是gcc的编译参数的前后顺序很讲究 +CFLAGS = $(LINK_STATIC) -lm -O3 -g -fPIC \ + -Wall $(STACK_MEM) -I$(shell realpath $(INC_DIR)) $(ARCH) $(FOMPFLAGS) # -fsanitize=address -lasan SRCS = $(wildcard $(SRC_DIR)/*.c) INCS = $(wildcard $(INC_DIR)/*.h) # 不同版本的目标文件目录 OBJS_FLOAT = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR_FLOAT)/%.o, $(SRCS)) -OBJS_FLOAT := $(filter-out $(BUILD_DIR_FLOAT)/grt_fmm.o, $(OBJS_FLOAT)) OBJS_DOUBLE = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR_DOUBLE)/%.o, $(SRCS)) -OBJS_DOUBLE := $(filter-out $(BUILD_DIR_DOUBLE)/grt_fmm.o, $(OBJS_DOUBLE)) - +DEPS_FLOAT := $(OBJS_FLOAT:.o=.d) +DEPS_DOUBLE := $(OBJS_DOUBLE:.o=.d) # 生成的库名称 TARGET_FLOAT = $(LIB_NAME)_float$(LIB_EXT) TARGET_DOUBLE = $(LIB_NAME)_double$(LIB_EXT) -.PHONY: all clean +.PHONY: all clean cleanbuild -all: $(TARGET_FLOAT) $(TARGET_DOUBLE) $(BIN_DIR)$(FILESEP)grt.fmm +all: $(BUILD_DIR) $(LIB_DIR) $(TARGET_FLOAT) $(TARGET_DOUBLE) -$(BIN_DIR)$(FILESEP)grt.fmm: $(OBJS_DOUBLE) $(BUILD_DIR_DOUBLE)/grt_fmm.o - @$(call MKDIR_FUNC, $(BIN_DIR)) - $(CC) -o $@ $^ $(CFLAGS) $(FOMPFLAGS) +$(BUILD_DIR): + @mkdir -p $@ -# 链接动态库，生成 float 版本 -$(TARGET_FLOAT): $(OBJS_FLOAT) - @$(call MKDIR_FUNC, $(LIB_DIR)) - $(CC) -o $@ $^ $(LDFLAGS) +$(LIB_DIR): + @mkdir -p $@ -# 链接动态库，生成 double 版本 -$(TARGET_DOUBLE): $(OBJS_DOUBLE) - @$(call MKDIR_FUNC, $(LIB_DIR)) - $(CC) -o $@ $^ $(LDFLAGS) +# ----------------------- Dependency generation ----------------------- +-include $(DEPS_FLOAT) +-include $(DEPS_DOUBLE) + +$(BUILD_DIR_FLOAT)/%.d: $(SRC_DIR)/%.c + @mkdir -p $(shell dirname $@) + @$(CC) $(CFLAGS) -MM $< > $@.$$$$; \ + sed 's,\($*\)\.o[ :]*,$(BUILD_DIR_FLOAT)/\1.o $@ : ,g' < $@.$$$$ > $@; \ + rm -f $@.$$$$ + +$(BUILD_DIR_DOUBLE)/%.d: $(SRC_DIR)/%.c + @mkdir -p $(shell dirname $@) + @$(CC) $(CFLAGS) -MM $< > $@.$$$$; \ + sed 's,\($*\)\.o[ :]*,$(BUILD_DIR_DOUBLE)/\1.o $@ : ,g' < $@.$$$$ > $@; \ + rm -f $@.$$$$ # 编译 float 版本的目标文件 -$(BUILD_DIR_FLOAT)/%.o: $(SRC_DIR)/%.c $(INCS) - @$(call MKDIR_FUNC, $(BUILD_DIR_FLOAT)) +$(BUILD_DIR_FLOAT)/%.o: $(SRC_DIR)/%.c $(CC) -o $@ -c $< $(CFLAGS) -DUSE_FLOAT # 编译 double 版本的目标文件 -$(BUILD_DIR_DOUBLE)/%.o: $(SRC_DIR)/%.c $(INCS) - @$(call MKDIR_FUNC, $(BUILD_DIR_DOUBLE)) +$(BUILD_DIR_DOUBLE)/%.o: $(SRC_DIR)/%.c $(CC) -o $@ -c $< $(CFLAGS) +# 链接动态库，生成 float 版本 +$(TARGET_FLOAT): $(OBJS_FLOAT) + $(CC) -shared -o $@ $^ $(LDFLAGS) + +# 链接动态库，生成 double 版本 +$(TARGET_DOUBLE): $(OBJS_DOUBLE) + $(CC) -shared -o $@ $^ $(LDFLAGS) + + +cleanbuild: + rm -rf $(BUILD_DIR) + clean: - $(call RMDIR_FUNC, $(BUILD_DIR)) - $(call RMDIR_FUNC, $(LIB_DIR)) + rm -rf $(BUILD_DIR) + rm -rf $(LIB_DIR) diff --git a/pyfmm/C_extension/include/sacio.h b/pyfmm/C_extension/include/sacio.h deleted file mode 100644 index afd3283..0000000 --- a/pyfmm/C_extension/include/sacio.h +++ /dev/null @@ -1,382 +0,0 @@ -/******************************************************************************* - Name: sacio.h - - Purpose: structure for header of a SAC (Seismic Analysis Code) - data file, and prototype for basic SAC I/O - - Notes: - Key to comment flags describing each field: - Column 1: - R required by SAC - (blank) optional - Column 2: - A = settable from a priori knowledge - D = available in data - F = available in or derivable from SEED fixed data header - T = available in SEED header tables - (blank) = not directly available from SEED data, header - tables, or elsewhere - - Problems: none known - - References: - * O'Neill, D. (1987). IRIS Interim Data Distribution Format - (SAC ASCII), Version 1.0 (12 November 1987). Incorporated - Research Institutions for Seismology, 1616 North Fort Myer - Drive, Suite 1440, Arlington, Virginia 22209. 11 pp. - * Tull, J. (1987). SAC User's Manual, Version 10.2, October 7, - 1987. Lawrence Livermore National Laboratory, L-205, - Livermore, California 94550. ??? pp. - - Language: C, hopefully ANSI standard - - Author: Dennis O'Neill - - Revisions: - 07/15/88 Dennis O'Neill Initial preliminary release 0.9 - 11/21/88 Dennis O'Neill Production release 1.0 - 01/27/91 Lorraine Hwang Header number is now version 6 - 07/06/93 Xiaoming Ding structure name sac -> sac_head - typedef structure to be SACHEAD - 12/06/96 Lupei Zhu prototype sacio functions - 05/04/13 Dongdong Tian modify headers according to SAC v101.5 -*******************************************************************************/ - -#ifndef SACIO_H -#define SACIO_H - -/******************************************************************************* - SAC header structure - - The SAC package is originally implemented in FORTRAN language, and the SAC - file format follows the conventions of FORTRAN language. Each character - string has a length of 8 bytes (16 bytes for kevnm). - - While reimplementing SAC in C, an extra byte is needed since C use '\0' to - mark the termination of a string to avoid wiping out the contents of the - last character. - - The header structure are needed, one for structure in memory and one for - reading/writing header structure from/to disk. - - While reading a file, first read in the numeric part to the header - structure, then read the string part to a temporary buffer, then map - strings from buffer to header structure. - - While writing a file, first write the numeric part from header structure - to disk, then map the strings from header structure to temporary buffer, - then write the buffer to disk. - -*******************************************************************************/ - -typedef struct sac_head { - float delta; /* RF increment between evenly spaced samples */ - float depmin; /* minimum value of dependent variable */ - float depmax; /* maximum value of dependent variable */ - float scale; /* amplitude scale factor (not used) */ - float odelta; /* Observed increment */ - float b; /* RD begining value of the independent variable */ - float e; /* RD ending value of the independent variable */ - float o; /* event origin time(seconds wrt referece time)*/ - float a; /* 1st arrival time (seconds wrt referece time)*/ - float internal1; /* internal use */ - float t0; /* user-defined time pick */ - float t1; /* user-defined time pick */ - float t2; /* user-defined time pick */ - float t3; /* user-defined time pick */ - float t4; /* user-defined time pick */ - float t5; /* user-defined time pick */ - float t6; /* user-defined time pick */ - float t7; /* user-defined time pick */ - float t8; /* user-defined time pick */ - float t9; /* user-defined time pick */ - float f; /* end time of event, sec > 0 */ - float resp0; /* instrument respnse parameter (not used) */ - float resp1; /* instrument respnse parameter (not used) */ - float resp2; /* instrument respnse parameter (not used) */ - float resp3; /* instrument respnse parameter (not used) */ - float resp4; /* instrument respnse parameter (not used) */ - float resp5; /* instrument respnse parameter (not used) */ - float resp6; /* instrument respnse parameter (not used) */ - float resp7; /* instrument respnse parameter (not used) */ - float resp8; /* instrument respnse parameter (not used) */ - float resp9; /* instrument respnse parameter (not used) */ - float stla; /* T station latititude (degree, north positive) */ - float stlo; /* T station longitude (degree, east positive) */ - float stel; /* T station elevation (meters, not used) */ - float stdp; /* T station depth (meters, not used) */ - float evla; /* event latitude (degree, north positive) */ - float evlo; /* event longitude (degree, east positive) */ - float evel; /* event elevation (meters, not used) */ - float evdp; /* event depth (kilometer, previously meters) */ - float mag; /* event magnitude */ - float user0; /* User defined variable storage area */ - float user1; /* User defined variable storage area */ - float user2; /* User defined variable storage area */ - float user3; /* User defined variable storage area */ - float user4; /* User defined variable storage area */ - float user5; /* User defined variable storage area */ - float user6; /* User defined variable storage area */ - float user7; /* User defined variable storage area */ - float user8; /* User defined variable storage area */ - float user9; /* User defined variable storage area */ - float dist; /* station-event distance (km) */ - float az; /* event-station azimuth */ - float baz; /* station-event azimuth */ - float gcarc; /* station-event great arc length (degrees) */ - float internal2; /* internal use */ - float internal3; /* internal use */ - float depmen; /* mean value of dependent variable */ - float cmpaz; /* T component azimuth (degree CW from north) */ - float cmpinc; /* T component inclination (degree from vertical)*/ - float xminimum; /* minimum value of X (spectral files only) */ - float xmaximum; /* maximum value of X (spectral files only) */ - float yminimum; /* minimum value of Y (spectral files only) */ - float ymaximun; /* maximum value of Y (spectral files only) */ - float unused1; /* reserved for future use */ - float unused2; /* reserved for future use */ - float unused3; /* reserved for future use */ - float unused4; /* reserved for future use */ - float unused5; /* reserved for future use */ - float unused6; /* reserved for future use */ - float unused7; /* reserved for future use */ - int nzyear; /* F GMT year corresponding to zero time of file */ - int nzjday; /* F GMT julia day */ - int nzhour; /* F GMT hour */ - int nzmin; /* F GMT minite */ - int nzsec; /* F GMT second */ - int nzmsec; /* F GMT millisecond */ - int nvhdr; /* R header version number (6) */ - int norid; /* origin ID (CSS 3.0) */ - int nevid; /* event ID (CSS 3.0) */ - int npts; /* RF number of points per data component */ - int internal4; /* internal use */ - int nwfid; /* waveform ID (CSS 3.0) */ - int nxsize; /* spectral length (spectral files only) */ - int nysize; /* spectral width (spectral files only) */ - int unused8; /* reserved for future use */ - int iftype; /* RA type of file */ - int idep; /* type of dependent variable */ - int iztype; /* reference time equivalence */ - int unused9; /* reserved for future use */ - int iinst; /* type of recording instrument (not used) */ - int istreg; /* station geographic region (not used) */ - int ievreg; /* event geographic region (not used) */ - int ievtyp; /* type of event */ - int iqual; /* quality of data (not used) */ - int isynth; /* synthetic data flag (not used) */ - int imagtyp; /* magnitude type */ - int imagsrc; /* source of magnitude information */ - int unused10; /* reserved for future use */ - int unused11; /* reserved for future use */ - int unused12; /* reserved for future use */ - int unused13; /* reserved for future use */ - int unused14; /* reserved for future use */ - int unused15; /* reserved for future use */ - int unused16; /* reserved for future use */ - int unused17; /* reserved for future use */ - int leven; /* RA TRUE if data is evenly spaced */ - int lpspol; /* station polarity flag (left hand rule) */ - int lovrok; /* overwrite permission */ - int lcalda; /* true if to calculate distance, azimuth */ - int unused18; /* reserved for future use */ - char kstnm[9]; /* F station name */ - char kevnm[18]; /* event name */ - char khole[9]; /* nuclear: hole id; Other: location id; */ - char ko[9]; /* event origin time id */ - char ka[9]; /* 1st arrival time id */ - char kt0[9]; /* time pick 0 id */ - char kt1[9]; /* time pick 1 id */ - char kt2[9]; /* time pick 2 id */ - char kt3[9]; /* time pick 3 id */ - char kt4[9]; /* time pick 4 id */ - char kt5[9]; /* time pick 5 id */ - char kt6[9]; /* time pick 6 id */ - char kt7[9]; /* time pick 7 id */ - char kt8[9]; /* time pick 8 id */ - char kt9[9]; /* time pick 9 id */ - char kf[9]; /* end of event id */ - char kuser0[9]; /* User defined variable storage area */ - char kuser1[9]; /* User defined variable storage area */ - char kuser2[9]; /* User defined variable storage area */ - char kcmpnm[9]; /* F channel name, three charaters */ - char knetwk[9]; /* name of seismic network */ - char kdatrd[9]; /* date data was read onto computer */ - char kinst[9]; /* generic name of recording instrument */ -} SACHEAD; - -/******************************************************************************* - SAC Enumerated type - - Definitions of constants for SAC enumerated data values. - - undocumented == couldn't find a definition for it (denio, 07/15/88) -*******************************************************************************/ -#define ITIME 1 /* file: time series data */ -#define IRLIM 2 /* file: real&imag spectrum */ -#define IAMPH 3 /* file: ampl&phas spectrum */ -#define IXY 4 /* file: gen'l x vs y data */ -#define IUNKN 5 /* x data: unknown type */ - /* zero time: unknown */ - /* event type: unknown */ -#define IDISP 6 /* x data: displacement (nm) */ -#define IVEL 7 /* x data: velocity (nm/sec) */ -#define IACC 8 /* x data: accel (cm/sec/sec)*/ -#define IB 9 /* zero time: start of file */ -#define IDAY 10 /* zero time: 0000 of GMT day*/ -#define IO 11 /* zero time: event origin */ -#define IA 12 /* zero time: 1st arrival */ -#define IT0 13 /* zero time: user timepick 0*/ -#define IT1 14 /* zero time: user timepick 1*/ -#define IT2 15 /* zero time: user timepick 2*/ -#define IT3 16 /* zero time: user timepick 3*/ -#define IT4 17 /* zero time: user timepick 4*/ -#define IT5 18 /* zero time: user timepick 5*/ -#define IT6 19 /* zero time: user timepick 6*/ -#define IT7 20 /* zero time: user timepick 7*/ -#define IT8 21 /* zero time: user timepick 8*/ -#define IT9 22 /* zero time: user timepick 9*/ -#define IRADNV 23 /* undocumented */ -#define ITANNV 24 /* undocumented */ -#define IRADEV 25 /* undocumented */ -#define ITANEV 26 /* undocumented */ -#define INORTH 27 /* undocumented */ -#define IEAST 28 /* undocumented */ -#define IHORZA 29 /* undocumented */ -#define IDOWN 30 /* undocumented */ -#define IUP 31 /* undocumented */ -#define ILLLBB 32 /* undocumented */ -#define IWWSN1 33 /* undocumented */ -#define IWWSN2 34 /* undocumented */ -#define IHGLP 35 /* undocumented */ -#define ISRO 36 /* undocumented */ -#define INUCL 37 /* event type: nuclear shot */ -#define IPREN 38 /* event type: nuke pre-shot */ -#define IPOSTN 39 /* event type: nuke post-shot*/ -#define IQUAKE 40 /* event type: earthquake */ -#define IPREQ 41 /* event type: foreshock */ -#define IPOSTQ 42 /* event type: aftershock */ -#define ICHEM 43 /* event type: chemical expl */ -#define IOTHER 44 /* event type: other source */ - /* data quality: other problm*/ -#define IGOOD 45 /* data quality: good */ -#define IGLCH 46 /* data quality: has glitches*/ -#define IDROP 47 /* data quality: has dropouts*/ -#define ILOWSN 48 /* data quality: low s/n */ -#define IRLDTA 49 /* data is real data */ -#define IVOLTS 50 /* file: velocity (volts) */ -#define IMB 52 /* undocumented */ -#define IMS 53 /* undocumented */ -#define IML 54 /* undocumented */ -#define IMW 55 /* undocumented */ -#define IMD 56 /* undocumented */ -#define IMX 57 /* undocumented */ -#define INEIC 58 /* undocumented */ -#define IPDEQ 59 /* undocumented */ -#define IPDEW 60 /* undocumented */ -#define IPDE 61 /* undocumented */ -#define IISC 62 /* undocumented */ -#define IREB 63 /* undocumented */ -#define IUSGS 64 /* undocumented */ -#define IBRK 65 /* undocumented */ -#define ICALTECH 66 /* undocumented */ -#define ILLNL 67 /* undocumented */ -#define IEVLOC 68 /* undocumented */ -#define IJSOP 69 /* undocumented */ -#define IUSER 70 /* undocumented */ -#define IUNKNOWN 71 /* undocumented */ -#define IQB 72 /* undocumented */ -#define IQB1 73 /* undocumented */ -#define IQB2 74 /* undocumented */ -#define IQBX 75 /* undocumented */ -#define IQMT 76 /* undocumented */ -#define IEQ 77 /* undocumented */ -#define IEQ1 78 /* undocumented */ -#define IEQ2 79 /* undocumented */ -#define IME 80 /* undocumented */ -#define IEX 81 /* undocumented */ -#define INU 82 /* undocumented */ -#define INC 83 /* undocumented */ -#define IO_ 84 /* undocumented */ -#define IL 85 /* undocumented */ -#define IR 86 /* undocumented */ -#define IT 87 /* undocumented */ -#define IU 88 /* undocumented */ -#define IEQ3 89 /* undocumented */ -#define IEQ0 90 /* undocumented */ -#define IEX0 91 /* undocumented */ -#define IQC 92 /* undocumented */ -#define IQB0 93 /* undocumented */ -#define IGEY 94 /* undocumented */ -#define ILIT 95 /* undocumented */ -#define IMET 96 /* undocumented */ -#define IODOR 97 /* undocumented */ -#define IOS 103 /* undocumented */ - -/* True/false definitions */ -#undef FALSE -#define FALSE 0 -#undef TRUE -#define TRUE 1 - -#define SAC_FLOAT_UNDEF (-12345.0) -#define SAC_INT_UNDEF (-12345) -#define SAC_CHAR8_UNDEF "-12345 " -#define SAC_CHAR16_UNDEF "-12345 " - -/* Format strings for writing headers for SAC ASCII files */ -#define FCS "%15.7f%15.7f%15.7f%15.7f%15.7f\n" /* for floats */ -#define ICS "%10d%10d%10d%10d%10d\n" /* for integers */ -#define CCS1 "%-8.8s%-8.8s%-8.8s\n" /* for strings */ -#define CCS2 "%-8.8s%-16.16s\n" /* for strings */ - -/* Number of floats in the SAC Header */ -#define SAC_HEADER_FLOATS 70 -/* Number of ints in the SAC Header */ -#define SAC_HEADER_INTS 40 /* 15 + 20 + 5 */ -/* Number of numeric values in the SAC Header (4 bytes) */ -#define SAC_HEADER_NUMBERS ( SAC_HEADER_FLOATS + SAC_HEADER_INTS ) -/* Number of strings in the SAC Header (8 or 9 bytes) */ -#define SAC_HEADER_STRINGS 24 /* 24 = 23 + 1 */ - -/* size of integers or floats on disk or in memory - * make sure sizeof(float) == 4 && sizeof(int)==4 - */ -#define SAC_DATA_SIZEOF 4 -/* Size of a character string stored on disk for a SAC header */ -#define SAC_HEADER_STRING_LENGTH_FILE 8 -/* Size of a character string stored in memory for a SAC header */ -#define SAC_HEADER_STRING_LENGTH 9 - -/* Size of floats in the SAC Header */ -#define SAC_HEADER_FLOATS_SIZE ( SAC_HEADER_FLOATS * SAC_DATA_SIZEOF ) -/* Size of ints in the SAC Header */ -#define SAC_HEADER_INTS_SIZE ( SAC_HEADER_INTS * SAC_DATA_SIZEOF ) -/* Size of numeric headers on disk */ -#define SAC_HEADER_NUMBERS_SIZE ( SAC_HEADER_FLOATS_SIZE + SAC_HEADER_INTS_SIZE ) -/* Size of string headers on disk */ -#define SAC_HEADER_STRINGS_SIZE ( SAC_HEADER_STRINGS * SAC_HEADER_STRING_LENGTH_FILE ) - -/* SAC Header Version Number */ -#define SAC_HEADER_MAJOR_VERSION 6 -/* offset of nvhdr relative to struct SACHEAD */ -#define SAC_VERSION_LOCATION 76 - -/* offset of T0 relative to pointer to struct SACHEAD */ -#define TMARK 10 - -/* offset of USER0 relative to pointer to struct SACHEAD */ -#define USERN 40 - -/* function prototype of basic SAC I/O */ -int read_sac_head(const char *name, SACHEAD *hd); -float *read_sac(const char *name, SACHEAD *hd); -int read_sac_xy(const char *name, SACHEAD *hd, float *xdata, float *ydata); -float *read_sac_pdw(const char *name, SACHEAD *hd, int tmark, float t1, float t2); -int write_sac(const char *name, SACHEAD hd, const float *ar); -int write_sac_xy(const char *name, SACHEAD hd, const float *xdata, const float *ydata); -SACHEAD new_sac_head(float dt, int ns, float b0); -int sac_head_index(const char *name); -int issac(const char *name); - -#endif /* sacio.h */ \ No newline at end of file diff --git a/pyfmm/C_extension/src/grt_fmm.c b/pyfmm/C_extension/src/grt_fmm.c deleted file mode 100644 index bdb16f5..0000000 --- a/pyfmm/C_extension/src/grt_fmm.c +++ /dev/null @@ -1,626 +0,0 @@ -/** - * @file grt_fmm.c - * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) - * @date 2024-11 - * - * 提供给GRT程序计算的辅助工具，计算对应全局初至波 - * - * GRT是计算1D层状半空间内的理论地震图的工具，详见 - * https://github.com/Dengda98/PyGRT - * -*/ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "const.h" -#include "interp.h" -#include "fmm.h" -#include "sacio.h" - - -extern char *optarg; -extern int optind; -extern int optopt; - -//****************** 在该文件以内的全局变量 ***********************// -// 命令名称 -static char *command = NULL; -// 模型路径，模型PYMODEL1D指针，全局最大最小速度 -static char *s_modelpath = NULL; -static char *s_modelname = NULL; -static MYINT nlay=0; -static double *Dep1d = NULL; -static double *Vp1d = NULL; -static double *Vs1d = NULL; -static double minh = 9e30; // 模型最细厚度 -// 三个维度的数组以及大小 -static double *xs=NULL, *zs=NULL; -static MYINT nx, nz, nxz; -static const double ys[1]={0.0}; -static const MYINT ny = 1; -// 震源坐标 -static double xsrc=0.0, zsrc=0.0; -// 场点坐标 -static double xrcv=0.0, zrcv=0.0; -// 慢度数组 -static MYREAL *SlwP = NULL, *SlwS = NULL; -// 走时数组 -static MYREAL *TTP = NULL, *TTS = NULL; -static MYINT nTT=0; -// grt命令的输出路径, 以及符合要求的待处理文件夹路径 -// 以及根据文件转化的震源深度、场点深度、震中距 -static char *s_grtoutdir = NULL; -static MYINT ndir=0; -static char **s_outsubdirs=NULL; -static double *depsrcs=NULL, *deprcvs=NULL, *epidists=NULL; -static double depmax=-1.0, epidistmax=-1.0; -// 在outsubdir中，要根据文件名中的非重复震源深度分配要计算多少个走时场， -// idx_nodup[i]反映的是depsrcs[i]第一次出现在depsrcs中的索引，长度为ndir -static MYINT *idx_nodup=NULL; -// 每个文件夹是否已计算对应的走时场，长度ndir -static bool *subdirs_hasTT=NULL; - -// 各选项的标志变量，初始化为0，定义了则为1 -static MYINT X_flag=0, Z_flag=0, M_flag=0, O_flag=0; - - -static void print_help(){ - -} - -/** - * 从路径字符串中找到用/或\\分隔的最后一项 - * - * @param path 路径字符串指针 - * - * @return 指向最后一项字符串的指针 - */ -static char* get_basename(char* path) { - // 找到最后一个 '/' - char* last_slash = strrchr(path, '/'); - -#ifdef _WIN32 - char* last_backslash = strrchr(path, '\\'); - if (last_backslash && (!last_slash || last_backslash > last_slash)) { - last_slash = last_backslash; - } -#endif - if (last_slash) { - // 返回最后一个 '/' 之后的部分 - return last_slash + 1; - } - // 如果没有 '/'，整个路径就是最后一项 - return path; -} - -/** - * 找出double数组中某数第一次出现的索引 - * - */ -static MYINT find_index(double *arr, MYINT n, double target){ - for(MYINT i=0; id_name, ".") == 0 || strcmp(entry->d_name, "..") == 0) { - continue; - } - - // 构建文件夹的完整路径 - snprintf(folder_path, sizeof(folder_path), "%s/%s", indir, entry->d_name); - - // 判断是否是文件夹 - if (opendir(folder_path) != NULL) { - // 检查文件夹名称是否符合模式 - double f1, f2, f3; - if (match_folder_name(entry->d_name, &f1, &f2, &f3)) { - printf("Find Matching folder: %s\n", entry->d_name); - // 写入数组 - s_outsubdirs = (char**)realloc(s_outsubdirs, sizeof(char*)*(ndir+1)); - s_outsubdirs[ndir] = NULL; - s_outsubdirs[ndir] = (char*)realloc(s_outsubdirs[ndir], sizeof(char)*lenmax); - strcpy(s_outsubdirs[ndir], folder_path); - - depsrcs = (double*)realloc(depsrcs, sizeof(double)*(ndir+1)); - depsrcs[ndir] = f1; - deprcvs = (double*)realloc(deprcvs, sizeof(double)*(ndir+1)); - deprcvs[ndir] = f2; - epidists = (double*)realloc(epidists, sizeof(double)*(ndir+1)); - epidists[ndir] = f3; - idx_nodup = (MYINT*)realloc(idx_nodup, sizeof(MYINT)*(ndir+1)); - idx_nodup[ndir] = find_index(depsrcs, ndir+1, f1); - subdirs_hasTT = (bool*)realloc(subdirs_hasTT, sizeof(bool)*(nlay+1)); - subdirs_hasTT[ndir] = false; - - if(depmax < f1) depmax=f1; - if(depmax < f2) depmax=f2; - if(epidistmax < f3) epidistmax = f3; - - ndir++; - - } - } - } - - closedir(dir); -} - - -/** - * 从1D模型文件中读取Vp、Vs、Depth - */ -static MYINT load_model1D(){ - FILE *fp; - if((fp = fopen(s_modelpath, "r")) == NULL){ - fprintf(stderr, "[%s] Model file open error.\n", command); - return -1; - } - - char line[1024]; - MYINT iline = 0; - double h, va, vb, rho, qa, qb; - h = va = vb = rho = qa = qb = -9.0; - nlay = 0; - double dep=0.0; - while(fgets(line, sizeof(line), fp)) { - iline++; - if(h > 0.0 && h < minh) minh = h; - - h = va = vb = rho = qa = qb = -9.0; - if(6 != sscanf(line, "%lf %lf %lf %lf %lf %lf\n", &h, &va, &vb, &rho, &qa, &qb)){ - fprintf(stderr, "[%s] Model file read error in line %d.\n", command, iline); - return -1; - }; - - if(h < 0.0 || va < 0.0 || vb < 0.0 || rho < 0.0 || qa < 0.0 || qb < 0.0){ - fprintf(stderr, "[%s] In line %d, negative value is not supported.\n", command, iline); - return -1; - } - - - Vp1d = (double*)realloc(Vp1d, sizeof(double)*(nlay+1)); - Vs1d = (double*)realloc(Vs1d, sizeof(double)*(nlay+1)); - Dep1d = (double*)realloc(Dep1d, sizeof(double)*(nlay+1)); - Vp1d[nlay] = va; - Vs1d[nlay] = vb; - Dep1d[nlay] = dep; - - nlay++; - dep += h; - } - - if(iline==0 || Vp1d==NULL){ - fprintf(stderr, "[%s] Model file read error.\n", command); - return -1; - } - - - - return 0; -} - -/** - * 根据设置的1D速度离散得到差分网格对应的速度 - * - * @param Vel1d 1D速度 - * @param Slw 2D差分网格慢度 - */ -static void interp1d_vel(const double *Vel1d, MYREAL *Slw){ - for(MYINT iz=0; izd_name, ".") == 0 || strcmp(entry->d_name, "..") == 0) { - continue; - } - - sprintf(fullpath, "%s/%s", subdir, entry->d_name); - - // 读sac文件 - float *arr = read_sac(fullpath, &hd); - // 写入走时 - hd.t0 = hd.o + travtP; - strcpy(hd.kt0, "P"); - hd.t1 = hd.o + travtS; - strcpy(hd.kt1, "S"); - // 写回sac - write_sac(fullpath, hd, arr); - - free(arr); - } - -} - - - -//======================================================================== -//======================================================================== -//======================================================================== -MYINT main(MYINT argc, char **argv){ - command = argv[0]; - getopt_from_command(argc, argv); - - // 提取符合条件的文件夹 - list_matching_folders(s_grtoutdir); - - // 读入一维模型 - if(load_model1D() < 0){ - exit(EXIT_FAILURE); - } - - // 自动设置Z网格 - if(Z_flag == 0){ - double ze = depmax, dz = minh/5.0; - - } - - // 自动设置X网格 - if(X_flag == 0){ - double xe = epidistmax; - } - - - // 计算维度 - nxz = nx*nz; - SlwP = (MYREAL*)malloc(sizeof(MYREAL)*nxz); - SlwS = (MYREAL*)malloc(sizeof(MYREAL)*nxz); - - // 插值速度场 - interp1d_vel(Vp1d, SlwP); - interp1d_vel(Vs1d, SlwS); - - // for(MYINT ix=0; ix -#include -#include -#include -#include -#include "sacio.h" - -/* function prototype for local use */ -static void byte_swap (char *pt, size_t n); -static int check_sac_nvhdr (const int nvhdr); -static void map_chdr_in (char *memar, char *buff); -static int read_head_in (const char *name, SACHEAD *hd, FILE *strm); -static void map_chdr_out (char *memar, char *buff); -static int write_head_out (const char *name, SACHEAD hd, FILE *strm); - -/* a SAC structure containing all null values */ -static SACHEAD sac_null = { - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345., -12345., -12345., -12345., -12345., - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - -12345 , -12345 , -12345 , -12345 , -12345 , - { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' }, - { '-','1','2','3','4','5',' ',' ' } -}; - -/* - * read_sac_head - * - * Description: Read binary SAC header from file. - * - * IN: - * const char *name : File name - * OUT: - * SACHEAD *hd : SAC header - * - * Return: 0 if success, -1 if failed - * - */ -int read_sac_head(const char *name, SACHEAD *hd) -{ - FILE *strm; - int lswap; - - if ((strm = fopen(name, "rb")) == NULL) { - fprintf(stderr, "Unable to open %s\n", name); - return -1; - } - - lswap = read_head_in(name, hd, strm); - - fclose(strm); - - return ((lswap == -1) ? -1 : 0); -} - -/* - * read_sac - * - * Description: Read binary SAC data from file. - * - * IN: - * const char *name : file name - * OUT: - * SACHEAD *hd : SAC header to be filled - * Return: float pointer to the data array, NULL if failed. - * - */ -float *read_sac(const char *name, SACHEAD *hd) -{ - FILE *strm; - float *ar; - int lswap; - size_t sz; - - if ((strm = fopen(name, "rb")) == NULL) { - fprintf(stderr, "Unable to open %s\n", name); - return NULL; - } - - lswap = read_head_in(name, hd, strm); - - if (lswap == -1) { - fclose(strm); - return NULL; - } - - sz = (size_t) hd->npts * SAC_DATA_SIZEOF; - if (hd->iftype == IXY) sz *= 2; - - if ((ar = (float *)malloc(sz)) == NULL) { - fprintf(stderr, "Error in allocating memory for reading %s\n", name); - fclose(strm); - return NULL; - } - - if (fread((char*)ar, sz, 1, strm) != 1) { - fprintf(stderr, "Error in reading SAC data %s\n", name); - free(ar); - fclose(strm); - return NULL; - } - fclose(strm); - - if (lswap == TRUE) byte_swap((char*)ar, sz); - - return ar; -} - -/* - * read_sac_xy - * - * Description: read SAC XY binary file - * - * IN: - * const char *name : file name - * OUT: - * SACHEAD *hd : SAC head to be filled - * float *xdata : pointer for X - * float *ydata : pointer for Y - * - * Return: 0 for success, -1 for fail - * - */ -int read_sac_xy(const char *name, SACHEAD *hd, float *xdata, float *ydata) -{ - float *data; - size_t npts; - - if ((data = read_sac(name, hd)) == NULL) return -1; - - npts = (size_t)hd->npts; - if ((xdata = (float *)malloc(npts*SAC_DATA_SIZEOF)) == NULL) { - fprintf(stderr, "Error in allocating memory for %s\n", name); - free(data); - return -1; - } - if ((ydata = (float *)malloc(npts*SAC_DATA_SIZEOF)) == NULL) { - fprintf(stderr, "Error in allocating memory for %s\n", name); - free(data); - free(xdata); - return -1; - } - - memcpy(xdata, data , npts*SAC_DATA_SIZEOF); - memcpy(ydata, data+npts, npts*SAC_DATA_SIZEOF); - - free(data); - return 0; -} - -/* - * write_sac - * - * Description: write binary SAC data - * - * IN: - * const char *name : file name - * SACHEAD hd : header - * const float *ar : float data array - * - * Return: - * -1 : fail - * 0 : succeed - * - */ -int write_sac(const char *name, SACHEAD hd, const float *ar) -{ - FILE *strm; - size_t sz; - - if ((strm = fopen(name, "wb")) == NULL) { - fprintf(stderr, "Error in opening file for writing %s\n", name); - return -1; - } - - if (write_head_out(name, hd, strm) == -1) { - fclose(strm); - return -1; - } - - sz = (size_t)hd.npts * SAC_DATA_SIZEOF; - if (hd.iftype == IXY) sz *= 2; - - if (fwrite(ar, sz, 1, strm) != 1) { - fprintf(stderr, "Error in writing SAC data for writing %s\n", name); - fclose(strm); - return -1; - } - fclose(strm); - return 0; -} - -/* - * write_sac_xy - * - * Description: write binary SAC XY data - * - * IN: - * const char *name : file name - * SACHEAD hd : header - * const float *xdata : float data array for X - * const float *ydata : float data array for Y - * - * Return: - * -1 : fail - * 0 : succeed - * - */ -int write_sac_xy(const char *name, SACHEAD hd, const float *xdata, const float *ydata) -{ - float *ar; - int npts; - int error; - size_t sz; - - npts = hd.npts; - sz = (size_t)npts * SAC_DATA_SIZEOF; - - if ((ar = (float *)malloc(sz*2)) == NULL) { - fprintf(stderr, "Error in allocating memory for file %s\n", name); - return -1; - } - memcpy(ar, xdata, sz); - memcpy(ar+npts, ydata, sz); - - /* needed for XY data */ - hd.iftype = IXY; - hd.leven = FALSE; - - error = write_sac(name, hd, ar); - - free(ar); - return error; -} - -/* - * read_sac_pdw - * - * Description: - * Read portion of data from file. - * - * Arguments: - * const char *name : file name - * SACHEAD *hd : SAC header to be filled - * int tmark : time mark in SAC header - * -5 -> b; - * -4 -> e; - * -3 -> o; - * -2 -> a; - * 0-9 -> Tn; - * others -> t=0; - * float t1 : begin time is tmark + t1 - * float t2 : end time is tmark + t2 - * - * Return: - * float pointer to the data array, NULL if failed. - * - */ -float *read_sac_pdw(const char *name, SACHEAD *hd, int tmark, float t1, float t2) -{ - FILE *strm; - int lswap; - float tref; - int nt1, nt2, npts, nn; - float *ar, *fpt; - - if ((strm = fopen(name, "rb")) == NULL) { - fprintf(stderr, "Error in opening %s\n", name); - return NULL; - } - - lswap = read_head_in(name, hd, strm); - - if (lswap == -1) { - fclose(strm); - return NULL; - } - - nn = (int)((t2-t1)/hd->delta); - if (nn<=0 || (ar = (float *)calloc((size_t)nn, SAC_DATA_SIZEOF)) == NULL) { - fprintf(stderr, "Errorin allocating memory for reading %s n=%d\n", name, nn); - fclose(strm); - return NULL; - } - - tref = 0.; - if (tmark>=-5 && tmark<=9 && tmark!=-1) { - tref = *((float *) hd + TMARK + tmark); - if (fabs(tref+12345.)<0.1) { - fprintf(stderr, "Time mark undefined in %s\n", name); - free(ar); - fclose(strm); - return NULL; - } - } - t1 += tref; - nt1 = (int)((t1 - hd->b) / hd->delta); - nt2 = nt1 + nn; - npts = hd->npts; - hd->npts = nn; - hd->b = t1; - hd->e = t1 + nn * hd->delta; - - if (nt1>npts || nt2 <0) return ar; /* return zero filled array */ - /* maybe warnings are needed! */ - - if (nt1<0) { - fpt = ar - nt1; - nt1 = 0; - } else { - if (fseek(strm, nt1*SAC_DATA_SIZEOF, SEEK_CUR) < 0) { - fprintf(stderr, "Error in seek %s\n", name); - free(ar); - fclose(strm); - return NULL; - } - fpt = ar; - } - if (nt2>npts) nt2 = npts; - nn = nt2 - nt1; - - if (fread((char *)fpt, (size_t)nn * SAC_DATA_SIZEOF, 1, strm) != 1) { - fprintf(stderr, "Error in reading SAC data %s\n", name); - free(ar); - fclose(strm); - return NULL; - } - fclose(strm); - - if (lswap == TRUE) byte_swap((char*)ar, (size_t)nn*SAC_DATA_SIZEOF); - - return ar; -} - -/* - * new_sac_head - * - * Description: create a new SAC header with required fields - * - * IN: - * float dt : sample interval - * int ns : number of points - * float b0 : starting time - */ -SACHEAD new_sac_head(float dt, int ns, float b0) -{ - SACHEAD hd = sac_null; - hd.delta = dt; - hd.npts = ns; - hd.b = b0; - hd.o = 0.; - hd.e = b0+(ns-1)*dt; - hd.iztype = IO; - hd.iftype = ITIME; - hd.leven = TRUE; - hd.nvhdr = SAC_HEADER_MAJOR_VERSION; - return hd; -} - -/* - * sac_head_index - * - * Description: return the index of a specified sac head field - * - * In: - * const char *name : name of sac head field - * Return: - * index of a specified field in sac head - * - */ -int sac_head_index(const char *name) -{ - const char fields[SAC_HEADER_NUMBERS+SAC_HEADER_STRINGS][10] = { - "delta", "depmin", "depmax", "scale", "odelta", - "b", "e", "o", "a", "internal1", - "t0", "t1", "t2", "t3", "t4", - "t5", "t6", "t7", "t8", "t9", - "f", "resp0", "resp1", "resp2", "resp3", - "resp4", "resp5", "resp6", "resp7", "resp8", - "resp9", "stla", "stlo", "stel", "stdp", - "evla", "evlo", "evel", "evdp", "mag", - "user0", "user1", "user2", "user3", "user4", - "user5", "user6", "user7", "user8", "user9", - "dist", "az", "baz", "gcarc", "internal2", - "internal3","depmen", "cmpaz", "cmpinc", "xminimum", - "xmaximum", "yminimum", "ymaximum", "unused1", "unused2", - "unused3", "unused4", "unused5", "unused6", "unused7", - "nzyear", "nzjday", "nzhour", "nzmin", "nzsec", - "nzmsec", "nvhdr", "norid", "nevid", "npts", - "internal4","nwfid", "nxsize", "nysize", "unused8", - "iftype", "idep", "iztype", "unused9", "iinst", - "istreg", "ievreg", "ievtyp", "iqual", "isynth", - "imagtyp", "imagsrc", "unused10", "unused11", "unused12", - "unused13", "unused14", "unused15", "unused16", "unused17", - "leven", "lpspol", "lovrok", "lcalda", "unused18", - "kstnm", "kevnm", "kevnmmore", - "khole", "ko", "ka", - "kt0", "kt1", "kt2", - "kt3", "kt4", "kt5", - "kt6", "kt7", "kt8", - "kt9", "kf", "kuser0", - "kuser1", "kuser2", "kcmpnm", - "knetwk", "kdatrd", "kinst", - }; - int i; - - for (i=0; invhdr); - if (lswap == -1) { - fprintf(stderr, "Warning: %s not in sac format.\n", name); - return -1; - } else if (lswap == TRUE) { - byte_swap((char *)hd, SAC_HEADER_NUMBERS_SIZE); - } - - /* read string parts of the SAC header */ - if ((buffer = (char *)malloc(SAC_HEADER_STRINGS_SIZE)) == NULL) { - fprintf(stderr, "Error in allocating memory %s\n", name); - return -1; - } - if (fread(buffer, SAC_HEADER_STRINGS_SIZE, 1, strm) != 1) { - fprintf(stderr, "Error in reading SAC header %s\n", name); - free(buffer); - return -1; - } - map_chdr_in((char *)(hd)+SAC_HEADER_NUMBERS_SIZE, buffer); - free(buffer); - - return lswap; -} - -/* - * map_chdr_out: - * map strings from memory to buffer - */ -static void map_chdr_out(char *memar, char *buff) -{ - char *ptr1; - char *ptr2; - int i; - - ptr1 = memar; - ptr2 = buff; - - memcpy(ptr2, ptr1, 8); - ptr1 += 9; - ptr2 += 8; - - memcpy(ptr2, ptr1, 16); - ptr1 += 18; - ptr2 += 16; - - for (i=0; i<21; i++) { - memcpy(ptr2, ptr1, 8); - ptr1 += 9; - ptr2 += 8; - } -} -/* - * - * write_head_out - * - * IN: - * const char *name : file name, only for debug - * SACHEAD hd : header to be written - * FILE *strm : file handler - * - * Return: - * -1 : failed. - * 0 : success. - */ -static int write_head_out(const char *name, SACHEAD hd, FILE *strm) -{ - char *buffer; - - if (sizeof(float) != SAC_DATA_SIZEOF || sizeof(int) != SAC_DATA_SIZEOF) { - fprintf(stderr, "Mismatch in size of basic data type!\n"); - return -1; - } - - if (fwrite(&hd, SAC_HEADER_NUMBERS_SIZE, 1, strm) != 1) { - fprintf(stderr, "Error in writing SAC data for writing %s\n", name); - return -1; - } - - if ((buffer = (char *)malloc(SAC_HEADER_STRINGS_SIZE)) == NULL){ - fprintf(stderr, "Error in allocating memory %s\n", name); - return -1; - } - map_chdr_out((char *)(&hd)+SAC_HEADER_NUMBERS_SIZE, buffer); - - if (fwrite(buffer, SAC_HEADER_STRINGS_SIZE, 1, strm) != 1) { - fprintf(stderr, "Error in writing SAC data for writing %s\n", name); - return -1; - } - free(buffer); - - return 0; -} \ No newline at end of file diff --git a/pyfmm/traveltime.py b/pyfmm/traveltime.py index 74dba3c..afdd97a 100644 --- a/pyfmm/traveltime.py +++ b/pyfmm/traveltime.py @@ -11,6 +11,7 @@ import numpy.ctypeslib as npct from ctypes import byref from scipy import interpolate +from typing import Union from . import c_interfaces from .logger import myLogger @@ -181,7 +182,7 @@ def travel_time_iniTT( def raytracing( TT:np.ndarray, srcloc:list, rcvloc:list, xarr:np.ndarray, yarr:np.ndarray, zarr:np.ndarray, - seglen:float, slw:np.ndarray|None=None, segfac:int=3, sphcoord:bool=False, maxdots:int=10000): + seglen:float, slw:Union[np.ndarray, None] = None, segfac:int=3, sphcoord:bool=False, maxdots:int=10000): r''' 根据给定源点坐标计算的走时场，使用梯度下降法做射线追踪 diff --git a/setup.py b/setup.py index 2d0eb1b..ed3af6e 100644 --- a/setup.py +++ b/setup.py @@ -1,34 +1,7 @@ from setuptools import setup, find_packages -from setuptools.command.build import build as build_orig -from setuptools.command.develop import develop as develop_orig -from setuptools.command.install import install as install_orig -import subprocess -import os, sys +import os -class BuildMake(build_orig): - def run(self): - # 执行 make 命令 - make_cmd = 'cd pyfmm/C_extension && make clean && make' - process = subprocess.Popen(make_cmd, stdout=sys.stdout, stderr=sys.stderr, shell=True) - process.wait() - if process.returncode != 0: - raise subprocess.CalledProcessError(process.returncode, 'make') - - super().run() - -# 强制install执行build -class Install(install_orig): - def run(self): - self.run_command('build') - super().run() - -# 强制develop执行build -class Develop(develop_orig): - def run(self): - self.run_command('build') - super().run() - # 读取版本号 def read_version(): @@ -58,16 +31,12 @@ def read_readme(): packages=find_packages(), package_data={'pyfmm': ['./C_extension/*']}, include_package_data=True, - cmdclass={ - 'build': BuildMake, - 'install': Install, - 'develop': Develop, # 添加 Develop 类到 cmdclass 中 - }, install_requires=[ 'numpy>=1.20, <2.0', 'scipy>=1.10', 'matplotlib>=3.5', 'jupyter', ], - python_requires='>=3.9', + python_requires='>=3.6', + zip_safe=False, # not compress the binary file )