python大地电磁正演

Question

    # ---------- 1. 导入（大幅减少单元数）----------
    import warnings
    from simpeg.utils import BreakingChangeWarning
    warnings.filterwarnings('ignore', category=BreakingChangeWarning)
    warnings.filterwarnings('ignore', category=UserWarning)
    warnings.filterwarnings('ignore', category=FutureWarning)

     import discretize
     from simpeg.electromagnetics import natural_source as NSEM
     from simpeg import utils
     import numpy as np
     import matplotlib.pyplot as plt

    # 正确导入求解器
     from pymatsolver import Solver

     def run(plotIt=True):
    print("=" * 50)
    print("快速测试模式（粗网格 + 少量频率，使用 pymatsolver.Solver）")
    print("=" * 50)

    # ---------- 2. 粗网格（大幅减少单元数）----------
    dx = 25           # 网格步长 50 m
    nx = int(600 / dx)   # 12
    ny = int(600 / dx)   # 12
    dz = 20           # 步长 25 m
    nz = int(1200 / dz)  # 48

    M = discretize.TensorMesh(
        [
            [(dx, nx)],
            [(dx, ny)],
            [(dz, nz)]
        ],
        x0=[-300, -300, -1200]
    )
    print(f"网格单元数: {M.n_cells} (nx={nx}, ny={ny}, nz={nz})")

    # ---------- 2. 电导率模型（基岩 + 立方体含水层）----------
    sig = np.ones(M.nC) * (1 / 2000)   # 基岩 2000 Ω·m
    cube = (
        (M.gridCC[:, 0] >= -35) & (M.gridCC[:, 0] <= 35) &
        (M.gridCC[:, 1] >= -35) & (M.gridCC[:, 1] <= 35) &
        (M.gridCC[:, 2] >= -150) & (M.gridCC[:, 2] <= -80)
    )
    sig[cube] = 1 / 50                 # 含水层 50 Ω·m
    sigBG = np.ones(M.nC) * (1 / 2000) # 背景仅为基岩

    # ---------- 3. 模型切片图 ----------
    if plotIt:
        fig, ax = plt.subplots(figsize=(7, 5))
        out = M.plot_slice(
            np.log10(sig),
            normal="X",
            ind=int(nx/2),
            ax=ax,
            grid=True,
            clim=(-4, 0),
            cmap="viridis"
        )[0]
        plt.colorbar(out, ax=ax, label='log10(σ)')
        ax.set_xlim(-300, 300)
        ax.set_ylim(-1200, 0)
        ax.set_title("Model (X=0 slice) - Quick")
        ax.set_xlabel("Y (m)")
        ax.set_ylabel("Z (m)")
        plt.show()
        print("模型切片图已显示")

    # ---------- 4. 测线（x=0，y 从 -300 到 300，间距 25 m）----------
    y_line = np.arange(-300, 301, 25)   # 25 个测点
    n_rx = len(y_line)
    rx_loc = np.c_[np.zeros_like(y_line), y_line, np.zeros_like(y_line)]

    # ---------- 5. 接收器（TM: yx, TE: xy）----------
    receiver_list = []
    for ori in ["xy", "yx"]:
        receiver_list.append(NSEM.Rx.Impedance(rx_loc, orientation=ori, component="real"))
        receiver_list.append(NSEM.Rx.Impedance(rx_loc, orientation=ori, component="imag"))

    # ---------- 6. 仅 3 个频率（10, 1000, 100000 Hz）----------
    freqs = np.logspace(1, 5, 3)
    n_freq = len(freqs)
    print(f"频率数: {n_freq} 个 ({freqs[0]:.0f}, {freqs[1]:.0f}, {freqs[2]:.0f} Hz)")

    source_list = [NSEM.Src.PlanewaveXYPrimary(receiver_list, f) for f in freqs]
    survey = NSEM.Survey(source_list)

    # ---------- 7. 正演（使用 pymatsolver.Solver，会尝试迭代求解器）----------
    print("开始正演计算（使用 pymatsolver.Solver，每个频率约 10~30 秒）...")
    problem = NSEM.Simulation3DPrimarySecondary(
        M,
        survey=survey,
        sigma=sig,
        sigmaPrimary=sigBG,
        solver=Solver,
        forward_only=True
    )
    data = problem.dpred()
    print("正演计算完成！")

    # ---------- 8. 提取 Zyx（TM模式）视电阻率 ----------
    print("提取视电阻率...")
    rho_app = np.zeros((n_freq, n_rx))
    mu0 = 4 * np.pi * 1e-7
    n_data_per_freq = len(receiver_list) * n_rx

    for i_f, f in enumerate(freqs):
        print(f"  处理频率 {f:.0f} Hz...")
        for i_rx in range(n_rx):
            base = i_f * n_data_per_freq
            idx_yx_real = base + (2 * 1) * n_rx + i_rx
            idx_yx_imag = base + (2 * 1 + 1) * n_rx + i_rx
            Zyx = data[idx_yx_real] + 1j * data[idx_yx_imag]
            rho_app[i_f, i_rx] = (np.abs(Zyx) ** 2) / (mu0 * 2 * np.pi * f)

    # ---------- 9. 拟断面图 ----------

    # =========================
    # 9. TM二维剖面图（绝对可靠版本）
    # =========================
    if plotIt:
        # 先诊断数据
        print("\n=== 绘图诊断 ===")
        print(f"频率数组: {freqs}")
        print(f"频率范围: {freqs.min():.2e} ~ {freqs.max():.2e} Hz")
        print(f"视电阻率形状: {rho_app.shape}")
        print(f"视电阻率有效范围: {np.nanmin(rho_app):.2e} ~ {np.nanmax(rho_app):.2e} Ω·m")
        print(f"NaN 数量: {np.isnan(rho_app).sum()}")

        fig2, ax2 = plt.subplots(figsize=(7, 5))

        # 方法1：使用 pcolormesh 且显式指定边缘（更准确）
        # 构建边缘坐标：测点边界和频率边界
        y_edges = np.linspace(-325, 325, len(y_line) + 1)  # 25个点 -> 26个边界
        freq_edges = np.logspace(np.log10(freqs.min()), np.log10(freqs.max()), len(freqs) + 1)

        X, Y = np.meshgrid(y_edges, freq_edges)
        # 注意 rho_app 的 shape 是 (n_freq, n_stations)，与 (n_freq, n_stations) 一致
        im = ax2.pcolormesh(X, Y, np.log10(rho_app), shading='flat', cmap='jet')

        ax2.set_yscale("log")
        ax2.set_xlabel("Distance (m)")
        ax2.set_ylabel("Frequency (Hz)")
        ax2.set_title("TM Apparent Resistivity Section (Zyx) - Fixed")
        cbar = plt.colorbar(im, ax=ax2, label='log10(ρa) [Ω·m]')

        # 可选：在图上标出实际频率点位置
        for f in freqs:
            ax2.axhline(y=f, color='white', linestyle=':', linewidth=0.5, alpha=0.5)

        plt.tight_layout()
        plt.show()
        print("拟断面图已显示\n")
    #if plotIt:
    #    fig2, ax2 = plt.subplots(figsize=(7, 5))
    #    Y, F = np.meshgrid(y_line, freqs)
    #    im = ax2.pcolormesh(Y, F, np.log10(rho_app), shading='auto', cmap='jet')
    #    ax2.set_yscale("log")
    #    ax2.set_xlabel("Distance (m)")
    #    ax2.set_ylabel("Frequency (Hz)")
    #    ax2.set_title("TM Apparent Resistivity - Quick (3 freqs)")
    #    plt.colorbar(im, ax=ax2, label='log10(ρa)')
    #    plt.tight_layout()
    #    plt.show()
    #    print("拟断面图已显示")

    print("\n快速测试完成！")

    if __name__ == "__main__":
     run(plotIt=True)

这是我的python立方体大地电磁正演代码,得出一个模型图和视电阻率拟断面图,若我将dx和dz改为25和20,那么一段时间后能运行出来,但是我最终目的是dx=25,dz=20,并且频率是freqs = np.logspace(1, 5, 27),但是现在仅有三个频率都运行不出来,请问我需要怎么改才能运行出来呢?或者有没有uu帮忙运行一下,谢谢!

Answer 1

首先，我们需要分析你的问题。你提到将dx和dz改为25和20后，代码可以运行，但当你尝试使用更高的频率分辨率（freqs = np.logspace(1, 5, 27)）时，代码运行失败。这可能是由于计算资源不足、内存溢出或者计算时间过长导致的。

以下是一些可能的解决方案：

1. 优化代码：

   • 确保你的代码中没有不必要的计算或者可以优化的部分。
   • 检查是否有重复的计算或者可以并行处理的部分。

2. 调整频率范围：

   • 减少频率的数量，比如使用freqs = np.logspace(1, 5, 11)或者更少，看看是否可以运行。
   • 如果减少频率数量后可以运行，再逐渐增加频率数量，直到找到无法运行的最大频率数量。

3. 增加计算资源：

   • 如果可能，尝试在具有更多内存和计算能力的机器上运行代码。
   • 如果你在服务器上运行，可以考虑增加服务器的资源。

4. 优化网格和模型：

   • 检查网格是否过于密集，如果可能，尝试使用更粗的网格。
   • 确保电导率模型不会导致计算不稳定。

5. 使用更高效的求解器：

   • 如果你使用的是pymatsolver.Solver，尝试使用其他求解器，比如scipy.sparse.linalg.spsolve。

下面是一个修改后的代码示例，其中减少了频率的数量，并尝试使用scipy的求解器：

from scipy.sparse.linalg import spsolve
import numpy as np

# ...（其他代码保持不变）

# 调整频率范围
freqs = np.logspace(1, 5, 11)

# ...（其他代码保持不变）

# 使用scipy的求解器
problem = NSEM.Simulation3DPrimarySecondary(
    M,
    survey=survey,
    sigma=sig,
    sigmaPrimary=sigBG,
    solver=spsolve,
    forward_only=True
)

# ...（其他代码保持不变）

请注意，这只是一个示例，你可能需要根据实际情况调整代码。

最后，如果你需要我帮你运行代码，请提供更多的信息，比如你的计算环境、代码的具体错误信息等。以上内容参考 GLM大模型编写提供！