本程序基于 FDTD（时域有限差分）方法对三维瞬变电磁（TEM）响应进行正演模拟。核心算法采用 Wang–Hohmann (1993)改进的 Du Fort–Frankel 方法：在 Yee 网格上对磁场 H 迭代求解，通过引入虚介电常数保证显式迭代的时间稳定性；采用共形网格技术处理起伏地形与任意形状异常体（以表面三角网格描述）；支持矩形回线源，可计算地面 TEM 、半航空（SATEM）全航空（ATEM）等模式。

程序的整体框架与三大核心技术分别源于以下工作：

整体框架与核心迭代算法（考虑关断时间的回线源激发 TEM 三维时域有限差分正演，Wang–Hohmann 改进的 Du Fort–Frankel 方法、虚介电常数、含关断时间的源波形），理论内容请见参考文献[1]；
CPML 吸收边界（瞬变电磁低频近似 Maxwell 方程的 CPML 吸收边界及施加方法），理论内容请见参考文献[2]；
共形网格技术（通过射线追踪方法将任意复杂形状的结构引入到Yee网格计算中），理论内容请见参考文献[3]。

各方法的具体原理、公式推导与实现细节详见文献 [1]–[3]，这些文献的作者都是对本开源项目做出突出贡献的人员。

代码结构：main.f90（主程序）、module/（全局参数与模块）、lib/（各功能子程序）。

1. 运行环境与编译

1.1 环境要求

项目	要求
操作系统	Windows 10/11(64 位)或 Linux(x86_64)，特别推荐使用国产操作系统 deepin 25
集成环境	Windows:Visual Studio 2019及以上;Linux/deepin:make；
Fortran 编译器	Intel oneAPI Fortran(ifort/ifx);
并行支持	OpenMP(多核 CPU 加速);额外的GPU支持在商业版中提供，请访问https://em3d.cn

1.2 用 VS2019 编译运行(推荐)

安装 Visual Studio 2019(勾选"C++ 桌面开发"工作负载)与 Intel oneAPI(安装时勾选 "Intel Fortran Compiler 的 Visual Studio 集成")。
双击打开工程文件 tem3dfdtd.sln(Intel Fortran 工程,对应 tem3dfdtd.vfproj)。工程已包含全部源文件,配置说明:
- Debug | x64 / Release | x64:使用 ifx 编译器(oneAPI 默认),推荐;
- Debug | Win32 / Release | Win32:使用 ifort 编译器。
选择 Release | x64 配置,点击 生成 → 生成解决方案。
运行前准备:程序在工作目录中查找 input.dat 及网格文件,因此请将 input.dat、Complex_Terrain.dat/.stl、Complex_anomalous.dat/.stl 放到 tem3dfdtd-open\tem3dfdtd目录中(或通过"项目属性 → 调试 → 工作目录"指定)。
直接运行 tem3dfdtd.exe,或按 F5 调试运行。

注意:程序使用 OpenMP,运行时需要 Intel 的运行时库 libiomp5md.dll (位于 Intel oneAPI 安装目录 bin/ 下)。若提示缺少该 DLL,可将它复制到 exe 同目录(本目录已放置一份),或将其所在目录加入系统 PATH。

1.3 命令行编译(可选，不建议在Windows下使用)

在"Intel oneAPI 命令行"环境(oneAPI Command Prompt)下:

ifx -c -O2 -Qopenmp module\*.f90
ifx -c -O2 -Qopenmp -Qopenmp lib\*.f90 main.f90
ifx -O2 -Qopenmp *.obj -o tem3dfdtd.exe

(链接时需要 MSVC 的 link.exe 与 Windows SDK 库,建议直接使用 VS 的"开发人员命令提示符", 并在 PATH 中加入 Intel oneAPI 的 bin 目录。)

1.4 Linux 环境编译(makefile)

在Linux下编译需要编写 makefile文件，请根据操作系统的配置和要求自行编写makefile，并使用make makefile进行编译。

说明与注意事项

与 Windows 版的差异:Linux 下可执行文件名为自定义，如果没有指定则默认为 main.exe,工作目录中同样需要 input.dat 与 Complex_Terrain.*、Complex_anomalous.* 网格文件。

2. 程序流程与模块结构

主程序执行顺序(main.f90):

GETDATA → CHECKPARAMETERS → MEMORY_USE_ESTIMATION → ALLOCATEMEMORY
→ GET_NON_UNIFORMGRID → ZERO → GET_COORDINATES
→ Get_Receiver_Gridlabel → RES_CONFIGURE → TIME_SERIOUS
→ Get_eps_r → (Logic_PML=1 时) Get_pml_parameters → Get_mstop
→ GetSourcePosition → Iteration → FREE_MEMORY

按功能划分为以下 6 个模块:

模块 1:程序控制与参数输入

文件	功能
`main.f90`	主程序,控制整个计算流程
`lib/getdata.f90`	读取参数控制文件 `input.dat`;检测地形/异常体网格文件的存在性并选择读取格式
`lib/checkparameters.f90`	将读入的计算参数回显到 `logfile.log`,便于人工检查
`lib/memory-use-estimation.f90`	根据网格规模估算所需内存并打印提示

模块 2:Yee 网格生成

文件	功能
`lib/allocatememory.f90`	根据输入参数动态分配所有全局数组(含 CPML 记忆变量数组,仅 `Logic_PML=1` 时)
`lib/get_non_uniformgrid.f90`	生成 x/y/z 三方向的非均匀网格(核心区均匀 + 外围按 1.3 倍递增扩展)
`lib/get_coordinates.f90`	计算各网格节点(含 Yee 节点)的坐标,坐标以源中心为原点
`lib/zero.f90`	将所有电磁场数组初始化为 0;den_* 置 1、c_h_zz 置 0、CPML 记忆变量清零

模块 3:电性参数构建

文件	功能
`lib/resistivity-configuration.f90`	构建模型电导率:无地形时按背景电导率 + 块状异常体赋值;有地形时调用共形网格;最后将电导率分配到 x/y/z 三个方向的棱边并写出 `conductivity.vtk`
`lib/Terrain_conformal.f90`	地形共形网格:从 `Complex_Terrain.dat/.stl` 读入地形三角网格,用射线–三角形求交(Möller–Trumbore 算法)沿 x/y/z 三个方向填充每个棱边的等效电导率,处理起伏地形与空气/地层的分界
`lib/Anomalous_conformal.f90`	异常体共形网格:从 `Complex_anomalous.dat/.stl` 读入异常体表面网格,采用同样的射线求交方法将异常体电导率(`tao_abnormal`)填充到棱边

说明:存在地形文件时,地形分支中空气电导率取 AIR_CONDUCTIVITY = 1e-6 S/m, 地层电导率取块状异常体参数中的第 2 个电导率 TAR_CONDUCTIVITY(2)(见 input.dat 第 9 行起的第二组数据);此时 input.dat 中块状异常体本身不直接生效,而是以 Complex_anomalous 网格文件描述异常体、以 tao_abnormal 赋予其电导率。

模块 4:激励源与时间序列

文件	功能
`lib/time-serious.f90`	生成整个计算的时间序列(含源波形),并根据 `MAX_OFF_TIME` 校正迭代步数 `NSTOP`;写出 `CTIME_TIXING_UPCOS.DAT`
`lib/tixing-source-upcos.f90`	源波形:梯形 + 余弦上升的关断电流波形(常用,`SOURCE_TYPE = 'TIXING_UPCOS'`)
`lib/tixing-source.f90`	纯梯形波形源(`TIXING_RAMP`)
`lib/sin-source.f90`	半正弦波形源(`HALF_SIN`)
`lib/triangle-source.f90`	三角波形源(`TRIANGLE`)
`lib/get-eps-r.f90`	计算虚介电常数 `EPS_R = 3·(Δt/Δx)²/μ₀` 及迭代系数,保证显式 FDTD 稳定
`lib/get-mstop.f90`	将总迭代过程切成若干"计算分段",每段独立分配缓存,便于内存管理

模块 5:FDTD 电磁场计算

文件	功能
`lib/GetSourcePosition.f90`	确定回线源在网格中的位置,标记源所在的棱边(电流赋值区域)
`lib/Iteration.f90`	核心迭代子程序:按分段循环推进时间步,更新 Ex/Ey/Ez、Hz 场,按源波形加载电流;在每个分段末尾对各接收点计算 Hz(由环绕该点的 8 个节点 Ex、Ey 差商加权得到)并写入结果文件。边界条件按 `Logic_PML` 开关切换:`1` 时在每个场更新中附加 CPML 记忆变量修正,`0` 时恢复原始 Dirichlet(零场)边界
`module/pml-parameters.f90`	CPML 吸收边界模块(Roden–Gedney 卷积 PML):σ/α/κ 多项式缩放参数(ma=3、mb=1)、26 个记忆变量 ψ 数组、be/c_e 卷积系数数组与 den*(=1/κ)缩放数组的声明
`lib/get-pml-paramters.f90`	构建 x/y/z 六个边界面的 σ/α/κ 分布(多项式从边界向内衰减,E/H 交错采样)及 den_* 缩放数组;仅在 `Logic_PML=1` 时由 main 调用

模块 6:输出

文件	功能
`lib/Iteration.f90`(输出部分)	写出各接收点响应文件 `dBzdt_1.txt`、`dBzdt_2.txt` …
`lib/resistivity-configuration.f90`(输出部分)	写出模型电导率分布 `conductivity.vtk`
`lib/time-serious.f90`(输出部分)	写出时间序列 `CTIME_TIXING_UPCOS.DAT`
`lib/free-memory.f90`	计算结束后释放所有动态内存

3. 输入文件格式

程序运行需要以下文件(全部放在 exe 的工作目录中):

文件	是否必须	说明
`input.dat`	必须	计算参数控制文件
`Complex_anomalous.dat` 或 `.stl`	可选(任选其一)	异常体表面三角网格;缺省时模型视为均匀背景
`Complex_Terrain.dat` 或 `.stl`	可选(任选其一)	地形表面三角网格;缺省时不考虑起伏地形

3.1 参数控制文件 input.dat

自由格式读取,按行顺序读取;数值后可加 ! 注释(可整行注释或行尾注释)。下面以本目录自带的 input.dat 为例逐行说明:

行号	示例	含义
1	`1`	计算模式 `CAL_TYPE`:`1` = 地面 TEM,`2` = 半航空(SATEM)
2	`500`	发射回线边长 `SourceLength`(m)
3	`101,101,100`	x、y、z 三个方向的网格数 `NX,NY,NZ`
4	`1`	边界条件开关 `Logic_PML`:`1` = CPML 吸收边界,`0` = 原始非均匀网格 Dirichlet(零场)边界
5	`10,10,10`	PML 层数 `PML_X,PML_Y,PML_Z`(x、y、z 方向,仅开关=1 时有效;建议 ≥ 5 层)
6	`25,25`	x 方向:核心均匀网格起始/结束区间编号 `UniGridNumX1,UniGridNumX2`
7	`25,25`	y 方向:核心均匀网格区间编号 `UniGridNumY1,UniGridNumY2`
8	`20,30`	z 方向:核心均匀网格区间编号 `UniGridNumZ1,UniGridNumZ2`
9	`20`	核心区均匀网格尺寸 `GridSize`(m)
10	`0.01`	背景介质电导率 `BACKGROUND_CONDUCTIVITY`(S/m)
11	`2`	块状异常体数量 `TEMP_II`(无地形时按棱柱体填充;设 0 表示均匀模型)
12–15	见下	第 1 个块状异常体参数,共 4 行
16–19	见下	第 2 个块状异常体参数,共 4 行
20	`4000000`	最大迭代次数 `NSTOP`
21	`90.101`	最大计算时间 `MAX_OFF_TIME`(单位 ms)
22	`1e-6,1e-9`	上升沿持续时间与时间步 `RAISETIME, RAISESTEP`(s)
23	`60000e-6`	平台阶段持续时间 `WAVE`(s,即 60 ms)
24	`1e-7,1e-9`	下降沿持续时间与时间步 `RAMP, RAMPSTEP`(s)
25	`1e-9`	初始时间步 `TIMESTEP`(s)
26	`1`	发射电流幅度 `AMP`(A)
27	`4.0`	异常体电导率 `tao_abnormal`(S/m,配合 `Complex_anomalous` 文件使用)
28	`TIXING_UPCOS`	源类型 `SOURCE_TYPE`:`TIXING_UPCOS` / `TIXING_RAMP` / `HALF_SIN` / `TRIANGLE`
29	`1`	接收点数 `Point_Num`
30	`1`	第 1 个接收点的编号
31	`0,0,0`	第 1 个接收点坐标(相对源中心,m)
32–33	`2` / `x,y,z`	(如有多余测点)第 2 个接收点(编号 + 坐标)

每个块状异常体由连续的 4 行组成:

行	示例	含义
①	`1,101`	x 方向网格起止编号 `TAR_X1, TAR_X2`
②	`1,101`	y 方向网格起止编号 `TAR_Y1, TAR_Y2`
③	`1,50`	z 方向网格起止编号 `TAR_Z1, TAR_Z2`
④	`1e-5`	该块电导率 `TAR_CONDUCTIVITY`(S/m)

本算例用两块"异常体"拼出半空间:第 1 块 z=1~~50(空气,1e-5 S/m)+ 第 2 块 z=51~~100(地层,1e-2 S/m)。接收点个数 Point_Num 后按每个测点 2 行排列 (编号 + 相对源中心坐标)。

3.2 地形网格文件 Complex_Terrain

地形由表面三角网格描述,支持两种格式,文件夹中只保留其中一个; 若两个同时存在,程序以 .dat 为优先并提示 .stl 被忽略。

格式 1:Complex_Terrain.dat(原始文本格式)

Number of Nodes and Elements:
10039               ← 节点总数 n_point
5426                ← 三角形单元总数 n_face
Nodes Coordinates:
1  -21000.0  -21000.0  224.08   ← n_point 行:节点编号, X, Y, Z
2  -21000.0  -20001.8  224.08
...                        (行中可带 ! 注释)
END Nodes Coordinates
NormalAreaElements:
1  1  2  10039       ← n_face 行:单元编号, 节点1, 节点2, 节点3
...
END NormalAreaElements

行	内容
第 1 行	标题行,可任意
第 2 行	节点总数 `n_point`
第 3 行	三角形单元总数 `n_face`
第 4 行	标题行,可任意
第 5 ~ 4+n_point 行	每个节点一行:`节点编号, X, Y, Z`
其后 1 行	区段结束标记 `END Nodes Coordinates`(程序按标题行跳过)
其后 1 行	面区标题 `NormalAreaElements:`(程序按标题行跳过)
其后 n_face 行	每个单元一行:`单元编号, 节点1编号, 节点2编号, 节点3编号`(节点按逆时针绕向)
末尾 1 行	结束标记 `END NormalAreaElements`(程序不读取)

节点区后的两个区段标记行(GiD 导出)与老版本格式的"2 行标题行"位置一致, 程序一律按标题行跳过,因此两种写法均兼容。

格式 2:Complex_Terrain.stl(ASCII STL 格式)

STL文件格式是一种用于描述三维物体表面几何形状的文件格式，广泛应用于快速成型、3D打印和计算机辅助制造（CAM）领域。 STL文件将物体表面细分为一系列小三角形，每个三角形由一个法线向量和三个顶点坐标来定义。 STL文件有两种格式：文本格式（ASCII）和二进制格式。

标准 ASCII STL,单元关键字为 facet/endfacet, 节点用 vertex 行表示，例如下面的格式:

facet normal nx ny nz
    outer loop
        vertex v1x v1y v1z
        vertex v2x v2y v2z
        vertex v3x v3y v3z
    endloop
endfacet

用户可以使用常用的CAD软件(如AutoCAD Blender FreeCAD MeshLab SketchUp Gid Maya、3ds Max等)创建和编辑STL文件。

程序读取时自动处理两点(无需用户操作):

顶点去重合并:STL 中每个面独立写顶点,重复顶点(容差 1e-5)自动合并为唯一节点表;
方向校正:比较每个面的叉积方向与文件中的 facet normal,若相反则交换该面第 2、3 个节点,保证法线方向约定与 .dat 格式一致。

3.3 异常体网格文件 Complex_anomalous

与地形文件完全相同:.dat / .stl 两种格式任选其一(同时存在时 .dat 优先), 读取方式、去重与方向校正规则均一致;.dat 的区段标记行(END Nodes Coordinates、 NormalAreaElements:、END NormalAreaElements)与地形文件一致,同样与程序兼容。

本目录的 Complex_anomalous.dat / .stl 描述的是起伏地形下的一个复杂三维异常体 (2663 个节点、5322 个三角形单元;范围 x ≈ -302 ~ 248 m、y ≈ -197 ~ 176 m、 z ≈ 51 ~ 285 m,嵌入地形面附近)。其电导率由 input.dat 第 27 行 tao_abnormal = 4.0 指定(低阻体)。

网格生成建议:用专业前处理软件建立地形面/异常体表面三角形网格后导出, 或选择"导出 → STL"生成 ASCII STL 文件。

3.4 建模注意事项

x、y 方向网格数建议设为奇数,使模型中心(源中心)恰好落在 Yee 网格面中心; 由于磁感应强度 B 定义于网格面中心,dBz/dt 测点应优先布置在 Yee 网格面中心位置, 以保证测点响应与场定义的对应关系。
发射回线边长应为网格尺寸的奇数倍,使回线中心落在 Yee 网格棱边位置,保证源电流棱边与网格棱边严格对位。
采用 Dirichlet 边界(Logic_PML=0)时,网格加密区域(核心均匀区)应覆盖发射源、接收测点与异常体范围,保证上述区域的计算精度,外围大网格用于扩展计算域、减弱零场边界对结果的影响。
发射源的 z 方向位置默认在 NZS+1 处(NZS=NZ/2,即网格中部),由程序自动确定,无需在输入文件中指定。
全航空(TEM)模拟:本程序同样支持全航空场景——将上半区设置为空气,并在源平面以下再多布置若干层空气网格,即可保证发射源与接收测点均处于空气中。

4. CPML 吸收边界

本程序在 Logic_PML=1 时采用 CPML(卷积完美匹配层, Roden & Gedney 2000) 作为吸收边界,在计算区域外围吸收向外传播的电磁场,模拟"无限大地层", 避免边界反射污染晚时响应。

4.1 实现位置

文件	作用
`module/pml-parameters.f90`	参数声明:σ/α/κ 最大值、PML 层数、26 个记忆变量 ψ 数组、卷积系数 `b_e/c_e`、缩放数组 `den_*(=1/κ)`、Hz 的 z 向递归系数 `c_h_zz`
`lib/get-pml-paramters.f90`	构建 x/y/z 六个边界面内 σ/α/κ 的空间分布(多项式由内向外递增,E/H 交错采样)、`den_*` 缩放数组与卷积系数;仅在 `Logic_PML=1` 时调用
`lib/Iteration.f90`(子程序 Iteration_cpml)	场更新主循环内内嵌记忆变量 ψ 的递推与修正项(Ex/Ey/Ez 与 Hx/Hy 共 24 个 ψ);Hz 的 z 方向采用基于 `c_h_zz` 的递归卷积(非 ψ)

4.2 参数与含义

参数	默认值	含义
`PML_X,PML_Y,PML_Z`	input.dat 第 5 行(建议 ≥5 层)	三个方向的 PML 层数
`ma`	3	σ 沿厚度方向的多项式阶数(由内向外幂律增长)
`mb`	1	α 沿厚度方向的多项式阶数
`sig_max`	1.0e2	PML 外侧最大电导率(决定吸收强度)
`alpha_max`	1.0e-1	复频移因子最大值
`kappa_max`	1.0	坐标拉伸系数最大值(1 表示不拉伸)

σ、α 沿厚度的空间分布(以 x 方向下层为例,其余边界对称):

σ(i) = sig_max · ((L−i)/(L−1))^ma
α(i) = alpha_max · ((i−1)/(L−1))^mb

E 场采样在整层、H 场采样在半层(交错),因此 H 方向的 σ/α/κ 按半层偏移 (i−0.5)构造,与 E 方向错开。

4.3 使用说明

input.dat 第 4 行 Logic_PML=1,第 5 行给出三个方向的 PML 层数 (如 15,15,15);PML 层内网格尺寸应与核心区一致(保持均匀)。
α_max 取 0.1 是关键调参:α(复频移因子)负责吸收低频扩散场。若取值过小(如 0.01),晚时段的低频反射场不能及时衰减,会在边界往返叠加,导致关断后约 10⁻⁵ s 量级出现指数发散(结果为 NaN)——这是 CPML 版最常见的不稳定来源,务必保持 α_max=1.0e-1。
切换回原始边界:第 4 行改为 0 即可,行为与旧版本完全一致,无需重新编译。

4.4 与 Dirichlet 边界的对比

	CPML(`Logic_PML=1`)	Dirichlet(`Logic_PML=0`)
边界处理	吸收层,模拟无界空间	边界处场直接为零
晚时精度	吸收反射,衰减曲线平直	边界反射可能污染晚时响应
计算量	每步多 24 个记忆变量递推(约 +30%)	无额外开销
网格要求	PML 层内需均匀网格	无特殊要求
稳定性	调参正确时稳定	稳定

一致性验证:均匀半空间算例(81×81×80 网格、PML 15 层、1 ms 平台、关断后对比)中,CPML 与 Dirichlet 两种边界在关断后早期(场尚未到达边界时)的响应曲线一致,差异 <0.2%(源自两版循环次序不同导致的浮点舍入累积, 非物理差异),证明 CPML 实现与主迭代等价、正确。

5. 输出文件说明

文件	内容
`dBzdt_1.txt`, `dBzdt_2.txt`, …	每个接收点一个文件。文件头 2 行为说明(测点编号、测点坐标),其后每行 3 列:迭代步数、关断后时间(s)、该时刻磁场响应值
`CTIME_TIXING_UPCOS.DAT`	计算时间序列,每行 3 列:累计时间、时间步长、源电流幅值(波形)
`conductivity.vtk`	模型电导率分布(规则网格 VTK 格式),可用 ParaView/Tecplot 等打开,检查模型是否正确构建
`logfile.log`	运行日志:参数回显、格式选择提示、运行错误等
`fort.5141`	共形网格计算过程的调试输出
`TEM_decay_curve.png`	衰减曲线图(由第 6 节的 `TEM_decay_plot.py` 生成)

运行结束时屏幕会打印各分段的迭代进度、总计算耗时;正常完成后 logfile.log 末尾出现 Computation finished!。

6. 衰减曲线快速成图(TEM_decay_plot.py)

程序目录下的 TEM_decay_plot.py 用于将正演结果 dBzdt_*.txt 快速绘制为 衰减电压曲线图(双对数坐标)。

6.1 使用方式

# 需要 numpy 与 matplotlib
pip install numpy matplotlib

# 在计算输出文件(dBzdt_*.txt)所在目录运行
python TEM_decay_plot.py

脚本自动搜索脚本同目录下所有 dBzdt_*.txt 文件,每个接收点画一条曲线, 测点编号与坐标自动标注在图例中;默认输出高分辨率图片 TEM_decay_curve.png (dpi=600)并弹窗显示。

6.2 主要可调参数(脚本顶部"User parameters"区)

参数	默认值	说明
`file_pattern`	`dBzdt_*.txt`	匹配的结果文件模式
`xmin, xmax`	`1e-6, 1e-1`	横轴(时间, s)显示范围
`ymin, ymax`	`None, None`	纵轴(响应)显示范围,`None` 表示自动
`use_abs`	`True`	`True` 画 \|dBz/dt\|(正响应),`False` 画带符号值;纵轴标签随之切换
`savefig`	`True`	是否保存图片
`save_name`	`TEM_decay_curve.png`	保存文件名
`dpi`	`600`	图片分辨率
`label_fontsize` 等	18/15/15	标签、刻度、图例字号
`linewidth`	`2.5`	曲线线宽

超过 8 条曲线时自动改用顺序蓝色渐变配色,不会循环重复颜色。

曲线取值规则:仅绘制 时间 > 0 且响应值 > 0 的点(双对数坐标下负值/零值无法显示)。

7. 快速上手(本目录自带算例)

本目录自带算例:起伏地形下的复杂三维异常体(模型 101×101×100 网格、网格尺寸 20 m、源边长 500 m、背景 0.01 S/m、异常体为低阻体 4.0 S/m、1 个测点(源中心正下方);异常体网格 2663 节点/5322 单元,地形网格 10039 节点/5426 单元)。

运行步骤:

确认目录下存在:input.dat、Complex_anomalous.dat(及 .stl,描述同一异常体)、 Complex_Terrain.dat(及 .stl,描述同一地形)。两类文件 .dat 与 .stl 同时存在时,程序以 .dat 优先并给出提示。
VS2019 打开 tem3dfdtd.sln → 选择 Release | x64 → 生成。
将生成的 exe 复制到本目录(或把输入文件放入 exe 目录)后运行。

观察屏幕输出,正常流程为:

Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.
Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.
The number of grids in the core area is odd
At least 320M memory is needed!
...
Conformal mesh of terrain is complete!
Ray tracing computation of terrain is complete!
Conformal mesh of terrain is finished
Start conformal processing of the anomalous body
...
Now computing fraction:  1
50 steps have just finished
...

计算完成后检查输出文件 dBzdt_1.txt 、conductivity.vtk。
(可选)运行 python TEM_decay_plot.py 生成衰减曲线图 TEM_decay_curve.png (见第 6 节)。

8. 常见问题

Q1:运行提示 libiomp5md.dll 找不到 OpenMP 运行时库缺失。将 Intel oneAPI 安装目录 bin/libiomp5md.dll 复制到 exe 旁 (或加入 PATH)。

Q2:提示 Both Complex_Terrain.dat and Complex_Terrain.stl exist! ... 两个格式文件都在。程序以 .dat 优先。若想用 STL,请将 .dat 文件移走或改名。

Q3:计算很慢 / 内存不足 减少 NX,NY,NZ 或增大 GridSize;控制 NSTOP;MAX_OFF_TIME 决定实际迭代步数, 程序会以两者中的较小者为准。运行前会打印所需内存估算。

Q4:如何只算均匀半空间(无异常体、无地形)? 半空间模型应包含"空气 + 大地"两部分。将 Complex_anomalous.* 与 Complex_Terrain.* 移走,并在 input.dat 中设 TEMP_II = 2:第 1 块设为上半部分 (空气,电导率如 1e-5),第 2 块设为下半部分(大地,电导率如 0.01),即构成均匀半空间。注意:TEMP_II = 0 时整个模型只填充背景电导率(全空间均匀介质, 不含空气层)。

Q5:接收点坐标怎么写? 坐标是相对回线源中心的局部坐标(单位 m),正负方向与坐标轴一致。

Q6:CPML 吸收边界与原始 Dirichlet 边界怎么选? input.dat 第 4 行开关 Logic_PML:1 启用 CPML 吸收边界(第 5 行 10,10,10 为三个方向的 PML 层数,可自行调整),能有效吸收边界反射,晚时 (大偏移/晚时间)衰减曲线更平直;0 使用原始非均匀网格 Dirichlet(零场) 边界。切换开关无需重新编译。

9. 参考文献

[1] 孙怀凤, 李貅, 李术才, 等. 考虑关断时间的回线源激发TEM三维时域有限差分正演[J]. 地球物理学报, 2013, 56(3): 1049-1064.

[2] 柳尚斌, 李雪峰, 蓝日彦, 等. 瞬变电磁低频近似Maxwell方程的CPML吸收边界及施加方法[J]. 地球物理学报, 2022, 65(4): 1472-1481.

[3] Li X, Zhao Q, Hu S, et al. Introducing complex geometries to Yee cells in FDTD for transient electromagnetic forward modeling[J]. Geophysics, 2025, 91(2): F1-F12.

10.贡献人员

整个项目是在山东大学孙怀凤教授的领导下开展的，主要贡献人员在代码注释、参考文献中已经写明。如需联系请访问 https://faculty.sdu.edu.cn/sun/

除了本开源代码库之外，我们还提供支持GPU高效计算的商业版本软件软件或专用求解器，如果需要，请访问 https://em3d.cn 获取更多信息。

11.声明

如需联系，请使用sunhuaifeng@email.sdu.edu.cn，不要继续使用代码中标注的gmail邮箱了，因为gmail邮箱经常会出现收发问题，谢谢。
代码中列出的tdem.org 由于精力原因暂时无法维护
所有代码均通过git仓库进行维护和发布https://git.em3d.cn/

README.md Unescape Escape

三维瞬变电磁 FDTD 正演程序 (tem3dfdtd)

目录