Support static displacements, strain and stress calculations in Python

.github/workflows/testbuild.yml

@@ -236,19 +236,20 @@ jobs:
           # don't use $(pwd) to get abspath, not valid on Windows
           echo "EXAMPLE_COPY_PATH=${{ env.PACK_NAME }}/test_tmp" >> $GITHUB_ENV
 
-      - name: Test 1 compare_results
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_results
+      - name: Test compare_c_py
+        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_c_py
         run: |
+          # 检查C程序和Python程序的计算结果是否一致
           chmod +x *.sh
-          ./run_milrow_grt.sh
-          python plot_cps_pygrt.py
+          ./run_grt.sh
+          python run_pygrt.py
 
-      - name: Test 2 site_effect
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/site_effect
-        run: |
-          chmod +x *.sh 
-          ./run1.sh
-          python run2.py
+      # - name: Test 2 site_effect
+      #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/site_effect
+      #   run: |
+      #     chmod +x *.sh 
+      #     ./run1.sh
+      #     python run2.py
 
       # - name: Test 3 multi_traces
       #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/multi_traces
@@ -258,12 +259,12 @@ jobs:
       #     ./run1.sh
       #   continue-on-error: true  # 即使失败，仍然标记为成功
 
-      - name: Test 4 lamb_problem
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/lamb_problem
-        run: |
-          chmod +x *.sh
-          ./run1.sh
-          python run2.py
+      # - name: Test 4 lamb_problem
+      #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/lamb_problem
+      #   run: |
+      #     chmod +x *.sh
+      #     ./run1.sh
+      #     python run2.py
 
       # - name: Test 5 far_field
       #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/far_field
@@ -274,24 +275,24 @@ jobs:
       #     python plot_compare_pygrt.py
       #   continue-on-error: true  # 即使失败，仍然标记为成功
 
-      - name: Test signal convolution
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals
-        run: |
-          chmod +x *.sh
-          ./run_grt.sh
-          python run_pygrt.py
+      # - name: Test signal convolution
+      #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals
+      #   run: |
+      #     chmod +x *.sh
+      #     ./run_grt.sh
+      #     python run_pygrt.py
 
-      - name: Test strain stress (dynamic)
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/dynamic
-        run: |
-          chmod +x *.sh
-          ./run_grt.sh
+      # - name: Test strain stress (dynamic)
+      #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/dynamic
+      #   run: |
+      #     chmod +x *.sh
+      #     ./run_grt.sh
 
-      - name: Test static displacement (static)
-        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/static
-        run: |
-          chmod +x *.sh
-          ./run_stgrt.sh
+      # - name: Test static displacement (static)
+      #   working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/static
+      #   run: |
+      #     chmod +x *.sh
+      #     ./run_stgrt.sh
 
       - name: Remove the test files
         run: |

example/compare_c_py/README

@@ -0,0 +1 @@
+Compare results from C API and Python API.

example/compare_c_py/milrow

@@ -0,0 +1,9 @@
+0.2  3.4  1.7  2.3  9e10  9e10
+0.6  3.7  1.9  2.4  9e10  9e10
+0.5  4.2  2.1  2.4  9e10  9e10
+0.5  4.6  2.3  2.5  9e10  9e10
+0.7  4.9    2.8  2.6  9e10  9e10
+0.5  5.1  2.9  2.7  9e10  9e10
+6.0  5.9  3.3  2.7  9e10  9e10
+28.  6.9  4.0  2.8  9e10  9e10
+-0.1   8.2  4.7  3.2  9e10  9e10

example/compare_c_py/run_grt.sh

@@ -0,0 +1,90 @@
+#!/bin/bash
+
+
+dist=10
+depsrc=2
+deprcv=3.3
+
+nt=1024
+dt=0.01
+
+modname="milrow"
+out="GRN"
+
+
+#-------------------------- Dynamic -----------------------------------------
+grt -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} -e
+
+# convolve different signals
+G=$out/${modname}_${depsrc}_${deprcv}_${dist}
+P=cout
+S=1e24
+az=39.2
+for N in "-N" "" ; do
+grt.syn -G$G -Osyn_exp -P$P -A$az -S$S -Dt/0.2/0.2/0.4 -e $N
+grt.strain syn_exp/$P
+grt.stress syn_exp/$P
+
+fn=2
+fe=-1
+fz=4
+grt.syn -G$G -Osyn_sf -P$P -A$az -S$S -F$fn/$fe/$fz -Dt/0.1/0.3/0.6 -e $N
+grt.strain syn_sf/$P
+grt.stress syn_sf/$P
+
+stk=77
+dip=88
+rak=99
+grt.syn -G$G -Osyn_dc -P$P -A$az -S$S -M$stk/$dip/$rak -Dp/0.6 -e $N
+grt.strain syn_dc/$P
+grt.stress syn_dc/$P
+
+M11=1
+M12=-2
+M13=-5
+M22=0.5
+M23=3
+M33=1.2
+grt.syn -G$G -Osyn_mt -P$P -A$az -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -Dr/3 -e $N
+grt.strain syn_mt/$P
+grt.stress syn_mt/$P
+done
+
+
+#-------------------------- Static -----------------------------------------
+x1=-3
+x2=3
+nx=11
+
+y1=-4
+y2=4
+ny=16
+
+mkdir -p static
+cd static
+stgrt -M../${modname} -D${depsrc}/${deprcv} -X$x1/$x2/$nx -Y$y1/$y2/$ny -e > grn
+
+for N in "-N" "" ; do
+stgrt.syn -S$S -e $N < grn > stsyn_exp$N
+stgrt.strain < stsyn_exp$N > strain_exp$N
+stgrt.stress < stsyn_exp$N > stress_exp$N
+
+stgrt.syn -S$S -F$fn/$fe/$fz -e $N < grn > stsyn_sf$N
+stgrt.strain < stsyn_sf$N > strain_sf$N
+stgrt.stress < stsyn_sf$N > stress_sf$N
+
+stgrt.syn -S$S -M$stk/$dip/$rak -e $N < grn > stsyn_dc$N
+stgrt.strain < stsyn_dc$N > strain_dc$N
+stgrt.stress < stsyn_dc$N > stress_dc$N
+
+stgrt.syn -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -e $N < grn > stsyn_mt$N
+stgrt.strain < stsyn_mt$N > strain_mt$N
+stgrt.stress < stsyn_mt$N > stress_mt$N
+done
+
+cd -
+
+for p in $(ls static/*); do
+    echo "----------------------------- $p -------------------------"
+    cat $p
+done

example/compare_c_py/run_pygrt.py

@@ -0,0 +1,201 @@
+import numpy as np
+import pygrt 
+from obspy import *
+
+
+def compare3(st_py:Stream, c_prefix:str, ZNE:bool=False, dim2:bool=False):
+    """return average relative error"""
+
+    if ZNE:
+        pattern = "[ZNE]"
+    else:
+        pattern = "[ZRT]"
+
+    if dim2:
+        pattern = pattern*2
+
+    st_c = read(f"{c_prefix}{pattern}.sac")
+
+    error = 0.0
+    nerr = 0
+    for tr_c in st_c:
+        tr_py = st_py.select(channel=tr_c.stats.channel)[0]
+        rerr = np.mean(np.abs(tr_c.data - tr_py.data) / np.abs(tr_c.data))
+        if np.isnan(rerr) or np.isinf(rerr):
+            rerr = 0.0
+        error += rerr
+        nerr += 1 
+
+    return error/nerr
+
+
+def static_compare3(resDct:dict, c_prefix:str, ZNE:bool=False, dim2:bool=False):
+    """return average relative error"""
+    midname = c_prefix.split("_")[-1][:2].upper()
+
+    if ZNE:
+        chs2 = ['Z', 'N', 'E']
+    else:
+        chs2 = ['Z', 'R', 'T']
+
+    chs1 = chs2.copy()
+
+    nx = len(resDct['_xarr'])
+    ny = len(resDct['_yarr'])
+
+    c_res = np.loadtxt(c_prefix)
+    c_resDct = {}
+    if not dim2:
+        for i, c in enumerate(chs2):
+            c_resDct[f"{midname}{c}"] = c_res[:,i+2].reshape((nx, ny), order='F')
+    else:
+        chs1 = []
+        for i1, c1 in enumerate(chs2):
+            for i2, c2 in enumerate(chs2):
+                if i2 < i1:
+                    continue
+                chs1.append(f"{c1}{c2}")
+
+        for i, cc in enumerate(chs1):
+            c_resDct[cc] = c_res[:,i+2].reshape((nx, ny), order='F')
+
+    error = 0.0
+    nerr = 0
+
+    for k in c_resDct.keys():
+        val1 = resDct[k]
+        val2 = c_resDct[k]
+        rerr = np.mean(np.abs(val1 - val2) / np.abs(val2))
+        if np.isnan(rerr) or np.isinf(rerr):
+            rerr = 0.0
+        error += rerr
+        nerr += 1 
+
+    return error/nerr
+
+
+
+dist=10
+depsrc=2
+deprcv=3.3
+
+nt=1024
+dt=0.01
+
+modname="milrow"
+
+modarr = np.loadtxt(modname)
+
+pymod = pygrt.PyModel1D(modarr, depsrc, deprcv)
+
+AVGRERR = []
+
+#-------------------------- Dynamic -----------------------------------------
+# compute green functions
+st_grn = pymod.compute_grn(dist, nt, dt, calc_upar=True)[0]
+
+S=1e24
+az=39.2
+
+fn=2
+fe=-1
+fz=4
+
+stk=77
+dip=88
+rak=99
+
+M11=1
+M12=-2
+M13=-5
+M22=0.5
+M23=3
+M33=1.2
+
+for ZNE in [False, True]:
+    # synthetic
+
+    st = pygrt.utils.gen_syn_from_gf_EXP(st_grn, S, az, ZNE=ZNE, calc_upar=True)
+    sigs = pygrt.sigs.gen_triangle_wave(0.4, dt)
+    pygrt.utils.stream_convolve(st, sigs)
+    AVGRERR.append(compare3(st, "syn_exp/cout", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(st)
+    ststress = pygrt.utils.compute_stress(st)
+    AVGRERR.append(compare3(ststrain, "syn_exp/cout.strain.", ZNE=ZNE, dim2=True))
+    AVGRERR.append(compare3(ststress, "syn_exp/cout.stress.", ZNE=ZNE, dim2=True))
+
+
+    st = pygrt.utils.gen_syn_from_gf_SF(st_grn, S, fn, fe, fz, az, ZNE=ZNE, calc_upar=True)
+    sigs = pygrt.sigs.gen_trap_wave(0.1, 0.3, 0.6, dt)
+    pygrt.utils.stream_convolve(st, sigs)
+    AVGRERR.append(compare3(st, "syn_sf/cout", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(st)
+    ststress = pygrt.utils.compute_stress(st)
+    AVGRERR.append(compare3(ststrain, "syn_sf/cout.strain.", ZNE=ZNE, dim2=True))
+    AVGRERR.append(compare3(ststress, "syn_sf/cout.stress.", ZNE=ZNE, dim2=True))
+
+    st = pygrt.utils.gen_syn_from_gf_DC(st_grn, S, stk, dip, rak, az, ZNE=ZNE, calc_upar=True)
+    sigs = pygrt.sigs.gen_parabola_wave(0.6, dt)
+    pygrt.utils.stream_convolve(st, sigs)
+    AVGRERR.append(compare3(st, "syn_dc/cout", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(st)
+    ststress = pygrt.utils.compute_stress(st)
+    AVGRERR.append(compare3(ststrain, "syn_dc/cout.strain.", ZNE=ZNE, dim2=True))
+    AVGRERR.append(compare3(ststress, "syn_dc/cout.stress.", ZNE=ZNE, dim2=True))
+
+    st = pygrt.utils.gen_syn_from_gf_MT(st_grn, S, [M11,M12,M13,M22,M23,M33], az, ZNE=ZNE, calc_upar=True)
+    sigs = pygrt.sigs.gen_ricker_wave(3, dt)
+    pygrt.utils.stream_convolve(st, sigs)
+    AVGRERR.append(compare3(st, "syn_mt/cout", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(st)
+    ststress = pygrt.utils.compute_stress(st)
+    AVGRERR.append(compare3(ststrain, "syn_mt/cout.strain.", ZNE=ZNE, dim2=True))
+    AVGRERR.append(compare3(ststress, "syn_mt/cout.stress.", ZNE=ZNE, dim2=True))
+
+
+#-------------------------- Static -----------------------------------------
+xarr = np.linspace(-3, 3, 11)
+yarr = np.linspace(-4, 4, 16)
+static_grn = pymod.compute_static_grn(xarr, yarr, calc_upar=True)
+
+for ZNE in [False, True]:
+    suffix = "-N" if ZNE else ""
+    static_syn = pygrt.utils.gen_syn_from_gf_EXP(static_grn, S, ZNE=ZNE, calc_upar=True)
+    AVGRERR.append(static_compare3(static_syn, f"static/stsyn_exp{suffix}", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(static_syn)
+    ststress = pygrt.utils.compute_stress(static_syn)
+    AVGRERR.append(static_compare3(ststrain, f"static/strain_exp{suffix}", ZNE=ZNE, dim2=True))
+    AVGRERR.append(static_compare3(ststress, f"static/stress_exp{suffix}", ZNE=ZNE, dim2=True))
+    print(static_syn, ststrain, ststress)
+
+    static_syn = pygrt.utils.gen_syn_from_gf_SF(static_grn, S, fn, fe, fz, ZNE=ZNE, calc_upar=True)
+    AVGRERR.append(static_compare3(static_syn, f"static/stsyn_sf{suffix}", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(static_syn)
+    ststress = pygrt.utils.compute_stress(static_syn)
+    AVGRERR.append(static_compare3(ststrain, f"static/strain_sf{suffix}", ZNE=ZNE, dim2=True))
+    AVGRERR.append(static_compare3(ststress, f"static/stress_sf{suffix}", ZNE=ZNE, dim2=True))
+    print(static_syn, ststrain, ststress)
+
+    static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, S, stk, dip, rak, ZNE=ZNE, calc_upar=True)
+    AVGRERR.append(static_compare3(static_syn, f"static/stsyn_dc{suffix}", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(static_syn)
+    ststress = pygrt.utils.compute_stress(static_syn)
+    AVGRERR.append(static_compare3(ststrain, f"static/strain_dc{suffix}", ZNE=ZNE, dim2=True))
+    AVGRERR.append(static_compare3(ststress, f"static/stress_dc{suffix}", ZNE=ZNE, dim2=True))
+    print(static_syn, ststrain, ststress)
+
+    static_syn = pygrt.utils.gen_syn_from_gf_MT(static_grn, S, [M11,M12,M13,M22,M23,M33], ZNE=ZNE, calc_upar=True)
+    AVGRERR.append(static_compare3(static_syn, f"static/stsyn_mt{suffix}", ZNE=ZNE))
+    ststrain = pygrt.utils.compute_strain(static_syn)
+    ststress = pygrt.utils.compute_stress(static_syn)
+    AVGRERR.append(static_compare3(ststrain, f"static/strain_mt{suffix}", ZNE=ZNE, dim2=True))
+    AVGRERR.append(static_compare3(ststress, f"static/stress_mt{suffix}", ZNE=ZNE, dim2=True))
+    print(static_syn, ststrain, ststress)
+
+
+AVGRERR = np.array(AVGRERR)
+print(AVGRERR)
+print(np.mean(AVGRERR), np.min(AVGRERR), np.max(AVGRERR))
+if np.mean(AVGRERR) > 0.01:
+    raise ValueError
+