546 行
49 KiB
Markdown
546 行
49 KiB
Markdown
# tem3dfdtd
|
||
|
||
[中文](./README.md) | [English](./README.en.md)
|
||
|
||
[!NOTE]
|
||
|
||
> All contents of this document are automatically translated by AI without human verification. Please refer to the original Chinese version if you need accurate information.
|
||
|
||
# 3D Transient Electromagnetic FDTD Forward Modeling Program (tem3dfdtd)
|
||
|
||
This program performs forward modeling of three-dimensional transient electromagnetic (TEM) responses based on the **FDTD (finite-difference time-domain) method**. The core algorithm adopts the Du Fort–Frankel method improved by Wang–Hohmann (1993): it iteratively solves the magnetic field H on a Yee grid, and ensures the time stability of the explicit iteration by introducing a **fictitious permittivity**. It uses **conformal grid techniques** to handle undulating terrain and arbitrarily shaped anomalous bodies (described by surface triangular meshes). It supports a rectangular loop source and can compute ground TEM, semi-airborne (SATEM), and airborne (ATEM) modes.
|
||
|
||
The overall framework and the three core techniques are derived from the following works:
|
||
|
||
- **Overall framework and core iterative algorithm** (loop-source TEM three-dimensional finite-difference time-domain forward modeling considering the turn-off time, the Wang–Hohmann improved Du Fort–Frankel method, the fictitious permittivity, and the source waveform including the turn-off time); for the theoretical details, see reference [1];
|
||
- **CPML absorbing boundary** (the CPML absorbing boundary for the low-frequency-approximated Maxwell equations in transient electromagnetics and its application method); for the theoretical details, see reference [2];
|
||
- **Conformal grid technique** (introducing arbitrarily complex-shaped structures into the Yee grid computation via a ray-tracing method); for the theoretical details, see reference [3].
|
||
|
||
The specific principles, formula derivations, and implementation details of each method are given in references [1]–[3]. The authors of these papers are all people who have made outstanding contributions to this open-source project.
|
||
|
||
Code structure: `main.f90` (main program), `module/` (global parameters and modules), `lib/` (functional subroutines).
|
||
|
||
---
|
||
|
||
## Table of Contents
|
||
|
||
- [tem3dfdtd](#tem3dfdtd)
|
||
- [3D Transient Electromagnetic FDTD Forward Modeling Program (tem3dfdtd)](#3d-transient-electromagnetic-fdtd-forward-modeling-program-tem3dfdtd)
|
||
- [Table of Contents](#table-of-contents)
|
||
- [1. Environment and Compilation](#1-environment-and-compilation)
|
||
- [1.1 Environment Requirements](#11-environment-requirements)
|
||
- [1.2 Build and Run with VS2019 (Recommended)](#12-build-and-run-with-vs2019-recommended)
|
||
- [1.3 Command-Line Compilation (Optional, Not Recommended on Windows)](#13-command-line-compilation-optional-not-recommended-on-windows)
|
||
- [1.4 Linux Environment Compilation (makefile)](#14-linux-environment-compilation-makefile)
|
||
- [2. Program Flow and Module Structure](#2-program-flow-and-module-structure)
|
||
- [Module 1: Program Control and Parameter Input](#module-1-program-control-and-parameter-input)
|
||
- [Module 2: Yee Grid Generation](#module-2-yee-grid-generation)
|
||
- [Module 3: Electrical Parameter Construction](#module-3-electrical-parameter-construction)
|
||
- [Module 4: Excitation Source and Time Series](#module-4-excitation-source-and-time-series)
|
||
- [Module 5: FDTD Electromagnetic Field Computation](#module-5-fdtd-electromagnetic-field-computation)
|
||
- [Module 6: Output](#module-6-output)
|
||
- [3. Input File Formats](#3-input-file-formats)
|
||
- [3.1 Parameter Control File input.dat](#31-parameter-control-file-inputdat)
|
||
- [3.2 Terrain Mesh File Complex_Terrain](#32-terrain-mesh-file-complex_terrain)
|
||
- [3.3 Anomalous Body Mesh File Complex_anomalous](#33-anomalous-body-mesh-file-complex_anomalous)
|
||
- [3.4 Modeling Notes](#34-modeling-notes)
|
||
- [4. CPML Absorbing Boundary](#4-cpml-absorbing-boundary)
|
||
- [4.1 Implementation Location](#41-implementation-location)
|
||
- [4.2 Parameters and Their Meanings](#42-parameters-and-their-meanings)
|
||
- [4.3 Usage Notes](#43-usage-notes)
|
||
- [4.4 Comparison with the Dirichlet Boundary](#44-comparison-with-the-dirichlet-boundary)
|
||
- [5. Output File Description](#5-output-file-description)
|
||
- [6. Quick Decay-Curve Plotting (TEM_decay_plot.py)](#6-quick-decay-curve-plotting-tem_decay_plotpy)
|
||
- [6.1 Usage](#61-usage)
|
||
- [6.2 Main Adjustable Parameters (the "User parameters" section at the top of the script)](#62-main-adjustable-parameters-the-user-parameters-section-at-the-top-of-the-script)
|
||
- [7. Quick Start (the example bundled in this directory)](#7-quick-start-the-example-bundled-in-this-directory)
|
||
- [8. FAQ](#8-faq)
|
||
- [9. References](#9-references)
|
||
- [10. Contributors](#10-contributors)
|
||
- [11. Statement](#11-statement)
|
||
|
||
---
|
||
|
||
## 1. Environment and Compilation
|
||
|
||
### 1.1 Environment Requirements
|
||
|
||
| Item | Requirement |
|
||
| ---------------- | -------------------------------------------------------------------------------------------------------------------------------- |
|
||
| Operating system | Windows 10/11 (64-bit) or Linux (x86_64); the domestic operating system deepin 25 is especially recommended |
|
||
| IDE | Windows: Visual Studio 2019 or later; Linux/deepin: make; |
|
||
| Fortran compiler | Intel oneAPI Fortran (ifort/ifx); |
|
||
| Parallel support | OpenMP (multi-core CPU acceleration); additional GPU support is provided in the commercial version, please visit https://em3d.cn |
|
||
|
||
### 1.2 Build and Run with VS2019 (Recommended)
|
||
|
||
1. Install Visual Studio 2019 (select the "Desktop development with C++" workload) and Intel oneAPI (select the "Visual Studio integration for Intel Fortran Compiler" during installation).
|
||
2. Double-click to open the project file **`tem3dfdtd.sln`** (an Intel Fortran project, corresponding to
|
||
`tem3dfdtd.vfproj`). The project already includes all source files. Configuration notes:
|
||
- `Debug | x64` / `Release | x64`: use the **ifx** compiler (oneAPI default), **recommended**;
|
||
- `Debug | Win32` / `Release | Win32`: use the ifort compiler.
|
||
3. Select the **`Release | x64`** configuration and click **Build → Build Solution**.
|
||
4. Preparation before running: the program looks for `input.dat` and the mesh files in the working directory, so place
|
||
`input.dat`, `Complex_Terrain.dat/.stl`, and `Complex_anomalous.dat/.stl` in the
|
||
**tem3dfdtd-open\tem3dfdtd** directory (or specify it via "Project Properties → Debugging → Working Directory").
|
||
5. Run `tem3dfdtd.exe` directly, or press F5 to debug-run.
|
||
|
||
> **Note**: the program uses OpenMP and requires Intel's runtime library `libiomp5md.dll` at run time
|
||
> (located in the Intel oneAPI installation directory `bin/`). If you get a missing-DLL error, copy it next to
|
||
> the exe (a copy is already placed in this directory), or add its directory to the system PATH.
|
||
|
||
### 1.3 Command-Line Compilation (Optional, Not Recommended on Windows)
|
||
|
||
In an "Intel oneAPI Command Prompt" environment:
|
||
|
||
```bat
|
||
ifx -c -O2 -Qopenmp module\*.f90
|
||
ifx -c -O2 -Qopenmp -Qopenmp lib\*.f90 main.f90
|
||
ifx -O2 -Qopenmp *.obj -o tem3dfdtd.exe
|
||
```
|
||
|
||
(Linking requires MSVC's link.exe and the Windows SDK libraries; it is recommended to use VS's "Developer Command Prompt"
|
||
and add Intel oneAPI's `bin` directory to PATH.)
|
||
|
||
### 1.4 Linux Environment Compilation (makefile)
|
||
|
||
To compile under Linux you need to write a `makefile`. Please write the makefile yourself according to your operating system's configuration and requirements, and compile with `make` (using the makefile).
|
||
|
||
**Notes and cautions**
|
||
|
||
- Differences from the Windows version: under Linux the executable name is custom; if not specified it defaults to `main.exe`. The working directory must likewise contain
|
||
`input.dat` and the `Complex_Terrain.*`, `Complex_anomalous.*` mesh files.
|
||
|
||
---
|
||
|
||
## 2. Program Flow and Module Structure
|
||
|
||
Main program execution order (`main.f90`):
|
||
|
||
```
|
||
GETDATA → CHECKPARAMETERS → MEMORY_USE_ESTIMATION → ALLOCATEMEMORY
|
||
→ GET_NON_UNIFORMGRID → ZERO → GET_COORDINATES
|
||
→ Get_Receiver_Gridlabel → RES_CONFIGURE → TIME_SERIOUS
|
||
→ Get_eps_r → (if Logic_PML=1) Get_pml_parameters → Get_mstop
|
||
→ GetSourcePosition → Iteration → FREE_MEMORY
|
||
```
|
||
|
||
By function, it is divided into the following 6 modules:
|
||
|
||
### Module 1: Program Control and Parameter Input
|
||
|
||
| File | Function |
|
||
| ------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------ |
|
||
| `main.f90` | Main program, controls the whole computation flow |
|
||
| `lib/getdata.f90` | Reads the parameter control file `input.dat`; detects the existence of terrain/anomalous-body mesh files and selects the read format |
|
||
| `lib/checkparameters.f90` | Echoes the read-in computation parameters to `logfile.log` for manual checking |
|
||
| `lib/memory-use-estimation.f90` | Estimates the required memory based on the grid size and prints a hint |
|
||
|
||
### Module 2: Yee Grid Generation
|
||
|
||
| File | Function |
|
||
| ----------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `lib/allocatememory.f90` | Dynamically allocates all global arrays according to the input parameters (including the CPML memory-variable arrays, only when `Logic_PML=1`) |
|
||
| `lib/get_non_uniformgrid.f90` | Generates the non-uniform grid in the x/y/z directions (uniform in the core area + expanded outward by a factor of 1.3) |
|
||
| `lib/get_coordinates.f90` | Computes the coordinates of all grid nodes (including Yee nodes); coordinates take the source center as origin |
|
||
| `lib/zero.f90` | Initializes all electromagnetic field arrays to 0; sets den\_\* to 1, c_h_zz to 0, and clears the CPML memory variables |
|
||
|
||
### Module 3: Electrical Parameter Construction
|
||
|
||
| File | Function |
|
||
| ----------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `lib/resistivity-configuration.f90` | Builds the model conductivity: without terrain, assigns background conductivity + block-shaped anomalous bodies; with terrain, calls the conformal mesh; finally distributes the conductivity onto edges in the x/y/z directions and writes `conductivity.vtk` |
|
||
| `lib/Terrain_conformal.f90` | **Terrain conformal mesh**: reads the terrain triangular mesh from `Complex_Terrain.dat/.stl`, fills the equivalent conductivity of each edge along the x/y/z directions using ray–triangle intersection (Möller–Trumbore algorithm), handling undulating terrain and the air/stratum interface |
|
||
| `lib/Anomalous_conformal.f90` | **Anomalous-body conformal mesh**: reads the anomalous-body surface mesh from `Complex_anomalous.dat/.stl`, and fills the anomalous-body conductivity (`tao_abnormal`) onto edges using the same ray-intersection method |
|
||
|
||
> **Note**: when a terrain file exists, in the terrain branch the air conductivity is `AIR_CONDUCTIVITY = 1e-6 S/m`,
|
||
> and the stratum conductivity is the **2nd conductivity in the block-shaped anomalous body parameters** `TAR_CONDUCTIVITY(2)` (see the second group of data from line 9 of input.dat). In this case the block-shaped anomalous bodies in input.dat do not take effect directly; instead the anomalous body is described by the `Complex_anomalous`
|
||
> mesh file and assigned its conductivity via `tao_abnormal`.
|
||
|
||
### Module 4: Excitation Source and Time Series
|
||
|
||
| File | Function |
|
||
| ----------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `lib/time-serious.f90` | Generates the time series for the whole computation (including the source waveform), and corrects the iteration step count `NSTOP` according to `MAX_OFF_TIME`; writes `CTIME_TIXING_UPCOS.DAT` |
|
||
| `lib/tixing-source-upcos.f90` | Source waveform: **trapezoidal + cosine-ramp** turn-off current waveform (commonly used, `SOURCE_TYPE = 'TIXING_UPCOS'`) |
|
||
| `lib/tixing-source.f90` | Pure trapezoidal waveform source (`TIXING_RAMP`) |
|
||
| `lib/sin-source.f90` | Half-sine waveform source (`HALF_SIN`) |
|
||
| `lib/triangle-source.f90` | Triangular waveform source (`TRIANGLE`) |
|
||
| `lib/get-eps-r.f90` | Computes the fictitious permittivity `EPS_R = 3·(Δt/Δx)²/μ₀` and the iteration coefficients, ensuring the stability of the explicit FDTD |
|
||
| `lib/get-mstop.f90` | Splits the total iteration into several "computation segments", each with its own cache allocation, for memory management |
|
||
|
||
### Module 5: FDTD Electromagnetic Field Computation
|
||
|
||
| File | Function |
|
||
| --------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `lib/GetSourcePosition.f90` | Determines the loop source's position in the grid and marks the source edges (the current assignment region) |
|
||
| `lib/Iteration.f90` | **Core iteration subroutine**: advances the time steps segment by segment, updates the Ex/Ey/Ez and Hz fields, and loads current according to the source waveform; at the end of each segment it computes Hz for each receiver point (obtained by weighted differencing of the Ex, Ey of the 8 nodes surrounding that point) and writes it to the result file. The boundary condition switches via `Logic_PML`: `1` adds the CPML memory-variable correction in each field update, `0` restores the original Dirichlet (zero-field) boundary |
|
||
| `module/pml-parameters.f90` | **CPML absorbing boundary module** (Roden–Gedney convolutional PML): declares the σ/α/κ polynomial scaling parameters (ma=3, mb=1), the 26 memory-variable ψ arrays, the b*e/c_e convolution-coefficient arrays, and the den*\*(=1/κ) scaling arrays |
|
||
| `lib/get-pml-paramters.f90` | Builds the σ/α/κ distributions on the six x/y/z boundary faces (polynomial decay from the boundary inward, E/H staggered sampling) and the den\_\* scaling arrays; called by main only when `Logic_PML=1` |
|
||
|
||
### Module 6: Output
|
||
|
||
| File | Function |
|
||
| ------------------------------------------------- | --------------------------------------------------------------------------------- |
|
||
| `lib/Iteration.f90` (output part) | Writes the response files `dBzdt_1.txt`, `dBzdt_2.txt`, … for each receiver point |
|
||
| `lib/resistivity-configuration.f90` (output part) | Writes the model conductivity distribution `conductivity.vtk` |
|
||
| `lib/time-serious.f90` (output part) | Writes the time series `CTIME_TIXING_UPCOS.DAT` |
|
||
| `lib/free-memory.f90` | Releases all dynamic memory after computation |
|
||
|
||
---
|
||
|
||
## 3. Input File Formats
|
||
|
||
The program requires the following files (all placed in the exe's working directory):
|
||
|
||
| File | Required? | Description |
|
||
| --------------------------------- | --------------------- | ----------------------------------------------------------------------------------------------- |
|
||
| `input.dat` | **Required** | Computation parameter control file |
|
||
| `Complex_anomalous.dat` or `.stl` | Optional (either one) | Anomalous-body surface triangular mesh; if absent, the model is treated as a uniform background |
|
||
| `Complex_Terrain.dat` or `.stl` | Optional (either one) | Terrain surface triangular mesh; if absent, no undulating terrain is considered |
|
||
|
||
### 3.1 Parameter Control File input.dat
|
||
|
||
**Free-format reading**, read line by line in order; a value may be followed by a `!` comment (whole-line or end-of-line comments). The following explains each line using the `input.dat` bundled in this directory as an example:
|
||
|
||
| Line | Example | Meaning |
|
||
| ----- | -------------- | ----------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| 1 | `1` | Computation mode `CAL_TYPE`: `1` = ground TEM, `2` = semi-airborne (SATEM) |
|
||
| 2 | `500` | Transmitting loop side length `SourceLength` (m) |
|
||
| 3 | `101,101,100` | Number of grids in the x, y, z directions `NX,NY,NZ` |
|
||
| 4 | `1` | **Boundary condition switch `Logic_PML`**: `1` = CPML absorbing boundary, `0` = original non-uniform-grid Dirichlet (zero-field) boundary |
|
||
| 5 | `10,10,10` | **PML layer count `PML_X,PML_Y,PML_Z`** (x, y, z directions; effective only when the switch = 1; ≥ 5 layers recommended) |
|
||
| 6 | `25,25` | x direction: core uniform-grid start/end interval indices `UniGridNumX1,UniGridNumX2` |
|
||
| 7 | `25,25` | y direction: core uniform-grid interval indices `UniGridNumY1,UniGridNumY2` |
|
||
| 8 | `20,30` | z direction: core uniform-grid interval indices `UniGridNumZ1,UniGridNumZ2` |
|
||
| 9 | `20` | Core-area uniform grid size `GridSize` (m) |
|
||
| 10 | `0.01` | Background medium conductivity `BACKGROUND_CONDUCTIVITY` (S/m) |
|
||
| 11 | `2` | Number of block-shaped anomalous bodies `TEMP_II` (filled as prisms when no terrain; set 0 for a uniform model) |
|
||
| 12–15 | see below | **Parameters of the 1st block-shaped anomalous body**, 4 lines in total |
|
||
| 16–19 | see below | **Parameters of the 2nd block-shaped anomalous body**, 4 lines in total |
|
||
| 20 | `4000000` | Maximum number of iterations `NSTOP` |
|
||
| 21 | `90.101` | Maximum computation time `MAX_OFF_TIME` (**unit: ms**) |
|
||
| 22 | `1e-6,1e-9` | Rising-edge duration and time step `RAISETIME, RAISESTEP` (s) |
|
||
| 23 | `60000e-6` | Plateau-stage duration `WAVE` (s, i.e. 60 ms) |
|
||
| 24 | `1e-7,1e-9` | Falling-edge duration and time step `RAMP, RAMPSTEP` (s) |
|
||
| 25 | `1e-9` | Initial time step `TIMESTEP` (s) |
|
||
| 26 | `1` | Transmitting current amplitude `AMP` (A) |
|
||
| 27 | `4.0` | Anomalous-body conductivity `tao_abnormal` (S/m, used with the `Complex_anomalous` file) |
|
||
| 28 | `TIXING_UPCOS` | Source type `SOURCE_TYPE`: `TIXING_UPCOS` / `TIXING_RAMP` / `HALF_SIN` / `TRIANGLE` |
|
||
| 29 | `1` | Number of receiver points `Point_Num` |
|
||
| 30 | `1` | Index of the 1st receiver point |
|
||
| 31 | `0,0,0` | Coordinates of the 1st receiver point (relative to the source center, m) |
|
||
| 32–33 | `2` / `x,y,z` | (if there are extra measurement points) the 2nd receiver point (index + coordinates) |
|
||
|
||
Each **block-shaped anomalous body** consists of 4 consecutive lines:
|
||
|
||
| Line | Example | Meaning |
|
||
| ---- | ------- | --------------------------------------------------- |
|
||
| ① | `1,101` | x-direction grid start/end indices `TAR_X1, TAR_X2` |
|
||
| ② | `1,101` | y-direction grid start/end indices `TAR_Y1, TAR_Y2` |
|
||
| ③ | `1,50` | z-direction grid start/end indices `TAR_Z1, TAR_Z2` |
|
||
| ④ | `1e-5` | Conductivity of this block `TAR_CONDUCTIVITY` (S/m) |
|
||
|
||
> This example uses two "anomalous bodies" to assemble a half-space: block 1 z=1~50 (air, 1e-5 S/m) + block 2
|
||
> z=51~100 (stratum, 1e-2 S/m). After the receiver count `Point_Num`, each measurement point occupies 2 lines
|
||
> (index + coordinates relative to the source center).
|
||
|
||
### 3.2 Terrain Mesh File Complex_Terrain
|
||
|
||
The terrain is described by a **surface triangular
|
||
mesh**, supported in two formats; keep **only one of them** in the folder;
|
||
if both exist at the same time, the program takes `.dat` as the priority and warns that `.stl` is ignored.
|
||
|
||
**Format 1: `Complex_Terrain.dat` (original text format)**
|
||
|
||
```
|
||
Number of Nodes and Elements:
|
||
10039 ← total number of nodes n_point
|
||
5426 ← total number of triangular elements n_face
|
||
Nodes Coordinates:
|
||
1 -21000.0 -21000.0 224.08 ← n_point lines: node index, X, Y, Z
|
||
2 -21000.0 -20001.8 224.08
|
||
... (a `!` comment may appear in a line)
|
||
END Nodes Coordinates
|
||
NormalAreaElements:
|
||
1 1 2 10039 ← n_face lines: element index, node1, node2, node3
|
||
...
|
||
END NormalAreaElements
|
||
```
|
||
|
||
| Line | Content |
|
||
| -------------------------- | -------------------------------------------------------------------------------------------------------------- |
|
||
| Line 1 | Title line, arbitrary |
|
||
| Line 2 | Total number of nodes `n_point` |
|
||
| Line 3 | Total number of triangular elements `n_face` |
|
||
| Line 4 | Title line, arbitrary |
|
||
| Lines 5 ~ 4+n_point | One line per node: `node index, X, Y, Z` |
|
||
| The following 1 line | Section end marker `END Nodes Coordinates` (skipped by the program as a title line) |
|
||
| The following 1 line | Face-section title `NormalAreaElements:` (skipped by the program as a title line) |
|
||
| The following n_face lines | One line per element: `element index, node1 index, node2 index, node3 index` (nodes in counterclockwise order) |
|
||
| The last 1 line | End marker `END NormalAreaElements` (not read by the program) |
|
||
|
||
> The two section-marker lines after the node section (GiD export) are in the same position as the "2 title lines" of the old format;
|
||
> the program always skips them as title lines, so both notations are compatible.
|
||
|
||
**Format 2: `Complex_Terrain.stl` (ASCII STL format)**
|
||
|
||
The STL file format is a file format used to describe the surface geometry of three-dimensional objects, widely used in rapid prototyping, 3D printing, and computer-aided manufacturing (CAM).
|
||
An STL file subdivides an object's surface into a series of small triangles, each defined by a normal vector and three vertex coordinates.
|
||
STL files have two formats: text format (ASCII) and binary format.
|
||
|
||
Standard ASCII STL uses the `facet`/`endfacet` element keywords,
|
||
with nodes represented by `vertex` lines, for example in the following format:
|
||
|
||
```
|
||
facet normal nx ny nz
|
||
outer loop
|
||
vertex v1x v1y v1z
|
||
vertex v2x v2y v2z
|
||
vertex v3x v3y v3z
|
||
endloop
|
||
endfacet
|
||
```
|
||
|
||
Users can create and edit STL files with common CAD software (such as AutoCAD, Blender, FreeCAD, MeshLab, SketchUp, Gid, Maya, 3ds Max, etc.).
|
||
|
||
The program automatically handles two things when reading (no user action needed):
|
||
|
||
1. **Vertex deduplication/merging**: in STL, each face writes its vertices independently; duplicate vertices (tolerance 1e-5) are automatically merged into a unique node table;
|
||
2. **Orientation correction**: compare each face's cross-product direction with the `facet normal` in the file; if opposite, swap the face's 2nd and 3rd
|
||
nodes to guarantee that the normal-direction convention is consistent with the `.dat` format.
|
||
|
||
### 3.3 Anomalous Body Mesh File Complex_anomalous
|
||
|
||
Exactly the same as the terrain file: either the `.dat` or `.stl` format (when both exist, `.dat` takes priority),
|
||
with the same reading, deduplication, and orientation-correction rules; the `.dat` section marker lines (`END Nodes Coordinates`,
|
||
`NormalAreaElements:`, `END NormalAreaElements`) are the same as the terrain file and are likewise compatible with the program.
|
||
|
||
The `Complex_anomalous.dat` / `.stl` in this directory describe **a complex three-dimensional anomalous body under undulating terrain**
|
||
(2663 nodes, 5322 triangular elements; extent x ≈ -302 ~ 248 m, y ≈ -197 ~ 176 m,
|
||
z ≈ 51 ~ 285 m, embedded near the terrain surface).
|
||
Its conductivity is specified by line 27 of `input.dat`,
|
||
`tao_abnormal = 4.0` (a low-resistivity body).
|
||
|
||
> Mesh generation suggestion: use professional preprocessing software to build the terrain surface / anomalous-body surface triangular mesh and export it,
|
||
> or choose "Export → STL" to generate an ASCII STL file.
|
||
|
||
### 3.4 Modeling Notes
|
||
|
||
1. **The number of grids in the x and y directions is recommended to be odd**, so that the model center (the source center) falls exactly at the center of a Yee grid face;
|
||
since the magnetic induction intensity B is defined at the grid-face center, the `dBz/dt` measurement point should **preferably be placed at the center of a Yee grid face**,
|
||
to guarantee the correspondence between the measurement-point response and the field definition.
|
||
2. **The transmitting loop side length should be an odd multiple of the grid size**, so that the loop center falls on a Yee grid edge, ensuring strict alignment between the source-current edges and the grid edges.
|
||
3. **When using the Dirichlet boundary (`Logic_PML=0`)**, the grid refinement region (core uniform area) should **cover the transmitting source,
|
||
the receiver points, and the anomalous body**, to guarantee the computational accuracy of the above region; the outer large grids are used to extend the computational domain and reduce
|
||
the influence of the zero-field boundary on the result.
|
||
4. **The z-direction position of the transmitting source is by default at `NZS+1`** (`NZS=NZ/2`, i.e. the middle of the grid), determined automatically by the program,
|
||
and does not need to be specified in the input file.
|
||
5. **Airborne (TEM) simulation**: the program also supports airborne scenarios — set the upper region to air, and
|
||
add several more layers of air grids below the source plane, so that both the transmitting source and the receiver points are in the air.
|
||
|
||
---
|
||
|
||
## 4. CPML Absorbing Boundary
|
||
|
||
When `Logic_PML=1`, the program uses **CPML (Convolutional Perfectly Matched Layer, Roden & Gedney 2000)**
|
||
as the absorbing boundary, absorbing the outward-propagating electromagnetic field at the periphery of the computational region to simulate an "infinite earth stratum"
|
||
and avoid boundary reflections contaminating the late-time response.
|
||
|
||
### 4.1 Implementation Location
|
||
|
||
| File | Function |
|
||
| ----------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `module/pml-parameters.f90` | Parameter declaration: σ/α/κ maxima, PML layer count, the 26 memory-variable ψ arrays, the convolution coefficients `b_e/c_e`, the scaling arrays `den_*(=1/κ)`, and the z-direction recursive coefficient `c_h_zz` for Hz |
|
||
| `lib/get-pml-paramters.f90` | Builds the spatial distributions of σ/α/κ on the six x/y/z boundary faces (polynomial increasing from inside to outside, E/H staggered sampling), the `den_*` scaling arrays, and the convolution coefficients; called only when `Logic_PML=1` |
|
||
| `lib/Iteration.f90` (subroutine Iteration_cpml) | In the field-update main loop, **embeds** the recurrence and correction terms of the memory variables ψ (24 ψ for Ex/Ey/Ez and Hx/Hy); the z direction of Hz uses a recursive convolution based on `c_h_zz` (not ψ) |
|
||
|
||
### 4.2 Parameters and Their Meanings
|
||
|
||
| Parameter | Default value | Meaning |
|
||
| ------------------- | ------------------------------------------- | --------------------------------------------------------------------------------------------- |
|
||
| `PML_X,PML_Y,PML_Z` | line 5 of input.dat (≥5 layers recommended) | PML layer count in the three directions |
|
||
| `ma` | 3 | Polynomial order of σ along the thickness direction (power-law growth from inside to outside) |
|
||
| `mb` | 1 | Polynomial order of α along the thickness direction |
|
||
| `sig_max` | 1.0e2 | Maximum conductivity at the outer side of the PML (determines absorption strength) |
|
||
| `alpha_max` | **1.0e-1** | Maximum value of the complex-frequency-shift factor |
|
||
| `kappa_max` | 1.0 | Maximum value of the coordinate-stretching factor (1 means no stretching) |
|
||
|
||
The spatial distributions of σ and α along the thickness (taking the lower side of the x direction as an example; the other boundaries are symmetric):
|
||
|
||
```
|
||
σ(i) = sig_max · ((L−i)/(L−1))^ma
|
||
α(i) = alpha_max · ((i−1)/(L−1))^mb
|
||
```
|
||
|
||
The E field is sampled over whole layers and the H field over half layers (staggered), so the H-direction σ/α/κ are constructed with a half-layer offset
|
||
(`i−0.5`), staggered with respect to the E direction.
|
||
|
||
### 4.3 Usage Notes
|
||
|
||
1. In `input.dat`, set line 4 `Logic_PML=1`, and give the PML layer counts for the three directions in line 5
|
||
(e.g. `15,15,15`); the grid size inside the PML layers should be consistent with the core area (kept uniform).
|
||
2. **Taking α_max = 0.1 is the key tuning parameter**: α (the complex-frequency-shift factor) is responsible for absorbing the low-frequency diffusion field.
|
||
If it is too small (e.g. 0.01), the low-frequency reflected field in the late-time stage cannot decay in time and will bounce back and forth across the boundary,
|
||
causing exponential divergence on the order of 10⁻⁵ s after turn-off (the result becomes NaN) — this is
|
||
the most common source of instability in the CPML version, so be sure to keep α_max = 1.0e-1.
|
||
3. To switch back to the original boundary: change line 4 to `0`; the behavior is exactly the same as the old version, no recompilation needed.
|
||
|
||
### 4.4 Comparison with the Dirichlet Boundary
|
||
|
||
| | CPML (`Logic_PML=1`) | Dirichlet (`Logic_PML=0`) |
|
||
| ------------------ | ---------------------------------------------------------- | ----------------------------------------------------------- |
|
||
| Boundary handling | Absorbing layer, simulates an unbounded space | Field is zero at the boundary |
|
||
| Late-time accuracy | Absorbs reflections, decay curve is flat | Boundary reflections may contaminate the late-time response |
|
||
| Computational cost | 24 extra memory-variable recurrences per step (about +30%) | No extra overhead |
|
||
| Grid requirement | Uniform grid required inside the PML layers | No special requirement |
|
||
| Stability | Stable when parameters are tuned correctly | Stable |
|
||
|
||
> **Consistency verification**: in a uniform half-space example (81×81×80 grid, 15 PML layers, 1 ms plateau,
|
||
> compared after turn-off), the CPML and Dirichlet boundaries give consistent response curves in the early stage after turn-off (when the field has not yet reached the boundary),
|
||
> with a difference < 0.2% (arising from floating-point rounding accumulation due to the different loop orderings in the two versions,
|
||
> not a physical difference), proving that the CPML implementation is equivalent to and correct with respect to the main iteration.
|
||
|
||
---
|
||
|
||
## 5. Output File Description
|
||
|
||
| File | Content |
|
||
| ------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||
| `dBzdt_1.txt`, `dBzdt_2.txt`, … | One file per receiver point. The first 2 lines of the header are descriptions (point index, point coordinates); each subsequent line has 3 columns: iteration step, time after turn-off (s), and the magnetic-field response value at that time |
|
||
| `CTIME_TIXING_UPCOS.DAT` | Computation time series; each line has 3 columns: cumulative time, time step, source current amplitude (waveform) |
|
||
| `conductivity.vtk` | Model conductivity distribution (regular-grid VTK format), can be opened with ParaView/Tecplot etc. to check that the model was built correctly |
|
||
| `logfile.log` | Run log: parameter echo, format-selection hints, run errors, etc. |
|
||
| `fort.5141` | Debug output of the conformal mesh computation process |
|
||
| `TEM_decay_curve.png` | Decay-curve figure (generated by `TEM_decay_plot.py` in [Section 6](#6-quick-decay-curve-plotting-tem_decay_plotpy)) |
|
||
|
||
At the end of the run the screen prints the iteration progress of each segment and the total computation time; on normal completion, the end of `logfile.log`
|
||
shows `Computation finished!`.
|
||
|
||
---
|
||
|
||
## 6. Quick Decay-Curve Plotting (TEM_decay_plot.py)
|
||
|
||
The `TEM_decay_plot.py` in the program directory is used to quickly plot the forward-modeling results `dBzdt_*.txt` as
|
||
**decay voltage curves** (log-log coordinates).
|
||
|
||
### 6.1 Usage
|
||
|
||
```bash
|
||
# requires numpy and matplotlib
|
||
pip install numpy matplotlib
|
||
|
||
# run in the directory where the computation output files (dBzdt_*.txt) are located
|
||
python TEM_decay_plot.py
|
||
```
|
||
|
||
The script automatically searches for all `dBzdt_*.txt` files in the same directory, draws one curve per receiver point,
|
||
and automatically annotates the point index and coordinates in the legend; by default it outputs a high-resolution image **`TEM_decay_curve.png`**
|
||
(dpi=600) and pops it up.
|
||
|
||
### 6.2 Main Adjustable Parameters (the "User parameters" section at the top of the script)
|
||
|
||
| Parameter | Default value | Description |
|
||
| ---------------------- | --------------------- | ------------------------------------------------------------------------------------------------------------------------- |
|
||
| `file_pattern` | `dBzdt_*.txt` | Pattern for matching result files |
|
||
| `xmin, xmax` | `1e-6, 1e-1` | Display range of the horizontal axis (time, s) |
|
||
| `ymin, ymax` | `None, None` | Display range of the vertical axis (response); `None` means automatic |
|
||
| `use_abs` | `True` | `True` plots \|dBz/dt\| (positive response), `False` plots the signed value; the vertical-axis label switches accordingly |
|
||
| `savefig` | `True` | Whether to save the figure |
|
||
| `save_name` | `TEM_decay_curve.png` | Saved file name |
|
||
| `dpi` | `600` | Image resolution |
|
||
| `label_fontsize`, etc. | 18/15/15 | Font sizes for labels, ticks, and legend |
|
||
| `linewidth` | `2.5` | Curve line width |
|
||
|
||
When there are more than 8 curves, a sequential blue-gradient colormap is automatically used, so colors are not recycled.
|
||
|
||
> Curve value rule: only points with `time > 0` and response value > 0 are plotted (negative/zero values cannot be displayed in log-log coordinates).
|
||
|
||
---
|
||
|
||
## 7. Quick Start (the example bundled in this directory)
|
||
|
||
This directory bundles an example: **a complex three-dimensional anomalous body under undulating terrain** (model 101×101×100 grids, grid size
|
||
20 m, source side length 500 m, background 0.01 S/m, the anomalous body is a low-resistivity body of 4.0 S/m, 1 measurement point (directly below
|
||
the source center); anomalous-body mesh 2663 nodes/5322 elements, terrain mesh 10039 nodes/5426 elements).
|
||
|
||
Run steps:
|
||
|
||
1. Confirm that the directory contains: `input.dat`, `Complex_anomalous.dat` (and `.stl`, describing the same anomalous body),
|
||
`Complex_Terrain.dat` (and `.stl`, describing the same terrain). When the `.dat` and `.stl` of the same kind
|
||
both exist, the program takes `.dat` as the priority and gives a hint.
|
||
2. Open `tem3dfdtd.sln` in VS2019 → select `Release | x64` → build.
|
||
3. Copy the generated exe to this directory (or put the input files in the exe directory) and run.
|
||
4. Observe the screen output; the normal flow is:
|
||
```
|
||
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
|
||
...
|
||
```
|
||
5. After the computation finishes, check the output files `dBzdt_1.txt` and `conductivity.vtk`.
|
||
6. (Optional) Run `python TEM_decay_plot.py` to generate the decay-curve figure `TEM_decay_curve.png`
|
||
(see [Section 6](#6-quick-decay-curve-plotting-tem_decay_plotpy)).
|
||
|
||
---
|
||
|
||
## 8. FAQ
|
||
|
||
**Q1: Running reports `libiomp5md.dll` not found**
|
||
The OpenMP runtime library is missing. Copy Intel oneAPI's `bin/libiomp5md.dll` next to the exe
|
||
(or add it to PATH).
|
||
|
||
**Q2: The prompt `Both Complex_Terrain.dat and Complex_Terrain.stl exist! ...` appears**
|
||
Both format files are present. The program gives `.dat` priority. If you want to use STL, move the `.dat` file away or rename it.
|
||
|
||
**Q3: The computation is very slow / out of memory**
|
||
Reduce `NX,NY,NZ` or increase `GridSize`; control `NSTOP`; `MAX_OFF_TIME` determines the actual number of iteration steps,
|
||
and the program uses the smaller of the two. Before running, the required memory estimate is printed.
|
||
|
||
**Q4: How to compute only a uniform half-space (no anomalous body, no terrain)?**
|
||
A half-space model should contain two parts: "air + earth". Move `Complex_anomalous.*` and
|
||
`Complex_Terrain.*` away, and set `TEMP_II = 2` in `input.dat`: set block 1 to the upper part
|
||
(air, conductivity e.g. `1e-5`), block 2 to the lower part (earth, conductivity e.g. `0.01`), which forms a
|
||
uniform half-space. Note: **when `TEMP_II = 0`, the whole model is filled only with the background conductivity (a whole-space uniform medium,
|
||
without an air layer)**.
|
||
|
||
**Q5: How to write the receiver coordinates?**
|
||
The coordinates are **local coordinates relative to the loop-source center** (unit: m), with positive/negative directions consistent with the coordinate axes.
|
||
|
||
**Q6: How to choose between the CPML absorbing boundary and the original Dirichlet boundary?**
|
||
Line 4 of `input.dat` is the switch `Logic_PML`: `1` enables the CPML absorbing boundary (line 5
|
||
`10,10,10` is the PML layer count in the three directions, which you can adjust), which effectively absorbs boundary reflections and gives a flatter late-time
|
||
(large-offset / late-time) decay curve; `0` uses the original non-uniform-grid Dirichlet (zero-field)
|
||
boundary. Switching does not require recompilation.
|
||
|
||
---
|
||
|
||
## 9. References
|
||
|
||
[1] Sun Huaifeng, Li Xiu, Li Shucai, et al. Three-dimensional FDTD forward modeling of loop-source TEM considering the turn-off time [J]. Chinese Journal of Geophysics, 2013, 56(3): 1049-1064.
|
||
|
||
[2] Liu Shangbin, Li Xuefeng, Lan Riyan, et al. CPML absorbing boundary for the low-frequency-approximated Maxwell equations in transient electromagnetics and its application method [J]. Chinese Journal of Geophysics, 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. Contributors
|
||
|
||
The entire project was carried out under the leadership of Professor Sun Huaifeng of Shandong University. The main contributors are already noted in the code comments and references. To contact us, please visit https://faculty.sdu.edu.cn/sun/
|
||
|
||
In addition to this open-source code repository, we also provide a commercial-version software or dedicated solver that supports efficient GPU computation. If needed, please visit https://em3d.cn for more information.
|
||
|
||
## 11. Statement
|
||
|
||
- For contact, please use sunhuaifeng@email.sdu.edu.cn, and do not continue using the gmail address noted in the code, because that gmail address often has sending/receiving problems. Thank you.
|
||
- The tdem.org listed in the code cannot be maintained for the time being due to limited energy.
|
||
- All code is maintained and released through the git repository https://git.em3d.cn/
|