Add single 3d field evaluation kernel

[jowezarek]

Adds the kernel eval_field_3d_once to psydac.core.field_evaluation_kernels.py and makes use of it in psydac.fem.tensor.py's method eval_field.

The purpose of this kernel is to accelerate the evaluation of a single 3d field at a single point. This is needed repeatedly, for example, when creating a Poincaré plot tracing a single fieldline.

To Do:

• Check for possibly faster implementation (@yguclu)
• Add 1D & 2D kernel
• Report timings

[jowezarek]

The current use of template does not work @yguclu:

@template(name='T1', types=['float[:,:,:]', 'complex[:,:,:]'])
@template(name='T2', types=['float[:]', 'complex[:]'])

I get the error message

ERROR at annotation (semantic) stage
pyccel:
 |error [semantic]: field_evaluation_kernels.py [28,16]| Variable has already been defined with type 'float64'. Type 'complex128' in assignment is incompatible (res=c3 * local_bases_0[i] * local_bases_1[j] * local_bases_2[k])

Locally I will switch back to the float-only case for now. Do you know how to properly make sure that dtype complex can also be handled?

[yguclu]

Related to issue #188

[yguclu]

The current use of template does not work @yguclu:

@template(name='T1', types=['float[:,:,:]', 'complex[:,:,:]'])
@template(name='T2', types=['float[:]', 'complex[:]'])

I get the error message

ERROR at annotation (semantic) stage
pyccel:
 |error [semantic]: field_evaluation_kernels.py [28,16]| Variable has already been defined with type 'float64'. Type 'complex128' in assignment is incompatible (res=c3 * local_bases_0[i] * local_bases_1[j] * local_bases_2[k])

Locally I will switch back to the float-only case for now. Do you know how to properly make sure that dtype complex can also be handled?

The templates should work now, provided that you upgrade Pyccel to version >= 1.11.2

[vcarlier]

The templates should work now, provided that you upgrade Pyccel to version >= 1.11.2

I just updated my Pyccel and the error is still there. Could you try on your machine?

[vcarlier]

@yguclu
I cannot reproduce the error in the tests on my laptop (even when running pytest in parallel), do you have any idea why?

[yguclu]

@yguclu I cannot reproduce the error in the tests on my laptop (even when running pytest in parallel), do you have any idea why?

You did not forget to run the script psydac_accelerate.py, did you?

[yguclu]

@yguclu I cannot reproduce the error in the tests on my laptop (even when running pytest in parallel), do you have any idea why?

You did not forget to run the script psydac_accelerate.py, did you?

I see the problem now. This branch has not been updated with devel, so it is still using the old version of Pyccel.

[yguclu]

@jowezarek Is there a reason for allowing the basis functions to be complex valued?

Instead of this code

@template(name='T1', types=['float[:,:,:]', 'complex[:,:,:]'])
@template(name='T2', types=['float[:]', 'complex[:]'])
def eval_field_3d_once(local_coeffs: 'T1', 
                       local_bases_0: 'T2', local_bases_1: 'T2', local_bases_2: 'T2'):

we can just write

@template(name='T', types=[float, complex])
def eval_field_3d_once(local_coeffs: 'T[:,:,:]', 
                       local_bases_0: 'float[:]', local_bases_1: 'float[:]', local_bases_2: 'float[:]'):

[jowezarek]

@jowezarek Is there a reason for allowing the basis functions to be complex valued?

Instead of this code

@template(name='T1', types=['float[:,:,:]', 'complex[:,:,:]'])
@template(name='T2', types=['float[:]', 'complex[:]'])
def eval_field_3d_once(local_coeffs: 'T1', 
                       local_bases_0: 'T2', local_bases_1: 'T2', local_bases_2: 'T2'):

we can just write

@template(name='T', types=[float, complex])
def eval_field_3d_once(local_coeffs: 'T[:,:,:]', 
                       local_bases_0: 'float[:]', local_bases_1: 'float[:]', local_bases_2: 'float[:]'):

@yguclu, is it possible that your proposed solution was not possible prior to your template-update PR?
Either way, that is exactly what I wanted in the first place. I never intended to define a template for complex coefficients and real basis function values / real coefficients and complex basis function values.

[yguclu]

@jowezarek Is there a reason for allowing the basis functions to be complex valued?
Instead of this code

@template(name='T1', types=['float[:,:,:]', 'complex[:,:,:]'])
@template(name='T2', types=['float[:]', 'complex[:]'])
def eval_field_3d_once(local_coeffs: 'T1', 
                       local_bases_0: 'T2', local_bases_1: 'T2', local_bases_2: 'T2'):

we can just write

@template(name='T', types=[float, complex])
def eval_field_3d_once(local_coeffs: 'T[:,:,:]', 
                       local_bases_0: 'float[:]', local_bases_1: 'float[:]', local_bases_2: 'float[:]'):

@yguclu, is it possible that your proposed solution was not possible prior to your template-update PR? Either way, that is exactly what I wanted in the first place. I never intended to define a template for complex coefficients and real basis function values / real coefficients and complex basis function values.

My PR related to templates makes this function easier to read, but even without it it was possible to write this:

@template(name='T', types=['float[:,:,:]', 'complex[:,:,:]'])
def eval_field_3d_once(local_coeffs: 'T', 
                       local_bases_0: 'float[:]', local_bases_1: 'float[:]', local_bases_2: 'float[:]'):

[vcarlier]

@yguclu everything seems to be fixed now