Benchmarks¶

Several different benchmarks are available to verify computed values of MTE code.

Verification aims of the benchmarks:

Assumption that the contribution of an element can be represented by a single dipole “far” away from source (see Governing equations).
Shape / orientation of the magnetic field outside source body, verifying rotation of vector components in translation from original facmag FORTRAN77 subroutine [Blakely, 1995], to current function embedded in MTE.
Assumption that volume elements are small enough and distance to P is large enough (see Computational approach).
Solving spatial problem with introduction of a small distance (\(\epsilon\)) suffices.

If the analytical solution takes the form of an equation, its solution will be generated using the function support.compute_analytical_solution(). All benchmarks outlined here do not require modification of the z-coordinate of nodes on the top and bottom surfaces and, consequently, do not call for a setup with additional triangles (see Computational approach). Therefore, unless otherwise specified, the function magnetostatics.compute_B_surface_integral_cuboid() computes the contribution of each element to the magnetic field strength \(\mathbf{B}\) for all benchmarks.

Benchmark 1: a single dipole¶

Analytical solution¶

Equation (3), provides us with a solution for the magnetic field produced by a single dipole [Blakely, 1995, Griffiths, 2013].
At very large distances, a sphere of uniform magnetization can be seen as a single dipole of volume \(V=\frac{4}{3}\pi r^3\), where \(r\) is the distance from \(\mathbf{P}\) to the center of the sphere \(z-z_s\) [Blakely, 1995].
If we align the magnetization of the sphere with one of the axes, \(\mathbf{M}=(0,0,M_{z0})\) and compute values on a line above the center of the sphere, \(\mathbf{r}=(0,0,z-z_s), r=z-z_s\),  eq.  (3) reduces to

(6)¶\[\begin{equation} \mathbf{B}_{dip}(\mathbf{r}) = \frac{\mu_0 V}{2\pi} \frac{M_{z0}}{r^3} \mathbf{e}_z \end{equation}\]

Model setup¶

The model setup is as follows: A spherical inclusion of radius \(ds=1m\) and \(\mathbf{M}= (0,0,7.5)\) in a domain of \(2\times2\times2m\) with a resolution of \(100\) elements in each direction.

The magnetic field values are computed along a vertical line positioned directly above the sphere’s center, where the distance \(r\) is progressively increased to exceed \(100m\).

Results¶

Figure 2 On the left, analytical solution for a single dipole and computed values at increasing distance from surface of a sphere. On the right, the difference between analytical solution for a single dipole and computed values at increasing distance from surface of a sphere. Both show an excellent match of the computed and analytical values, even at close proximity to the surface. At \(0.25m\) the error remains approximately \(\sim \lvert 0.01 \rvert \mu T\).¶

As illustrated in Figure 2, the discrepancies between the analytical solution and computed values are minimal.

Even at a height of \(0.25m\), the smallest height above the topography measured in the Etna case study [Meyer and de Groot, 2024], the error remains approximately \(\sim \lvert 0.01 \rvert \mu T\).

Reproduce¶

Steps to reproduce the results and figures

Please note basic setup in Installation

In MTE.py, modify benchmark attribution to 1:
Code 6 /main/MTE.py¶
49
benchmark = '1'

Run “zoomed” setup & rename/move files

Code 7 /main/MTE.py¶

if benchmark == '1':
    # General settings
    remove_zerotopo = False
    compute_analytical = True
    do_spiral_measurements = False
    do_path_measurements = False

    # Domain settings
    Lx, Ly, Lz = 2, 2, 2  # Length of the domain in each direction.
    nelx, nely, nelz = 100, 100, 100  # Amount of elements in each direction.
    Mx0, My0, Mz0 = 0, 0, 7.5  # Magnetization [A/m], do not change Mx0 and My0.

    # Sphere settings
    sphere_R = 1  # Radius of the sphere.
    sphere_xc, sphere_yc, sphere_zc = Lx / 2, Ly / 2, -Lz / 2  # Center position of the sphere.

    # Line measurement settings
    do_line_measurements = True  # Do computations along a observation line (path).
    line_nmeas = 100  # Amount of observation points.
    xstart, xend = Lx / 2, Lx / 2  # x-coordinates of start and end of observation path.
    ystart, yend = Ly / 2, Ly / 2  # y-coordinates "".
    zstart, zend = 0.01, 2 # Zoomed setup.
    #zstart, zend = 0.01, 100  # Non-zoomed setup.

    # Plane measurement settings
    do_plane_measurements = False  # Do computations on a observation plane.
    plane_nnx, plane_nny = 3, 3  # Amount of observation points in each direction.
    plane_x0, plane_y0, plane_z0 = -Lx / 2, -Ly / 2, 1
    plane_Lx, plane_Ly = 2 * Lx, 2 * Ly  # Length of observation plane in each direction.

Code 8 /main/ (runtime: ~15 min)¶

python3 -u MTE.py | tee log.txt

Code 9 /main/¶

mkdir -p benchmarks/benchmark_1/results_zoom && mv log.txt *.vtu *.ascii $_

Run regular setup & move files

Code 10 /main/MTE.py¶

# Line measurement settings
do_line_measurements = True  # Do computations along a observation line (path).
line_nmeas = 100  # Amount of observation points.
xstart, xend = Lx / 2, Lx / 2  # x-coordinates of start and end of observation path.
ystart, yend = Ly / 2, Ly / 2  # y-coordinates "".
#zstart, zend = 0.01, 2 # Zoomed setup.
zstart, zend = 0.01, 100  # Non-zoomed setup.

Code 11 /main/ (runtime: ~15 min)¶

python3 -u MTE.py | tee log.txt

Code 12 /main/¶

mv log.txt *.vtu *.ascii benchmarks/benchmark_1/

Go to directory & plot
Code 13 /main/¶
```
cd benchmarks/benchmark_1
```
Code 14 /main/benchmarks/benchmark_1/¶

gnuplot plot_script_B1.p
Code 15 /main/benchmarks/benchmark_1/¶

python3 plot_script_B1.py

Benchmark 2: internal cancellation¶

Analytical solution¶

According to theory, all internal magnetic forces, or contributions, on the surfaces within the magnetized object should cancel out [Jackson, 1999]. Hence, regardless of variations on the internal surfaces of elements in our domain, the computed values at any point above the surface should be consistent.

Model setup¶

Figure 3 Visualization of different model setups by cross sectional (z-y) planes trough middle of each domain. On the left the undeformed base domain, in the middle deformation setup (1), where perturbations within a range of \([-0.1:0.1m]\) are introduced to the z-coordinates of the internal nodes of the domain. On the right deformation setup (2), where the same perturbations are repeated, additionally, the aspect ratio of elements are increased.¶

To verify this, a domain of \(10\times10\times10m\), with an initial element size of \(2\times2\times2m\) and \(\mathbf{M}= (0,0,7.5)\), is deformed in two ways:

A random value between \(-0.1\) and \(0.1m\) is added to the z coordinates of internal nodes.
In addition to situation (1), the aspect ratio of elements is significantly increased, with each element’s dimension now \(5\times1\times0.2m\).

The original domain and results will be referred to as the base setup, see Figure 3 for visualization of the deformation setups.

Results¶

Figure 4 The magnetic field strength \(\mathbf{B}\) on a plane \(1m\) above the surface of a box with a spatial extent that is twice as large as the (undeformed) domain beneath. Black lines denote the domain edges, while the grey lines connect the observation points on the plane. The observed pattern reveals the extent of the magnetization of the cuboid domain, rapidly decreasing with distance from the domain. These results provided the base for further testing in this benchmark.¶

Figure 5 Difference between Figure 4 and results from the deformed domain setups. The left side depicts the variance between deformation setup (1) and the base, while the right side shows the difference of deformation setup (2) relative to the base. Notably, the errors in both deformation setups are minimal.¶

As expected and observable in Figure 5, the computed values of the magnetic field on the observation plane, located one meter above the domain, remained consistent (up to machine precision) across these scenarios. However, an increased error is noticeable in deformation setup (1). This might be due to the high aspect ratio, which results in fewer observation points corresponding with the domain’s spatial instabilities, see Computational approach. However, the errors in both experiments are minimal, indicating that internal modifications, which are theoretically expected to neutralize each other, do not significantly impact the results when employing the MTE code.

Reproduce¶

Steps to reproduce the results and figures

Please note basic setup in Installation

In MTE.py, modify benchmark attribution to 2a:
Code 16 /main/MTE.py¶
49
benchmark = '2a'

Run base setup & rename/move files

Code 17 /main/MTE.py¶

if benchmark == '2a':
   # General settings
   remove_zerotopo = False
   compute_analytical = False
   do_spiral_measurements = False
   do_path_measurements = False

   # Domain settings
   Lx, Ly, Lz = 10, 10, 10
   nelx, nely, nelz = 5, 5, 5
   Mx0, My0, Mz0 = 0, 0, 7.5  # Magnetization [A/m].

   # Plane measurement settings
   do_plane_measurements = True
   plane_nnx, plane_nny = 11, 11
   plane_x0, plane_y0, plane_z0 = -Lx / 2, -Ly / 2, 1
   plane_Lx, plane_Ly = 2 * Lx, 2 * Ly

   # Line measurement settings
   do_line_measurements = False
   line_nmeas = 47
   xstart, xend = 0.23 + ((Lx - 50) / 2), 49.19 + ((Ly - 50) / 2)
   ystart, yend = Ly / 2 - 0.221, Ly / 2 - 0.221
   zstart, zend = 1, 1  # 1m above surface.

Code 18 /main/ (runtime: ~3 s)¶

python3 -u MTE.py | tee log.txt

Code 19 /main/¶

mkdir -p benchmarks/benchmark_2/d0 && mv log.txt *.vtu *.ascii $_

In MTE.py, modify benchmark attribution to 2b:
Code 20 /main/MTE.py¶
49
benchmark = '2b'

Run deformation setup (1) & move files

Code 21 /main/MTE.py¶

if benchmark == '2b':
   # General settings
   remove_zerotopo = False
   compute_analytical = False
   do_spiral_measurements = False
   do_path_measurements = False
   do_line_measurements = False

   # Domain settings
   Lx, Ly, Lz = 10, 10, 10
   nelx, nely, nelz = 5, 5, 5
   #nelx, nely, nelz = 2, 10, 50
   Mx0, My0, Mz0 = 0, 0, 7.5
   dz = 0.1  # Amplitude random for perturbations in domain.

   # Plane measurement settings
   do_plane_measurements = True
   plane_nnx, plane_nny = 11, 11
   plane_x0, plane_y0, plane_z0 = -Lx / 2, -Ly / 2, 1
   plane_Lx, plane_Ly = 2 * Lx, 2 * Ly

Code 22 /main/ (runtime: ~5 s)¶

python3 -u MTE.py | tee log.txt

Code 23 /main/¶

mkdir -p benchmarks/benchmark_2/d0_1 && mv log.txt *.vtu *.ascii $_

Run deformation setup (2) & move files

Code 24 /main/MTE.py¶

# Domain settings
Lx, Ly, Lz = 10, 10, 10
#nelx, nely, nelz = 5, 5, 5
nelx, nely, nelz = 2, 10, 50
Mx0, My0, Mz0 = 0, 0, 7.5
dz = 0.1  # Amplitude random for perturbations in domain.

Code 25 /main/ (runtime: ~5 s)¶

python3 -u MTE.py | tee log.txt

Code 26 /main/¶

mkdir -p benchmarks/benchmark_2/d0_1_2_10_50 && mv log.txt *.vtu *.ascii $_

Go to directory & use paraview or plotting to visualize
Code 27 /main/¶
```
cd benchmarks/benchmark_2
```
Code 28 /main/benchmarks/benchmark_2/¶
```
paraview --state=plot_result_b2_final.pvsm
```
Code 29 /main/benchmarks/benchmark_2/¶

gnuplot plot_script_B2.p
Code 30 /main/benchmarks/benchmark_2/¶

python3 plot_script_B2.py

Benchmark 3: a magnetized sphere¶

Analytical solution¶

Using the boundary conditions of a magnetized sphere present in a magnetic field \(\mathbf{B_0}\), equation (5) can be simplified (see appendix Complete derivation of the equations).
This is applicable if the sphere is uniformly magnetized with \(\mathbf{M}\) parallel to \(\hat{k}\), the polar direction and if the origin of the coordinate system is placed at the center of the sphere (see Figure 26 in Complete derivation of the equations for visualization).
Then, the magnetic field outside this sphere is defined as [Reitz and Milford, 1960]

(7)¶\[\begin{equation} \mathbf{B_t(r)} = B_0\mathbf{\hat{k}} + \frac{\mu_{0}}{3}M\left(\frac{a^3}{r^3}\right) \left(2\mathbf{\hat{r}}\cos{\theta}+\mathbf{{\hat{\theta}}}\sin{\theta}\right) \end{equation}\]

where \(r\) is the distance from the center of the sphere to the observation point, \(a\) is the radius of the sphere, \(\mathbf{\hat{r}}\) is the unit vector in the direction of \(r\), \(\mathbf{\hat{\theta}}\) is the unit vector in the direction of \(\theta\), \(\theta\) is the angle between \(\mathbf{\hat{r}}\) and \(\mathbf{\hat{k}}\) increasing clockwise from \(\mathbf{\hat{k}}\) and both \(\mathbf{M}\) and \(\mathbf{B_0}\) are in the direction of \(\mathbf{\hat{k}}\).

Figure 6 Visualization of the model setup with numbering along Fibonacci spiral for uniform distribution above the tessellated sphere. The numbering sequence begins at the top of the sphere and proceeds downward in a counterclockwise spiral. The magnetization is assigned to any elements within the spherical domain, and aligned with the z-axis.¶

Model setup¶

The model setup is as follows, see Figure 6: A spherical inclusion similar to the first benchmark, but now featuring a radius of \(a=10m\) is placed in a domain of \(20\times20\times20m\) with a magnetization of \(\mathbf{M}= (0,0,7.5)\) assigned exclusively to elements within the spherical inclusion. Given that a sphere is a complex shape to accurately represent using hexahedron elements, a large number of elements are anticipated to be necessary to produce adequate results. A Fibonacci spiral is used to uniformly distribute \(101\) computation points at \(0.25m\) and \(0.5m\) above the surface of a sphere with a domain resolution of either \(3\) or \(6\) elements per meter.

Results¶

_images/B3sphere_dif_mp_splitcase_all.png

Figure 7 Comparative analysis of the analytical solution versus computed values at \(101\) distinct observation points, located either \(0.25m\) or \(0.5m\) above a sphere’s surface. The resolutions tested are \(3\) or \(6\) elements per meter (\(el/m\)). It is demonstrated that increasing the number of elements substantially reduces the discrepancy between analytical and computed values.¶

The results are shown in the Figure 7. The findings confirm that proximity to the surface necessitates a higher resolution. A notable feature of this relationship is its scalability; for instance, at twice the altitude, the required resolution can be halved while maintaining accuracy.

At a distance of \(0.5m\) above the sphere’s surface, a resolution of \(3\) elements per meter is found to be adequate, whilst at a distance of \(0.25m\), a resolution of \(6\) elements per meter is required.

Reproduce¶

Steps to reproduce the results and figures

Please note basic setup in Installation

In MTE.py, modify benchmark attribution to 3:
Code 31 /main/MTE.py¶
49
benchmark = '3'

Run 25cm above setup & rename/move files

Code 32 /main/MTE.py¶

if benchmark == '3':
   # General settings
   remove_zerotopo = False
   compute_analytical = True
   do_line_measurements = False
   do_path_measurements = False

   # Domain settings
   Lx, Ly, Lz = 20, 20, 20
   nelx, nely, nelz = 60, 60, 60  # 3 el/m.
   #nelx, nely, nelz = 120, 120, 120  # 6 el/m.
   Mx0, My0, Mz0 = 0, 0, 7.5

   # Sphere settings
   sphere_R = 10  # Do not change, or change radius_spiral as well.
   sphere_xc, sphere_yc, sphere_zc = Lx / 2, Ly / 2, -Lz / 2

   # Spiral measurement settings
   do_spiral_measurements = True
   npts_spiral = 101  # keep odd
   radius_spiral = 1.025 * sphere_R  # 25 cm above surface sphere.
   #radius_spiral = 1.05 * sphere_R  # 50 cm above surface sphere.

   # Plane measurement settings
   do_plane_measurements = False
   plane_nnx, plane_nny = 30, 30
   plane_x0, plane_y0, plane_z0 = -Lx / 2, -Ly / 2, 0.5
   plane_Lx, plane_Ly = 2 * Lx, 2 * Ly

Code 33 /main/ (runtime: ~3 min)¶

python3 -u MTE.py | tee log.txt

Code 34 /main/¶

mkdir -p benchmarks/benchmark_3/0_25_above && mv log.txt *.vtu *.ascii $_

Run 25cm above setup with double amount of elements & rename/move files

Code 35 /main/MTE.py¶

# Domain settings
Lx, Ly, Lz = 20, 20, 20
#nelx, nely, nelz = 60, 60, 60  # 3 el/m.
nelx, nely, nelz = 120, 120, 120  # 6 el/m.
Mx0, My0, Mz0 = 0, 0, 7.5

# Sphere settings
sphere_R = 10  # Do not change, or change radius_spiral as well.
sphere_xc, sphere_yc, sphere_zc = Lx / 2, Ly / 2, -Lz / 2

# Spiral measurement settings
do_spiral_measurements = True
npts_spiral = 101  # keep odd
radius_spiral = 1.025 * sphere_R  # 25 cm above surface sphere.
#radius_spiral = 1.05 * sphere_R  # 50 cm above surface sphere.

Code 36 /main/ (runtime: ~25 min)¶

python3 -u MTE.py | tee log.txt

Code 37 /main/¶

mkdir -p benchmarks/benchmark_3/0_25_2_above && mv log.txt *.vtu *.ascii $_

Run 50cm above setup & rename/move files

Code 38 /main/MTE.py¶

# Domain settings
Lx, Ly, Lz = 20, 20, 20
nelx, nely, nelz = 60, 60, 60  # 3 el/m.
#nelx, nely, nelz = 120, 120, 120  # 6 el/m.
Mx0, My0, Mz0 = 0, 0, 7.5

# Sphere settings
sphere_R = 10  # Do not change, or change radius_spiral as well.
sphere_xc, sphere_yc, sphere_zc = Lx / 2, Ly / 2, -Lz / 2

# Spiral measurement settings
do_spiral_measurements = True
npts_spiral = 101  # keep odd
#radius_spiral = 1.025 * sphere_R  # 25 cm above surface sphere.
radius_spiral = 1.05 * sphere_R  # 50 cm above surface sphere.

Code 39 /main/ (runtime: ~3 min)¶

python3 -u MTE.py | tee log.txt

Code 40 /main/¶

mkdir -p benchmarks/benchmark_3/0_5_above && mv log.txt *.vtu *.ascii $_

Run 50cm above setup with double amount of elements & rename/move files

Code 41 /main/MTE.py¶

# Domain settings
Lx, Ly, Lz = 20, 20, 20
#nelx, nely, nelz = 60, 60, 60  # 3 el/m.
nelx, nely, nelz = 120, 120, 120  # 6 el/m.
Mx0, My0, Mz0 = 0, 0, 7.5

# Sphere settings
sphere_R = 10  # Do not change, or change radius_spiral as well.
sphere_xc, sphere_yc, sphere_zc = Lx / 2, Ly / 2, -Lz / 2

# Spiral measurement settings
do_spiral_measurements = True
npts_spiral = 101  # keep odd
#radius_spiral = 1.025 * sphere_R  # 25 cm above surface sphere.
radius_spiral = 1.05 * sphere_R  # 50 cm above surface sphere.

Code 42 /main/ (runtime: ~25 min)¶

python3 -u MTE.py | tee log.txt

Code 43 /main/¶

mkdir -p benchmarks/benchmark_3/0_5_2_above && mv log.txt *.vtu *.ascii $_

Go to directory & plot
Code 44 /main/¶
```
cd benchmarks/benchmark_3
```
Code 45 /main/benchmarks/benchmark_3/¶

gnuplot plot_script_B3.p
Code 46 /main/benchmarks/benchmark_3/¶

python3 plot_script_B3.py

(OPTIONAL) Use paraview to visualize model setups

Code 47 /main/benchmarks/benchmark_3/¶

tee ./0_5_above/model_setup.pvsm ./0_5_2_above/model_setup.pvsm ./0_25_2_above/model_setup.pvsm ./0_25_above/model_setup.pvsm < ./model_setup.pvsm >/dev/null

Code 48 /main/benchmarks/benchmark_3/¶

paraview --state=0_5_2_above/model_setup.pvsm

paraview --state=0_5_above/model_setup.pvsm

paraview --state=0_25_above/model_setup.pvsm

paraview --state=0_25_2_above/model_setup.pvsm

Benchmark 4: a prismatic body¶

Analytical solution¶

In this benchmark, rather than relying on theoretically derived analytical solutions for specific scenarios, we will assess our model against the numerical outcomes of another study presenting a (similar) closed-form analytical solution for a homogeneous polyhedral magnetic target [Ren et al., 2017, Ren et al., 2019]. The findings of this study have been corroborated against [Philip J. Heath and Direen, 2005] using a prismatic body model setup. Their paper only presented results for the gradient tensor components. However, it also involved calculations for the magnetic field components.

Fortunately, full access to the original code and the magnetic field strengths dataset is provided, (github), enabling us to perform the comparison for this numerical configuration.

Model setup¶

Our model setup is designed to mirror the configuration delineated by [Ren et al., 2019]. Nevertheless, certain parameters are adjusted to align optimally with our coding framework.
The evaluations are conducted on an equilateral cuboid with side lengths of \(L=10m\), centrally positioned at \((0,0,-5m)\). The top face of this cuboid intersects the plane at a height of \(z=0m\). Each element within the cuboid is assigned a homogeneous magnetization of \(\mathbf{M}= (0,0,200) A/m\). \(21\) observation points are equally distributed along a line on the top surface of the prism at \(z = 0m\), with an x-coordinate of \(x = 6m\) and y-coordinates ranging from \([-25:25m]\).
[Ren et al., 2019]’s model setup section described that the prism is discretized into eight tetrahedral elements. Contrarily, their code stated a division into just six elements. Given this disparity, our modeling approach segmented the domain into ten elements, erring on the side of caution. Nonetheless, theoretical understanding posits that the resolution in this configuration is irrelevant, as shown in the parameter section.

Another key distinction to note: the analytical solution proposed by [Ren et al., 2019] is engineered to counteract numerical instabilities occurring when the observation point aligns with an element’s edge. Our model does exhibit this instability (see Computational approach). To closely replicate their observation points, minor adjustments to the positioning of our study’s observation points are done.

It should be noted, that while this model setup does not necessitate modifications to the top or bottom surface (absence of topography), employing magnetostatics.compute_B_surface_integral_cuboid() would suffice. Nonetheless, to validate our proposed solution (see support.shift_observation_points_edge()) for additional singularities on diagonals of domain elements, function magnetostatics.compute_B_surface_integral_wtopo() is still utilized.

Results¶

Figure 8 Comparison of magnetic field components \(\mathbf{B_x}\), \(\mathbf{B_y}\), \(\mathbf{B_z}\) for the prismatic body. As observation site location are displaced from [Ren et al., 2017], the x-axis now refers to index relating to the observation point rather than distance. The computed values of our study match those of [Ren et al., 2019].¶

Figure 8 shows a near-perfect alignment of the two numerical outputs. This observation underscores that either approach is viable. One can either utilize an analytical solution, acknowledging the instability at the domain edges and minimally adjusting observation points or adopt numerically stable solutions tailored for logarithmic singularities as highlighted in [Ren et al., 2017].

Discussion¶

Figure 9 Comparison of magnetic field components \(\mathbf{B_x}\), \(\mathbf{B_y}\), \(\mathbf{B_z}\) for the prismatic body. As observation site location are displaced from [Ren et al., 2017], the x-axis now refers to index relating to the observation point rather than distance. The differences are small compared to the absolute values depicted in the previous graph.¶

When reproducing the benchmark’s results, discrepancies may arise between the computed values from our model and the numerical results from [Ren et al., 2019], as shown in Figure 9. These discrepancies mainly stem from our method of handling observation points at singularities. Our code shifts problematic observation points by introducing a small, random artificial distance to allow computation (see Computational approach). As a result, when observation points fall at these singularities, the comparison between our model’s results and Ren’s numerical solutions can show notable and variable deviations in each run.

While such variations might initially seem to indicate instability in our method, it is important to consider the specific aims of Ren’s numerical benchmark. Ren’s work was specifically designed to test observation points at singularities, presenting a code with significant advancement over previous methods. In the context of our study, however, such singularities are seldom encountered in practice. Our observation points are usually derived from real field measurements or are selected for synthetic topography simulations. Therefore, despite potential variability in these particular cases, our method remains stable and effective for the types of applications and scenarios we focus on.

Most significantly, it is worth noting that these discrepancies are relatively minor when compared to the absolute values, as can be observed in Figure 8. This underscores the overall reliability and applicability of our approach in a broad range of practical scenarios, despite the occasional variations encountered under specific conditions.

Reproduce¶

Steps to reproduce the results and figures

Please note basic setup in Installation

In MTE.py, modify benchmark attribution to 4:
Code 49 /main/MTE.py¶
49
benchmark = '4'

Run setup & rename/move files

Code 50 /main/MTE.py¶

if benchmark == '4':
   # General settings
   remove_zerotopo = False
   compute_analytical = False
   do_plane_measurements = False
   do_spiral_measurements = False
   do_path_measurements = False

   # Domain settings
   Lx, Ly, Lz = 10, 10, 10
   nelx, nely, nelz = int(Lx), int(Ly), 10
   Mx0, My0, Mz0 = 0, 0, 200

   # Line measurement settings
   do_line_measurements = True
   line_nmeas = 21
   xstart, xend = 6, 6
   ystart, yend = -25, 25
   zstart, zend = 0, 0

   # Reading in values from Ren.
   pathfile = 'sites/B.dat'
   with open(pathfile, 'r') as path:
        lines_path = path.readlines()
   BxB4, ByB4, BzB4  = np.zeros((3, len(lines_path)), dtype=np.float64)  # Bx, By, Bz from Ren.
   data = np.array([list(map(float, line.split())) for line in lines_path])
   BxB4, ByB4, BzB4 = data[:, 0], data[:, 1], data[:, 2]

Code 51 /main/ (runtime: ~5 s)¶

python3 -u MTE.py | tee log.txt

Code 52 /main/¶

mv log.txt *.vtu *.ascii benchmarks/benchmark_4/

Go to directory & plot
Code 53 /main/¶
```
cd benchmarks/benchmark_4
```
Code 54 /main/benchmarks/benchmark_4/¶

gnuplot plot_script_B4.p
Code 55 /main/benchmarks/benchmark_4/¶

python3 plot_script_B4.py
Note

If the resulted difference is unsatisfactory, see discussion section and repeat run to randomly generate new observation point for the model.