Custom magnetic fields ignored by the propagation modules

Describe the bug A custom magnetic field defined in a Python class is accepted, but apparently not used by the propagation modules PropagationCK and PropagationBP, resulting in a constant zero field being effectively employed. Though the problem affects both PropagationCK and PropagationBP, it seems easier to demonstrate using the latter as it provides access to the field it is using.

Is the bug related to the physics in CRPropa? No.

To Reproduce

custom_mf.py script where a simple custom magnetic field is defined:

#!/usr/bin/env python
# coding=utf-8

import time
import crpropa
from crpropa import *

class MyUniform(MagneticField):
    def __init__(self, val):
        super().__init__()

        self.val = val

    def getField(self, position):
        return Vector3d(self.val)

if __name__ == "__main__":
    field = MyUniform(1*gauss)

    source = Source()
    source.add(SourcePosition(Vector3d(0, 0, 0)))
    source.add(SourceParticleType(nucleusId(1, 1)))
    source.add(SourceEnergy(1*EeV))

    observer = Observer()
    observer.add(ObserverSurface(Sphere(Vector3d(0, 0, 0), 1*Mpc)))

    prop = PropagationBP(field, 1e-4, 1*pc, 1*Mpc)

    simulation = ModuleList()
    simulation.add(prop)

    simulation.add(observer)
    simulation.setShowProgress(True)

    tstart = time.time()
    simulation.run(source, 10, True)
    tstop = time.time()
    print(f'Run time: {tstop - tstart:.2f} s')

    print(f'Field value (direct): {field.getField(Vector3d(-1, 0, 0))}')
    print(f'Field value (propagation): {prop.getFieldAtPosition(Vector3d(-1, 0, 0), 0)}')

Upon execution it prints:

python3 custom_mf.py 
crpropa::ModuleList: Number of Threads: 8
Run time: 0.00 s
Field value (direct): Vector(0.0001, 0.0001, 0.0001)
Field value (propagation): Vector(0, 0, 0)

The value retrieved from the propagation module is different from the value the field actually has. The run time also indicates the zero field value is actually used. Replacing PropagationBP with PropagationCK does not change the outcome (the last two lines should be commented out in this case as there's no getFieldAtPosition() method in this case).

If the custom field is replaced with a built-in analogue - i.e. L.20 is replaced with field = UniformMagneticField(Vector3d(1*gauss)), - the two field value estimates match and run time become noticeable, indicating the propagation works properly:

python3 custom_mf.py 
crpropa::ModuleList: Number of Threads: 8
Run time: 2.87 s
Field value (direct): Vector(0.0001, 0.0001, 0.0001)
Field value (propagation): Vector(0.0001, 0.0001, 0.0001)

Expected behavior The supplied non-zero custom field should be used by the propagation module.

System (please complete the following information):

Operating system: Ubuntu Linux 20.04
Laptop or computing cluster? Singularity and Docker containers with aforementioned OS running on a laptop
Python version: 3.8
Swig version: 4.0.1
Cmake version: 3.16.3

Additional context The same issue encountered also if the custom field is built as a CRpropa plugin (the swig "-builtin" option needs to be disabled in this case, following #305). For brevity the corresponding code is not attached here, but can provided upon request, if needed. The issue was tested with CRPropa v3.1.7.

Hi @IevgenVovk, thanks for reporting the error and the detailed description. I was able to reproduce the error.

Problem As you have a Python class that inherits from the wrapped C++ MagneticField class containing virtual methods, cross-language polymorphism via directors is needed. Currently, your child Python class MyUniform should override the getField base class function. However, as we don't have enabled directors for the MagneticField class, your overwritten getField function will not be called.

Possible solution By inserting %feature("director") crpropa::MagneticField; above %include "crpropa/magneticField/MagneticField.h" in the2_headers.i file, the director of the MagneticField class is enabled and the bug should be fixed. I have made the change in the branch https://github.com/reichherzerp/CRPropa3/tree/365_custom_bfield. If you want to test the change, you might have to delete the build folder first and recreate it according to the installation instructions, otherwise, the changes might not be applied.

With the following tiny changes of your custom_mf.py script:

#!/usr/bin/env python
# coding=utf-8

import time
import crpropa
from crpropa import *

class MyUniform(MagneticField):
    def __init__(self, val):
        super().__init__()

        self.val = val

    def getField(self, position):
        return Vector3d(self.val)

    def getField(self, position, z):
        return Vector3d(self.val)

if __name__ == "__main__":
    field = MyUniform(1*gauss)

    source = Source()
    source.add(SourcePosition(Vector3d(0, 0, 0)))
    source.add(SourceParticleType(nucleusId(1, 1)))
    source.add(SourceEnergy(1*EeV))

    observer = Observer()
    observer.add(ObserverSurface(Sphere(Vector3d(0, 0, 0), 1*Mpc)))

    prop = PropagationBP(field, 1e-4, 1*pc, 1*Mpc)

    simulation = ModuleList()
    simulation.add(prop)

    simulation.add(observer)
    simulation.setShowProgress(True)

    tstart = time.time()
    simulation.run(source, 10, True)
    tstop = time.time()
    print(f'Run time: {tstop - tstart:.2f} s')

    print(f'Field value (direct): {field.getField(Vector3d(-1, 0, 0), 0)}')
    print(f'Field value (propagation): {prop.getFieldAtPosition(Vector3d(-1, 0, 0), 0)}')

I get the excepted result. I'll check if this fix solves the bug completely before I prepare the pull request, but I wanted to share the current status already now.

Hello @reichherzerp , thanks a lot for checking! Indeed, CRPropa version from your branch does not suffer from this issue. I will then use your solution till this issue lands in the stable CRPropa release.

Great, thanks for testing. I prepared the PR and took the opportunity to introduce the getFieldAtPosition functions for the other propagation modules (PropagationCK and DiffusionSDE). Thanks again for the suggestions!

CRPropa / CRPropa3

Custom magnetic fields ignored by the propagation modules #365