thliebig / openEMS

openEMS is a free and open-source electromagnetic field solver using the EC-FDTD method.
http://openEMS.de
GNU General Public License v3.0
413 stars 146 forks source link

Semi-Automatic & Fully-Automatic Correctness Tests Needed #106

Open biergaizi opened 1 year ago

biergaizi commented 1 year ago

Recently some developers are interested in exploring different ways to speedup the FDTD engine, examples include my #100 proposal and #105 patch, and @MircoController's experimental CUDA kernel at https://github.com/thliebig/openEMS-Project/discussions/36. To allow making significant changes to the FDTD kernel with confidence, there's a need for semi-automatic and fully-automatic correctness tests.

Also, when creating a package for a system, it's often customary to run the built-in tests that come within a project to ensure the program is running correctly. This is why FreeBSD developer @yurivict created Issue #78 on "please explain how to run tests." However, currently there is no test.

I'm seeing the possibility of using two types of checks:

  1. Comparison against analytical results: The easiest check is to create simulations for some structures with well-known analytical results, such as the cutoff frequency of waveguide, the characteristic impedance of transmission lines, or the radiation pattern of antennas. However, only simple structure can be compared with this method.
  2. Comparison against previous versions: Another type of check is comparing the result of a new result against the result of a known-good historical version. This allows the use of arbitrary simulations as checks without referring to analytical results.

The checks should eventually be fully-automatic, but for now a semi-automatic check is enough - for example, plot the S-parameter / radiation pattern of a device and the known-good result on the same chart on the screen. This allows developers to easily spot a deviation that indicates the existence of bugs. I think many examples from the Tutorial directory can already serve as tests. Eventually, tests should be fully automatic - each data point of the result would be automatically compared with a known-good reference. An error is generated if a deviation beyond reasonable tolerance is detected.

JohannesLau commented 1 year ago

I'm thinking of taking a shot at this. Any preference on python or octave? I'm thinking of running the examples from the tutorial and saving the results as a reference. Then making a script that will run the examples and compare them to the saved reference. I guess there should be some adjustable tolerance on the what count as an identical result.

0xCoto commented 1 year ago

Any preference on python or octave?

I'm personally relying heavily on Python, so that would be much more preferable for me, but that's probably more of a question for folks making engine/performance improvements etc. who wish to test performance updates. I.e., if they've already tested it in Python or Octave and validated the results are virtually the same, I probably wouldn't worry too much about re-running an extensive set of tests.

I guess there should be some adjustable tolerance on the what count as an identical result.

I agree, that had also been my thought process a while back.

biergaizi commented 1 year ago

During my optimization effort, I'm using a simple field dump test for correctness verification. Just find a simulation example, and turn the field dump on (for throughout testing I probably need to dump the entire simulation domain, but for now I'm dumping a 2D plane only as I'm using a PCB example). Then I compare the VTK output of both the upstream reference run and my own fork.

Since my code is not supposed to touch the floating-point arithmetic, only memory access, I'm checking floating-point for exact equivalence, which is okay if you're running the same code on the same system...

import os
import pyvista

def check_shape(dataset_a, dataset_b):
    if dataset_a.n_cells != dataset_b.n_cells:
        return False
    elif dataset_a.n_points != dataset_b.n_points:
        return False
    elif dataset_a.bounds != dataset_b.bounds:
        return False
    else:
        return True 

def check_timestep(upstream_dataset, dev_dataset):
    if not check_shape(upstream_dataset, dev_dataset):
        print("dataset do not have the same shape!")
        exit(1)

    points = 0 
    bad_points = 0 
    for i, val_i in enumerate(upstream_dataset.point_data['E-Field']):
        for j, val_j in enumerate(val_i):
            if upstream_dataset.point_data['E-Field'][i][j] != dev_dataset.point_data['E-Field'][i][j]:
                bad_points += 1
            points += 1

    if bad_points > 0:
        print("%d points in the dataset are not equal!" % bad_points)

    print("... %d points" % points)

files = list(os.scandir("./upstream/sim/dump_0/"))
files.sort(key=lambda file: file.name)
for dataset in files:
    print("checking...", dataset.name)

    upstream_dataset = pyvista.read(dataset.path)
    dev_dataset = pyvista.read("./dev/sim/dump_0/" + dataset.name)
    check_timestep(upstream_dataset, dev_dataset)