|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +""" |
| 4 | +This analysis script checks that there is no spurious charge build-up when a particle is absorbed by an embedded boundary. |
| 5 | +
|
| 6 | +More specifically, this test simulates two particles of oppposite charge that are initialized at |
| 7 | +the same position and then move in opposite directions. The particles are surrounded by a cylindrical |
| 8 | +embedded boundary, and are absorbed when their trajectory intersects this boundary. With an |
| 9 | +electromagnetic solver, this can lead to spurious charge build-up (i.e., div(E)!= rho/epsion_0) |
| 10 | +that remains at the position where particle was absorbed. |
| 11 | +
|
| 12 | +Note that, in this test, there will also be a (non-spurious) component of div(E) that propagates |
| 13 | +along the embedded boundary, due to electromagnetic waves reflecting on this boundary. |
| 14 | +When checking for static, spurious charge build-up, we average div(E) in time to remove this component. |
| 15 | +
|
| 16 | +The test is performed in 2D, 3D and RZ. |
| 17 | +(In 2D, the cylindrical embedded boundary becomes two parallel plates) |
| 18 | +""" |
| 19 | + |
| 20 | +from openpmd_viewer import OpenPMDTimeSeries |
| 21 | + |
| 22 | +ts = OpenPMDTimeSeries("./diags/diag1/") |
| 23 | + |
| 24 | +divE_stacked = ts.iterate( |
| 25 | + lambda iteration: ts.get_field("divE", iteration=iteration)[0] |
| 26 | +) |
| 27 | +start_avg_iter = 25 |
| 28 | +end_avg_iter = 100 |
| 29 | +divE_avg = divE_stacked[start_avg_iter:end_avg_iter].mean(axis=0) |
| 30 | + |
| 31 | +# Adjust the tolerance so that the remaining error due to the propagating |
| 32 | +# div(E) (after averaging) is below this tolerance, but so that any typical |
| 33 | +# spurious charge build-up is above this tolerance. This is dimension-dependent. |
| 34 | +dim = ts.fields_metadata["divE"]["geometry"] |
| 35 | +if dim == "3dcartesian": |
| 36 | + tolerance = 7e-11 |
| 37 | +elif dim == "2dcartesian": |
| 38 | + tolerance = 3.5e-10 |
| 39 | +elif dim == "thetaMode": |
| 40 | + # In RZ: there are issues with divE on axis |
| 41 | + # Set the few cells around the axis to 0 for this test |
| 42 | + divE_avg[13:19] = 0 |
| 43 | + tolerance = 4e-12 |
| 44 | + |
| 45 | + |
| 46 | +def check_tolerance(array, tolerance): |
| 47 | + assert abs(array).max() <= tolerance, ( |
| 48 | + f"Test did not pass: the max error {abs(array).max()} exceeded the tolerance of {tolerance}." |
| 49 | + ) |
| 50 | + print("All elements of are within the tolerance.") |
| 51 | + |
| 52 | + |
| 53 | +check_tolerance(divE_avg, tolerance) |
0 commit comments