ems 0.0.2

Pure Rust FDTD electromagnetic field solver, compatible with openEMS models
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

//! csx plane-wave excitations (`ExcitationType::PlaneWave`, type 10) →
//! [`PlaneWaveSource`]s, installed by [`crate::Simulation`] as
//! [`PlaneWave`](crate::engine::PlaneWave) operators.
//!
//! Mapping rules:
//!
//! * Every primitive of an enabled plane-wave excitation is snapped to the
//!   mesh like a probe box ([`super`] on snapping) and becomes one
//!   total-field box; it must have a volume on the mesh (otherwise it is
//!   ignored with a warning). Several primitives give several independent
//!   sources, so they should not overlap.
//! * `PropDir` is the propagation direction (normalised; a zero vector
//!   ignores the excitation with a warning), the excitation vector
//!   `Excite` the incident E field in V/m at amplitude 1 of the waveform.
//!   A component of `Excite` along `PropDir` is removed (warning); weight
//!   functions are not used for plane waves (warning if set).
//! * The incident wave is `Excite·a(t − Delay − ξ/c)` with `ξ` the distance
//!   along `PropDir` from the box corner the wave enters first (see
//!   [`crate::engine::planewave`]).
//! * The background medium is the material on the box surface
//!   (after rasterisation): the effective edge
//!   materials of the faces that do not lie on the domain boundary must be
//!   one homogeneous, isotropic medium (relative tolerance 1e-4, PEC edges
//!   ignored); its `ε_r`, `μ_r`, `κ` and `σ_m` go into the spec (lossy
//!   backgrounds included). Otherwise compiling fails with
//!   [`Error::Invalid`] naming the plane wave and
//!   the offending edges. Dispersive materials on the surface are not
//!   detected (only their static part is seen).
//! * The auxiliary grid's dispersion is matched at the centre of the
//!   waveform: `f0` of a Gaussian (`fc/2` for `f0 = 0`), `f0` of a
//!   sinusoid, otherwise half the maximum frequency of the settings, or
//!   the operator's default (20 cells per wavelength) without one.

use csx::properties::ExcitationType;
use csx::{ContinuousStructure, PropertyKind};

use super::snap_primitive;
use crate::engine::medium::{interior_faces, surface_medium};
use crate::engine::{EdgeMaterials, Grid, PlaneWaveSpec};
use crate::fdtd::Fdtd;
use crate::waveform::Waveform;
use crate::{Error, Result};

/// A plane-wave excitation of the model (see the [module docs](self)).
#[derive(Debug, Clone, PartialEq)]
pub struct PlaneWaveSource {
    /// Property name.
    pub name: String,
    /// Excitation number.
    pub number: u32,
    /// Delay \[s\].
    pub delay: f64,
    /// Geometry and physics of the engine operator (engine node indices).
    pub spec: PlaneWaveSpec,
}

/// Whether `w` is a csx discretisation warning about plane wave `name`
/// (csx does not discretise type-10 excitations; the compile step does).
pub(super) fn is_csx_skip_warning(w: &str, names: &[&str]) -> bool {
    w.contains("is not discretised")
        && names
            .iter()
            .any(|n| w.starts_with(&format!("property {n:?} primitive #")))
}

/// The matching frequency for the settings' waveform (0: operator default).
fn match_frequency(fdtd: &Fdtd) -> f64 {
    match fdtd.get_waveform() {
        Some(&Waveform::Gaussian { f0, .. }) if f0 > 0.0 => f0,
        Some(&Waveform::Gaussian { fc, .. }) => 0.5 * fc,
        Some(&Waveform::Sinusoid { f0 }) => f0,
        _ => fdtd.get_max_frequency().map_or(0.0, |f| 0.5 * f),
    }
}

/// The plane-wave sources of `model` (see the [module docs](self)).
pub(super) fn plane_wave_sources(
    model: &ContinuousStructure,
    fdtd: &Fdtd,
    warnings: &mut Vec<String>,
) -> Vec<PlaneWaveSource> {
    let mut out = Vec::new();
    for p in model.properties() {
        let PropertyKind::Excitation(ex) = &p.kind else {
            continue;
        };
        if ex.excitation_type != ExcitationType::PlaneWave || !ex.enabled {
            continue;
        }
        let name = &p.name;
        let dir = ex.prop_dir.each_ref().map(|s| s.value());
        let len = (dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]).sqrt();
        if !(len.is_finite() && len > 0.0) {
            warnings.push(format!(
                "plane wave {name:?}: no propagation direction (PropDir {dir:?}); ignored"
            ));
            continue;
        }
        let k = dir.map(|d| d / len);
        let e = ex.excite.each_ref().map(|s| s.value());
        let along = e[0] * k[0] + e[1] * k[1] + e[2] * k[2];
        let pol = [0, 1, 2].map(|a| e[a] - along * k[a]);
        let e_len = (e[0] * e[0] + e[1] * e[1] + e[2] * e[2]).sqrt();
        let p_len = (pol[0] * pol[0] + pol[1] * pol[1] + pol[2] * pol[2]).sqrt();
        if !(p_len.is_finite() && p_len > 1e-12 * e_len.max(f64::MIN_POSITIVE)) {
            warnings.push(format!(
                "plane wave {name:?}: excitation vector {e:?} has no component across \
                 the propagation direction; ignored"
            ));
            continue;
        }
        if along.abs() > 1e-9 * e_len {
            warnings.push(format!(
                "plane wave {name:?}: excitation vector component along the propagation \
                 direction removed ({e:?} → {pol:?})"
            ));
        }
        if ex.weight.iter().any(|w| !w.is_unity()) {
            warnings.push(format!(
                "plane wave {name:?}: weighting functions are not used for plane waves"
            ));
        }
        if p.primitives.is_empty() {
            warnings.push(format!("plane wave {name:?}: no primitive; ignored"));
        }
        for prim in &p.primitives {
            let (a, b) = snap_primitive(&model.grid, prim, "plane wave", name, warnings);
            let start = [0, 1, 2].map(|d| a[d].min(b[d]));
            let stop = [0, 1, 2].map(|d| a[d].max(b[d]));
            if (0..3).any(|d| start[d] == stop[d]) {
                warnings.push(format!(
                    "plane wave {name:?}: primitive has no volume on the mesh; ignored"
                ));
                continue;
            }
            out.push(PlaneWaveSource {
                name: name.clone(),
                number: ex.number,
                delay: ex.delay.value(),
                spec: PlaneWaveSpec::new(start, stop, k, pol)
                    .with_match_frequency(match_frequency(fdtd)),
            });
        }
    }
    out
}

/// Sets the background of every source to the medium on its box surface
/// (see the [module docs](self)).
///
/// # Errors
/// [`Error::Invalid`] if the surface of a box is not homogeneous.
pub(super) fn set_backgrounds(
    sources: &mut [PlaneWaveSource],
    grid: &Grid,
    materials: &EdgeMaterials,
) -> Result<()> {
    for s in sources {
        let (lo, hi) = (s.spec.start, s.spec.stop);
        let faces = interior_faces(grid, lo, hi);
        match surface_medium(grid, materials, lo, hi, faces) {
            Ok(Some(m)) => s.spec = s.spec.clone().with_medium(m),
            Ok(None) => {}
            Err(e) => {
                return Err(Error::Invalid(format!(
                    "plane wave {:?} (box {lo:?}..={hi:?}): {e}; the total-field box must \
                     lie in a homogeneous background, with scatterers strictly inside",
                    s.name
                )));
            }
        }
    }
    Ok(())
}