From a8dda508282c5d0d6912a845bad688b762bc1ce4 Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Thu, 27 Jun 2024 15:41:37 +0200 Subject: [PATCH 01/26] proper initialization of Dnom_tunnel --- src/Ionization/IonizationTunnel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Ionization/IonizationTunnel.cpp b/src/Ionization/IonizationTunnel.cpp index 6d396a980..abadcc5bd 100755 --- a/src/Ionization/IonizationTunnel.cpp +++ b/src/Ionization/IonizationTunnel.cpp @@ -53,7 +53,7 @@ void IonizationTunnel::operator()( Particles *particles, unsigned int ipart_min, unsigned int Z, Zp1, newZ, k_times; double TotalIonizPot, E, invE, factorJion, delta, ran_p, Mult, D_sum, P_sum, Pint_tunnel; - vector IonizRate_tunnel( atomic_number_ ), Dnom_tunnel( atomic_number_ ); + vector IonizRate_tunnel( atomic_number_ ), Dnom_tunnel( atomic_number_ , 0.); LocalFields Jion; double factorJion_0 = au_to_mec2 * EC_to_au*EC_to_au * invdt; From 08985b0a8c13176cf5b8f0dbe4e1e6532e0ed07f Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Thu, 27 Jun 2024 17:05:06 +0200 Subject: [PATCH 02/26] Fix bug for species specific diagnostics with AM vecto --- doc/Sphinx/Overview/releases.rst | 9 +++++++ src/Projector/ProjectorAM2OrderV.cpp | 40 ++++++++++++++++++++++------ src/Projector/ProjectorAM2OrderV.h | 4 +-- 3 files changed, 43 insertions(+), 10 deletions(-) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index e271b32c5..8ad9acd9a 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -40,6 +40,15 @@ Ongoing projects * Spectral solvers +---- + +New release +^^^^^^^^^^^^^^^^^^^^^ + +* **Bug fixes** : + + * Fixed a bug which made tunnel ionization wrong in some cases for high atomic numbers. + * Fixed species specific diagnostics in AM geometry with Vectorization. ---- diff --git a/src/Projector/ProjectorAM2OrderV.cpp b/src/Projector/ProjectorAM2OrderV.cpp index 890d37332..37cd5c668 100755 --- a/src/Projector/ProjectorAM2OrderV.cpp +++ b/src/Projector/ProjectorAM2OrderV.cpp @@ -69,10 +69,11 @@ void ProjectorAM2OrderV::currentsAndDensity( ElectroMagnAM *emAM, double * __restrict__ deltaold, std::complex * __restrict__ array_eitheta_old, int npart_total, - int ipart_ref ) + int ipart_ref, + int ispec ) { - currents( emAM, particles, istart, iend, invgf, iold, deltaold, array_eitheta_old, npart_total, ipart_ref ); + currents( emAM, particles, istart, iend, invgf, iold, deltaold, array_eitheta_old, npart_total, ipart_ref, ispec+1 ); //ispec+1 is passed as a marker of diag int ipo = iold[0]; int jpo = iold[1]; @@ -131,6 +132,9 @@ void ProjectorAM2OrderV::currentsAndDensity( ElectroMagnAM *emAM, int iloc0 = ipom2*nprimr_+jpom2; for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { rho = &( *emAM->rho_AM_[imode] )( 0 ); + unsigned int n_species = emAM->rho_AM_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec; + rho = emAM->rho_AM_s [ifield] ? &( * ( emAM->rho_AM_s [ifield] ) )( 0 ) : &( *emAM->rho_AM_ [imode] )( 0 ) ; int iloc = iloc0; for( unsigned int i=0 ; i<5 ; i++ ) { #pragma omp simd @@ -451,7 +455,8 @@ void ProjectorAM2OrderV::currents( ElectroMagnAM *emAM, double * __restrict__ deltaold, std::complex * __restrict__ array_eitheta_old, int npart_total, - int ipart_ref ) + int ipart_ref, + int ispec ) { // ------------------------------------- // Variable declaration & initialization @@ -533,7 +538,13 @@ void ProjectorAM2OrderV::currents( ElectroMagnAM *emAM, int iloc0 = ipom2*nprimr_+jpom2; for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { - Jl = &( *emAM->Jl_[imode] )( 0 ); + if (ispec == 0) { // When diags are not needed ispec is always 0. + Jl = &( *emAM->Jl_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jl_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jl = emAM->Jl_s [ifield] ? &( * ( emAM->Jl_s [ifield] ) )( 0 ) : &( *emAM->Jl_ [imode] )( 0 ) ; + } int iloc = iloc0; for( unsigned int i=1 ; i<5 ; i++ ) { iloc += nprimr_; @@ -552,7 +563,14 @@ void ProjectorAM2OrderV::currents( ElectroMagnAM *emAM, for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { - Jr = &( *emAM->Jr_[imode] )( 0 ); + if (ispec == 0) { // When diags are not needed ispec is always 0. + Jr = &( *emAM->Jr_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jr_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jr = emAM->Jr_s [ifield] ? &( * ( emAM->Jr_s [ifield] ) )( 0 ) : &( *emAM->Jr_ [imode] )( 0 ) ; + } + int iloc = iloc0 + ipom2 + 1; for( unsigned int i=0 ; i<5 ; i++ ) { #pragma omp simd @@ -570,7 +588,13 @@ void ProjectorAM2OrderV::currents( ElectroMagnAM *emAM, } for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { - Jt = &( *emAM->Jt_[imode] )( 0 ); + if (ispec == 0) { // When diags are not needed ispec is always 0. + Jt = &( *emAM->Jt_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jt_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jt = emAM->Jt_s [ifield] ? &( * ( emAM->Jt_s [ifield] ) )( 0 ) : &( *emAM->Jt_ [imode] )( 0 ) ; + } int iloc = iloc0; for( unsigned int i=0 ; i<5 ; i++ ) { #pragma omp simd @@ -592,7 +616,7 @@ void ProjectorAM2OrderV::currents( ElectroMagnAM *emAM, // --------------------------------------------------------------------------------------------------------------------- //! Wrapper for projection // --------------------------------------------------------------------------------------------------------------------- -void ProjectorAM2OrderV::currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool is_spectral, int /*ispec*/, int scell, int ipart_ref ) +void ProjectorAM2OrderV::currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool is_spectral, int ispec, int scell, int ipart_ref ) { if( istart == iend ) { return; //Don't treat empty cells. @@ -626,7 +650,7 @@ void ProjectorAM2OrderV::currentsAndDensityWrapper( ElectroMagn *EMfields, Parti //double *b_Jy = EMfields->Jy_s [ispec] ? &( *EMfields->Jy_s [ispec] )( 0 ) : &( *EMfields->Jy_ )( 0 ) ; //double *b_Jz = EMfields->Jz_s [ispec] ? &( *EMfields->Jz_s [ispec] )( 0 ) : &( *EMfields->Jz_ )( 0 ) ; //double *b_rho = EMfields->rho_s[ispec] ? &( *EMfields->rho_s[ispec] )( 0 ) : &( *EMfields->rho_ )( 0 ) ; - currentsAndDensity( emAM, particles, istart, iend, invgf->data(), iold, delta->data(), array_eitheta_old->data(), invgf->size(), ipart_ref ); + currentsAndDensity( emAM, particles, istart, iend, invgf->data(), iold, delta->data(), array_eitheta_old->data(), invgf->size(), ipart_ref, ispec ); } } diff --git a/src/Projector/ProjectorAM2OrderV.h b/src/Projector/ProjectorAM2OrderV.h index 5320ac33b..7a3f92b87 100755 --- a/src/Projector/ProjectorAM2OrderV.h +++ b/src/Projector/ProjectorAM2OrderV.h @@ -17,10 +17,10 @@ class ProjectorAM2OrderV : public ProjectorAM ~ProjectorAM2OrderV(); //! Project global current densities (EMfields->Jl_/Jr_/Jt_) - void currents(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0 ); + void currents(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0, int ispec =0 ); //! Project global current densities (EMfields->Jl_/Jr_/Jt_/rho), diagFields timestep - void currentsAndDensity(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0 ); + void currentsAndDensity(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0, int ispec = 0 ); //! Project global current charge (EMfields->rho_), frozen & diagFields timestep void basicForComplex( std::complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode ) override final; From 922cc856f259569b7467214386f6cd541f20de7b Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Fri, 28 Jun 2024 10:08:44 +0200 Subject: [PATCH 03/26] fix diag custom functions with python 3.12 --- doc/Sphinx/Overview/releases.rst | 25 ++++++------- src/Diagnostic/Histogram.h | 64 ++++++++++++++++---------------- src/Patch/VectorPatch.cpp | 2 + src/Tools/PyTools.h | 22 +++++++++++ 4 files changed, 66 insertions(+), 47 deletions(-) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 8ad9acd9a..51911dd91 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -16,13 +16,19 @@ Get Smilei You can find older, `unsupported versions here `_ -.. -.. ---- +---- + +.. _latestVersion: + +Changes made in the repository (not released) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +* **Bug fixes** : -.. .. _latestVersion: + * Tunnel ionization was wrong in some cases for high atomic numbers. + * Species-specific diagnostics in AM geometry with vectorization. + * Custom functions in ``ParticleBinning`` with python 3.12 -.. Changes made in the repository (not released) -.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ---- @@ -40,15 +46,6 @@ Ongoing projects * Spectral solvers ----- - -New release -^^^^^^^^^^^^^^^^^^^^^ - -* **Bug fixes** : - - * Fixed a bug which made tunnel ionization wrong in some cases for high atomic numbers. - * Fixed species specific diagnostics in AM geometry with Vectorization. ---- diff --git a/src/Diagnostic/Histogram.h b/src/Diagnostic/Histogram.h index 2b65c4793..49e970869 100755 --- a/src/Diagnostic/Histogram.h +++ b/src/Diagnostic/Histogram.h @@ -578,25 +578,24 @@ class HistogramAxis_user_function : public HistogramAxis private: void calculate_locations( Species *s, double *array, int *index, unsigned int npart, SimWindow * ) { - #pragma omp critical + SMILEI_PY_ACQUIRE_GIL + // Expose particle data as numpy arrays + particleData.resize( npart ); + particleData.set( s->particles ); + // run the function + PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); + particleData.clear(); + // Copy the result to "array" + double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); + for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { - // Expose particle data as numpy arrays - particleData.resize( npart ); - particleData.set( s->particles ); - // run the function - PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); - particleData.clear(); - // Copy the result to "array" - double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); - for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) - { - if( index[ipart]<0 ) { - continue; - } - array[ipart] = arr[ipart]; + if( index[ipart]<0 ) { + continue; } - Py_DECREF( ret ); + array[ipart] = arr[ipart]; } + Py_DECREF( ret ); + SMILEI_PY_RELEASE_GIL }; PyObject *function; @@ -1167,26 +1166,25 @@ class Histogram_user_function : public Histogram void valuate( Species *s, double *array, int *index ) { unsigned int npart = s->getNbrOfParticles(); - #pragma omp critical - { - // Expose particle data as numpy arrays - particleData.resize( npart ); - particleData.set( s->particles ); - // run the function - PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); - particleData.clear(); - // Copy the result to "array" - double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); - for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { - if( index[ipart]<0 ) { - continue; - } - array[ipart] = arr[ipart]; + SMILEI_PY_ACQUIRE_GIL + // Expose particle data as numpy arrays + particleData.resize( npart ); + particleData.set( s->particles ); + // run the function + PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); + particleData.clear(); + // Copy the result to "array" + double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); + for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { + if( index[ipart]<0 ) { + continue; } - Py_DECREF( ret ); + array[ipart] = arr[ipart]; } + Py_DECREF( ret ); + SMILEI_PY_RELEASE_GIL }; - + PyObject *function; ParticleData particleData; }; diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 42f4dd3d8..5639a4656 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -1443,10 +1443,12 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int if( globalDiags[idiag]->theTimeIsNow_ ) { // All patches run + SMILEI_PY_SAVE_MASTER_THREAD #pragma omp for schedule(runtime) for( unsigned int ipatch=0 ; ipatchrun( ( *this )( ipatch ), itime, simWindow ); } + SMILEI_PY_RESTORE_MASTER_THREAD // MPI procs gather the data and compute #pragma omp single smpi->computeGlobalDiags( globalDiags[idiag], itime ); diff --git a/src/Tools/PyTools.h b/src/Tools/PyTools.h index 055875718..1059e0130 100755 --- a/src/Tools/PyTools.h +++ b/src/Tools/PyTools.h @@ -32,6 +32,28 @@ #include #endif +// The following definitions are required for some interaction between python and openmp +// Specifically, it is needed for an openmp loop that includes calls to python +// Usage: +// SMILEI_PY_SAVE_MASTER_THREAD +// #pragma omp for +// for(....) { +// SMILEI_PY_ACQUIRE_GIL +// python calls +// SMILEI_PY_RELEASE_GIL +// } +// SMILEI_PY_RESTORE_MASTER_THREAD +#define SMILEI_PY_SAVE_MASTER_THREAD \ + PyThreadState *_save = NULL; \ + _Pragma("omp master") \ + _save = PyEval_SaveThread(); +#define SMILEI_PY_RESTORE_MASTER_THREAD \ + _Pragma("omp master") \ + PyEval_RestoreThread( _save ); +#define SMILEI_PY_ACQUIRE_GIL PyGILState_STATE _state = PyGILState_Ensure(); +#define SMILEI_PY_RELEASE_GIL PyGILState_Release( _state ); + + //! tools to query python nemlist and get back C++ values and vectors class PyTools From 049f169180acba2960d23ace3e799439e76f0ba7 Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Fri, 28 Jun 2024 16:38:58 +0200 Subject: [PATCH 04/26] Handle AM prescribed fields as Cartesian ones --- doc/Sphinx/Overview/releases.rst | 3 ++- doc/Sphinx/Use/namelist.rst | 5 +++- src/ElectroMagn/ElectroMagnAM.cpp | 38 ------------------------------- src/ElectroMagn/ElectroMagnAM.h | 2 -- 4 files changed, 6 insertions(+), 42 deletions(-) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 51911dd91..f30088fd5 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -27,7 +27,8 @@ Changes made in the repository (not released) * Tunnel ionization was wrong in some cases for high atomic numbers. * Species-specific diagnostics in AM geometry with vectorization. - * Custom functions in ``ParticleBinning`` with python 3.12 + * Custom functions in ``ParticleBinning`` with python 3.12. + * Handle prescribed fields in AM. ---- diff --git a/doc/Sphinx/Use/namelist.rst b/doc/Sphinx/Use/namelist.rst index a07f19005..e29830115 100755 --- a/doc/Sphinx/Use/namelist.rst +++ b/doc/Sphinx/Use/namelist.rst @@ -2056,11 +2056,14 @@ This feature is accessible using the ``PrescribedField`` block:: .. py:data:: field Field name: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``. + AM field name: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"`` or ``"Bt_m_mode_m"``. .. warning:: When prescribing a magnetic field, always use the time-centered fields ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``. These fields are those used in the particle pusher, and are defined at integer time-steps. + When prescribing a field in AM geometry, the mode "m" must be specified explicitly in the name of the field and the profile + must return a complex value. .. py:data:: profile @@ -3641,4 +3644,4 @@ namelist. They should not be re-defined by the user! .. note:: These variables can be access during ``happi`` post-processing, e.g. - ``S.namelist.smilei_mpi_size``. \ No newline at end of file + ``S.namelist.smilei_mpi_size``. diff --git a/src/ElectroMagn/ElectroMagnAM.cpp b/src/ElectroMagn/ElectroMagnAM.cpp index 700c0b901..286136746 100755 --- a/src/ElectroMagn/ElectroMagnAM.cpp +++ b/src/ElectroMagn/ElectroMagnAM.cpp @@ -1961,44 +1961,6 @@ void ElectroMagnAM::compute_B_m_fromEB() } - - -void ElectroMagnAM::applyPrescribedFields( Patch *patch, double time ) -{ - -#ifdef _TODO_AM -#endif - int Nmodes = El_.size(); - - Field *field; - - for (int imode=0;imode::iterator extfield=prescribedFields.begin(); extfield!=prescribedFields.end(); extfield++ ) { - string name = LowerCase( extfield->savedField->name ); - if( El_[imode] && name==LowerCase( El_[imode]->name ) ) { - field = El_[imode]; - } else if( Er_[imode] && name==LowerCase( Er_[imode]->name ) ) { - field = Er_[imode]; - } else if( Et_[imode] && name==LowerCase( Et_[imode]->name ) ) { - field = Et_[imode]; - } else if( Bl_[imode] && name==LowerCase( Bl_[imode]->name ) ) { - field = Bl_[imode]; - } else if( Br_[imode] && name==LowerCase( Br_[imode]->name ) ) { - field = Br_[imode]; - } else if( Bt_[imode] && name==LowerCase( Bt_[imode]->name ) ) { - field = Bt_[imode]; - } else { - field = NULL; - } - - if( field ){ - applyPrescribedField( field, extfield->profile, patch, time ); - } - } - } - -} - void ElectroMagnAM::applyExternalField( Field *my_field, Profile *profile, Patch *patch ) { cField2D *field2D=static_cast( my_field ); diff --git a/src/ElectroMagn/ElectroMagnAM.h b/src/ElectroMagn/ElectroMagnAM.h index cd3063113..1da60f42e 100755 --- a/src/ElectroMagn/ElectroMagnAM.h +++ b/src/ElectroMagn/ElectroMagnAM.h @@ -196,8 +196,6 @@ class ElectroMagnAM : public ElectroMagn void computePoynting( unsigned int axis, unsigned int side ) override; //! Method used to impose external fields void applyExternalFields( Patch *patch ) override; - //! Method used to impose external time fields - void applyPrescribedFields( Patch *patch, double time ) override; //! Method used to impose one external field void applyExternalField( Field *, Profile *, Patch * ) override; //! Method used to impose one external time field From e108d241698e75f97c81844bb7cc78ffc449b208 Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Fri, 28 Jun 2024 16:51:11 +0200 Subject: [PATCH 05/26] Small doc improvement on external and prescribed fields --- doc/Sphinx/Use/namelist.rst | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/doc/Sphinx/Use/namelist.rst b/doc/Sphinx/Use/namelist.rst index e29830115..67fd1a680 100755 --- a/doc/Sphinx/Use/namelist.rst +++ b/doc/Sphinx/Use/namelist.rst @@ -2009,8 +2009,8 @@ at the beginning of the simulation using the ``ExternalField`` block:: .. py:data:: field - Field name in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx"``, ``"By"``, ``"Bz"``, ``"Bx_m"``, ``"By_m"``, ``"Bz_m"`` - Field name in AM geometry: ``"El"``, ``"Er"``, ``"Et"``, ``"Bl"``, ``"Br"``, ``"Bt"``, ``"Bl_m"``, ``"Br_m"``, ``"Bt_m"``, ``"A"``, ``"A0"`` . + Field names in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx"``, ``"By"``, ``"Bz"``, ``"Bx_m"``, ``"By_m"``, ``"Bz_m"``. + Field names in AM geometry: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_mode_m"``, ``"Br_mode_m"``, ``"Bt_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"``, ``"Bt_m_mode_m"``, ``"A_mode_1"``, ``"A0_mode_1"`` . .. py:data:: profile @@ -2055,13 +2055,16 @@ This feature is accessible using the ``PrescribedField`` block:: .. py:data:: field - Field name: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``. - AM field name: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"`` or ``"Bt_m_mode_m"``. + Field names in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``. + Field names in AM geometry: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"`` or ``"Bt_m_mode_m"``. .. warning:: When prescribing a magnetic field, always use the time-centered fields ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``. These fields are those used in the particle pusher, and are defined at integer time-steps. + +.. warning:: + When prescribing a field in AM geometry, the mode "m" must be specified explicitly in the name of the field and the profile must return a complex value. From 78c57bfba88f7dd1a8301177c6ad75e53762b012 Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Fri, 28 Jun 2024 17:17:39 +0200 Subject: [PATCH 06/26] fix python changes for tornado --- src/Tools/PyTools.h | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/Tools/PyTools.h b/src/Tools/PyTools.h index 1059e0130..3a165f92b 100755 --- a/src/Tools/PyTools.h +++ b/src/Tools/PyTools.h @@ -43,16 +43,17 @@ // SMILEI_PY_RELEASE_GIL // } // SMILEI_PY_RESTORE_MASTER_THREAD -#define SMILEI_PY_SAVE_MASTER_THREAD \ - PyThreadState *_save = NULL; \ - _Pragma("omp master") \ - _save = PyEval_SaveThread(); -#define SMILEI_PY_RESTORE_MASTER_THREAD \ - _Pragma("omp master") \ - PyEval_RestoreThread( _save ); -#define SMILEI_PY_ACQUIRE_GIL PyGILState_STATE _state = PyGILState_Ensure(); -#define SMILEI_PY_RELEASE_GIL PyGILState_Release( _state ); - +#if PY_MAJOR_VERSION > 2 && PY_MINOR_VERSION > 11 + #define SMILEI_PY_SAVE_MASTER_THREAD PyThreadState *_save = NULL; _Pragma("omp master") { _save = PyEval_SaveThread(); } + #define SMILEI_PY_RESTORE_MASTER_THREAD _Pragma("omp master") PyEval_RestoreThread( _save ); + #define SMILEI_PY_ACQUIRE_GIL PyGILState_STATE _state = PyGILState_Ensure(); + #define SMILEI_PY_RELEASE_GIL PyGILState_Release( _state ); +#else + #define SMILEI_PY_SAVE_MASTER_THREAD + #define SMILEI_PY_RESTORE_MASTER_THREAD + #define SMILEI_PY_ACQUIRE_GIL _Pragma("omp critical") { + #define SMILEI_PY_RELEASE_GIL } +#endif //! tools to query python nemlist and get back C++ values and vectors From 20c472efc6625247543029ff8c25e1767434f028 Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Mon, 1 Jul 2024 17:23:45 +0200 Subject: [PATCH 07/26] Handle species reflective boundary conditions on Rmax in AM. --- doc/Sphinx/Overview/releases.rst | 1 + src/ParticleBC/BoundaryConditionType.cpp | 54 ++++++++++++++---------- src/ParticleBC/PartBoundCond.cpp | 2 +- src/Species/SpeciesFactory.h | 2 +- 4 files changed, 34 insertions(+), 25 deletions(-) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index f30088fd5..b6a53afee 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -29,6 +29,7 @@ Changes made in the repository (not released) * Species-specific diagnostics in AM geometry with vectorization. * Custom functions in ``ParticleBinning`` with python 3.12. * Handle prescribed fields in AM. + * Handle species reflective boundary conditions on Rmax in AM. ---- diff --git a/src/ParticleBC/BoundaryConditionType.cpp b/src/ParticleBC/BoundaryConditionType.cpp index 304656eca..55579a2c7 100755 --- a/src/ParticleBC/BoundaryConditionType.cpp +++ b/src/ParticleBC/BoundaryConditionType.cpp @@ -143,7 +143,7 @@ void reflect_particle_wall( Species *species, int imin, int imax, int direction, } // direction not used below, direction is "r" -void refl_particle_AM( Species *species, int imin, int imax, int /*direction*/, double limit_sup, double /*dt*/, std::vector &/*invgf*/, Random * /*rand*/, double &energy_change ) +void refl_particle_AM( Species *species, int imin, int imax, int /*direction*/, double limit_sup, double /*dt*/, std::vector &invgf, Random * /*rand*/, double &energy_change ) { energy_change = 0.; // no energy loss during reflection @@ -151,36 +151,44 @@ void refl_particle_AM( Species *species, int imin, int imax, int /*direction*/, double* position_z = species->particles->getPtrPosition(2); double* momentum_y = species->particles->getPtrMomentum(1); double* momentum_z = species->particles->getPtrMomentum(2); - //limite_sup = 2*Rmax. - //We look for the coordiunate of the point at which the particle crossed the boundary - //We need to fine the parameter t at which (y + py*t)+(z+pz*t) = Rmax^2 + //We look for the coordinate of the point at which the particle crossed the boundary + //We need to fine the parameter t at which (y - vy*t)^2+(z-vz*t)^2 = Rmax^2 for (int ipart=imin ; ipart= limit_sup*limit_sup ) { - limit_sup *= 2.; - double b = 2*( position_y[ ipart ]*momentum_y[ ipart ] + position_z[ ipart ]*momentum_z[ ipart ] ); - double r2 = ( position_y[ ipart ]*position_y[ ipart ] + position_z[ ipart ]*position_z[ ipart ] ); - double pr2 = ( momentum_y[ ipart ]*momentum_y[ ipart ] + momentum_z[ ipart ]*momentum_z[ ipart ] ); - double delta = b*b - pr2*( 4*r2 - limit_sup*limit_sup ); + double r2 = position_y[ipart]*position_y[ipart]+position_z[ipart]*position_z[ipart]; + if ( r2 >= limit_sup*limit_sup ) { + //solving v2*t2 + b*t + r2-rmax^2 = 0 + double b = -2.*( position_y[ ipart ]*momentum_y[ ipart ] + position_z[ ipart ]*momentum_z[ ipart ] )*invgf[ipart]; + double p2 = ( momentum_y[ ipart ]*momentum_y[ ipart ] + momentum_z[ ipart ]*momentum_z[ ipart ] ); + double v2 = p2 * invgf[ipart]*invgf[ipart]; + double delta = b*b - 4.*v2*( r2 - limit_sup*limit_sup ); - //b and delta are neceseraliy >=0 otherwise there are no solution which means that something unsual happened - if( b < 0 || delta < 0 ) { + //delta is neceseraliy >=0 otherwise there are no solution which means that something unsual happened + if( delta < 0 ) { ERROR( "There are no solution to reflexion. This should never happen" ); } - double t = ( -b + sqrt( delta ) )/pr2*0.5; + double t = ( -b - sqrt( delta ) )/v2*0.5; double y0, z0; //Coordinates of the crossing point 0 - y0 = position_y[ ipart ] + momentum_y[ ipart ]*t ; - z0 = position_z[ ipart ] + momentum_z[ ipart ]*t ; + y0 = position_y[ ipart ] - momentum_y[ ipart ]*t*invgf[ipart] ; + z0 = position_z[ ipart ] - momentum_z[ ipart ]*t*invgf[ipart] ; + + //Update momentum after reflexion at crossing point + //pr = p.er + //pt = p.et + //er = (y0/limit_sup, z0/limit_sup) + //et = (-z0/limit_sup, y0/limit_sup) + double pr = ( momentum_y[ ipart ]*y0 + momentum_z[ ipart ]*z0) / limit_sup; + double pt = (-momentum_y[ ipart ]*z0 + momentum_z[ ipart ]*y0) / limit_sup; + + // New p = -pr.er + pt.et + momentum_y[ipart] = (-pr*y0 - pt*z0)/limit_sup; + momentum_z[ipart] = (-pr*z0 + pt*y0)/limit_sup; + + //Update new particle position as a movement from (y0,z0) with new momentum during time t: - //Update new particle position as a reflexion to the plane tangent to the circle at 0. - - position_y[ ipart ] -= 4*( position_y[ ipart ]-y0 )*y0/limit_sup ; - position_z[ ipart ] -= 4*( position_z[ ipart ]-z0 )*z0/limit_sup ; - - momentum_y[ ipart ] *= 1 - 2*y0/limit_sup ; - momentum_z[ ipart ] *= 1 - 2*z0/limit_sup ; + position_y[ ipart ] = y0 + momentum_y[ipart]*invgf[ipart]*t ; + position_z[ ipart ] = z0 + momentum_z[ipart]*invgf[ipart]*t ; } } } diff --git a/src/ParticleBC/PartBoundCond.cpp b/src/ParticleBC/PartBoundCond.cpp index ca90d6358..e7377f345 100755 --- a/src/ParticleBC/PartBoundCond.cpp +++ b/src/ParticleBC/PartBoundCond.cpp @@ -246,7 +246,7 @@ PartBoundCond::PartBoundCond( Params ¶ms, Species *species, Patch *patch ) }//nDim_particle>1 else if( isAM ) { - + // Ymax bc_ymin = &internal_inf_AM; bc_ymax = &internal_sup_AM; diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index b250182cd..edef5e334 100755 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -669,7 +669,7 @@ class SpeciesFactory t << " " << this_species->boundary_conditions_[iDim][ii]; } } - if( (params.geometry=="AMcylindrical") && ( this_species->boundary_conditions_[1][1] != "remove" ) && ( this_species->boundary_conditions_[1][1] != "stop" ) ) { + if( (params.geometry=="AMcylindrical") && ( this_species->boundary_conditions_[1][1] != "remove" ) && ( this_species->boundary_conditions_[1][1] != "stop" ) && ( this_species->boundary_conditions_[1][1] != "reflective" ) ) { ERROR_NAMELIST( " In AM geometry particle boundary conditions supported in Rmax are 'remove' and 'stop' ", LINK_NAMELIST + std::string("#species") From 4ba471a83cb984e70a6b95290f6bb3f77da40957 Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Tue, 2 Jul 2024 17:14:33 +0200 Subject: [PATCH 08/26] fix happi average (last bin) --- doc/Sphinx/Overview/releases.rst | 12 ++++++++---- happi/_Diagnostics/Diagnostic.py | 4 ++-- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index b6a53afee..3f4d92fd4 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -23,13 +23,17 @@ You can find older, `unsupported versions here Date: Fri, 5 Jul 2024 18:02:04 +0200 Subject: [PATCH 09/26] clarification on the Prescribed Fields in the doc --- doc/Sphinx/Use/namelist.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/doc/Sphinx/Use/namelist.rst b/doc/Sphinx/Use/namelist.rst index 67fd1a680..3018fb956 100755 --- a/doc/Sphinx/Use/namelist.rst +++ b/doc/Sphinx/Use/namelist.rst @@ -2068,6 +2068,12 @@ This feature is accessible using the ``PrescribedField`` block:: When prescribing a field in AM geometry, the mode "m" must be specified explicitly in the name of the field and the profile must return a complex value. +.. warning:: + + ``PrescribedFields`` are not visible in the ``Field`` diagnostic, + but can be visualised through ``Probes`` and with the fields attributes of ``TrackParticles`` + (since they sample the total field acting on the macro-particles). + .. py:data:: profile :type: float or :doc:`profile ` From 9b42649f20b4b25411eb965cc0c860561e04efa4 Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Thu, 11 Jul 2024 16:31:32 +0200 Subject: [PATCH 10/26] move initExchParticles outside of dynamics --- src/Patch/VectorPatch.cpp | 35 ++++++++++++----------------------- src/Patch/VectorPatch.h | 5 +++-- src/Smilei.cpp | 7 ++++--- 3 files changed, 19 insertions(+), 28 deletions(-) diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 5639a4656..bbce82778 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -398,19 +398,6 @@ void VectorPatch::dynamics( Params ¶ms, timers.multiphoton_Breit_Wheeler_timer.updateThreaded( *this, params.printNow( itime ) ); #endif - timers.syncPart.restart(); - for( unsigned int ispec=0 ; ispec<( *this )( 0 )->vecSpecies.size(); ispec++ ) { - Species *spec = species( 0, ispec ); - if ( (!params.Laser_Envelope_model) && (spec->isProj( time_dual, simWindow )) ){ - SyncVectorPatch::initExchParticles( ( *this ), ispec, params, smpi ); // Included sortParticles - } // end condition on Species and on envelope model - } // end loop on species - //MESSAGE("exchange particles"); - timers.syncPart.update( params.printNow( itime ) ); - -#ifdef __DETAILED_TIMERS - timers.sorting.update( *this, params.printNow( itime ) ); -#endif } // END dynamics // --------------------------------------------------------------------------------------------------------------------- @@ -457,6 +444,18 @@ void VectorPatch::projectionForDiags( Params ¶ms, } // END projection for diags +void VectorPatch::initExchParticles( Params ¶ms, SmileiMPI *smpi, SimWindow *simWindow, + double time_dual, Timers &timers, int itime ) +{ + timers.syncPart.restart(); + for( unsigned int ispec=0 ; ispec<( *this )( 0 )->vecSpecies.size(); ispec++ ) { + if( species( 0, ispec )->isProj( time_dual, simWindow ) ) { + SyncVectorPatch::initExchParticles( *this, ispec, params, smpi ); + } + } + timers.syncPart.update( params.printNow( itime ) ); +} + // --------------------------------------------------------------------------------------------------------------------- //! For all patches, exchange particles and sort them. // --------------------------------------------------------------------------------------------------------------------- @@ -4588,16 +4587,6 @@ void VectorPatch::ponderomotiveUpdatePositionAndCurrents( Params ¶ms, timers.cell_keys.update( *this, params.printNow( itime ) ); #endif - timers.syncPart.restart(); - for( unsigned int ispec=0 ; ispec<( *this )( 0 )->vecSpecies.size(); ispec++ ) { - if( ( *this )( 0 )->vecSpecies[ispec]->isProj( time_dual, simWindow ) ) { - SyncVectorPatch::initExchParticles( ( *this ), ispec, params, smpi ); // Included sortParticles - } // end condition on species - } // end loop on species - timers.syncPart.update( params.printNow( itime ) ); - - - } // END ponderomotiveUpdatePositionAndCurrents diff --git a/src/Patch/VectorPatch.h b/src/Patch/VectorPatch.h index 051d78276..d2e12e545 100755 --- a/src/Patch/VectorPatch.h +++ b/src/Patch/VectorPatch.h @@ -157,9 +157,10 @@ public : Timers &timers, int itime ); //! For all patches, exchange particles and sort them. + void initExchParticles( Params ¶ms, SmileiMPI *smpi, SimWindow *simWindow, + double time_dual, Timers &timers, int itime ); void finalizeExchParticlesAndSort( Params ¶ms, SmileiMPI *smpi, SimWindow *simWindow, - double time_dual, - Timers &timers, int itime ); + double time_dual, Timers &timers, int itime ); void finalizeSyncAndBCFields( Params ¶ms, SmileiMPI *smpi, SimWindow *simWindow, double time_dual, Timers &timers, int itime ); diff --git a/src/Smilei.cpp b/src/Smilei.cpp index 81ba6c258..11e43976e 100755 --- a/src/Smilei.cpp +++ b/src/Smilei.cpp @@ -521,7 +521,9 @@ int main( int argc, char *argv[] ) if( params.Laser_Envelope_model ) { vecPatches.runEnvelopeModule( params, &smpi, simWindow, time_dual, timers, itime ); } // end condition if Laser Envelope Model is used - + + vecPatches.initExchParticles( params, &smpi, simWindow, time_dual, timers, itime ); + // Sum densities vecPatches.sumDensities( params, time_dual, timers, itime, simWindow, &smpi ); @@ -629,8 +631,7 @@ int main( int argc, char *argv[] ) #pragma omp parallel shared (time_dual,smpi,params, vecPatches, region, simWindow, checkpoint, itime) { // finalize particle exchanges and sort particles - vecPatches.finalizeExchParticlesAndSort( params, &smpi, simWindow, - time_dual, timers, itime ); + vecPatches.finalizeExchParticlesAndSort( params, &smpi, simWindow, time_dual, timers, itime ); // Particle merging vecPatches.mergeParticles(params, time_dual,timers, itime ); From 42d097b8d4173a584d3e80ea9c13d27039a72c1c Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Thu, 11 Jul 2024 17:00:16 +0200 Subject: [PATCH 11/26] remove useless test in exch particles --- src/Patch/Patch.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/Patch/Patch.cpp b/src/Patch/Patch.cpp index ca76c6ece..1a96ab654 100755 --- a/src/Patch/Patch.cpp +++ b/src/Patch/Patch.cpp @@ -548,10 +548,6 @@ void Patch::copyExchParticlesToBuffers( int ispec, Params ¶ms ) sendBuffer[2*iDim+0] = buffer.partSend[iDim][0]; sendBuffer[2*iDim+1] = buffer.partSend[iDim][1]; } - if( params.geometry == "AMcylindrical" ) { - copy[0] = copy[0] && ( Pcoordinates[0]!=0 || vecSpecies[ispec]->boundary_conditions_[0][0]=="periodic" ); - copy[1] = copy[1] && ( Pcoordinates[0]!=params.number_of_patches[0]-1 || vecSpecies[ispec]->boundary_conditions_[0][1]=="periodic" ); - } part.copyLeavingParticlesToBuffers( copy, sendBuffer ); From da16dc037f0c6f73ac086e29f1c586e410b7d7c0 Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Fri, 12 Jul 2024 15:45:47 +0000 Subject: [PATCH 12/26] New am shape functions --- .../tstAM_04d_laser_propagation_Order1.py | 98 ++ doc/Sphinx/Use/namelist.rst | 4 + .../InterpolatorAM1OrderRuyten.cpp | 886 +++++++++++++ src/Interpolator/InterpolatorAM1OrderRuyten.h | 118 ++ .../InterpolatorAM1OrderRuytenV.cpp | 323 +++++ .../InterpolatorAM1OrderRuytenV.h | 40 + src/Interpolator/InterpolatorFactory.h | 16 +- src/Params/Params.cpp | 12 +- src/Projector/ProjectorAM1OrderRuyten.cpp | 1133 +++++++++++++++++ src/Projector/ProjectorAM1OrderRuyten.h | 58 + src/Projector/ProjectorAM1OrderRuytenV.cpp | 844 ++++++++++++ src/Projector/ProjectorAM1OrderRuytenV.h | 350 +++++ src/Projector/ProjectorFactory.h | 14 +- ...date_tstAM_04d_laser_propagation_Order1.py | 29 + .../tstAM_04d_laser_propagation_Order1.py.txt | Bin 0 -> 33947 bytes 15 files changed, 3914 insertions(+), 11 deletions(-) create mode 100755 benchmarks/tstAM_04d_laser_propagation_Order1.py create mode 100755 src/Interpolator/InterpolatorAM1OrderRuyten.cpp create mode 100755 src/Interpolator/InterpolatorAM1OrderRuyten.h create mode 100755 src/Interpolator/InterpolatorAM1OrderRuytenV.cpp create mode 100755 src/Interpolator/InterpolatorAM1OrderRuytenV.h create mode 100755 src/Projector/ProjectorAM1OrderRuyten.cpp create mode 100755 src/Projector/ProjectorAM1OrderRuyten.h create mode 100755 src/Projector/ProjectorAM1OrderRuytenV.cpp create mode 100755 src/Projector/ProjectorAM1OrderRuytenV.h create mode 100755 validation/analyses/validate_tstAM_04d_laser_propagation_Order1.py create mode 100644 validation/references/tstAM_04d_laser_propagation_Order1.py.txt diff --git a/benchmarks/tstAM_04d_laser_propagation_Order1.py b/benchmarks/tstAM_04d_laser_propagation_Order1.py new file mode 100755 index 000000000..ffd8d2792 --- /dev/null +++ b/benchmarks/tstAM_04d_laser_propagation_Order1.py @@ -0,0 +1,98 @@ +# ---------------------------------------------------------------------------------------- +# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI +# ---------------------------------------------------------------------------------------- + +import math +dx = 0.251327 +dtrans = 1.96349 +dt = 0.96 * dx +nx = 960 +ntrans = 256 +Lx = nx * dx +Ltrans = ntrans * dtrans +npatch_x = 64 +npatch_trans =32 +Nit = 2000 + + +Main( + geometry = "AMcylindrical", + number_of_AM=2, + interpolation_order = 1, + timestep = dt, + simulation_time = dt*Nit, + cell_length = [dx, dtrans], + grid_length = [ Lx, Ltrans], + number_of_patches = [npatch_x, npatch_trans], + cluster_width = 5, + EM_boundary_conditions = [ + ["silver-muller","silver-muller"], + ["PML","PML"], + ], + number_of_pml_cells = [[0,0], [10,10]], + solve_poisson = False, + print_every = 100 +) + +MovingWindow( + time_start = Main.grid_length[0] - 50*dx, #Leaves 2 patches untouched, in front of the laser. + velocity_x = 0.996995486 +) + +ne = 0.0045 +begin_upramp = 10. #Density is 0 before that and up ramp starts. +Lupramp = 100. #Length of the upramp +Lplateau = 15707. #Length of the plateau +Ldownramp = 2356.19 #Length of the down ramp +xplateau = begin_upramp + Lupramp # Start of the plateau +begin_downramp = xplateau + Lplateau # Beginning of the output ramp. +finish = begin_downramp + Ldownramp # End of plasma + +g = polygonal(xpoints=[begin_upramp, xplateau, begin_downramp, finish], xvalues=[0, ne, ne, 0.]) + +def my_profile(x,y): + return g(x,y) + +Species( + name = "electron", + position_initialization = "regular", + momentum_initialization = "cold", + ionization_model = "none", + particles_per_cell = 30, + c_part_max = 1.0, + mass = 1.0, + charge = -1.0, + charge_density = my_profile, # Here absolute value of the charge is 1 so charge_density = nb_density + mean_velocity = [0., 0., 0.], + time_frozen = 0.0, + boundary_conditions = [ + ["remove", "remove"], + ["reflective", "remove"], + ], +) + +laser_fwhm = 82. +LaserGaussianAM( + box_side = "xmin", + a0 = 2., + focus = [10.], + waist = 120., + time_envelope = tgaussian(center=2**0.5*laser_fwhm, fwhm=laser_fwhm) +) + +DiagProbe( + every = 1000, + origin = [0., 2*dtrans, 0.], + corners = [ + [Main.grid_length[0], 2*dtrans, 0.] + ], + number = [nx], +) + +DiagPerformances( + every = 1000, +) + +DiagScalar( + every = 20 +) diff --git a/doc/Sphinx/Use/namelist.rst b/doc/Sphinx/Use/namelist.rst index 3018fb956..163aa8654 100755 --- a/doc/Sphinx/Use/namelist.rst +++ b/doc/Sphinx/Use/namelist.rst @@ -143,9 +143,13 @@ The block ``Main`` is **mandatory** and has the following syntax:: Interpolation order, defines particle shape function: + * ``1`` : 2 points stencil in r with Ruyten correction, 3 points stencil in x. Supported only in AM geometry. * ``2`` : 3 points stencil, supported in all configurations. * ``4`` : 5 points stencil, not supported in vectorized 2D geometry. + The Ruyten correction is the scheme described bu equation 4.2 in `this paper `_ . + It allows for a more accurate description on axis at the cost of a higher statistic noise so it often requires the use of more macro-particles. + .. py:data:: interpolator :default: ``"momentum-conserving"`` diff --git a/src/Interpolator/InterpolatorAM1OrderRuyten.cpp b/src/Interpolator/InterpolatorAM1OrderRuyten.cpp new file mode 100755 index 000000000..28313aa8a --- /dev/null +++ b/src/Interpolator/InterpolatorAM1OrderRuyten.cpp @@ -0,0 +1,886 @@ +#include "InterpolatorAM1OrderRuyten.h" + +#include +#include +#include +#include "ElectroMagn.h" +#include "ElectroMagnAM.h" +#include "cField2D.h" +#include "Particles.h" +#include "LaserEnvelope.h" +#include +#include "dcomplex.h" + +using namespace std; + + +// --------------------------------------------------------------------------------------------------------------------- +// Creator for InterpolatorAM1OrderRuyten +// --------------------------------------------------------------------------------------------------------------------- +InterpolatorAM1OrderRuyten::InterpolatorAM1OrderRuyten( Params ¶ms, Patch *patch ) : InterpolatorAM( patch ) +{ + + D_inv_[0] = 1.0/params.cell_length[0]; + D_inv_[1] = 1.0/params.cell_length[1]; + nmodes_ = params.nmodes; + +} + +// --------------------------------------------------------------------------------------------------------------------- +// 1st Order Interpolation of the fields with Ruyten correction at a the particle position (2 nodes are used) +// --------------------------------------------------------------------------------------------------------------------- +void InterpolatorAM1OrderRuyten::fields( ElectroMagn *EMfields, Particles &particles, int ipart, int nparts, double *ELoc, double *BLoc ) +{ + + //Treat mode 0 first + + // Static cast of the electromagnetic fields + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + + + // Normalized particle position + double xpn = particles.position( 0, ipart ) * D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + exp_m_theta_ = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; //exp(-i theta) + complex exp_mm_theta = 1. ; //exp(-i m theta) + // Compute coeffs + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + //Here we assume that mode 0 is real !! + // Interpolation of El^(d,p) + *( ELoc+0*nparts ) = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + *( ELoc+1*nparts ) = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + *( ELoc+2*nparts ) = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + *( BLoc+0*nparts ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + *( BLoc+1*nparts ) = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + *( BLoc+2*nparts ) = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + + exp_mm_theta *= exp_m_theta_ ; + + *( ELoc+0*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[1], El, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( ELoc+1*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[1], Er, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( ELoc+2*nparts ) += std::real( compute( &coeffxp[1], &coeffyp[1], Et, idx_p[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[1], Bl, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( BLoc+1*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[1], Br, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+2*nparts ) += std::real( compute( &coeffxd[1], &coeffyd[1], Bt, idx_d[0], idx_d[1] )* exp_mm_theta ) ; + } + + //Translate field into the cartesian y,z coordinates + double delta2 = std::real( exp_m_theta_ ) * *( ELoc+1*nparts ) + std::imag( exp_m_theta_ ) * *( ELoc+2*nparts ); + *( ELoc+2*nparts ) = -std::imag( exp_m_theta_ ) * *( ELoc+1*nparts ) + std::real( exp_m_theta_ ) * *( ELoc+2*nparts ); + *( ELoc+1*nparts ) = delta2 ; + delta2 = std::real( exp_m_theta_ ) * *( BLoc+1*nparts ) + std::imag( exp_m_theta_ ) * *( BLoc+2*nparts ); + *( BLoc+2*nparts ) = -std::imag( exp_m_theta_ ) * *( BLoc+1*nparts ) + std::real( exp_m_theta_ ) * *( BLoc+2*nparts ); + *( BLoc+1*nparts ) = delta2 ; + +} // END InterpolatorAM1OrderRuyten + +void InterpolatorAM1OrderRuyten::fieldsAndCurrents( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *, int ithread, LocalFields *JLoc, double *RhoLoc ) +{ + int ipart = *istart; + + double *ELoc = &( smpi->dynamics_Epart[ithread][ipart] ); + double *BLoc = &( smpi->dynamics_Bpart[ithread][ipart] ); + double *BLocyBTIS3; + double *BLoczBTIS3; + if (smpi->use_BTIS3){ + BLocyBTIS3 = &( smpi->dynamics_Bpart_yBTIS3[ithread][ipart] ); + BLoczBTIS3 = &( smpi->dynamics_Bpart_zBTIS3[ithread][ipart] ); + } + + // Interpolate E, B + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + cField2D *Jl = ( static_cast( EMfields ) )->Jl_[0]; + cField2D *Jr = ( static_cast( EMfields ) )->Jr_[0]; + cField2D *Jt = ( static_cast( EMfields ) )->Jt_[0]; + cField2D *Rho= ( static_cast( EMfields ) )->rho_AM_[0]; + cField2D *Br_BTIS3; + cField2D *Bt_BTIS3; + if (smpi->use_BTIS3){ + Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[0]; + Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[0]; + } + + // Normalized particle position + double xpn = particles.position( 0, ipart ) * D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + complex exp_mm_theta = 1. ; + + // Calculate coeffs + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + int nparts( particles.numberOfParticles() ); + + *( ELoc+0*nparts ) = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + *( ELoc+1*nparts ) = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + *( ELoc+2*nparts ) = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + *( BLoc+0*nparts ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + *( BLoc+1*nparts ) = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + *( BLoc+2*nparts ) = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + + // Interpolation of Jl^(d,p,p) + JLoc->x = std::real( compute( &coeffxd[1], &coeffyp[0], Jl, idx_d[0], idx_p[1] ) ); + // Interpolation of Jr^(p,d,p) + JLoc->y = std::real( compute( &coeffxp[1], &coeffyd[0], Jr, idx_p[0], idx_d[1]) ); + // Interpolation of Jt^(p,p,d) + JLoc->z = std::real( compute( &coeffxp[1], &coeffyp[0], Jt, idx_p[0], idx_p[1]) ); + // Interpolation of Rho^(p,p,p) + ( *RhoLoc ) = std::real( compute( &coeffxp[1], &coeffyp[0], Rho, idx_p[0], idx_p[1] ) ); + + if (smpi->use_BTIS3){ + // BTIS fields, in the x direction they are centered as Ey and Ez + // Interpolation of Br^(p,p) for BTIS + *( BLocyBTIS3+0*nparts ) = std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] ) ); + // Interpolation of Bt^(p,d) for BTIS + *( BLoczBTIS3+0*nparts ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] ) ); + } + + if (r > 0){ + exp_m_theta_ = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; + } else { + exp_m_theta_ = 1. ; + } + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + Jl = ( static_cast( EMfields ) )->Jl_[imode]; + Jr = ( static_cast( EMfields ) )->Jr_[imode]; + Jt = ( static_cast( EMfields ) )->Jt_[imode]; + Rho= ( static_cast( EMfields ) )->rho_AM_[imode]; + + if (smpi->use_BTIS3){ + Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[imode]; + Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[imode]; + } + + exp_mm_theta *= exp_m_theta_ ; + + *( ELoc+0*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( ELoc+1*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( ELoc+2*nparts ) += std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( BLoc+1*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+2*nparts ) += std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] )* exp_mm_theta ) ; + JLoc->x += std::real( compute( &coeffxd[1], &coeffyp[0], Jl, idx_d[0], idx_p[1] ) * exp_mm_theta ) ; + JLoc->y += std::real( compute( &coeffxp[1], &coeffyd[0], Jr, idx_p[0], idx_d[1] ) * exp_mm_theta ) ; + JLoc->z += std::real( compute( &coeffxp[1], &coeffyp[0], Jt, idx_p[0], idx_p[1] ) * exp_mm_theta ) ; + ( *RhoLoc ) += std::real( compute( &coeffxp[1], &coeffyp[0], Rho, idx_p[0], idx_p[1] )* exp_mm_theta ) ; + + if (smpi->use_BTIS3){ + // BTIS3 fields, in the x direction they are centered as Ey and Ez + *( BLocyBTIS3+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] )* exp_mm_theta ); + *( BLoczBTIS3+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] )* exp_mm_theta ); + } + } + double delta2 = std::real( exp_m_theta_ ) * *( ELoc+1*nparts ) + std::imag( exp_m_theta_ ) * *( ELoc+2*nparts ); + *( ELoc+2*nparts ) = -std::imag( exp_m_theta_ ) * *( ELoc+1*nparts ) + std::real( exp_m_theta_ ) * *( ELoc+2*nparts ); + *( ELoc+1*nparts ) = delta2 ; + delta2 = std::real( exp_m_theta_ ) * *( BLoc+1*nparts ) + std::imag( exp_m_theta_ ) * *( BLoc+2*nparts ); + *( BLoc+2*nparts ) = -std::imag( exp_m_theta_ ) * *( BLoc+1*nparts ) + std::real( exp_m_theta_ ) * *( BLoc+2*nparts ); + *( BLoc+1*nparts ) = delta2 ; + delta2 = std::real( exp_m_theta_ ) * JLoc->y + std::imag( exp_m_theta_ ) * JLoc->z; + JLoc->z = -std::imag( exp_m_theta_ ) * JLoc->y + std::real( exp_m_theta_ ) * JLoc->z; + JLoc->y = delta2 ; + if (smpi->use_BTIS3){ + delta2 = std::real( exp_m_theta_ ) * *( BLocyBTIS3+0*nparts ) + std::imag( exp_m_theta_ ) * *( BLoczBTIS3+0*nparts ); + *( BLoczBTIS3+0*nparts ) = -std::imag( exp_m_theta_ ) * *( BLocyBTIS3+0*nparts ) + std::real( exp_m_theta_ ) * *( BLoczBTIS3+0*nparts ); + *( BLocyBTIS3+0*nparts ) = delta2 ; + } + +} + +// Interpolator on another field than the basic ones +void InterpolatorAM1OrderRuyten::oneField( Field **field, Particles &particles, int *istart, int *iend, double *Jxloc, double *Jyloc, double *Jzloc, double *Rholoc ) +{ + + // **field points to the first field of the species of interest in EM->allFields + // They are ordered as Jx0, Jy0, Jz0, Rho0, Jx1, Jy1, Jz1, Rho1, etc. + + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + double xpn = particles.position( 0, ipart )*D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + complex exp_m_theta_ = 1., exp_mm_theta = 1. ; + if (r > 0) { + exp_m_theta_ = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; + } + + double Jx_ = 0., Jy_ = 0., Jz_ = 0., Rho_ = 0.; + for( unsigned int imode = 0; imode < nmodes_ ; imode++ ) { + cField2D *Jl = static_cast( *(field+4*imode+0) ); + cField2D *Jr = static_cast( *(field+4*imode+1) ); + cField2D *Jt = static_cast( *(field+4*imode+2) ); + cField2D *Rho = static_cast( *(field+4*imode+3) ); + Jx_ += std::real( compute( &coeffxd[1], &coeffyp[0], Jl , idx_d[0], idx_p[1] ) * exp_mm_theta ); + Jy_ += std::real( compute( &coeffxp[1], &coeffyd[0], Jr , idx_p[0], idx_d[1] ) * exp_mm_theta ); + Jz_ += std::real( compute( &coeffxp[1], &coeffyp[0], Jt , idx_p[0], idx_p[1] ) * exp_mm_theta ); + Rho_ += std::real( compute( &coeffxp[1], &coeffyp[0], Rho, idx_p[0], idx_p[1] ) * exp_mm_theta ); + + exp_mm_theta *= exp_m_theta_; + } + Jxloc [ipart] = Jx_; + Jyloc [ipart] = std::real( exp_m_theta_ ) * Jy_ + std::imag( exp_m_theta_ ) * Jz_; + Jzloc [ipart] = -std::imag( exp_m_theta_ ) * Jy_ + std::real( exp_m_theta_ ) * Jz_; + Rholoc[ipart] = Rho_; + } +} + +void InterpolatorAM1OrderRuyten::fieldsWrapper( ElectroMagn *EMfields, + Particles &particles, + SmileiMPI *smpi, + int *istart, + int *iend, + int ithread, + unsigned int, + int ) +{ + + double *Epart = &( smpi->dynamics_Epart[ithread][0] ); + double *Bpart = &( smpi->dynamics_Bpart[ithread][0] ); + int *iold = &( smpi->dynamics_iold[ithread][0] ); + double *delta = &( smpi->dynamics_deltaold[ithread][0] ); + std::complex *eitheta_old = &( smpi->dynamics_eithetaold[ithread][0] ); + + //Loop on bin particles + int nparts( particles.numberOfParticles() ); + + if (!smpi->use_BTIS3){ // without B-TIS3 interpolation + + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + complex exp_m_theta_local ; + complex exp_mm_theta_local = 1. ; + + // Normalized particle position + double xpn = particles.position( 0, ipart ) * D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + exp_m_theta_local = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; //exp(-i theta) + //exp(-i m theta) + + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + // Static cast of the electromagnetic fields, mode 0 + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + + //Here we assume that mode 0 is real !! + // Interpolation of El^(d,p) + *( Epart+0*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + *( Epart+1*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + *( Epart+2*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + *( Bpart+0*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + *( Bpart+1*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + *( Bpart+2*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + + exp_mm_theta_local *= exp_m_theta_local ; + + *( Epart+0*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Epart+1*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + *( Epart+2*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Bpart+0*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + *( Bpart+1*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Bpart+2*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] )* exp_mm_theta_local ) ; + } + + //Translate field into the cartesian y,z coordinates + double delta2 = std::real( exp_m_theta_local ) * *( Epart+1*nparts+ipart ) + std::imag( exp_m_theta_local ) * *( Epart+2*nparts+ipart ); + *( Epart+2*nparts+ipart ) = -std::imag( exp_m_theta_local ) * *( Epart+1*nparts+ipart ) + std::real( exp_m_theta_local ) * *( Epart+2*nparts+ipart ); + *( Epart+1*nparts+ipart ) = delta2 ; + delta2 = std::real( exp_m_theta_local ) * *( Bpart+1*nparts+ipart ) + std::imag( exp_m_theta_local ) * *( Bpart+2*nparts+ipart ); + *( Bpart+2*nparts+ipart ) = -std::imag( exp_m_theta_local ) * *( Bpart+1*nparts+ipart ) + std::real( exp_m_theta_local ) * *( Bpart+2*nparts+ipart ); + *( Bpart+1*nparts+ipart ) = delta2 ; + + // store indices and delta + *( iold+0*nparts+ipart) = idx_p[0]; + *( iold+1*nparts+ipart) = idx_p[1]; + *( delta+0*nparts+ipart) = delta_p[0]; + *( delta+1*nparts+ipart) = delta_p[1]; + *( eitheta_old+ipart ) = 2.*std::real(exp_m_theta_local) - exp_m_theta_local ; //exp(i theta) + + } // end ipart loop + + } else { // with B-TIS3 interpolation + + double *BLocyBTIS3 = &( smpi->dynamics_Bpart_yBTIS3[ithread][0] ); + double *BLoczBTIS3 = &( smpi->dynamics_Bpart_zBTIS3[ithread][0] ); + + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + complex exp_m_theta_local ; + complex exp_mm_theta_local = 1. ; + + // Normalized particle position + double xpn = particles.position( 0, ipart ) * D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + exp_m_theta_local = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; //exp(-i theta) + //exp(-i m theta) + + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + // Static cast of the electromagnetic fields, mode 0 + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + cField2D *Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[0]; + cField2D *Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[0]; + + //Here we assume that mode 0 is real !! + // Interpolation of El^(d,p) + *( Epart+0*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + *( Epart+1*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + *( Epart+2*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + *( Bpart+0*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + *( Bpart+1*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + *( Bpart+2*nparts+ipart ) = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + + // BTIS3 fields, in the x direction they are centered as Ey and Ez + // Interpolation of Br^(p,p) for BTIS + *( BLocyBTIS3+0*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] ) ); + // Interpolation of Bt^(p,d) for BTIS + *( BLoczBTIS3+0*nparts+ipart ) = std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] ) ); + + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[imode]; + Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[imode]; + + exp_mm_theta_local *= exp_m_theta_local ; + + *( Epart+0*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyp[0], El , idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Epart+1*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyd[0], Er , idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + *( Epart+2*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyp[0], Et , idx_p[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Bpart+0*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bl , idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + *( Bpart+1*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyp[0], Br , idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + *( Bpart+2*nparts+ipart ) += std::real( compute( &coeffxd[1], &coeffyd[0], Bt , idx_d[0], idx_d[1] )* exp_mm_theta_local ) ; + // BTIS3 fields, in the x direction they are centered as Ey and Ez + *( BLocyBTIS3+0*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] )* exp_mm_theta_local ); + *( BLoczBTIS3+0*nparts+ipart ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] )* exp_mm_theta_local ); + } + + //Translate field into the cartesian y,z coordinates + double delta2 = std::real( exp_m_theta_local ) * *( Epart+1*nparts+ipart ) + std::imag( exp_m_theta_local ) * *( Epart+2*nparts+ipart ); + *( Epart+2*nparts+ipart ) = -std::imag( exp_m_theta_local ) * *( Epart+1*nparts+ipart ) + std::real( exp_m_theta_local ) * *( Epart+2*nparts+ipart ); + *( Epart+1*nparts+ipart ) = delta2 ; + delta2 = std::real( exp_m_theta_local ) * *( Bpart+1*nparts+ipart ) + std::imag( exp_m_theta_local ) * *( Bpart+2*nparts+ipart ); + *( Bpart+2*nparts+ipart ) = -std::imag( exp_m_theta_local ) * *( Bpart+1*nparts+ipart ) + std::real( exp_m_theta_local ) * *( Bpart+2*nparts+ipart ); + *( Bpart+1*nparts+ipart ) = delta2 ; + delta2 = std::real( exp_m_theta_local ) * *( BLocyBTIS3+0*nparts+ipart ) + std::imag( exp_m_theta_local ) * *( BLoczBTIS3+0*nparts+ipart ); + *( BLoczBTIS3+0*nparts+ipart ) = -std::imag( exp_m_theta_local ) * *( BLocyBTIS3+0*nparts+ipart ) + std::real( exp_m_theta_local ) * *( BLoczBTIS3+0*nparts+ipart ); + *( BLocyBTIS3+0*nparts+ipart ) = delta2 ; + + // store indices and delta + *( iold+0*nparts+ipart) = idx_p[0]; + *( iold+1*nparts+ipart) = idx_p[1]; + *( delta+0*nparts+ipart) = delta_p[0]; + *( delta+1*nparts+ipart) = delta_p[1]; + *( eitheta_old+ipart ) = 2.*std::real(exp_m_theta_local) - exp_m_theta_local ; //exp(i theta) + + } // end ipart loop + + } // end with B-TIS3 interpolation + +} + +// Interpolator specific to tracked particles. A selection of particles may be provided +void InterpolatorAM1OrderRuyten::fieldsSelection( ElectroMagn *EMfields, Particles &particles, double *buffer, int offset, vector *selection ) +{ + if( selection ) { + + int nsel_tot = selection->size(); + for( int isel=0 ; isel *Epart = &( smpi->dynamics_Epart[ithread] ); + std::vector *Bpart = &( smpi->dynamics_Bpart[ithread] ); + + std::vector *PHIpart = &( smpi->dynamics_PHIpart[ithread] ); + std::vector *GradPHIpart = &( smpi->dynamics_GradPHIpart[ithread] ); + + std::vector *iold = &( smpi->dynamics_iold[ithread] ); + std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); + std::vector> *eitheta_old = &( smpi->dynamics_eithetaold[ithread] ); + + // Static cast of the envelope fields + Field2D *Phi = static_cast( EMfields->envelope->Phi_ ); + Field2D *GradPhil = static_cast( EMfields->envelope->GradPhil_ ); + Field2D *GradPhir = static_cast( EMfields->envelope->GradPhir_ ); + + // auxiliary quantities + double delta2, xpn, r, rpn; + int nparts = particles.numberOfParticles() ; + + if (!smpi->use_BTIS3){ // without B-TIS3 interpolation + + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + complex exp_m_theta_local ; + complex exp_mm_theta_local = 1. ; + + // Normalized particle position + xpn = particles.position( 0, ipart ) * D_inv_[0]; + r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + rpn = r * D_inv_[1]; + + // Compute coefficients + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + // mode 0 is treated first + + // Interpolate E, B + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + + // Interpolation of El^(d,p) + ( *Epart ) [ 0*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + ( *Epart ) [ 1*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + ( *Epart ) [ 2*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + ( *Bpart ) [ 0*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + ( *Bpart ) [ 1*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + ( *Bpart ) [ 2*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + // Interpolation of Phi^(p,p) + ( *PHIpart ) [ 0*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], Phi, idx_p[0], idx_p[1] ) ; + // Interpolation of GradPhil^(p,p) + ( *GradPHIpart ) [ 0*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], GradPhil, idx_p[0], idx_p[1] ) ; + // Interpolation of GradPhir^(p,p) + ( *GradPHIpart ) [ 1*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], GradPhir, idx_p[0], idx_p[1] ) ; + // GradPhit = 0 in cylindrical symmetry + ( *GradPHIpart ) [ 2*nparts+ipart ] = 0.; + + if (r > 0){ + exp_m_theta_local = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; + } else { + exp_m_theta_local = 1. ; + } + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + + exp_mm_theta_local *= exp_m_theta_local ; + + ( *Epart ) [ 0*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Epart ) [ 1*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + ( *Epart ) [ 2*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 0*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 1*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 2*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] )* exp_mm_theta_local ) ; + } + + // project on x,y,z, remember that GradPhit = 0 in cylindrical symmetry + delta2 = std::real( exp_m_theta_local ) * ( *Epart ) [ 1*nparts+ipart ] + std::imag( exp_m_theta_local ) * ( *Epart ) [ 2*nparts+ipart ]; + ( *Epart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *Epart ) [ 1*nparts+ipart ] + std::real( exp_m_theta_local ) * ( *Epart ) [ 2*nparts+ipart ]; + ( *Epart ) [ 1*nparts+ipart ] = delta2 ; + delta2 = std::real( exp_m_theta_local ) * ( *Bpart ) [ 1*nparts+ipart ] + std::imag( exp_m_theta_local ) * ( *Bpart ) [ 2*nparts+ipart ]; + ( *Bpart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *Bpart ) [ 1*nparts+ipart ] + std::real( exp_m_theta_local ) * ( *Bpart ) [ 2*nparts+ipart ]; + ( *Bpart ) [ 1*nparts+ipart ] = delta2 ; + + delta2 = std::real( exp_m_theta_local ) * ( *GradPHIpart ) [ 1*nparts+ipart ] ; + ( *GradPHIpart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *GradPHIpart ) [ 1*nparts+ipart ] ; + ( *GradPHIpart ) [ 1*nparts+ipart ] = delta2 ; + + //Buffering of iold and delta + ( *iold )[ipart+0*nparts] = idx_p[0]; + ( *iold )[ipart+1*nparts] = idx_p[1]; + ( *delta )[ipart+0*nparts] = delta_p[0]; + ( *delta )[ipart+1*nparts] = delta_p[1]; + ( *eitheta_old)[ipart] = std::conj(exp_m_theta_local); //exp(i theta) + + } // end ipart loop + + } else { // with B-TIS3 interpolation + + std::vector *BLocyBTIS3 = &( smpi->dynamics_Bpart_yBTIS3[ithread] ); + std::vector *BLoczBTIS3 = &( smpi->dynamics_Bpart_zBTIS3[ithread] ); + + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + int idx_p[2], idx_d[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + double coeffxd[3], coeffyd[2]; + complex exp_m_theta_local ; + complex exp_mm_theta_local = 1. ; + + // Normalized particle position + xpn = particles.position( 0, ipart ) * D_inv_[0]; + r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + rpn = r * D_inv_[1]; + + // Compute coefficients + coeffs( xpn, rpn, idx_p, idx_d, coeffxp, coeffyp, coeffxd, coeffyd, delta_p ); + + // mode 0 is treated first + + // Interpolate E, B + cField2D *El = ( static_cast( EMfields ) )->El_[0]; + cField2D *Er = ( static_cast( EMfields ) )->Er_[0]; + cField2D *Et = ( static_cast( EMfields ) )->Et_[0]; + cField2D *Bl = ( static_cast( EMfields ) )->Bl_m[0]; + cField2D *Br = ( static_cast( EMfields ) )->Br_m[0]; + cField2D *Bt = ( static_cast( EMfields ) )->Bt_m[0]; + cField2D *Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[0]; + cField2D *Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[0]; + + // Interpolation of El^(d,p) + ( *Epart ) [ 0*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] ) ); + // Interpolation of Er^(p,d) + ( *Epart ) [ 1*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] ) ); + // Interpolation of Et^(p,p) + ( *Epart ) [ 2*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] ) ); + // Interpolation of Bl^(p,d) + ( *Bpart ) [ 0*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] ) ); + // Interpolation of Br^(d,p) + ( *Bpart ) [ 1*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] ) ); + // Interpolation of Bt^(d,d) + ( *Bpart ) [ 2*nparts+ipart ] = std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] ) ); + // Interpolation of Phi^(p,p) + ( *PHIpart ) [ 0*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], Phi, idx_p[0], idx_p[1] ) ; + // Interpolation of GradPhil^(p,p) + ( *GradPHIpart ) [ 0*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], GradPhil, idx_p[0], idx_p[1] ) ; + // Interpolation of GradPhir^(p,p) + ( *GradPHIpart ) [ 1*nparts+ipart ] = compute( &coeffxp[1], &coeffyp[0], GradPhir, idx_p[0], idx_p[1] ) ; + // GradPhit = 0 in cylindrical symmetry + ( *GradPHIpart ) [ 2*nparts+ipart ] = 0.; + // BTIS3 fields, in the x direction they are centered as Ey and Ez + // Interpolation of Br^(p,p) for BTIS + ( *BLocyBTIS3)[ 0*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] ) ); + // Interpolation of Bt^(p,d) for BTIS + ( *BLoczBTIS3)[ 0*nparts+ipart ] = std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] ) ); + + if (r > 0){ + exp_m_theta_local = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; + } else { + exp_m_theta_local = 1. ; + } + for( unsigned int imode = 1; imode < nmodes_ ; imode++ ) { + El = ( static_cast( EMfields ) )->El_[imode]; + Er = ( static_cast( EMfields ) )->Er_[imode]; + Et = ( static_cast( EMfields ) )->Et_[imode]; + Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + Br = ( static_cast( EMfields ) )->Br_m[imode]; + Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + Br_BTIS3 = ( static_cast( EMfields ) )->Br_mBTIS3[imode]; + Bt_BTIS3 = ( static_cast( EMfields ) )->Bt_mBTIS3[imode]; + + exp_mm_theta_local *= exp_m_theta_local ; + + ( *Epart ) [ 0*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyp[0], El , idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Epart ) [ 1*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyd[0], Er , idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + ( *Epart ) [ 2*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyp[0], Et , idx_p[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 0*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyd[0], Bl , idx_p[0], idx_d[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 1*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyp[0], Br , idx_d[0], idx_p[1] )* exp_mm_theta_local ) ; + ( *Bpart ) [ 2*nparts+ipart ] += std::real( compute( &coeffxd[1], &coeffyd[0], Bt , idx_d[0], idx_d[1] )* exp_mm_theta_local ) ; + // BTIS3 fields, in the x direction they are centered as Ey and Ez + ( *BLocyBTIS3)[ 0*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyp[0], Br_BTIS3, idx_p[0], idx_p[1] )* exp_mm_theta_local ); + ( *BLoczBTIS3)[ 0*nparts+ipart ] += std::real( compute( &coeffxp[1], &coeffyd[0], Bt_BTIS3, idx_p[0], idx_d[1] )* exp_mm_theta_local ); + } + + // project on x,y,z, remember that GradPhit = 0 in cylindrical symmetry + delta2 = std::real( exp_m_theta_local ) * ( *Epart ) [ 1*nparts+ipart ] + std::imag( exp_m_theta_local ) * ( *Epart ) [ 2*nparts+ipart ]; + ( *Epart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *Epart ) [ 1*nparts+ipart ] + std::real( exp_m_theta_local ) * ( *Epart ) [ 2*nparts+ipart ]; + ( *Epart ) [ 1*nparts+ipart ] = delta2 ; + delta2 = std::real( exp_m_theta_local ) * ( *Bpart ) [ 1*nparts+ipart ] + std::imag( exp_m_theta_local ) * ( *Bpart ) [ 2*nparts+ipart ]; + ( *Bpart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *Bpart ) [ 1*nparts+ipart ] + std::real( exp_m_theta_local ) * ( *Bpart ) [ 2*nparts+ipart ]; + ( *Bpart ) [ 1*nparts+ipart ] = delta2 ; + + delta2 = std::real( exp_m_theta_local ) * ( *GradPHIpart ) [ 1*nparts+ipart ] ; + ( *GradPHIpart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *GradPHIpart ) [ 1*nparts+ipart ] ; + ( *GradPHIpart ) [ 1*nparts+ipart ] = delta2 ; + + delta2 = std::real( exp_m_theta_local ) * ( *BLocyBTIS3) [ 0*nparts+ipart ] + std::imag( exp_m_theta_local ) * ( *BLoczBTIS3) [ 0*nparts+ipart ]; + ( *BLoczBTIS3) [0*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *BLocyBTIS3)[ 0*nparts+ipart ] + std::real( exp_m_theta_local ) * ( *BLoczBTIS3) [ 0*nparts+ipart ]; + ( *BLocyBTIS3) [0*nparts+ipart ] = delta2 ; + + //Buffering of iold and delta + ( *iold )[ipart+0*nparts] = idx_p[0]; + ( *iold )[ipart+1*nparts] = idx_p[1]; + ( *delta )[ipart+0*nparts] = delta_p[0]; + ( *delta )[ipart+1*nparts] = delta_p[1]; + ( *eitheta_old)[ipart] = std::conj(exp_m_theta_local); //exp(i theta) + + } // end ipart loop + + } // end with B-TIS3 interpolation + + +} // END InterpolatorAM1OrderRuyten + +void InterpolatorAM1OrderRuyten::timeCenteredEnvelope( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, int ) +{ + // // Static cast of the envelope fields + // Static cast of the envelope fields + Field2D *Phi_m2Dcyl = static_cast( EMfields->envelope->Phi_m ); + Field2D *GradPhil_m2Dcyl = static_cast( EMfields->envelope->GradPhil_m ); + Field2D *GradPhir_m2Dcyl = static_cast( EMfields->envelope->GradPhir_m ); + + std::vector *PHI_mpart = &( smpi->dynamics_PHI_mpart[ithread] ); + std::vector *GradPHI_mpart = &( smpi->dynamics_GradPHI_mpart[ithread] ); + + std::vector *iold = &( smpi->dynamics_iold[ithread] ); + std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); + std::vector> *eitheta_old = &( smpi->dynamics_eithetaold[ithread] ); + + double r, delta2, xpn, rpn; + + complex exp_m_theta_local ; + //Loop on bin particles + int nparts = particles.numberOfParticles() ; + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + // Normalized particle position + xpn = particles.position( 0, ipart ) * D_inv_[0]; + r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + rpn = r * D_inv_[1]; + + int idx_p[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + + // Compute coefficients + coeffs( xpn, rpn, idx_p, NULL, coeffxp, coeffyp, NULL, NULL, delta_p ); + + // only mode 0 is used + + // ------------------------- + // Interpolation of Phi_m^(p,p) + // ------------------------- + ( *PHI_mpart )[ipart] = compute( &coeffxp[1], &coeffyp[0], Phi_m2Dcyl, idx_p[0], idx_p[1] ); + + // ------------------------- + // Interpolation of GradPhi_m^(p,p), l component + // ------------------------- + ( *GradPHI_mpart )[ipart+0*nparts] = compute( &coeffxp[1], &coeffyp[0], GradPhil_m2Dcyl, idx_p[0], idx_p[1] ); + + // ------------------------- + // Interpolation of GradPhi_m^(p,p), r component + // ------------------------- + ( *GradPHI_mpart )[ipart+1*nparts] = compute( &coeffxp[1], &coeffyp[0], GradPhir_m2Dcyl, idx_p[0], idx_p[1] ); + + // ------------------------- + // Interpolation of GradPhi_m^(p,p), theta component + // ------------------------- + ( *GradPHI_mpart )[ipart+2*nparts] = 0.; // zero with cylindrical symmetry + + if (r > 0){ + exp_m_theta_local = ( particles.position( 1, ipart ) - Icpx * particles.position( 2, ipart ) ) / r ; + } else { + exp_m_theta_local = 1. ; + } + + // project on x,y,z, remember that GradPhit = 0 in cylindrical symmetry + delta2 = std::real( exp_m_theta_local ) * ( *GradPHI_mpart ) [ 1*nparts+ipart ] ; + ( *GradPHI_mpart ) [ 2*nparts+ipart ] = -std::imag( exp_m_theta_local ) * ( *GradPHI_mpart ) [ 1*nparts+ipart ] ; + ( *GradPHI_mpart ) [ 1*nparts+ipart ] = delta2 ; + + //Buffering of iold and delta + // store indices and delta + ( *iold )[0*nparts+ipart] = idx_p[0]; + ( *iold )[1*nparts+ipart] = idx_p[1]; + ( *delta )[0*nparts+ipart] = delta_p[0]; + ( *delta )[1*nparts+ipart] = delta_p[1]; + ( *eitheta_old)[ipart] = std::conj(exp_m_theta_local); //exp(i theta) + + } + +} // END InterpolatorAM1OrderRuyten + + +void InterpolatorAM1OrderRuyten::envelopeAndSusceptibility( ElectroMagn *EMfields, Particles &particles, int ipart, double *Env_A_abs_Loc, double *Env_Chi_Loc, double *Env_E_abs_Loc, double *Env_Ex_abs_Loc ) +{ + // Static cast of the electromagnetic fields + Field2D *Env_A_abs_2Dcyl = static_cast( EMfields->Env_A_abs_ ); + Field2D *Env_Chi_2Dcyl = static_cast( EMfields->Env_Chi_ ); + Field2D *Env_E_abs_2Dcyl = static_cast( EMfields->Env_E_abs_ ); + Field2D *Env_Ex_abs_2Dcyl = static_cast( EMfields->Env_Ex_abs_ ); + + // Normalized particle position + double xpn = particles.position( 0, ipart ) * D_inv_[0]; + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + double rpn = r * D_inv_[1]; + + // Compute coefficients + int idx_p[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + coeffs( xpn, rpn, idx_p, NULL, coeffxp, coeffyp, NULL, NULL, delta_p ); + + // ------------------------- + // Interpolation of Env_A_abs_^(p,p) + // ------------------------- + *( Env_A_abs_Loc ) = compute( &coeffxp[1], &coeffyp[0], Env_A_abs_2Dcyl, idx_p[0], idx_p[1]); + + // ------------------------- + // Interpolation of Env_Chi_^(p,p) + // ------------------------- + *( Env_Chi_Loc ) = compute( &coeffxp[1], &coeffyp[0], Env_Chi_2Dcyl, idx_p[0], idx_p[1]); + + // ------------------------- + // Interpolation of Env_E_abs_^(p,p) + // ------------------------- + *( Env_E_abs_Loc ) = compute( &coeffxp[1], &coeffyp[0], Env_E_abs_2Dcyl, idx_p[0], idx_p[1]); + + // ------------------------- + // Interpolation of Env_Ex_abs_^(p,p) + // ------------------------- + *( Env_Ex_abs_Loc ) = compute( &coeffxp[1], &coeffyp[0], Env_Ex_abs_2Dcyl, idx_p[0], idx_p[1]); + +} // END InterpolatorAM1OrderRuyten + +void InterpolatorAM1OrderRuyten::envelopeFieldForIonization( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, int ) +{ + // Static cast of the envelope fields + Field2D *EnvEabs = static_cast( EMfields->Env_E_abs_ ); + Field2D *EnvExabs = static_cast( EMfields->Env_Ex_abs_ ); + + std::vector *EnvEabs_part = &( smpi->dynamics_EnvEabs_part[ithread] ); + std::vector *EnvExabs_part = &( smpi->dynamics_EnvExabs_part[ithread] ); + + double xpn,rpn,r; + + //Loop on bin particles + for( int ipart=*istart ; ipart<*iend; ipart++ ) { + + // Normalized particle position + xpn = particles.position( 0, ipart ) * D_inv_[0]; + r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ) ; + rpn = r * D_inv_[1]; + + // Compute coefficients + int idx_p[2]; + double delta_p[2]; + double coeffxp[3], coeffyp[2]; + coeffs( xpn, rpn, idx_p, NULL, coeffxp, coeffyp, NULL, NULL, delta_p ); + + // only mode 0 is used + + // --------------------------------- + // Interpolation of Env_E_abs^(p,p) + // --------------------------------- + ( *EnvEabs_part )[ipart] = compute( &coeffxp[1], &coeffyp[0], EnvEabs, idx_p[0], idx_p[1] ); + + // --------------------------------- + // Interpolation of Env_Ex_abs^(p,p) + // --------------------------------- + ( *EnvExabs_part )[ipart] = compute( &coeffxp[1], &coeffyp[0], EnvExabs, idx_p[0], idx_p[1] ); + } + +} // END InterpolatorAM1OrderRuyten diff --git a/src/Interpolator/InterpolatorAM1OrderRuyten.h b/src/Interpolator/InterpolatorAM1OrderRuyten.h new file mode 100755 index 000000000..69740dd31 --- /dev/null +++ b/src/Interpolator/InterpolatorAM1OrderRuyten.h @@ -0,0 +1,118 @@ +#ifndef INTERPOLATORAM1ORDERRUYTEN_H +#define INTERPOLATORAM1ORDERRUYTEN_H + + +#include "InterpolatorAM.h" +#include "cField2D.h" +#include "Field2D.h" + + +// -------------------------------------------------------------------------------------------------------------------- +//! Class for 1st order interpolator with Ruyten correction for AM simulations based on the article from W. Ruyten: Density-Conserving Shape Factors for Particle Simulations in Cylindrical and Spherical Coordinates J. Comp. Phys. 105, 224-232 (1993) +// -------------------------------------------------------------------------------------------------------------------- +class InterpolatorAM1OrderRuyten final : public InterpolatorAM +{ +public: + InterpolatorAM1OrderRuyten( Params &, Patch * ); + ~InterpolatorAM1OrderRuyten() override final {}; + + inline void __attribute__((always_inline)) fields( ElectroMagn *EMfields, Particles &particles, int ipart, int nparts, double *ELoc, double *BLoc ); + void fieldsAndCurrents( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, LocalFields *JLoc, double *RhoLoc ) override final ; + void fieldsWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, unsigned int scell = 0, int ipart_ref = 0 ) override final ; + void fieldsSelection( ElectroMagn *EMfields, Particles &particles, double *buffer, int offset, std::vector *selection ) override final; + void oneField( Field **field, Particles &particles, int *istart, int *iend, double *FieldLoc, double *l1=NULL, double *l2=NULL, double *l3=NULL ) override final; + + void fieldsAndEnvelope( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, int ipart_ref = 0 ) override final; + void timeCenteredEnvelope( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, int ipart_ref = 0 ) override final; + void envelopeAndSusceptibility( ElectroMagn *EMfields, Particles &particles, int ipart, double *Env_A_abs_Loc, double *Env_Chi_Loc, double *Env_E_abs_Loc, double *Env_Ex_abs_Loc ) override final; + void envelopeFieldForIonization( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, int ipart_ref = 0 ) override final; + + inline std::complex __attribute__((always_inline)) compute( double *coeffx, double *coeffy, cField2D *f, int idx, int idy ) + { + std::complex interp_res( 0. ); + for( int iloc=-1 ; iloc<2 ; iloc++ ) { + for( int jloc=0 ; jloc<2 ; jloc++ ) { + //std::cout << "Er standard idx = " << idx << " idy = " << idy << " iloc = " << iloc << " jloc = " << jloc << " " << *( coeffx+iloc ) << " " << *( coeffy+jloc ) << " " << ( *f )( idx+iloc, idy+jloc ) << std::endl; + interp_res += *( coeffx+iloc ) * *( coeffy+jloc ) * ( ( *f )( idx+iloc, idy+jloc ) ) ; + + //std::cout<<"f "< 0.) ? ( double )idx_d[1] - 0.5 : 0. ); + //Ruyten scheme + coeffyd[0] = 0.5*(1. - delta)*(5*x_n + 2 - ypn)/(x_np1*x_np1-x_n*x_n); + coeffyd[1] = 1. - coeffyd[0]; //conservation de la charge + + idx_d[0] = idx_d[0] - i_domain_begin_; + idx_d[1] = idx_d[1] - j_domain_begin_; + } + } + + // exp m theta + std::complex exp_m_theta_; + //! Number of modes; + unsigned int nmodes_; + +};//END class + +#endif diff --git a/src/Interpolator/InterpolatorAM1OrderRuytenV.cpp b/src/Interpolator/InterpolatorAM1OrderRuytenV.cpp new file mode 100755 index 000000000..108ef27d8 --- /dev/null +++ b/src/Interpolator/InterpolatorAM1OrderRuytenV.cpp @@ -0,0 +1,323 @@ +#include "InterpolatorAM1OrderRuytenV.h" + +#include +#include +#include +#include "ElectroMagn.h" +#include "ElectroMagnAM.h" +#include "cField2D.h" +#include "Particles.h" +#include "LaserEnvelope.h" +#include +#include "dcomplex.h" + +using namespace std; + + +// --------------------------------------------------------------------------------------------------------------------- +// Creator for InterpolatorAM1OrderRuytenV +// --------------------------------------------------------------------------------------------------------------------- +InterpolatorAM1OrderRuytenV::InterpolatorAM1OrderRuytenV( Params ¶ms, Patch *patch ) : InterpolatorAM( patch ) +{ + + nmodes_ = params.nmodes; + D_inv_[0] = 1.0/params.cell_length[0]; + D_inv_[1] = 1.0/params.cell_length[1]; + nscellr_ = params.patch_size_[1] + 1; + oversize_[0] = params.oversize[0]; + oversize_[1] = params.oversize[1]; +} + +void InterpolatorAM1OrderRuytenV::fieldsWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int *istart, int *iend, int ithread, unsigned int scell, int ipart_ref ) +{ + if( istart[0] == iend[0] ) { + return; //Don't treat empty cells. + } + + int nparts( ( smpi->dynamics_invgf[ithread] ).size() ); + + double * __restrict__ Epart[3]; + double * __restrict__ Bpart[3]; + + double * __restrict__ position_x = particles.getPtrPosition(0); + double * __restrict__ position_y = particles.getPtrPosition(1); + double * __restrict__ position_z = particles.getPtrPosition(2); + + + double * __restrict__ deltaO[2]; //Delta is the distance of the particle from its primal node in cell size. Delta is in [-0.5, +0.5[ + std::complex * __restrict__ eitheta_old; //eithetaold stores exp(i theta) of the particle before pusher. + + int idx[2], idxO[2]; + //Primal indices are constant over the all cell + idx[0] = scell/nscellr_+oversize_[0]+i_domain_begin_; + idxO[0] = idx[0] - i_domain_begin_ -1 ; + idx[1] = ( scell%nscellr_ )+oversize_[1]+j_domain_begin_; + idxO[1] = idx[1] - j_domain_begin_ -1 ; + + double coeff[2][2][3][32]; + double dual[1][32]; // Size ndim. Boolean converted into double indicating if the part has a dual indice equal to the primal one (dual=0) or if it is +1 (dual=1). + + int vecSize = 32; + double delta, delta2; + + + int cell_nparts( ( int )iend[0]-( int )istart[0] ); + + std::vector> exp_m_theta_( vecSize), exp_mm_theta( vecSize) ; //exp(-i theta), exp(-i m theta) + + //Loop on groups of vecSize particles + for( int ivect=0 ; ivect < cell_nparts; ivect += vecSize ) { + + int np_computed( min( cell_nparts-ivect, vecSize ) ); + deltaO[0] = &( smpi->dynamics_deltaold[ithread][0 + ivect + istart[0] - ipart_ref] ); + deltaO[1] = &( smpi->dynamics_deltaold[ithread][nparts + ivect + istart[0] - ipart_ref] ); + eitheta_old = &( smpi->dynamics_eithetaold[ithread][ ivect + istart[0] - ipart_ref] ); + + + #pragma omp simd private(delta2, delta) + for( int ipart=0 ; ipart X + // j=0 primal + delta = position_x[ipart2]*D_inv_[0] - (double)idx[0]; + delta2 = delta*delta; + coeff[0][0][0][ipart] = 0.5 * ( delta2-delta+0.25 ); + coeff[0][0][1][ipart] = ( 0.75 - delta2 ); + coeff[0][0][2][ipart] = 0.5 * ( delta2+delta+0.25 ); + deltaO[0][ipart] = delta; + + // j=1 dual + dual [0][ipart] = ( delta >= 0. ); + //delta dual = distance to dual node + delta = delta - dual[0][ipart] + 0.5 ; + delta2 = delta*delta; + coeff[0][1][0][ipart] = 0.5 * ( delta2-delta+0.25 ); + coeff[0][1][1][ipart] = ( 0.75 - delta2 ); + coeff[0][1][2][ipart] = 0.5 * ( delta2+delta+0.25 ); + + // i= 1 ==> Y + // j=0 primal + + + double x_n, ypn; + + ypn = r * D_inv_[1]; + delta = ypn - (double)idx[1]; + int rshift = (delta >= 0); + x_n = (double)idx[1]-1. + (double)rshift; + + coeff[1][0][rshift][ipart] = (x_n+1-ypn)*(5*x_n + 2 - ypn)/(4.*x_n + 2.); + coeff[1][0][rshift+1][ipart] = 1. - coeff[1][0][rshift][ipart]; + coeff[1][0][(rshift+2)%3][ipart] = 0.; + deltaO[1][ipart] = delta; + + // j=1 dual + double x_np1 = ( double )idx[1] + 0.5 ; // This is known thanks to the cell sorting + x_n = ( (( double )idx[1] - 0.5 > 0.) ? ( double )idx[1] - 0.5 : 0. ); + + coeff[1][1][0][ipart] = (x_np1-ypn)*(5*x_n + 2 - ypn)/(4.*x_n + 2.); + coeff[1][1][1][ipart] = 1. - coeff[1][0][0][ipart]; + //coeff[1][0][2][ipart] = 0.; // This coefficient is not supposed to be used + } + + double interp_res; + double * __restrict__ coeffld = &( coeff[0][1][1][0] ); + double * __restrict__ coefflp = &( coeff[0][0][1][0] ); + double * __restrict__ coeffrp = &( coeff[1][0][1][0] ); + double * __restrict__ coeffrd = &( coeff[1][1][1][0] ); + + for( unsigned int k=0; k<3; k++ ) { + Epart[k]= &( smpi->dynamics_Epart[ithread][k*nparts-ipart_ref+ivect+istart[0]] ); + Bpart[k]= &( smpi->dynamics_Bpart[ithread][k*nparts-ipart_ref+ivect+istart[0]] ); + #pragma omp simd + for( int ipart=0 ; ipart field_buffer[4][4]; + + for( unsigned int imode = 0; imode < nmodes_ ; imode++ ) { + // Static cast of the electromagnetic fields + cField2D * __restrict__ El = ( static_cast( EMfields ) )->El_[imode]; + cField2D * __restrict__ Er = ( static_cast( EMfields ) )->Er_[imode]; + cField2D * __restrict__ Et = ( static_cast( EMfields ) )->Et_[imode]; + cField2D * __restrict__ Bl = ( static_cast( EMfields ) )->Bl_m[imode]; + cField2D * __restrict__ Br = ( static_cast( EMfields ) )->Br_m[imode]; + cField2D * __restrict__ Bt = ( static_cast( EMfields ) )->Bt_m[imode]; + + + + // Field buffers for vectorization (required on A64FX) + for( int iloc=-1 ; iloc<3 ; iloc++ ) { + for( int jloc=-1 ; jloc<2 ; jloc++ ) { + field_buffer[iloc+1][jloc+1] = ( *El )( idxO[0]+1+iloc, idxO[1]+1+jloc ); + } + } + + #pragma omp simd private(interp_res) + for( int ipart=0 ; ipart * ) override final {}; + void oneField( Field **, Particles &, int *, int *, double *, double * = NULL, double * = NULL, double * = NULL ) override final {}; + + void fieldsAndEnvelope( ElectroMagn *, Particles &, SmileiMPI *, int *, int *, int, int = 0 ) override final {}; + void timeCenteredEnvelope( ElectroMagn *, Particles &, SmileiMPI *, int *, int *, int, int = 0 ) override final {}; + void envelopeAndSusceptibility( ElectroMagn *, Particles &, int , double *, double *, double *, double * ) override final {}; + void envelopeFieldForIonization( ElectroMagn *, Particles &, SmileiMPI *, int *, int *, int , int = 0 ) override final {}; + + +private: + + //! Number of modes; + unsigned int nmodes_; + + +};//END class + +#endif diff --git a/src/Interpolator/InterpolatorFactory.h b/src/Interpolator/InterpolatorFactory.h index 37e1042fb..e3c4b17d8 100755 --- a/src/Interpolator/InterpolatorFactory.h +++ b/src/Interpolator/InterpolatorFactory.h @@ -10,6 +10,7 @@ #include "Interpolator3D2Order.h" #include "Interpolator3D4Order.h" #include "InterpolatorAM1Order.h" +#include "InterpolatorAM1OrderRuyten.h" #include "InterpolatorAM2Order.h" #include "Interpolator1DWT2Order.h" #include "Interpolator1DWT4Order.h" @@ -24,6 +25,7 @@ #include "Interpolator3D2OrderV.h" #include "Interpolator3D4OrderV.h" #include "InterpolatorAM2OrderV.h" +#include "InterpolatorAM1OrderRuytenV.h" #include "Interpolator1DWT2OrderV.h" #include "Interpolator2DWT2OrderV.h" @@ -154,12 +156,20 @@ class InterpolatorFactory else if( params.geometry == "AMcylindrical" ) { if ( !params.is_spectral){ if( !vectorization ) { - Interp = new InterpolatorAM2Order( params, patch ); + if ( params.interpolation_order == 1 ) { + Interp = new InterpolatorAM1OrderRuyten( params, patch ); + } else { + Interp = new InterpolatorAM2Order( params, patch ); // default + } } else { - Interp = new InterpolatorAM2OrderV( params, patch ); + if ( params.interpolation_order == 1 ) { + Interp = new InterpolatorAM1OrderRuytenV( params, patch ); + } else { + Interp = new InterpolatorAM2OrderV( params, patch ); + } } - } else { + } else { // Spectral Interp = new InterpolatorAM1Order( params, patch ); } } diff --git a/src/Params/Params.cpp b/src/Params/Params.cpp index 69973d104..27e743fd1 100755 --- a/src/Params/Params.cpp +++ b/src/Params/Params.cpp @@ -282,10 +282,6 @@ Params::Params( SmileiMPI *smpi, std::vector namelistsFiles ) : ERROR_NAMELIST( "Main.interpolation_order " << interpolation_order << " should be 1 for PSATD solver", LINK_NAMELIST + std::string("#main-variables") ); } - if( interpolation_order != 2 && !is_spectral ){ - ERROR_NAMELIST( "Main.interpolation_order " << interpolation_order << " should be 2 for FDTD solver.", - LINK_NAMELIST + std::string("#main-variables")); - } } else if( interpolation_order!=2 && interpolation_order!=4 && !is_spectral ) { ERROR_NAMELIST( "Main.interpolation_order " << interpolation_order << " should be 2 or 4", LINK_NAMELIST + std::string("#main-variables")); @@ -554,8 +550,12 @@ Params::Params( SmileiMPI *smpi, std::vector namelistsFiles ) : // This method is detailed in P.-L. Bourgeois and X. Davoine (2023) https://doi.org/10.1017/S0022377823000223 use_BTIS3 = false; PyTools::extract( "use_BTIS3_interpolation", use_BTIS3, "Main" ); - if (use_BTIS3 && interpolation_order != 2 ){ - ERROR("B-TIS3 interpolation implemented only at order 2."); + if (use_BTIS3 && interpolation_order != 2 && (interpolation_order != 1 || geometry != "AMcylindrical" || is_spectral==true )){ + if (geometry=="AMcylindrical"){ + ERROR("B-TIS3 interpolation is not implemented for PSATD solver."); + } else { + ERROR("B-TIS3 interpolation is implemented only at order 2 for Cartesian geometries."); + } } // Current filter properties diff --git a/src/Projector/ProjectorAM1OrderRuyten.cpp b/src/Projector/ProjectorAM1OrderRuyten.cpp new file mode 100755 index 000000000..b3db66ba8 --- /dev/null +++ b/src/Projector/ProjectorAM1OrderRuyten.cpp @@ -0,0 +1,1133 @@ +#include "ProjectorAM1OrderRuyten.h" + +#include +#include +#include +#include "dcomplex.h" +#include "ElectroMagnAM.h" +#include "cField2D.h" +#include "Particles.h" +#include "Tools.h" +#include "Patch.h" +#include "PatchAM.h" + +using namespace std; + + +// --------------------------------------------------------------------------------------------------------------------- +// Constructor for ProjectorAM1OrderRuyten +// --------------------------------------------------------------------------------------------------------------------- +ProjectorAM1OrderRuyten::ProjectorAM1OrderRuyten( Params ¶ms, Patch *patch ) : ProjectorAM( params, patch ) +{ + dt = params.timestep; + dr = params.cell_length[1]; + dl_inv_ = 1.0/params.cell_length[0]; + dl_ov_dt_ = params.cell_length[0] / params.timestep; + dr_ov_dt_ = params.cell_length[1] / params.timestep; + dr_inv_ = 1.0 / dr; + one_ov_dt = 1.0 / params.timestep; + Nmode_=params.nmodes; + i_domain_begin_ = patch->getCellStartingGlobalIndex( 0 ); + j_domain_begin_ = patch->getCellStartingGlobalIndex( 1 ); + + nprimr_ = params.patch_size_[1] + 2*params.oversize[1] + 1; + npriml_ = params.patch_size_[0] + 2*params.oversize[0] + 1; + + invR_ = &((static_cast( patch )->invR)[0]); + invRd_ = &((static_cast( patch )->invRd)[0]); + + dts2 = params.timestep/2.; + dts4 = params.timestep/4.; +} + + +// --------------------------------------------------------------------------------------------------------------------- +// Destructor for ProjectorAM1OrderRuyten +// --------------------------------------------------------------------------------------------------------------------- +ProjectorAM1OrderRuyten::~ProjectorAM1OrderRuyten() +{ +} + +// --------------------------------------------------------------------------------------------------------------------- +//! Project local currents for all modes +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuyten::currents( ElectroMagnAM *emAM, + Particles &particles, + unsigned int ipart, + double invgf, + int *iold, + double *deltaold, + std::complex *array_eitheta_old, + bool diag_flag, int ispec) +{ + + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- int iloc, + int nparts= particles.size(); + int iloc, jloc, linindex; + // (x,y,z) components of the current density for the macro-particle + double charge_weight = inv_cell_volume * ( double )( particles.charge( ipart ) )*particles.weight( ipart ); + double crl_p = charge_weight*dl_ov_dt_; + double crr_p = charge_weight*one_ov_dt; + + // variable declaration + double xpn, ypn; + double delta, delta2; + // arrays used for the Esirkepov projection method + double Sl0[5], Sl1[5], Sr0[4], Sr1[4], DSl[5], DSr[4]; + complex Jl_p[5], Jr_p[4]; + complex e_delta, e_delta_m1, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; + complex *Jl, *Jr, *Jt, *rho; + + for( unsigned int i=0; i<5; i++ ) { + Sl1[i] = 0.; + Sr1[i] = 0.; + } + Sl0[0] = 0.; + Sl0[4] = 0.; + Sr0[0] = 0.; + Sr0[3] = 0.; + // -------------------------------------------------------- + // Locate particles & Calculate Esirkepov coef. S, DS and W + // -------------------------------------------------------- + + // locate the particle on the primal grid at former time-step & calculate coeff. S0 + delta = deltaold[0*nparts]; + delta2 = delta*delta; + Sl0[1] = 0.5 * ( delta2-delta+0.25 ); + Sl0[2] = 0.75-delta2; + Sl0[3] = 0.5 * ( delta2+delta+0.25 ); + + delta = deltaold[1*nparts]; + + double x_n = iold[1*nparts]+j_domain_begin_ ; + // Ruyten scheme + Sr0[1] = (1. - delta)*(4*x_n + 2 - delta)/(4.*x_n + 2.); + Sr0[2] = 1. - Sr0[1]; //conservation de la charge + + //calculate exponential coefficients + + double rp = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); + std::complex eitheta_old = array_eitheta_old[0]; + std::complex eitheta = ( particles.position( 1, ipart ) + Icpx * particles.position( 2, ipart ) ) / rp ; //exp(i theta) + e_bar = 1.; + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn ); + int ipo = iold[0*nparts]; + int ip_m_ipo = ip-ipo-i_domain_begin_; + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[ip_m_ipo+1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[ip_m_ipo+2] = 0.75-delta2; + Sl1[ip_m_ipo+3] = 0.5 * ( delta2+delta+0.25 ); + + ypn = rp *dr_inv_ ; + int jp = floor( ypn ); + int jpo = iold[1*nparts]; + int jp_m_jpo = jp-jpo-j_domain_begin_; + delta = ypn - ( double )jp; + x_n = jp ; + // Ruyten scheme + Sr1[jp_m_jpo+1] = (1. - delta)*(5*x_n + 2 - ypn)/(4.*x_n + 2.); + Sr1[jp_m_jpo+2] = 1. - Sr1[jp_m_jpo+1]; //conservation de la charge + + for( unsigned int i=0; i < 5; i++ ) { + DSl[i] = Sl1[i] - Sl0[i]; + } + for( unsigned int i=0; i < 4; i++ ) { + DSr[i] = Sr1[i] - Sr0[i]; + } + + double r_bar = ((jpo + j_domain_begin_ + deltaold[1*nparts])*dr + rp) * 0.5; // r at t = t0 - dt/2 + + e_delta_m1 = std::sqrt(eitheta * std::conj(eitheta_old)); //ei^(theta-theta_old)/2. std::sqrt keeps the root with positive real part which is what we need here. + e_bar_m1 = eitheta_old * e_delta_m1; //ei^(theta+theta_old)/2. + + ipo -= 2; //This minus 2 come from the order 2 scheme, based on a 5 points stencil from -2 to +2. + // i/j/kpo stored with - i/j/k_domain_begin in Interpolator + jpo -= 1; // Order 1 in R: 4 points stencil from -1 to +2. + + double *invR__local = &(invR_[jpo]); + + // ------------------------------------------------ + // Local current created by the particle + // calculate using the charge conservation equation + // ------------------------------------------------ + + // --------------------------- + // Calculate the total current + // --------------------------- + + //initial value of crt_p for imode = 0. + complex crt_p= charge_weight*( particles.momentum( 2, ipart )* real(e_bar_m1) - particles.momentum( 1, ipart )*imag(e_bar_m1) ) * invgf; + + // Compute everything independent of theta + double tmpJl[4]; + for( unsigned int j=0 ; j<4 ; j++ ) { + tmpJl[j] = crl_p * ( Sr0[j] + 0.5*DSr[j] )* invR__local[j]; + } + Jl_p[0]= 0.; + for( unsigned int i=1 ; i<5 ; i++ ) { + Jl_p[i]= Jl_p[i-1] - DSl[i-1]; + } + + + double Vd[3]; + double tmpJr[3]; + for( int j=2 ; j>=0 ; j-- ) { + jloc = j+jpo+1; + Vd[j] = abs( jloc + j_domain_begin_ + 0.5 )* invRd_[jloc]*dr ; + tmpJr[j] = crr_p * DSr[j+1] * invRd_[jpo+j+1]*dr; + } + Jr_p[3]= 0.; + for( int j=2 ; j>=0 ; j-- ) { + Jr_p[j] = Jr_p[j+1] * Vd[j] + tmpJr[j]; + } + + e_delta = 0.5; + + //Compute division by R in advance for Jt and rho evaluation. + for( unsigned int j=0 ; j<4 ; j++ ) { + Sr0[j] *= invR__local[j]; + Sr1[j] *= invR__local[j]; + } + + for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { + + if (imode > 0){ + e_delta *= e_delta_m1; + e_bar *= e_bar_m1; + C_m = 2. * e_bar ; //multiply modes > 0 by 2 and C_m = 1 otherwise. + crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.*r_bar; + } + + // Add contribution J_p to global array + if (!diag_flag){ + Jl = &( *emAM->Jl_[imode] )( 0 ); + Jr = &( *emAM->Jr_[imode] )( 0 ); + Jt = &( *emAM->Jt_[imode] )( 0 ); + } else { + unsigned int n_species = emAM->Jl_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec; + Jl = emAM->Jl_s [ifield] ? &( * ( emAM->Jl_s [ifield] ) )( 0 ) : &( *emAM->Jl_ [imode] )( 0 ) ; + Jr = emAM->Jr_s [ifield] ? &( * ( emAM->Jr_s [ifield] ) )( 0 ) : &( *emAM->Jr_ [imode] )( 0 ) ; + Jt = emAM->Jt_s [ifield] ? &( * ( emAM->Jt_s [ifield] ) )( 0 ) : &( *emAM->Jt_ [imode] )( 0 ) ; + rho = emAM->rho_AM_s[ifield] ? &( * ( emAM->rho_AM_s[ifield] ) )( 0 ) : &( *emAM->rho_AM_[imode] )( 0 ) ; + + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_; + for( unsigned int j=0 ; j<4 ; j++ ) { + jloc = j+jpo; + linindex = iloc+jloc; + rho [linindex] += C_m*charge_weight* Sl1[i]*Sr1[j]; + } + }//i + } + + // Jl^(d,p) + for( unsigned int i=1 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_+jpo; + for( unsigned int j=0 ; j<4 ; j++ ) { + linindex = iloc+j; + Jl [linindex] += C_m * Jl_p[i]*tmpJl[j] ; + } + }//i + + // Jr^(p,d) + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*( nprimr_+1 )+jpo+1; + for( unsigned int j=0 ; j<3 ; j++ ) { + linindex = iloc+j; + Jr [linindex] += C_m * ( Sl0[i] + 0.5*DSl[i] ) * Jr_p[j] ; + } + }//i + + // Jt^(p,p) + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_ + jpo; + for( unsigned int j=0 ; j<4 ; j++ ) { + linindex = iloc+j; + Jt [linindex] -= crt_p*(Sr1[j]*Sl1[i]*( e_delta-1. ) - Sr0[j]*Sl0[i]*(std::conj(e_delta) - 1.)); + } + } + + if (imode == 0) e_delta = 1. ; //Restore e_delta correct initial value. + }// end loop on modes + +} // END Project local current densities (Jl, Jr, Jt, sort) + +// --------------------------------------------------------------------------------------------------------------------- +//! Project for diags and frozen species - +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuyten::basicForComplex( complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode ) +{ + //Warning : this function is not charge conserving. + // This function also assumes that particles position is evaluated at the same time as currents which is usually not true (half time-step difference). + // It will therefore fail to evaluate the current accurately at t=0 if a plasma is already in the box. + + + + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + int iloc, nr( nprimr_ ); + double charge_weight = inv_cell_volume * ( double )( particles.charge( ipart ) )*particles.weight( ipart ); + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); + + if( type > 0 ) { //if current density + charge_weight *= 1./sqrt( 1.0 + particles.momentum( 0, ipart )*particles.momentum( 0, ipart ) + + particles.momentum( 1, ipart )*particles.momentum( 1, ipart ) + + particles.momentum( 2, ipart )*particles.momentum( 2, ipart ) ); + if( type == 1 ) { //if Jl + charge_weight *= particles.momentum( 0, ipart ); + } else if( type == 2 ) { //if Jr + charge_weight *= ( particles.momentum( 1, ipart )*particles.position( 1, ipart ) + particles.momentum( 2, ipart )*particles.position( 2, ipart ) )/ r ; + nr++; + } else { //if Jt + charge_weight *= ( -particles.momentum( 1, ipart )*particles.position( 2, ipart ) + particles.momentum( 2, ipart )*particles.position( 1, ipart ) ) / r ; + } + } + + complex e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; + complex C_m = 1.; + if( imode > 0 ) { + C_m = 2.; + } + for( unsigned int i=0; i<( unsigned int )imode; i++ ) { + C_m *= e_theta; + } + + double xpn, ypn; + double delta, delta2; + double Sl1[5], Sr1[5]; + + // -------------------------------------------------------- + // Locate particles & Calculate Esirkepov coef. S, DS and W + // -------------------------------------------------------- + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn + 0.5 * ( type==1 ) ); + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[2] = 0.75-delta2; + Sl1[3] = 0.5 * ( delta2+delta+0.25 ); + ypn = r * dr_inv_ ; + int jp = round( ypn + 0.5*( type==2 ) ); + delta = ypn - ( double )jp; + delta2 = delta*delta; + Sr1[1] = 0.5 * ( delta2-delta+0.25 ); + Sr1[2] = 0.75-delta2; + Sr1[3] = 0.5 * ( delta2+delta+0.25 ); + + // --------------------------- + // Calculate the total charge + // --------------------------- + ip -= i_domain_begin_ + 2; + jp -= j_domain_begin_ + 2; + + if( type != 2 ) { + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + rhoj [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invR_[j+jp]; + } + }//i + } else { + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + rhoj [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invRd_[j+jp]; + } + }//i + } +} // END Project for diags local current densities + +// --------------------------------------------------------------------------------------------------------------------- +//! Project for diags and frozen species - +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuyten::basicForComplexOnBuffer( complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode, int bdim0, int bin_shift ) +{ + //Warning : this function is not charge conserving. + // This function also assumes that particles position is evaluated at the same time as currents which is usually not true (half time-step difference). + // It will therefore fail to evaluate the current accurately at t=0 if a plasma is already in the box. + + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + int iloc, nr( nprimr_ ); + double charge_weight = inv_cell_volume * ( double )( particles.charge( ipart ) )*particles.weight( ipart ); + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); + + if( type > 0 ) { //if current density + ERROR("This projector can be used only for charge density at the moment."); + } + + complex e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; + complex C_m = 1.; + if( imode > 0 ) { + C_m = 2.; + } + for( unsigned int i=0; i<( unsigned int )imode; i++ ) { + C_m *= e_theta; + } + + double xpn, ypn; + double delta, delta2; + double Sl1[5], Sr1[5]; + + // -------------------------------------------------------- + // Locate particles & Calculate Esirkepov coef. S, DS and W + // -------------------------------------------------------- + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn + 0.5 * ( type==1 ) ); + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[2] = 0.75-delta2; + Sl1[3] = 0.5 * ( delta2+delta+0.25 ); + ypn = r * dr_inv_ ; + int jp = round( ypn + 0.5*( type==2 ) ); + delta = ypn - ( double )jp; + delta2 = delta*delta; + Sr1[1] = 0.5 * ( delta2-delta+0.25 ); + Sr1[2] = 0.75-delta2; + Sr1[3] = 0.5 * ( delta2+delta+0.25 ); + + // --------------------------- + // Calculate the total charge + // --------------------------- + ip -= i_domain_begin_ + 2 + bin_shift; + jp -= j_domain_begin_ + 2; + + // complex *rho = &(rhoj[ imode*(bdim0*nr ) ] ) ; + + int mode_shift = imode*(bdim0*nr ); + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + rhoj [mode_shift+iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invR_[j+jp]; + } + }//i +} // END Project for diags local current densities + +// Apply boundary conditions on axis for currents and densities +void ProjectorAM1OrderRuyten::axisBC(ElectroMagnAM *emAM, bool diag_flag ) +{ + + for (unsigned int imode=0; imode < Nmode_; imode++){ + + std::complex *rhoj = &( *emAM->rho_AM_[imode] )( 0 ); + std::complex *Jl = &( *emAM->Jl_[imode] )( 0 ); + std::complex *Jr = &( *emAM->Jr_[imode] )( 0 ); + std::complex *Jt = &( *emAM->Jt_[imode] )( 0 ); + + apply_axisBC(rhoj, Jl, Jr, Jt, imode, diag_flag); + } + + if (diag_flag){ + unsigned int n_species = emAM->Jl_s.size() / Nmode_; + for( unsigned int imode = 0 ; imode < emAM->Jl_.size() ; imode++ ) { + for( unsigned int ispec = 0 ; ispec < n_species ; ispec++ ) { + unsigned int ifield = imode*n_species+ispec; + complex *Jl = emAM->Jl_s [ifield] ? &( * ( emAM->Jl_s [ifield] ) )( 0 ) : NULL ; + complex *Jr = emAM->Jr_s [ifield] ? &( * ( emAM->Jr_s [ifield] ) )( 0 ) : NULL ; + complex *Jt = emAM->Jt_s [ifield] ? &( * ( emAM->Jt_s [ifield] ) )( 0 ) : NULL ; + complex *rho = emAM->rho_AM_s[ifield] ? &( * ( emAM->rho_AM_s[ifield] ) )( 0 ) : NULL ; + apply_axisBC( rho , Jl, Jr, Jt, imode, diag_flag ); + } + } + } +} + +void ProjectorAM1OrderRuyten::apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ) +{ + // Mode 0 contribution "below axis" is added. + // Non zero modes are substracted because a particle sitting exactly on axis has a non defined theta and can not contribute to a theta dependent mode. + double sign = (imode == 0) ? 1 : -1 ; + + if (diag_flag && rhoj) { + for( unsigned int i=2 ; i 0){ + rhoj[i] = 0.; + rhoj[i-1] = - rhoj[i+1]; // Zero Jl mode > 0 on axis. + } else { + rhoj[i] = rhoj[i+1]; //This smoothing is just for cosmetics on the picture, rho has no influence on the results. + rhoj[i-1] = rhoj[i+1]; // Non zero Jl mode > 0 on axis. + } + } + } + + if (Jl) { + for( unsigned int i=2 ; i<(npriml_+1)*nprimr_+2; i+=nprimr_ ) { + if (imode > 0){ + Jl [i] = 0. ; + Jl[i-1] = -Jl[i+1]; // Zero Jl mode > 0 on axis. + } else { + //Jl mode 0 on axis should be left as is. It looks over estimated but it might be necessary to conserve a correct divergence and a proper evaluation on the field on axis. + Jl [i-1] = Jl [i+1] ; // Non zero Jl mode 0 on axis. + } + } + } + + if (Jt && Jr) { + for( unsigned int i=0 ; i( Jl ); + cField2D *JrAM = static_cast( Jr ); + cField2D *JtAM = static_cast( Jt ); + + + //Declaration of local variables + int ip, id, jp, jd; + double xpn, xpmxip, xpmxip2, xpmxid, xpmxid2; + double ypn, ypmyjp, ypmyjp2, ypmyjd, ypmyjd2; + double Slp[3], Sld[3], Srp[3], Srd[3]; + + // weighted currents + double weight = inv_cell_volume * particles.weight( ipart ); + double Jl_ion = Jion.x * weight; + double Jr_ion = Jion.y * weight; + double Jt_ion = Jion.z * weight; + + //Locate particle on the grid + xpn = particles.position( 0, ipart ) * dl_inv_; // normalized distance to the first node + ypn = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) )*dr_inv_ ; + // x-primal index + ip = round( xpn ); // x-index of the central node + xpmxip = xpn - ( double )ip; // normalized distance to the nearest grid point + xpmxip2 = xpmxip*xpmxip; // square of the normalized distance to the nearest grid point + + // x-dual index + id = round( xpn+0.5 ); // x-index of the central node + xpmxid = xpn - ( double )id + 0.5; // normalized distance to the nearest grid point + xpmxid2 = xpmxid*xpmxid; // square of the normalized distance to the nearest grid point + + // y-primal index + jp = round( ypn ); // y-index of the central node + ypmyjp = ypn - ( double )jp; // normalized distance to the nearest grid point + ypmyjp2 = ypmyjp*ypmyjp; // square of the normalized distance to the nearest grid point + + // y-dual index + jd = round( ypn+0.5 ); // y-index of the central node + ypmyjd = ypn - ( double )jd + 0.5; // normalized distance to the nearest grid point + ypmyjd2 = ypmyjd*ypmyjd; // square of the normalized distance to the nearest grid point + + Slp[0] = 0.5 * ( xpmxip2-xpmxip+0.25 ); + Slp[1] = ( 0.75-xpmxip2 ); + Slp[2] = 0.5 * ( xpmxip2+xpmxip+0.25 ); + + Sld[0] = 0.5 * ( xpmxid2-xpmxid+0.25 ); + Sld[1] = ( 0.75-xpmxid2 ); + Sld[2] = 0.5 * ( xpmxid2+xpmxid+0.25 ); + + Srp[0] = 0.5 * ( ypmyjp2-ypmyjp+0.25 ); + Srp[1] = ( 0.75-ypmyjp2 ); + Srp[2] = 0.5 * ( ypmyjp2+ypmyjp+0.25 ); + + Srd[0] = 0.5 * ( ypmyjd2-ypmyjd+0.25 ); + Srd[1] = ( 0.75-ypmyjd2 ); + Srd[2] = 0.5 * ( ypmyjd2+ypmyjd+0.25 ); + + ip -= i_domain_begin_; + id -= i_domain_begin_; + jp -= j_domain_begin_; + jd -= j_domain_begin_; + + for( unsigned int i=0 ; i<3 ; i++ ) { + //int iploc=ip+i-1; + int idloc=id+i-1; + for( unsigned int j=0 ; j<3 ; j++ ) { + int jploc=jp+j-1; + //int jdloc=jd+j-1; + if( jploc+ j_domain_begin_ ==0 ) { + // Jl^(d,p) + ( *JlAM )( idloc, jploc ) += Jl_ion*8. /dr * Sld[i]*Srp[j]; + ( *JrAM )( idloc, jploc ) += Jr_ion*8. /dr * Slp[i]*Srd[j]; + ( *JtAM )( idloc, jploc ) += Jt_ion*8. /dr * Slp[i]*Srp[j]; //A corriger dualite et repliement + } else { + ( *JlAM )( idloc, jploc ) += Jl_ion /( ( jploc+ j_domain_begin_ )*dr ) * Sld[i]*Srp[j]; + ( *JrAM )( idloc, jploc ) += Jr_ion /( ( jploc+ j_domain_begin_ )*dr ) * Slp[i]*Srd[j]; + ( *JtAM )( idloc, jploc ) += Jt_ion /( ( jploc+ j_domain_begin_ )*dr ) * Slp[i]*Srp[j]; + } + + } + }//i + + +} // END Project global current densities (ionize) + +// --------------------------------------------------------------------------------------------------------------------- +//! Project global current densities : ionization for tasks NOT DONE YET +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuyten::ionizationCurrentsForTasks( double */*b_Jx*/, double */*b_Jy*/, double */*b_Jz*/, Particles &/*particles*/, int /*ipart*/, LocalFields /*Jion*/, int /*bin_shift*/ ) +{ + // cField2D *JlAM = static_cast( Jl ); + // cField2D *JrAM = static_cast( Jr ); + // cField2D *JtAM = static_cast( Jt ); + // + // + // //Declaration of local variables + // int ip, id, jp, jd; + // double xpn, xpmxip, xpmxip2, xpmxid, xpmxid2; + // double ypn, ypmyjp, ypmyjp2, ypmyjd, ypmyjd2; + // double Slp[3], Sld[3], Srp[3], Srd[3]; + // + // // weighted currents + // double weight = inv_cell_volume * particles.weight( ipart ); + // double Jl_ion = Jion.x * weight; + // double Jr_ion = Jion.y * weight; + // double Jt_ion = Jion.z * weight; + // + // //Locate particle on the grid + // xpn = particles.position( 0, ipart ) * dl_inv_; // normalized distance to the first node + // ypn = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) )*dr_inv_ ; + // // x-primal index + // ip = round( xpn ); // x-index of the central node + // xpmxip = xpn - ( double )ip; // normalized distance to the nearest grid point + // xpmxip2 = xpmxip*xpmxip; // square of the normalized distance to the nearest grid point + // + // // x-dual index + // id = round( xpn+0.5 ); // x-index of the central node + // xpmxid = xpn - ( double )id + 0.5; // normalized distance to the nearest grid point + // xpmxid2 = xpmxid*xpmxid; // square of the normalized distance to the nearest grid point + // + // // y-primal index + // jp = round( ypn ); // y-index of the central node + // ypmyjp = ypn - ( double )jp; // normalized distance to the nearest grid point + // ypmyjp2 = ypmyjp*ypmyjp; // square of the normalized distance to the nearest grid point + // + // // y-dual index + // jd = round( ypn+0.5 ); // y-index of the central node + // ypmyjd = ypn - ( double )jd + 0.5; // normalized distance to the nearest grid point + // ypmyjd2 = ypmyjd*ypmyjd; // square of the normalized distance to the nearest grid point + // + // Slp[0] = 0.5 * ( xpmxip2-xpmxip+0.25 ); + // Slp[1] = ( 0.75-xpmxip2 ); + // Slp[2] = 0.5 * ( xpmxip2+xpmxip+0.25 ); + // + // Sld[0] = 0.5 * ( xpmxid2-xpmxid+0.25 ); + // Sld[1] = ( 0.75-xpmxid2 ); + // Sld[2] = 0.5 * ( xpmxid2+xpmxid+0.25 ); + // + // Srp[0] = 0.5 * ( ypmyjp2-ypmyjp+0.25 ); + // Srp[1] = ( 0.75-ypmyjp2 ); + // Srp[2] = 0.5 * ( ypmyjp2+ypmyjp+0.25 ); + // + // Srd[0] = 0.5 * ( ypmyjd2-ypmyjd+0.25 ); + // Srd[1] = ( 0.75-ypmyjd2 ); + // Srd[2] = 0.5 * ( ypmyjd2+ypmyjd+0.25 ); + // + // ip -= i_domain_begin_; + // id -= i_domain_begin_; + // jp -= j_domain_begin_; + // jd -= j_domain_begin_; + // + // for( unsigned int i=0 ; i<3 ; i++ ) { + // //int iploc=ip+i-1; + // int idloc=id+i-1; + // for( unsigned int j=0 ; j<3 ; j++ ) { + // int jploc=jp+j-1; + // //int jdloc=jd+j-1; + // if( jploc+ j_domain_begin ==0 ) { + // // Jl^(d,p) + // ( *JlAM )( idloc, jploc ) += Jl_ion*8. /dr * Sld[i]*Srp[j]; + // ( *JrAM )( idloc, jploc ) += Jr_ion*8. /dr * Slp[i]*Srd[j]; + // ( *JtAM )( idloc, jploc ) += Jt_ion*8. /dr * Slp[i]*Srp[j]; //A corriger dualite et repliement + // } else { + // ( *JlAM )( idloc, jploc ) += Jl_ion /( ( jploc+ j_domain_begin )*dr ) * Sld[i]*Srp[j]; + // ( *JrAM )( idloc, jploc ) += Jr_ion /( ( jploc+ j_domain_begin )*dr ) * Slp[i]*Srd[j]; + // ( *JtAM )( idloc, jploc ) += Jt_ion /( ( jploc+ j_domain_begin )*dr ) * Slp[i]*Srp[j]; + // } + // + // } + // }//i + // + +} // END Project global current densities (ionize) for tasks + +//------------------------------------// +//Wrapper for projection +void ProjectorAM1OrderRuyten::currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool /*is_spectral*/, int ispec, int /*icell*/, int /*ipart_ref*/ ) +{ + + std::vector *iold = &( smpi->dynamics_iold[ithread] ); + std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); + std::vector *invgf = &( smpi->dynamics_invgf[ithread] ); + std::vector> *array_eitheta_old = &( smpi->dynamics_eithetaold[ithread] ); + ElectroMagnAM *emAM = static_cast( EMfields ); + + for( int ipart=istart ; ipart *b_Jl, std::complex *b_Jr, std::complex *b_Jt, std::complex *b_rhoAM, int bin_shift, int bdim0, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, int /*ipart_ref*/ ) +{ + std::vector *iold = &( smpi->dynamics_iold[ithread] ); + std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); + std::vector *invgf = &( smpi->dynamics_invgf[ithread] ); + std::vector> *array_eitheta_old = &( smpi->dynamics_eithetaold[ithread] ); + ElectroMagnAM *emAM = static_cast( EMfields ); + + + + for( unsigned int ipart= (unsigned int) istart ; ipart< (unsigned int ) iend; ipart++ ) { + // cerr << ipart << endl; + // cerr << ( *iold )[ipart] << endl; + + // Jx is used for Jl + // Jy is used fot Jr + // Jt is used for Jt + currentsForTasks( emAM, b_Jl, b_Jr, b_Jt, b_rhoAM, particles, ipart, ( *invgf )[ipart], &( *iold )[ipart], &( *delta )[ipart], &( *array_eitheta_old )[ipart], bin_shift, bdim0, diag_flag ); + } + +} + +// --------------------------------------------------------------------------------------------------------------------- +//! Project local currents for all modes on bin sub-grid +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuyten::currentsForTasks( ElectroMagnAM */*emAM*/, std::complex *b_Jl, std::complex *b_Jr, std::complex *b_Jt, std::complex *b_rhoAM, Particles &particles, unsigned int ipart, double invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int bin_shift, int bdim0, bool diag_flag ) +{ + + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- int iloc, + int nparts= particles.size(); + int iloc, jloc, linindex; + // (x,y,z) components of the current density for the macro-particle + double charge_weight = inv_cell_volume * ( double )( particles.charge( ipart ) )*particles.weight( ipart ); + double crl_p = charge_weight*dl_ov_dt_; + double crr_p = charge_weight*one_ov_dt; + + // variable declaration + double xpn, ypn; + double delta, delta2; + // arrays used for the Esirkepov projection method + double Sl0[5], Sl1[5], Sr0[5], Sr1[5], DSl[5], DSr[5]; + complex Jl_p[5][5], Jr_p[5][5]; + complex e_delta, e_delta_m1, e_delta_inv, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; + + for( unsigned int i=0; i<5; i++ ) { + Sl1[i] = 0.; + Sr1[i] = 0.; + } + Sl0[0] = 0.; + Sl0[4] = 0.; + Sr0[0] = 0.; + Sr0[4] = 0.; + // -------------------------------------------------------- + // Locate particles & Calculate Esirkepov coef. S, DS and W + // -------------------------------------------------------- + + // locate the particle on the primal grid at former time-step & calculate coeff. S0 + delta = deltaold[0*nparts]; + delta2 = delta*delta; + Sl0[1] = 0.5 * ( delta2-delta+0.25 ); + Sl0[2] = 0.75-delta2; + Sl0[3] = 0.5 * ( delta2+delta+0.25 ); + + delta = deltaold[1*nparts]; + delta2 = delta*delta; + Sr0[1] = 0.5 * ( delta2-delta+0.25 ); + Sr0[2] = 0.75-delta2; + Sr0[3] = 0.5 * ( delta2+delta+0.25 ); + //calculate exponential coefficients + + // double yp = particles.position( 1, ipart ); + // double zp = particles.position( 2, ipart ); + double rp = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); + std::complex theta_old = array_eitheta_old[0]; + std::complex eitheta = ( particles.position( 1, ipart ) + Icpx * particles.position( 2, ipart ) ) / rp ; //exp(i theta) + e_bar = 1.; + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn ); + int ipo = iold[0*nparts]; + int ip_m_ipo = ip-ipo-i_domain_begin_; + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[ip_m_ipo+1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[ip_m_ipo+2] = 0.75-delta2; + Sl1[ip_m_ipo+3] = 0.5 * ( delta2+delta+0.25 ); + + ypn = rp *dr_inv_ ; + int jp = round( ypn ); + int jpo = iold[1*nparts]; + int jp_m_jpo = jp-jpo-j_domain_begin_; + delta = ypn - ( double )jp; + delta2 = delta*delta; + Sr1[jp_m_jpo+1] = 0.5 * ( delta2-delta+0.25 ); + Sr1[jp_m_jpo+2] = 0.75-delta2; + Sr1[jp_m_jpo+3] = 0.5 * ( delta2+delta+0.25 ); + + for( unsigned int i=0; i < 5; i++ ) { + DSl[i] = Sl1[i] - Sl0[i]; + DSr[i] = Sr1[i] - Sr0[i]; + } + + double r_bar = ((jpo + j_domain_begin_ + deltaold[1*nparts])*dr + rp) * 0.5; // r at t = t0 - dt/2 + e_delta_m1 = std::sqrt(eitheta * (2.*std::real(theta_old) - theta_old)); // std::sqrt keeps the root with positive real part which is what we need here. + e_bar_m1 = theta_old * e_delta_m1; + + ipo -= 2 + bin_shift; //This minus 2 come from the order 2 scheme, based on a 5 points stencil from -2 to +2. + // i/j/kpo stored with - i/j/k_domain_begin in Interpolator + jpo -= 2; + + double *invR_local = &(invR_[jpo]); + + //initial value of crt_p for imode = 0. + complex crt_p= charge_weight*( particles.momentum( 2, ipart )* real(e_bar_m1) - particles.momentum( 1, ipart )*imag(e_bar_m1) ) * invgf; + + // ------------------------------------------------ + // Local current created by the particle + // calculate using the charge conservation equation + // ------------------------------------------------ + for( unsigned int j=0 ; j<5 ; j++ ) { + Jl_p[0][j]= 0.; + } + for( unsigned int i=0 ; i<5 ; i++ ) { + Jr_p[i][4]= 0.; + } + + // --------------------------- + // Calculate the total current + // --------------------------- + + // Compute everything independent of theta + for( unsigned int j=0 ; j<5 ; j++ ) { + double tmp = crl_p * ( Sr0[j] + 0.5*DSr[j] )* invR_local[j]; + for( unsigned int i=1 ; i<5 ; i++ ) { + Jl_p[i][j]= Jl_p[i-1][j] - DSl[i-1] * tmp; + } + } + + for( int j=3 ; j>=0 ; j-- ) { + jloc = j+jpo+1; + double Vd = abs( jloc + j_domain_begin_ + 0.5 )* invRd_[jloc]*dr ; + double tmp = crr_p * DSr[j+1] * invRd_[jpo+j+1]*dr; + for( unsigned int i=0 ; i<5 ; i++ ) { + Jr_p[i][j] = Jr_p[i][j+1] * Vd + ( Sl0[i] + 0.5*DSl[i] ) * tmp; + } + } + + e_delta = 0.5; + e_delta_inv = 0.5; + + //Compute division by R in advance for Jt and rho evaluation. + for( unsigned int j=0 ; j<5 ; j++ ) { + Sr0[j] *= invR_local[j]; + Sr1[j] *= invR_local[j]; + } + + int factor_shift_JlJtRho = (bdim0*nprimr_ ); + int factor_shift_Jr = (bdim0*(nprimr_+1)); + + for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { + + if (imode > 0){ + e_delta *= e_delta_m1; + e_bar *= e_bar_m1; + C_m = 2. * e_bar ; //multiply modes > 0 by 2 and C_m = 1 otherwise. + e_delta_inv =1./e_delta - 1.; + crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.*r_bar; + } + + int mode_shift_JlJtRho = imode*factor_shift_JlJtRho; + int mode_shift_Jr = imode*factor_shift_Jr; + // Add contribution J_p to global array + if (diag_flag){ + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_; + for( unsigned int j=0 ; j<5 ; j++ ) { + jloc = j+jpo; + linindex = iloc+jloc; + b_rhoAM [mode_shift_JlJtRho+linindex] += C_m*charge_weight* Sl1[i]*Sr1[j]; + } + }//i + } + + // Jl^(d,p) + for( unsigned int i=1 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_+jpo; + for( unsigned int j=0 ; j<5 ; j++ ) { + linindex = iloc+j; + b_Jl [mode_shift_JlJtRho+linindex] += C_m * Jl_p[i][j] ; + } + }//i + + // Jr^(p,d) + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*( nprimr_+1 )+jpo+1; + for( unsigned int j=0 ; j<4 ; j++ ) { + linindex = iloc+j; + b_Jr [mode_shift_Jr+linindex] += C_m * Jr_p[i][j] ; + } + }//i + + // Jt^(p,p) + for( unsigned int i=0 ; i<5 ; i++ ) { + iloc = ( i+ipo )*nprimr_ + jpo; + for( unsigned int j=0 ; j<5 ; j++ ) { + linindex = iloc+j; + b_Jt [mode_shift_JlJtRho+linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. )); + } + } + + if (imode == 0) e_delta = 1. ; //Restore e_delta correct initial value. + }// end loop on modes + +} // END Project local current densities (Jl, Jr, Jt, sort) on bin subgrid + +// Projector for susceptibility used as source term in envelope equation +void ProjectorAM1OrderRuyten::susceptibility( ElectroMagn *EMfields, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int /*icell*/, int /*ipart_ref*/ ) + +{ + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + double *Chi_envelope = &( *EMfields->Env_Chi_ )( 0 ); + + std::vector *Epart = &( smpi->dynamics_Epart[ithread] ); + std::vector *Phipart = &( smpi->dynamics_PHIpart[ithread] ); + std::vector *GradPhipart = &( smpi->dynamics_GradPHIpart[ithread] ); + std::vector *inv_gamma_ponderomotive = &( smpi->dynamics_inv_gamma_ponderomotive[ithread] ); + + + double gamma_ponderomotive, gamma0, gamma0_sq; + double charge_over_mass_dts2, charge_sq_over_mass_sq_dts4, charge_sq_over_mass_sq; + double pxsm, pysm, pzsm; + double one_over_mass=1./species_mass; + double momentum[3]; + + int nparts = particles.size(); + double *Ex = &( ( *Epart )[0*nparts] ); + double *Ey = &( ( *Epart )[1*nparts] ); + double *Ez = &( ( *Epart )[2*nparts] ); + double *Phi = &( ( *Phipart )[0*nparts] ); + double *GradPhix = &( ( *GradPhipart )[0*nparts] ); + double *GradPhiy = &( ( *GradPhipart )[1*nparts] ); + double *GradPhiz = &( ( *GradPhipart )[2*nparts] ); + + for( int ipart=istart ; ipart e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; + double C_m = 1.; // only mode 0 + + + double xpn, ypn; + double delta, delta2; + double Sl1[3], Sr1[2]; + + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn ); + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[0] = 0.5 * ( delta2-delta+0.25 ); + Sl1[1] = 0.75-delta2; + Sl1[2] = 0.5 * ( delta2+delta+0.25 ); + + ypn = r * dr_inv_ ; + int jp = floor( ypn ); + delta = ypn - ( double )jp; + double x_n = jp ; + Sr1[0] = (1. - delta)*(5*x_n + 2 - ypn)/(4.*x_n + 2.); + Sr1[1] = 1. - Sr1[0]; //conservation de la charge + + // --------------------------- + // Calculate the total charge + // --------------------------- + ip -= i_domain_begin_ + 1; + jp -= j_domain_begin_ ; + + for( unsigned int i=0 ; i<3 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=0 ; j<2 ; j++ ) { + Chi_envelope [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invR_[j+jp]; + } + }//i + } +} + +// Projector for susceptibility used as source term in envelope equation +void ProjectorAM1OrderRuyten::susceptibilityOnBuffer( ElectroMagn */*EMfields*/, double *b_ChiAM, int bin_shift, int /*bdim0*/, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int /*icell*/, int /*ipart_ref*/ ) +{ + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + // double *Chi_envelope = &( *EMfields->Env_Chi_ )( 0 ); + + std::vector *Epart = &( smpi->dynamics_Epart[ithread] ); + std::vector *Phipart = &( smpi->dynamics_PHIpart[ithread] ); + std::vector *GradPhipart = &( smpi->dynamics_GradPHIpart[ithread] ); + std::vector *inv_gamma_ponderomotive = &( smpi->dynamics_inv_gamma_ponderomotive[ithread] ); + + + double gamma_ponderomotive, gamma0, gamma0_sq; + double charge_over_mass_dts2, charge_sq_over_mass_sq_dts4, charge_sq_over_mass_sq; + double pxsm, pysm, pzsm; + double one_over_mass=1./species_mass; + double momentum[3]; + + int nparts = particles.size(); + double *Ex = &( ( *Epart )[0*nparts] ); + double *Ey = &( ( *Epart )[1*nparts] ); + double *Ez = &( ( *Epart )[2*nparts] ); + double *Phi = &( ( *Phipart )[0*nparts] ); + double *GradPhix = &( ( *GradPhipart )[0*nparts] ); + double *GradPhiy = &( ( *GradPhipart )[1*nparts] ); + double *GradPhiz = &( ( *GradPhipart )[2*nparts] ); + + for( int ipart=istart ; ipart e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; + double C_m = 1.; // only mode 0 + + + double xpn, ypn; + double delta, delta2; + double Sl1[5], Sr1[5]; + + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn ); + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[2] = 0.75-delta2; + Sl1[3] = 0.5 * ( delta2+delta+0.25 ); + ypn = r * dr_inv_ ; + int jp = round( ypn ); + delta = ypn - ( double )jp; + delta2 = delta*delta; + Sr1[1] = 0.5 * ( delta2-delta+0.25 ); + Sr1[2] = 0.75-delta2; + Sr1[3] = 0.5 * ( delta2+delta+0.25 ); + + // --------------------------- + // Calculate the total charge + // --------------------------- + ip -= i_domain_begin_ + 2 + bin_shift; + jp -= j_domain_begin_ + 2; + + + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + b_ChiAM [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invR_[j+jp]; + } + }//i + + + + } + +} + + +void ProjectorAM1OrderRuyten::axisBCEnvChi( double *EnvChi ) +{ + +return; +} diff --git a/src/Projector/ProjectorAM1OrderRuyten.h b/src/Projector/ProjectorAM1OrderRuyten.h new file mode 100755 index 000000000..7eea6e257 --- /dev/null +++ b/src/Projector/ProjectorAM1OrderRuyten.h @@ -0,0 +1,58 @@ +#ifndef PROJECTORAM1ORDERRUYTEN_H +#define PROJECTORAM1ORDERRUYTEN_H + +#include + +#include "ProjectorAM.h" +#include "ElectroMagnAM.h" + + +class ProjectorAM1OrderRuyten : public ProjectorAM +{ +public: + ProjectorAM1OrderRuyten( Params &, Patch *patch ); + ~ProjectorAM1OrderRuyten(); + + inline void currents( ElectroMagnAM *emAM, Particles &particles, unsigned int ipart, double invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, bool diag_flag, int ispec); + + inline void currentsForTasks( ElectroMagnAM *emAM, std::complex *b_Jl, std::complex *b_Jr, std::complex *b_Jt, std::complex *b_rhoAM, Particles &particles, unsigned int ipart, double invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int bin_shift, int bdim0, bool diag_flag ); + + //! Project global current charge (EMfields->rho_), frozen & diagFields timestep + void basicForComplex( std::complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode ) override final; + + //! Project global current charge (EMfields->rho_), frozen & diagFields timestep + void basicForComplexOnBuffer( std::complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode, int bdim0, int bin_shift ) override final; + + //! Apply boundary conditions on Rho and J + void axisBC( ElectroMagnAM *emAM, bool diag_flag ) override final; + void apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ); + + //! Apply boundary conditions on Env_Chi + void axisBCEnvChi( double *EnvChi ) override final; + + //! Project global current densities if Ionization in Species::dynamics, + void ionizationCurrents( Field *Jl, Field *Jr, Field *Jt, Particles &particles, int ipart, LocalFields Jion ) override final; + + //! Project global current densities if Ionization in Species::dynamicsTasks + void ionizationCurrentsForTasks( double *b_Jx, double *b_Jy, double *b_Jz, Particles &particles, int ipart, LocalFields Jion, int bin_shift ) override final; + + //!Wrapper + void currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool is_spectral, int ispec, int icell = 0, int ipart_ref = 0 ) override final; + + //!Wrapper for projection on buffers + void currentsAndDensityWrapperOnBuffers( double *, double *, double *, double *, int, Particles &, SmileiMPI *, int, int, int, bool, bool, int, int = 0, int = 0 ) override final {}; + + //!Wrapper for projection on AM buffers + void currentsAndDensityWrapperOnAMBuffers( ElectroMagn *EMfields, std::complex *b_Jl, std::complex *b_Jr, std::complex *b_Jt, std::complex *b_rhoAM, int bin_shift, int bdim0, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, int ipart_ref = 0 ) override final; + + void susceptibility( ElectroMagn *EMfields, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int icell = 0, int ipart_ref = 0 ) override final; + + void susceptibilityOnBuffer( ElectroMagn *EMfields, double *b_ChiAM, int bin_shift, int bdim0, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int icell = 0, int ipart_ref = 0 ) override final; + + + +private: + double dt, dts2, dts4; +}; + +#endif diff --git a/src/Projector/ProjectorAM1OrderRuytenV.cpp b/src/Projector/ProjectorAM1OrderRuytenV.cpp new file mode 100755 index 000000000..631987ac8 --- /dev/null +++ b/src/Projector/ProjectorAM1OrderRuytenV.cpp @@ -0,0 +1,844 @@ +#include "ProjectorAM1OrderRuytenV.h" + +#include +#include +#include +#include "dcomplex.h" +#include "ElectroMagnAM.h" +#include "cField2D.h" +#include "Particles.h" +#include "Tools.h" +#include "Patch.h" +#include "PatchAM.h" + +using namespace std; + + +// --------------------------------------------------------------------------------------------------------------------- +// Constructor for ProjectorAM1OrderRuytenV +// --------------------------------------------------------------------------------------------------------------------- +ProjectorAM1OrderRuytenV::ProjectorAM1OrderRuytenV( Params ¶ms, Patch *patch ) : ProjectorAM( params, patch ) +{ + dt = params.timestep; + dr = params.cell_length[1]; + dl_inv_ = 1.0/params.cell_length[0]; + dl_ov_dt_ = params.cell_length[0] / params.timestep; + one_ov_dt = 1.0 / params.timestep; + dr_inv_ = 1.0/dr; + dr_ov_dt_ = dr / dt; + + i_domain_begin_ = patch->getCellStartingGlobalIndex( 0 ); + j_domain_begin_ = patch->getCellStartingGlobalIndex( 1 ); + + nscellr_ = params.patch_size_[1] + 1; + oversize_[0] = params.oversize[0]; + oversize_[1] = params.oversize[1]; + nprimr_ = nscellr_ + 2*oversize_[1]; + npriml_ = params.patch_size_[0] + 1 + 2*oversize_[0]; + + Nmode_=params.nmodes; + dq_inv_[0] = dl_inv_; + dq_inv_[1] = dr_inv_; + + invR_ = &((static_cast( patch )->invR)[0]); + invRd_ = &((static_cast( patch )->invRd)[0]); + + DEBUG( "cell_length "<< params.cell_length[0] ); + +} + + +// --------------------------------------------------------------------------------------------------------------------- +// Destructor for ProjectorAM1OrderRuytenV +// --------------------------------------------------------------------------------------------------------------------- +ProjectorAM1OrderRuytenV::~ProjectorAM1OrderRuytenV() +{ +} + + + +// --------------------------------------------------------------------------------------------------------------------- +//! Project current densities & charge : diagFields timstep (not vectorized) +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuytenV::currentsAndDensity( ElectroMagnAM *emAM, + Particles &particles, + unsigned int istart, + unsigned int iend, + double * __restrict__ invgf, + int * __restrict__ iold, + double * __restrict__ deltaold, + std::complex * __restrict__ array_eitheta_old, + int npart_total, + int ipart_ref, + int ispec ) +{ + + currents( emAM, particles, istart, iend, invgf, iold, deltaold, array_eitheta_old, npart_total, ipart_ref, ispec+1 ); //ispec+1 is passed as a marker of diag + + int ipo = iold[0]; + int jpo = iold[1]; + int ipom2 = ipo-2; + int jpom2 = jpo-2; + + int vecSize = 8; + int bsize = 5*5*vecSize*Nmode_; + + std::complex brho[bsize] __attribute__( ( aligned( 64 ) ) ); + + double Sl0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + double Sr0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + double DSl[40] __attribute__( ( aligned( 64 ) ) ); + double DSr[40] __attribute__( ( aligned( 64 ) ) ); + double charge_weight[8] __attribute__( ( aligned( 64 ) ) ); + double r_bar[8] __attribute__( ( aligned( 64 ) ) ); + complex * __restrict__ rho; + + double *invR_local = &(invR_[jpom2]); + + // Pointer for GPU and vectorization on ARM processors + double * __restrict__ position_x = particles.getPtrPosition(0); + double * __restrict__ position_y = particles.getPtrPosition(1); + double * __restrict__ position_z = particles.getPtrPosition(2); + double * __restrict__ weight = particles.getPtrWeight(); + short * __restrict__ charge = particles.getPtrCharge(); + int * __restrict__ cell_keys = particles.getPtrCellKeys(); + + #pragma omp simd + for( unsigned int j=0; j<200*Nmode_; j++ ) { + brho[j] = 0.; + } + + // Closest multiple of 8 higher or equal than npart = iend-istart. + int cell_nparts( ( int )iend-( int )istart ); + + for( int ivect=0 ; ivect < cell_nparts; ivect += vecSize ) { + + int np_computed = min( cell_nparts-ivect, vecSize ); + int istart0 = ( int )istart + ivect; + complex e_bar[8], e_delta_m1[8]; + + #pragma omp simd + for( int ipart=0 ; ipartrho_AM_[imode] )( 0 ); + unsigned int n_species = emAM->rho_AM_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec; + rho = emAM->rho_AM_s [ifield] ? &( * ( emAM->rho_AM_s [ifield] ) )( 0 ) : &( *emAM->rho_AM_ [imode] )( 0 ) ; + int iloc = iloc0; + for( unsigned int i=0 ; i<5 ; i++ ) { + #pragma omp simd + for( unsigned int j=0 ; j<5 ; j++ ) { + complex tmprho( 0. ); + int ilocal = ( i*5+j )*vecSize; + UNROLL(8) + for( int ipart=0 ; ipart<8; ipart++ ) { + tmprho += brho [200*imode + ilocal+ipart]; + } + rho[iloc+j] += tmprho; + } + iloc += nprimr_; + } + } + +} // END Project local current densities at dag timestep. + +// --------------------------------------------------------------------------------------------------------------------- +//! Project for diags and frozen species - +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuytenV::basicForComplex( complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode ) +{ + //Warning : this function is not charge conserving. + // This function also assumes that particles position is evaluated at the same time as currents which is usually not true (half time-step difference). + // It will therefore fail to evaluate the current accurately at t=0 if a plasma is already in the box. + + + + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + int iloc, nr( nprimr_ ); + double charge_weight = inv_cell_volume * ( double )( particles.charge( ipart ) )*particles.weight( ipart ); + double r = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); + + if( type > 0 ) { //if current density + charge_weight *= 1./sqrt( 1.0 + particles.momentum( 0, ipart )*particles.momentum( 0, ipart ) + + particles.momentum( 1, ipart )*particles.momentum( 1, ipart ) + + particles.momentum( 2, ipart )*particles.momentum( 2, ipart ) ); + if( type == 1 ) { //if Jl + charge_weight *= particles.momentum( 0, ipart ); + } else if( type == 2 ) { //if Jr + charge_weight *= ( particles.momentum( 1, ipart )*particles.position( 1, ipart ) + particles.momentum( 2, ipart )*particles.position( 2, ipart ) )/ r ; + nr++; + } else { //if Jt + charge_weight *= ( -particles.momentum( 1, ipart )*particles.position( 2, ipart ) + particles.momentum( 2, ipart )*particles.position( 1, ipart ) ) / r ; + } + } + + complex e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; + complex C_m = 1.; + if( imode > 0 ) { + C_m = 2.; + } + for( unsigned int i=0; i<( unsigned int )imode; i++ ) { + C_m *= e_theta; + } + + double xpn, ypn; + double delta, delta2; + double Sl1[5], Sr1[5]; + + // -------------------------------------------------------- + // Locate particles & Calculate Esirkepov coef. S, DS and W + // -------------------------------------------------------- + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + xpn = particles.position( 0, ipart ) * dl_inv_; + int ip = round( xpn + 0.5 * ( type==1 ) ); + delta = xpn - ( double )ip; + delta2 = delta*delta; + Sl1[1] = 0.5 * ( delta2-delta+0.25 ); + Sl1[2] = 0.75-delta2; + Sl1[3] = 0.5 * ( delta2+delta+0.25 ); + ypn = r * dr_inv_ ; + int jp = round( ypn + 0.5*( type==2 ) ); + delta = ypn - ( double )jp; + delta2 = delta*delta; + Sr1[1] = 0.5 * ( delta2-delta+0.25 ); + Sr1[2] = 0.75-delta2; + Sr1[3] = 0.5 * ( delta2+delta+0.25 ); + + // --------------------------- + // Calculate the total charge + // --------------------------- + ip -= i_domain_begin_ + 2; + jp -= j_domain_begin_ + 2; + + if( type != 2 ) { + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + rhoj [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invR_[j+jp]; + } + }//i + } else { + for( unsigned int i=1 ; i<4 ; i++ ) { + iloc = ( i+ip )*nr+jp; + for( unsigned int j=1 ; j<4 ; j++ ) { + rhoj [iloc+j] += C_m*charge_weight* Sl1[i]*Sr1[j] * invRd_[j+jp]; + } + }//i + } +} // END Project for diags local current densities + +// Apply boundary conditions on axis for currents and densities +void ProjectorAM1OrderRuytenV::axisBC(ElectroMagnAM *emAM, bool diag_flag ) +{ + + for (unsigned int imode=0; imode < Nmode_; imode++){ + + std::complex *rhoj = &( *emAM->rho_AM_[imode] )( 0 ); + std::complex *Jl = &( *emAM->Jl_[imode] )( 0 ); + std::complex *Jr = &( *emAM->Jr_[imode] )( 0 ); + std::complex *Jt = &( *emAM->Jt_[imode] )( 0 ); + + apply_axisBC(rhoj, Jl, Jr, Jt, imode, diag_flag); + } + + if (diag_flag){ + unsigned int n_species = emAM->Jl_s.size() / Nmode_; + for( unsigned int imode = 0 ; imode < emAM->Jl_.size() ; imode++ ) { + for( unsigned int ispec = 0 ; ispec < n_species ; ispec++ ) { + unsigned int ifield = imode*n_species+ispec; + complex *Jl = emAM->Jl_s [ifield] ? &( * ( emAM->Jl_s [ifield] ) )( 0 ) : NULL ; + complex *Jr = emAM->Jr_s [ifield] ? &( * ( emAM->Jr_s [ifield] ) )( 0 ) : NULL ; + complex *Jt = emAM->Jt_s [ifield] ? &( * ( emAM->Jt_s [ifield] ) )( 0 ) : NULL ; + complex *rho = emAM->rho_AM_s[ifield] ? &( * ( emAM->rho_AM_s[ifield] ) )( 0 ) : NULL ; + apply_axisBC( rho , Jl, Jr, Jt, imode, diag_flag ); + } + } + } +} + +void ProjectorAM1OrderRuytenV::apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ) +{ + // Mode 0 contribution "below axis" is added. + // Non zero modes are substracted because a particle sitting exactly on axis has a non defined theta and can not contribute to a theta dependent mode. + double sign = (imode == 0) ? 1 : -1 ; + + if (diag_flag && rhoj) { + for( unsigned int i=2 ; i 0){ + rhoj[i] = 0.; + rhoj[i-1] = - rhoj[i+1]; // Zero Jl mode > 0 on axis. + } else { + rhoj[i] = rhoj[i+1]; //This smoothing is just for cosmetics on the picture, rho has no influence on the results. + rhoj[i-1] = rhoj[i+1]; // Non zero Jl mode > 0 on axis. + } + } + } + + if (Jl) { + for( unsigned int i=2 ; i<(npriml_+1)*nprimr_+2; i+=nprimr_ ) { + //Fold Jl + for( unsigned int j=1 ; j<3; j++ ) { + Jl [i+j] += sign * Jl[i-j]; //Add even modes, substract odd modes since el(theta=0 = el(theta=pi) at all r. + } + if (imode > 0){ + Jl [i] = 0. ; + Jl[i-1] = -Jl[i+1]; // Zero Jl mode > 0 on axis. + } else { + //Jl mode 0 on axis should be left as is. It looks over estimated but it might be necessary to conserve a correct divergence and a proper evaluation on the field on axis. + Jl [i-1] = Jl [i+1] ; // Non zero Jl mode 0 on axis. + } + } + } + + if (Jt && Jr) { + for( unsigned int i=0 ; i( Jx ); + Field2D *Jy2D = static_cast( Jy ); + Field2D *Jz2D = static_cast( Jz ); + + + //Declaration of local variables + int ip, id, jp, jd; + double xpn, xpmxip, xpmxip2, xpmxid, xpmxid2; + double ypn, ypmyjp, ypmyjp2, ypmyjd, ypmyjd2; + double Sxp[3], Sxd[3], Syp[3], Syd[3]; + + // weighted currents + double weight = inv_cell_volume * particles.weight( ipart ); + double Jx_ion = Jion.x * weight; + double Jy_ion = Jion.y * weight; + double Jz_ion = Jion.z * weight; + + //Locate particle on the grid + xpn = particles.position( 0, ipart ) * dx_inv_; // normalized distance to the first node + ypn = particles.position( 1, ipart ) * dy_inv_; // normalized distance to the first node + + // x-primal index + ip = round( xpn ); // x-index of the central node + xpmxip = xpn - ( double )ip; // normalized distance to the nearest grid point + xpmxip2 = xpmxip*xpmxip; // square of the normalized distance to the nearest grid point + + // x-dual index + id = round( xpn+0.5 ); // x-index of the central node + xpmxid = xpn - ( double )id + 0.5; // normalized distance to the nearest grid point + xpmxid2 = xpmxid*xpmxid; // square of the normalized distance to the nearest grid point + + // y-primal index + jp = round( ypn ); // y-index of the central node + ypmyjp = ypn - ( double )jp; // normalized distance to the nearest grid point + ypmyjp2 = ypmyjp*ypmyjp; // square of the normalized distance to the nearest grid point + + // y-dual index + jd = round( ypn+0.5 ); // y-index of the central node + ypmyjd = ypn - ( double )jd + 0.5; // normalized distance to the nearest grid point + ypmyjd2 = ypmyjd*ypmyjd; // square of the normalized distance to the nearest grid point + + Sxp[0] = 0.5 * ( xpmxip2-xpmxip+0.25 ); + Sxp[1] = ( 0.75-xpmxip2 ); + Sxp[2] = 0.5 * ( xpmxip2+xpmxip+0.25 ); + + Sxd[0] = 0.5 * ( xpmxid2-xpmxid+0.25 ); + Sxd[1] = ( 0.75-xpmxid2 ); + Sxd[2] = 0.5 * ( xpmxid2+xpmxid+0.25 ); + + Syp[0] = 0.5 * ( ypmyjp2-ypmyjp+0.25 ); + Syp[1] = ( 0.75-ypmyjp2 ); + Syp[2] = 0.5 * ( ypmyjp2+ypmyjp+0.25 ); + + Syd[0] = 0.5 * ( ypmyjd2-ypmyjd+0.25 ); + Syd[1] = ( 0.75-ypmyjd2 ); + Syd[2] = 0.5 * ( ypmyjd2+ypmyjd+0.25 ); + + ip -= i_domain_begin_; + id -= i_domain_begin_; + jp -= j_domain_begin_; + jd -= j_domain_begin_; + + for( unsigned int i=0 ; i<3 ; i++ ) { + int iploc=ip+i-1; + int idloc=id+i-1; + for( unsigned int j=0 ; j<3 ; j++ ) { + int jploc=jp+j-1; + int jdloc=jd+j-1; + // Jx^(d,p) + ( *Jx2D )( idloc, jploc ) += Jx_ion * Sxd[i]*Syp[j]; + // Jy^(p,d) + ( *Jy2D )( iploc, jdloc ) += Jy_ion * Sxp[i]*Syd[j]; + // Jz^(p,p) + ( *Jz2D )( iploc, jploc ) += Jz_ion * Sxp[i]*Syp[j]; + } + }//i*/ + + +} // END Project global current densities (ionize) + + +// --------------------------------------------------------------------------------------------------------------------- +//! Project current densities : main projector vectorized +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuytenV::currents( ElectroMagnAM *emAM, + Particles &particles, + unsigned int istart, + unsigned int iend, + double * __restrict__ invgf, + int * __restrict__ iold, + double * __restrict__ deltaold, + std::complex * __restrict__ array_eitheta_old, + int npart_total, + int ipart_ref, + int ispec ) +{ + // ------------------------------------- + // Variable declaration & initialization + // ------------------------------------- + + int ipo = iold[0]; + int jpo = iold[1]; + int ipom2 = ipo-2; + int jpom2 = jpo-2; + + int vecSize = 8; + int bsize = 5*5*vecSize*Nmode_; + + std::complex bJl[bsize] __attribute__( ( aligned( 64 ) ) ); + std::complex bJr[bsize] __attribute__( ( aligned( 64 ) ) ); + std::complex bJt[bsize] __attribute__( ( aligned( 64 ) ) ); + + double Sl0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + double Sr0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + double DSl[40] __attribute__( ( aligned( 64 ) ) ); + double DSr[40] __attribute__( ( aligned( 64 ) ) ); + double charge_weight[8] __attribute__( ( aligned( 64 ) ) ); + double r_bar[8] __attribute__( ( aligned( 64 ) ) ); + complex * __restrict__ Jl; + complex * __restrict__ Jr; + complex * __restrict__ Jt; + + double *invR_local = &(invR_[jpom2]); + double *invRd_local = &(invRd_[jpom2]); + + // Pointer for GPU and vectorization on ARM processors + double * __restrict__ position_x = particles.getPtrPosition(0); + double * __restrict__ position_y = particles.getPtrPosition(1); + double * __restrict__ position_z = particles.getPtrPosition(2); + double * __restrict__ momentum_y = particles.getPtrMomentum(1); + double * __restrict__ momentum_z = particles.getPtrMomentum(2); + double * __restrict__ weight = particles.getPtrWeight(); + short * __restrict__ charge = particles.getPtrCharge(); + int * __restrict__ cell_keys = particles.getPtrCellKeys(); + + #pragma omp simd + for( unsigned int j=0; j<200*Nmode_; j++ ) { + bJl[j] = 0.; + bJr[j] = 0.; + bJt[j] = 0.; + } + + // Closest multiple of 8 higher or equal than npart = iend-istart. + int cell_nparts( ( int )iend-( int )istart ); + + for( int ivect=0 ; ivect < cell_nparts; ivect += vecSize ) { + + int np_computed = min( cell_nparts-ivect, vecSize ); + int istart0 = ( int )istart + ivect; + complex e_bar[8], e_delta_m1[8]; + + #pragma omp simd + for( int ipart=0 ; ipartJl_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jl_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jl = emAM->Jl_s [ifield] ? &( * ( emAM->Jl_s [ifield] ) )( 0 ) : &( *emAM->Jl_ [imode] )( 0 ) ; + } + int iloc = iloc0; + for( unsigned int i=1 ; i<5 ; i++ ) { + iloc += nprimr_; + #pragma omp simd + for( unsigned int j=0 ; j<5 ; j++ ) { + complex tmpJl( 0. ); + int ilocal = ( i*5+j )*vecSize; + UNROLL(8) + for( int ipart=0 ; ipart<8; ipart++ ) { + tmpJl += bJl [200*imode + ilocal+ipart]; + } + Jl[iloc+j] += tmpJl; + } + } + } + + + for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { + if (ispec == 0) { // When diags are not needed ispec is always 0. + Jr = &( *emAM->Jr_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jr_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jr = emAM->Jr_s [ifield] ? &( * ( emAM->Jr_s [ifield] ) )( 0 ) : &( *emAM->Jr_ [imode] )( 0 ) ; + } + + int iloc = iloc0 + ipom2 + 1; + for( unsigned int i=0 ; i<5 ; i++ ) { + #pragma omp simd + for( unsigned int j=0 ; j<4 ; j++ ) { + complex tmpJr( 0. ); + int ilocal = ( i*5+j+1 )*vecSize; + UNROLL(8) + for( int ipart=0 ; ipart<8; ipart++ ) { + tmpJr += bJr [200*imode + ilocal+ipart]; + } + Jr[iloc+j] += tmpJr; + } + iloc += nprimr_+1; + } + } + + for( unsigned int imode=0; imode<( unsigned int )Nmode_; imode++ ) { + if (ispec == 0) { // When diags are not needed ispec is always 0. + Jt = &( *emAM->Jt_[imode] )( 0 ); + } else { // When diags are needed ispec+1 is passed. + unsigned int n_species = emAM->Jt_s.size() / Nmode_; + unsigned int ifield = imode*n_species+ispec-1; + Jt = emAM->Jt_s [ifield] ? &( * ( emAM->Jt_s [ifield] ) )( 0 ) : &( *emAM->Jt_ [imode] )( 0 ) ; + } + int iloc = iloc0; + for( unsigned int i=0 ; i<5 ; i++ ) { + #pragma omp simd + for( unsigned int j=0 ; j<5 ; j++ ) { + complex tmpJt( 0. ); + int ilocal = ( i*5+j )*vecSize; + UNROLL(8) + for( int ipart=0 ; ipart<8; ipart++ ) { + tmpJt += bJt [200*imode + ilocal+ipart]; + } + Jt[iloc+j] += tmpJt; + } + iloc += nprimr_; + } + } +} // END Projection currents vectorized + + +// --------------------------------------------------------------------------------------------------------------------- +//! Wrapper for projection +// --------------------------------------------------------------------------------------------------------------------- +void ProjectorAM1OrderRuytenV::currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool is_spectral, int ispec, int scell, int ipart_ref ) +{ + if( istart == iend ) { + return; //Don't treat empty cells. + } + + //Independent of cell. Should not be here + //{ + std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); + std::vector *invgf = &( smpi->dynamics_invgf[ithread] ); + std::vector> *array_eitheta_old = &( smpi->dynamics_eithetaold[ithread] ); + ElectroMagnAM *emAM = static_cast( EMfields ); + + + //} + int iold[2]; + iold[0] = scell/nscellr_+oversize_[0]; + iold[1] = ( scell%nscellr_ )+oversize_[1]; + + + // If no field diagnostics this timestep, then the projection is done directly on the total arrays + if( !diag_flag ) { + if( !is_spectral ) { + currents( emAM, particles, istart, iend, invgf->data(), iold, delta->data(), array_eitheta_old->data(), invgf->size(), ipart_ref ); + } else { + ERROR( "Vectorized projection is not supported in spectral AM" ); + } + + // Otherwise, the projection may apply to the species-specific arrays + } else { + //double *b_Jx = EMfields->Jx_s [ispec] ? &( *EMfields->Jx_s [ispec] )( 0 ) : &( *EMfields->Jx_ )( 0 ) ; + //double *b_Jy = EMfields->Jy_s [ispec] ? &( *EMfields->Jy_s [ispec] )( 0 ) : &( *EMfields->Jy_ )( 0 ) ; + //double *b_Jz = EMfields->Jz_s [ispec] ? &( *EMfields->Jz_s [ispec] )( 0 ) : &( *EMfields->Jz_ )( 0 ) ; + //double *b_rho = EMfields->rho_s[ispec] ? &( *EMfields->rho_s[ispec] )( 0 ) : &( *EMfields->rho_ )( 0 ) ; + currentsAndDensity( emAM, particles, istart, iend, invgf->data(), iold, delta->data(), array_eitheta_old->data(), invgf->size(), ipart_ref, ispec ); + } +} + +// Project susceptibility +void ProjectorAM1OrderRuytenV::susceptibility( ElectroMagn *EMfields, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int icell, int ipart_ref ) +{ + double dts2 = dt/2.; + double dts4 = dt/4.; + double * __restrict__ Chi_envelope = &( *EMfields->Env_Chi_ )( 0 ) ; + + int iold[2]; + iold[0] = icell/nscellr_+oversize_[0]; + iold[1] = ( icell%nscellr_ )+oversize_[1]; + + int ipom2 = iold[0]-2; + int jpom2 = iold[1]-2; + + int vecSize = 8; + int bsize = 5*5*vecSize; // Chi has only one mode //*Nmode_; + + std::vector *Epart = &( smpi->dynamics_Epart[ithread] ); + std::vector *Phipart = &( smpi->dynamics_PHIpart[ithread] ); + std::vector *GradPhipart = &( smpi->dynamics_GradPHIpart[ithread] ); + double * __restrict__ inv_gamma_ponderomotive = &( smpi->dynamics_inv_gamma_ponderomotive[ithread][0] ); + + + int nparts = smpi->dynamics_invgf[ithread].size(); + double * __restrict__ Ex = &( ( *Epart )[0*nparts] ); + double * __restrict__ Ey = &( ( *Epart )[1*nparts] ); + double * __restrict__ Ez = &( ( *Epart )[2*nparts] ); + double * __restrict__ Phi = &( ( *Phipart )[0*nparts] ); + double * __restrict__ GradPhix = &( ( *GradPhipart )[0*nparts] ); + double * __restrict__ GradPhiy = &( ( *GradPhipart )[1*nparts] ); + double * __restrict__ GradPhiz = &( ( *GradPhipart )[2*nparts] ); + + double bChi[bsize] __attribute__( ( aligned( 64 ) ) ); + + double Sl1[32] __attribute__( ( aligned( 64 ) ) ); + double Sr1[32] __attribute__( ( aligned( 64 ) ) ); + // double Sl0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + // double Sr0_buff_vect[32] __attribute__( ( aligned( 64 ) ) ); + // double DSl[40] __attribute__( ( aligned( 64 ) ) ); + // double DSr[40] __attribute__( ( aligned( 64 ) ) ); + double charge_weight[8] __attribute__( ( aligned( 64 ) ) ); + // double r_bar[8] __attribute__( ( aligned( 64 ) ) ); + + // Pointer for GPU and vectorization on ARM processors + double * __restrict__ position_x = particles.getPtrPosition(0); + double * __restrict__ position_y = particles.getPtrPosition(1); + double * __restrict__ position_z = particles.getPtrPosition(2); + double * __restrict__ momentum_x = particles.getPtrMomentum(0); + double * __restrict__ momentum_y = particles.getPtrMomentum(1); + double * __restrict__ momentum_z = particles.getPtrMomentum(2); + double * __restrict__ weight = particles.getPtrWeight(); + short * __restrict__ charge = particles.getPtrCharge(); + //int * __restrict__ cell_keys = particles.getPtrCellKeys(); + + // Closest multiple of 8 higher or equal than npart = iend-istart. + int cell_nparts( ( int )iend-( int )istart ); + + double one_over_mass=1./species_mass; + + for( int ivect=0 ; ivect < cell_nparts; ivect += vecSize ) { + + int np_computed = min( cell_nparts-ivect, vecSize ); + int istart0 = ( int )istart + ivect; + + #pragma omp simd + for( unsigned int j=0; j<200*Nmode_; j++ ) { + bChi[j] = 0.; + } + + #pragma omp simd + for( int ipart=0 ; ipart= 0.); + Sr1[1*vecSize+ipart] = 1. - Sr1[0*vecSize+ipart] - Sr1[2*vecSize+ipart]; + + Sr1[0*vecSize+ipart] *= invR_[1+jp]; + Sr1[1*vecSize+ipart] *= invR_[2+jp]; + Sr1[2*vecSize+ipart] *= invR_[3+jp]; + + } // end ipart loop + + #pragma omp simd + for( int ipart=0 ; ipart +#include "ProjectorAM.h" +#include +#include "dcomplex.h" +#include "Pragma.h" +#include + +using namespace std; + +class ProjectorAM1OrderRuytenV : public ProjectorAM +{ +public: + ProjectorAM1OrderRuytenV( Params &, Patch *patch ); + ~ProjectorAM1OrderRuytenV(); + + //! Project global current densities (EMfields->Jl_/Jr_/Jt_) + void currents(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0, int ispec =0 ); + + //! Project global current densities (EMfields->Jl_/Jr_/Jt_/rho), diagFields timestep + void currentsAndDensity(ElectroMagnAM *emAM, Particles &particles, unsigned int istart, unsigned int iend, double *invgf, int *iold, double *deltaold, std::complex *array_eitheta_old, int npart_total, int ipart_ref = 0, int ispec = 0 ); + + //! Project global current charge (EMfields->rho_), frozen & diagFields timestep + void basicForComplex( std::complex *rhoj, Particles &particles, unsigned int ipart, unsigned int type, int imode ) override final; + + //! Apply boundary conditions on Rho and J + void axisBC( ElectroMagnAM *emAM, bool diag_flag ) override final; + void apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ); + + //! Apply boundary conditions on Env_Chi + void axisBCEnvChi( double *EnvChi ) override final; + + //! Project global current densities if Ionization in Species::dynamics, + void ionizationCurrents( Field *Jl, Field *Jr, Field *Jt, Particles &particles, int ipart, LocalFields Jion ) override final; + // + //!Wrapper + void currentsAndDensityWrapper( ElectroMagn *EMfields, Particles &particles, SmileiMPI *smpi, int istart, int iend, int ithread, bool diag_flag, bool is_spectral, int ispec, int icell, int ipart_ref ) override final; + + //!Wrapper for projection on buffers + void currentsAndDensityWrapperOnBuffers( double *, double *, double *, double *, int, Particles &, SmileiMPI *, int, int, int, bool, bool, int, int = 0, int = 0 ) override final {}; + + //!Wrapper for projection on AM buffers + void currentsAndDensityWrapperOnAMBuffers( ElectroMagn *, std::complex *, std::complex *, std::complex *, std::complex *, int, int, Particles &, SmileiMPI *, int, int, int, bool, int = 0 ) override final {}; + + // Project susceptibility + void susceptibility( ElectroMagn *EMfields, Particles &particles, double species_mass, SmileiMPI *smpi, int istart, int iend, int ithread, int icell, int ipart_ref ) override final; + +private: + + inline void __attribute__((always_inline)) compute_distances( double * __restrict__ position_x, + double * __restrict__ position_y, + double * __restrict__ position_z, + int * __restrict__, + int npart_total, int ipart, int istart, int ipart_ref, + double *deltaold, std::complex *array_eitheta_old, int *iold, + double *Sl0, double *Sr0, double *DSl, double *DSr, + double *r_bar, std::complex *e_bar, std::complex *e_delta_m1) + { + + int ipo = iold[0]; + int jpo = iold[1]; + int vecSize = 8; + + // locate the particle on the primal grid at former time-step & calculate coeff. S0 + // L // + double delta = deltaold[istart+ipart-ipart_ref]; + double delta2 = delta*delta; + Sl0[ ipart] = 0.5 * ( delta2-delta+0.25 ); + Sl0[ vecSize+ipart] = 0.75-delta2; + Sl0[2*vecSize+ipart] = 0.5 * ( delta2+delta+0.25 ); + Sl0[3*vecSize+ipart] = 0.; + // R // + delta = deltaold[istart+ipart-ipart_ref+npart_total]; // delta = x - x_n si delta > 0 et delta = x - x_n -1 si delta < 0. + double pos = (double)jpo + j_domain_begin_ + delta; + double x_n = floor(pos); + double coeff0 = (x_n+1-pos)*(5*x_n + 2 - pos)/(4.*x_n + 2.); + + Sr0[ ipart] = coeff0 * (delta < 0. ); + Sr0[2*vecSize+ipart] = (1.-coeff0) * (delta >= 0. ); + Sr0[ vecSize+ipart] = 1. - Sr0[ipart] - Sr0[2*vecSize+ipart]; + Sr0[3*vecSize+ipart] = 0.; + + + // locate the particle on the primal grid at current time-step & calculate coeff. S1 + // L // + pos = position_x[istart + ipart] * dl_inv_; + int cell = round( pos ); + int cell_shift = cell-ipo-i_domain_begin_; + delta = pos - ( double )cell; + delta2 = delta*delta; + double deltam = 0.5 * ( delta2-delta+0.25 ); + double deltap = 0.5 * ( delta2+delta+0.25 ); + delta2 = 0.75 - delta2; + double m1 = ( cell_shift == -1 ); + double c0 = ( cell_shift == 0 ); + double p1 = ( cell_shift == 1 ); + DSl [ ipart] = m1 * deltam ; + DSl [ vecSize+ipart] = c0 * deltam + m1 * delta2 - Sl0[ ipart]; + DSl [2*vecSize+ipart] = p1 * deltam + c0 * delta2 + m1* deltap - Sl0[ vecSize+ipart]; + DSl [3*vecSize+ipart] = p1 * delta2 + c0* deltap - Sl0[2*vecSize+ipart]; + DSl [4*vecSize+ipart] = p1* deltap ; + + double rp = sqrt( position_y[istart+ipart]*position_y[istart+ipart] + position_z[istart+ipart]*position_z[istart+ipart] ); + pos = rp * dr_inv_; + + x_n = floor(pos); + coeff0 = (x_n+1-pos)*(5*x_n + 2 - pos)/(4.*x_n + 2.); + + cell = round( pos ); + cell_shift = cell-jpo-j_domain_begin_; + delta = pos - ( double )cell; + deltam = coeff0 * (double)(delta < 0 ); + deltap = (1.-coeff0) * (double)(delta >= 0 ); + delta2 = 1. - deltam - deltap; + m1 = ( cell_shift == -1 ); + c0 = ( cell_shift == 0 ); + p1 = ( cell_shift == 1 ); + DSr [ ipart] = m1 * deltam ; + DSr [ vecSize+ipart] = c0 * deltam + m1 * delta2 - Sr0[ ipart] ; + DSr [2*vecSize+ipart] = p1 * deltam + c0 * delta2 + m1* deltap - Sr0[ vecSize+ipart] ; + DSr [3*vecSize+ipart] = p1 * delta2 + c0* deltap - Sr0[2*vecSize+ipart] ; + DSr [4*vecSize+ipart] = p1* deltap ; + + r_bar[ipart] = ((jpo + j_domain_begin_ + deltaold[istart+ipart-ipart_ref+npart_total])*dr + rp) * 0.5; // r at t = t0 - dt/2 + std::complex eitheta = ( position_y[istart+ipart] + Icpx * position_z[istart+ipart] ) / rp ; //exp(i theta) + e_delta_m1[ipart] = std::sqrt(eitheta * (2.*std::real(array_eitheta_old[istart+ipart-ipart_ref]) - array_eitheta_old[istart+ipart-ipart_ref])); + e_bar[ipart] = array_eitheta_old[istart+ipart-ipart_ref] * e_delta_m1[ipart]; + + } + + inline void __attribute__((always_inline)) computeJl( int ipart, double *charge_weight, double *DSl, double *DSr, double *Sr0, std::complex *bJ, double, double *invR_local, std::complex *e_bar ) + { + + int vecSize = 8; + double sum[5]; + std::complex C_m =2.; + double crl_p = charge_weight[ipart]*dl_ov_dt_; + + sum[0] = 0.; + UNROLL_S(4) + for( unsigned int k=1 ; k<5 ; k++ ) { + sum[k] = sum[k-1]-DSl[( k-1 )*vecSize+ipart]; + } + + //mode 0 + double tmp = crl_p * ( 0.5*DSr[ipart] ) * invR_local[0] ; + UNROLL_S(4) + for( unsigned int i=1 ; i<5 ; i++ ) { + bJ [( i*5 )*vecSize+ipart] += sum[i] * tmp; + } + UNROLL_S(4) + for ( unsigned int j=1; j<5 ; j++ ) { + tmp = crl_p * ( Sr0[(j-1)*vecSize+ipart] + 0.5*DSr[j*vecSize+ipart] ) * invR_local[j]; + UNROLL_S(4) + for( unsigned int i=1 ; i<5 ; i++ ) { + bJ [(i*5+j )*vecSize+ipart] += sum[i] * tmp; + } + } + //mode > 0 + for (unsigned int imode=1; imode *bJ, double one_ov_dt, double *invRd_local, std::complex *e_bar, int jpo ) + { + + int vecSize = 8; + double sum[5]; + std::complex C_m =2.; + double crr_p = charge_weight[ipart]*one_ov_dt; + + sum[4] = 0.; + UNROLL_S(4) + for( int j=3 ; j>=0 ; j-- ) { + sum[j] = sum[j+1] * abs( jpo+j+1 + j_domain_begin_ + 0.5 )*dr * invRd_local[j+1] + crr_p * DSr[(j+1)*vecSize+ipart] * invRd_local[j+1]*dr; + } + + //mode 0 + double tmp = 0.5*DSl[ipart]; + UNROLL_S(4) + for( unsigned int j=0 ; j<4 ; j++ ) { + bJ [( j+1 )*vecSize+ipart] += sum[j] * tmp; + } + UNROLL_S(4) + for ( unsigned int i=1; i<5 ; i++ ) { + tmp = Sl0[(i-1)*vecSize + ipart] + 0.5*DSl[i*vecSize + ipart]; + UNROLL_S(4) + for( unsigned int j=0 ; j<4 ; j++ ) { + bJ [(i*5+j+1 )*vecSize + ipart] += sum[j] * tmp; + } + } + + //mode > 0 + for (unsigned int imode=1; imode *bJ, double *invR_local, double *r_bar, std::complex *e_bar_m1, std::complex *e_delta_m1, + double one_ov_dt) + { + + int vecSize = 8; + //mode 0 + //i=0 and i=4: Sl0=Sr0=0 + //Sr0 and Sr1 need to be divided by inv_R + // S1 = DS + S0 + double Sl1[5], Sr1[5]; + + Sl1[0] = DSl[ ipart] ; + Sl1[1] = DSl[1*vecSize + ipart] + Sl0_buff_vect[ ipart]; + Sl1[2] = DSl[2*vecSize + ipart] + Sl0_buff_vect[1*vecSize + ipart]; + Sl1[3] = DSl[3*vecSize + ipart] + Sl0_buff_vect[2*vecSize + ipart]; + Sl1[4] = DSl[4*vecSize + ipart] ; + + Sr1[0] = (DSr[ ipart] ) * invR_local[0]; + Sr1[1] = (DSr[1*vecSize + ipart] + Sr0_buff_vect[ ipart]) * invR_local[1]; + Sr1[2] = (DSr[2*vecSize + ipart] + Sr0_buff_vect[1*vecSize + ipart]) * invR_local[2]; + Sr1[3] = (DSr[3*vecSize + ipart] + Sr0_buff_vect[2*vecSize + ipart]) * invR_local[3]; + Sr1[4] = (DSr[4*vecSize + ipart] ) * invR_local[4]; + + Sr0_buff_vect[ ipart] *= invR_local[1]; + Sr0_buff_vect[1*vecSize + ipart] *= invR_local[2]; + Sr0_buff_vect[2*vecSize + ipart] *= invR_local[3]; + + //mode 0 + std::complex crt_p= charge_weight[ipart]*( momentum_z[ipart]* real(e_bar_m1[ipart]) - momentum_y[ipart]*imag(e_bar_m1[ipart]) ) * invgf[ipart]; + std::complex e_delta = 0.5; + std::complex e_bar = 1.; + + + for (unsigned int imode=0; imode 0){ + e_delta *= e_delta_m1[ipart]; + e_bar *= e_bar_m1[ipart]; + crt_p = charge_weight[ipart]*Icpx*e_bar * one_ov_dt * 2. * r_bar[ipart] /( double )imode ; + } + + //j=0 case + UNROLL_S(5) + for( unsigned int i=0 ; i<5 ; i++ ) { + bJ [200*imode + (i*5 )*vecSize + ipart] -= crt_p*(Sr1[0]*Sl1[i]*( e_delta-1. ) ); + } + //i=0 case + UNROLL_S(4) + for( unsigned int j=1 ; j<5 ; j++ ) { + bJ [200*imode + (j)*vecSize + ipart] -= crt_p*(Sr1[j]*Sl1[0]*( e_delta-1. ) ); + } + + + UNROLL_S(4) + for( unsigned int i=1 ; i<5 ; i++ ) { + UNROLL_S(4) + for ( unsigned int j=1; j<5 ; j++ ) { + bJ [200*imode + (i*5+j )*vecSize + ipart] -= crt_p*(Sr1[j]*Sl1[i]*( e_delta-1. ) - Sr0_buff_vect[(j-1)*vecSize + ipart]*Sl0_buff_vect[(i-1)*vecSize + ipart]*(std::conj(e_delta) - 1.)); + } + } + + if (imode == 0) e_delta = 1. ; //Restore e_delta correct initial value. + } + } + + inline void __attribute__((always_inline)) computeRho( int ipart, + double *charge_weight, + double *DSl, double *DSr, double *Sl0_buff_vect, double *Sr0_buff_vect, + std::complex *brho, double *invR_local, std::complex *e_bar_m1) + { + + int vecSize = 8; + //mode 0 + //i=0 and i=4: Sl0=Sr0=0 + //Sr0 and Sr1 need to be divided by inv_R + // S1 = DS + S0 + double Sl1[5], Sr1[5]; + + Sl1[0] = DSl[ ipart] ; + Sl1[1] = DSl[1*vecSize + ipart] + Sl0_buff_vect[ ipart]; + Sl1[2] = DSl[2*vecSize + ipart] + Sl0_buff_vect[1*vecSize + ipart]; + Sl1[3] = DSl[3*vecSize + ipart] + Sl0_buff_vect[2*vecSize + ipart]; + Sl1[4] = DSl[4*vecSize + ipart] ; + + Sr1[0] = (DSr[ ipart] ) * invR_local[0]; + Sr1[1] = (DSr[1*vecSize + ipart] + Sr0_buff_vect[ ipart]) * invR_local[1]; + Sr1[2] = (DSr[2*vecSize + ipart] + Sr0_buff_vect[1*vecSize + ipart]) * invR_local[2]; + Sr1[3] = (DSr[3*vecSize + ipart] + Sr0_buff_vect[2*vecSize + ipart]) * invR_local[3]; + Sr1[4] = (DSr[4*vecSize + ipart] ) * invR_local[4]; + + Sr0_buff_vect[ ipart] *= invR_local[1]; + Sr0_buff_vect[1*vecSize + ipart] *= invR_local[2]; + Sr0_buff_vect[2*vecSize + ipart] *= invR_local[3]; + + //mode 0 + std::complex C_m = 1.; + std::complex e_bar = 1.; + + for (unsigned int imode=0; imode 0){ + e_bar *= e_bar_m1[ipart]; + C_m = 2. * e_bar; + } + + UNROLL_S(5) + for( unsigned int i=0 ; i<5 ; i++ ) { + UNROLL_S(5) + for ( unsigned int j=0; j<5 ; j++ ) { + brho [200*imode + (i*5+j )*vecSize + ipart] += C_m * charge_weight[ipart]*(Sr1[j]*Sl1[i]); + } + } + } + } + +}; + +#endif + diff --git a/src/Projector/ProjectorFactory.h b/src/Projector/ProjectorFactory.h index 5b1f50e37..dfe5281c1 100755 --- a/src/Projector/ProjectorFactory.h +++ b/src/Projector/ProjectorFactory.h @@ -12,6 +12,7 @@ #include "Projector3D2OrderGPU.h" #include "Projector3D4Order.h" #include "ProjectorAM1Order.h" +#include "ProjectorAM1OrderRuyten.h" #include "ProjectorAM2Order.h" #include "Projector2D2OrderV.h" @@ -19,6 +20,7 @@ #include "Projector3D2OrderV.h" #include "Projector3D4OrderV.h" #include "ProjectorAM2OrderV.h" +#include "ProjectorAM1OrderRuytenV.h" #include "Params.h" #include "Patch.h" @@ -94,10 +96,18 @@ class ProjectorFactory Proj = new ProjectorAM1Order( params, patch ); } else { if( !vectorization ) { - Proj = new ProjectorAM2Order( params, patch ); + if ( params.interpolation_order == 1 ) { + Proj = new ProjectorAM1OrderRuyten( params, patch ); + } else { + Proj = new ProjectorAM2Order( params, patch ); + } } else { - Proj = new ProjectorAM2OrderV( params, patch ); + if ( params.interpolation_order == 1 ) { + Proj = new ProjectorAM1OrderRuytenV( params, patch ); + } else { + Proj = new ProjectorAM2OrderV( params, patch ); + } } } } else { diff --git a/validation/analyses/validate_tstAM_04d_laser_propagation_Order1.py b/validation/analyses/validate_tstAM_04d_laser_propagation_Order1.py new file mode 100755 index 000000000..ab71ed39e --- /dev/null +++ b/validation/analyses/validate_tstAM_04d_laser_propagation_Order1.py @@ -0,0 +1,29 @@ +import os, re, numpy as np, math, h5py +import happi + +S = happi.Open(["./restart*"], verbose=False) + + +# COMPARE THE Ey FIELD in polarization direction +Ey = S.Probe(0, "Ey", timesteps=2000.).getData()[0] +Validate("Ey field at iteration 2000", Ey, 0.01) + +# COMPARE THE Jy FIELD in polarization direction +Jy = S.Probe(0, "Jy", timesteps=2000.).getData()[0] +Validate("Jy field at iteration 2000", Jy, 0.0005) + +for sc in ["Ubal","Uelm","Uexp","Uelm_bnd"]: + Validate("Scalar "+sc, S.Scalar(sc).getData(), 1e4) +for sc in ["Uelm_inj_mvw","Uelm_out_mvw"]: + Validate("Scalar "+sc, S.Scalar(sc).getData(), 100.) + +## Performances non regression +#timer_particle = S.Performances(raw="timer_particles").getData()[-1].mean() +#Validate("Mean time spent in particles", timer_particle, 2.) + +#timer_mw = S.Performances(raw="timer_movWindow").getData()[-1].mean() +#Validate("Mean time spent in moving windows", timer_mw, 0.5) +# +#timer_syncdens = S.Performances(raw="timer_syncDens").getData()[-1].mean() +#Validate("Mean time spent in sync densities", timer_syncdens, 2.) + diff --git a/validation/references/tstAM_04d_laser_propagation_Order1.py.txt b/validation/references/tstAM_04d_laser_propagation_Order1.py.txt new file mode 100644 index 0000000000000000000000000000000000000000..4b2dfccf07f1d9a1abaf545b05db0f6e54cd8ac4 GIT binary patch literal 33947 zcmd4ZcOaGjA3uDQ1|_Q$l9g0Qk?p*Q$mV2=?7c3el$LRjG9poFQlccQVZ6x5&MMhk znc1Vnt&iXD```WdeV;$hACKqTbnAoU9ZB5XI5Sz%vWVt)dOV#Fqav($>e%(~geVYD}~C z`1Vo+od%T#bqGOECbkJ`Y%r!?Q#oS0wzjtV>e}tUI}Iuyg26iM1Qiz*mDsuD(dZY3%&+~pM`9j*nW<8GkTVC=)YD)DH;(=b;u@E8;CzF-LX zl7wC2Y}_B(UtrGAFYi+qU*b{Ap@DRtm)OCykGEen1v{LUqIU8pV?;ao+4$PuEj(`F z5Hq}sw3WpU^<_cGoSgr)n#KhRI``HwZj1&gC;qoR&mO`c^Wuiu8_AG05Np-$mI~i@ zk8-56zJ|P(i0FfS8Sr?B{H>_&4an}{r0x0g8XA}E{Nr3wp-U*~e0fL`EZ@!I5$?MO z;k9lx9m~PM>#{>N;>S6h(|2Xv6A_R0jNRci*-7Z^^Kh2^^eZ%3os_EydW$Os#o?Tt zSvY&dzB!^l8+((vat;n=V?#wq6Ng|nmShh@^K*h`0%~5BhR6YyH@@wX9h@5^KvD~XTvMG#ZY~(cd+fD zbf}R02gr+PuiB0ufX3#N%IoI$P%ZUC(@#AYngvhKpBQ}$8+T~F91D2?` zFGi1x=RyiPiqP)#HXh>#h4|rCVF)$r2W%KW#r7usEpl!7cH9wBEQw&D!Al9q*mDFg)POvEu&eK9p9JvrSd)_<#0L8gOUnN2;9R~&QO68fydRHZMP-U{{#H^{@$6?*uo{f?)2YSTsK;ikW%c-h zd&kFwhz7hCn|3HmqyfvE^G(Zs*Wq%n5|yZXCF*?BeD$*Z69&(p4NJe7jb_C&-Aew> z$Z+a$AiGKyq@>4}9Stpk-;43@=rgOKNmA>WFINL(Jh>o$#jOcmSL|;5s^1Ju_jpIz zPc}n`Q6jy^RO8xyBNE*%P!9ob-I-R?D&Q;I`7%}20;m#HJ9(w)0c;l75mV{$1{Y6` zHW_J`Aiwvq7+1MkWK+(vENE)PPUeWRO$WYV(mS25y<;u-rD5plJ^t?)aQiupNZU88 z;n?eyy000_&fY(Uf4-nw+^g#sK9u6v^^dI@QW+Q|t90Gk^BJhJ);bvre}aXCnCc6k z>fpewi^=Al&A=VXLdV133Kh49^EE!TLyAfd|H|VIh&FJ&;uY2o-7Ona4069iDtysrG-G(-nP4IhdylOyLCB53vq%v{hMD)9kLaQ1Q#)U!*7>9 zJ_(m^$B)7rv&MV7@Um-D_eYl=V6Ym zYFg!&24L861a$sv}7m{;q&afo(VJzc4o9k}}aMFjx zMk#L))xRG78)ZI-Up~0Ljq>QngRPUnSsFcP!=gtmZSoy$r1#Tiu2v)U(9^27?oZL| z$nm`i*FHg}fq_idiDnq)dt;@a)&=|g1BG5C_k)PRm-h#9h9JOXyuUN}C;VQJ66YxY z38p75J(N-U0oB%_b3x(*Ky~1a{OoWy+-fl%$-MFvto&9JT#tW(UxAY~RDUz@$v`S? zn_2_fTa`?BHz9WnmKioe&S~qbD^B#5fpHRzJ`ubjCfOjGKy;yH?r(w zKkD)mpHH>bWoQlHcCT2IpoLCU(v%o~;aHCuR|i7g72LtSnu59AN2|bn-02L@-!{-q z>LBin%U6`jcCpT@oQUi2KJ^JaesTWKt%Zrd{SO>WXA0Kc@C!eGdHY%9!Zgy$@{Mm* zn#H7<7GwRQS!7_VIC6Jv26xmSkW?L=!XnoY%krdgOwGQnvw8O*PJU6ozWDkpUJMgb zxAwe3sQXczUd&y(%Lf*;=p2@qv;I{6S19y)NfXw^l--jJ)Ks0t|;LhEZ`2Or|tq(l?$a%?2(RA-G z{8OIU-Iq3t&#igD)pikc6nBjnx%t4sm?=;w zw`AXUVhCEQEN-wE*F(f^iM?DqDp5l?X5)rO1K2O8Fo{h2Hhw z_*6pgIp>DI7))lf7pVJ-bQ+exmA#A~cUW3(d$EK93ig?+^Rw6+Z%-q?X9D?z&+lMY z=|R@0cPDp;CZURX1N}v3q^~51o36k_^b=_fvt<}J>PvB1d%cz{LT8nyr-3`JzOP{IxOQ2z_qY7jCZGzt z%BOkcGrG98-G~ty#JJ4S;mGJIjO8n-w0By-QoU_meFJ~+xV&#?ThU(S){J{mYvg~tL=TRkbE(kRy@pJu+*4Z6>ILAIAGQgOP(q6H_yasJB z>#kj%sW1jlQpitiT4sUv+Ag8{qQBvUdgPyP{eOYqlurFj$;S;d&o7j`a)Eu*n7 zj0JiuV05EA=*v%GyCqg!S?I@&kLVL}@{2Iq=2U3TqZWARTCFBcKMIR`cru0cXJP5O zxT50V5{%wpGgE5$3xTEH=hOyPf#+D(CC>T8*)Ik;c*an+DVA4=WtPB?yX5|c&uhDW`f$88zKGoOS1 zVA;SmUGC%+w76SPsq=LeZK#UmsmueG z9}c2>{F|&u{!09IZ$AEeX(JeCWe4tb{s|8+lzQ%5oQ6qC4#h2e5#s6+EBx;)14+TM zSD%!H$TNSi|$MK$yq}24lEIQY2P1skugePQW-;-$n;-(ZQ4rcom z6pYiq(bxMIGp>0@3f=mH-NIwe_t*CK#Q z_VbPXM5om!M~;|IW}sh2`u4B03y{W17@Kx2!HwCov~KyoK}zF`hsXF5>y{PWsM zmhWevcHx1jkIy&+v`38xeCY+>EdryX^QB;)q{K5ST!+#j;i`8(_2cBtYwebACNM5x zGr!T-8I-NOtFwh+0p}LKb4k@LVu(smuFL8o7BrBCFFs$um~X;mQysJTsZ0A$QNtv@ zHve=l_Q4MnH~;oa_5OFHihn<7>q0^DiO36KJ6pgkH=~v96u$xfC-NTbZK%_OT8 ziT!G=l~;4HD4}2eFGn2&{<>?QAjHcR7gQ&r|LozGBpuU}Y^Dnhwu&@2()FILfzl0L`crELJbl7B9 zpwu&z&=4+tQ__q+B-!CVC4IkndviI)V{K7FA)m2XSDg3~~yxoCr8r=kx zl!_=*_%5C9smahU{QB(T)puV;kS1cn`xfs2x)hw}yqNbLx9e|9OJpj-iC40U!ezyf zY3gz3B2z0A#L$xi3j4u>)B4Q)2g7ih_uJ02nsMlEvy6MZeG)>x4;_^`J_({Nny)9) zCqO0gch|sfj2ATWeS@8c@%g7Clf%w4DI!KuRT4)g0s*rDQO9;#4-r+s}3 z0_Z!iGEGIy@ozsKxoW}9p86BT2Tv(%C>lj~+0+|ybmMqM=6l`LwK1$d*hqdLGlG9O z#SK)Jrz~YYvFnSGX8#6i~ zb2G15T24Pa5mym(R{Q}{8<*0=G>2jDl#Y#D(J)9(h;R$Wt)1`QzHNMGU=ZY(HZ%Vn z?FINQ&TKN-3IeIcw}a$s!BhOJm)-d!crZ!&-1qN#)ag>uyYOb@dmW{bd9Vu=3pb?k z_VnRpZg2ho*FoH|8k2oDehBA`kL0YF4Pm_L;cS+|0puVQ8p-XwNK`hv{Ay~)nTKI9 z&m9}_#mixaPm)EbHS7QU%B^Q0k>h-L`D7(5C?&AXwtj^><(;e@H#=d&n^*jYl6oQM z(7um?M*T38t-9%5c|Qbk2H(81&PKpD z)6?~0+oS#JyO3-`}^7~hS4M-E26KiiITo9^kW-u#L(DuOR88EVn0kWQIhFb^dh zL`R;7-G*p4_vDGRVz}gKkZ6173$V)Vp)q>&6)yjcFZEk)1)mcv>4y?KAaQ)bGRvwH zwrxz5wb1PVU!70YGNG+d-OLqsQ0XghIrwQ;Cwu`0dO)Wb$;~nKG< zN#Dw`w1>EP&jE2{X{*Sb)End(rI|SO)fUc)MGscHW&kgn&*$?W3gK0iUi2=*GWb|> z=B80|CFqOA>`RQVfxM6Q9~-XJ0^xF?$U341Y&pZlxAm>ng(Q{c71GM!kd^FUHWa~T zv2ET%zq6oy#|L{q*?Vx|oj-H3^L-RpHdzpG$-*?BkD* z%JH=4`Nb8%3bg1d|Jlu7fiXRIf8{+W$8MDlr^JF1e66HiJ>-YD@>w{mElNXl?d$^xDCXC zx~GOr(I}R7ohg9)62Cj`fBJzo8+Vtg**g9BfagqG{8=XpaOZdL%iC8FgU=RL3{e;1 z1Ci5IhiD4%;2r%Rp5h;oHRSr|k;!ZnJ@b=KWyfp0da21n^+h5kILl4vkVt69d%32x z{655O&y)OxsnF(M_JA`l6QnH_9<+VTg~Jy=T<((j0Ef&QLPl6V!tmrTnnCygEsf`D zHjL$hkT=#Sny2^~E(MblWANm}8{4ec+L(H2`z14xXm~eH6{zx(&Y*6P=m7F^p65sx-V8RY~*(@f+2jdL;Za8Cz!nB)UEkuj06Q}zkG@(`b6 zmwyI(EUUfshMz$p>&MPZGS6Tsuduko?FqQv?^rOjOoSnyq(gUj<;6+Oo zi=Og=+NngP>v%_9zCM>Kf-4EE4c(tACJ}NFtJMI`8%ueT4PH zh1*1|A0ZaJsY(t@#0Sc9wO1?F#)UCQN#@4lP{&qNgScDBK5#9)^7J*lqncZBwNnpU znyyzyKlFw4JrCSPIYPkq$ovX;M?r^1aby8G7QBQn_85MP10nV)bMr@cp`OJ)M87l+ zl9V|Z^54Wl0>^U9LVPq-b@KLH9k~HQLLPJT-Phpe_1F}XpN_C|rTyOg&a*h&!+B$}Tf ze#~`$)KU<1yNS+l4Tk|O565HI&YRFPJ$IQ1xCLEPUD1JX6N=bq2JJtG!3>FGC)EWC zJdB+$3>6K4_x`o-U-pFf$-#a&312F#v#bTx~0M%Y3)!#jQ9jQjxFLdnm$Le#= z4_AG>aigbcW=xR{lG%!uW9>yy%abRiQ^*(+Ue!PRK5h#hHXpYMzjB0|9h%WO`i@Xq zIiq~f#15oddlmfntYGWUgF~W`CgA*IU(nNTJrGInoy(0^g;qb?xBlVspuXeV&Cz>D zfqN*$rD795cqv*oN>%R!lc=Qt+q_Mnd{^AA`_nvWyQE^F(9c1V^_@``4fSSH?{HAQ zF#9JG=gbE0YQ8rlewIJ)8Cw!a?U%O=%Vh4LH_w%f@o7X&(XR0n$ zwJW2@bIh*h2%9^}mmT-HoNKNj|H*SN?~nUJ@?wnfHSW~g<@0+i&nDa+lT%2e)%G>4 zmP;!Vx92u`DL1e)z$xMbMJ{35^NC-sja=V^iQ}&@1-VUb#qz@6*yXJKnZovMpOej8 z@=|y*_C=On#_>nG?F-r88~A13L|&71Iy6WsN;Z^@U+%s5bVyj1S?;M~&bAG*Keq1V zP*8C~FEf|@6*!Fj!YvJKzDVlSI*77&pc{!&GJKSWq&>k5p z6iC~&b%kxFuyG(IZ#9@2c~pC}%>6f`^kK1@I&MbPc{VWN$i$C9jsvZ~5`qkT+{h*IW5b2pJjfocl;Pmci`J2%Pu@lBL$WPjch@i{KDMf` z`sTxq7M`E#T}OB0mT2^zj|HHfF`lhi{pSj)3MX05ZM$E% zcx>_WgP1&jLGIsK6 z8XfN~jNAHgKi8Ax!b7^;NwlZwG3i!pg6c90TF)K4RUXBOieIm0zxL-tAF@#4@WF%V zz|o{|q+S@6pANqKt}cun%wEr5bPC~-&um;_wn8}1Mej1;FM_}8>-IgvwOp4s^4p&S zsBkdE{Rk^BE|$L)*It`nR^+J!jw5uaD<)BwzM+})bN1X*o>yCeC9f&&B0n2s4$_?7 zeTNsQs%bW;x(Y&T(IB%<^#Mo`7Gc#$-VY5a%x}b2`C;#?U9aW5c|l`#wQKeCPO#+S zso>LGB1vx*O7FI%#SX&R>97zNdRBQ>%LyFDU~9qCo!#PiucLM<@R|&6GQLs}N+*v& z-FgPeFBGtlXM8LDUq$rYT4=W=RT1~R9y0B*S44@c3#CE1B#dX|xF7IX5?4}#7d!#@W$(`a>xXy=`S@Zp^2hXX*ZP`&ecsQnlEbL zDL6J`6Q+fe$zsClceIc#xGd$jx&}Vn9w;64Rs~fW?`#yBAmRO>?z)3Dhp>e0Vb#+a zHfYHw=Zl|{0hYbjDfN=75bsYT<|(fYI=^3sJEZBsB#Y;f1A(XE;5!bX+nM^XkylD* zepDZvL#l*L#Pva&qn`a}h8_q*p8J#ins8pI;e@4<5^TYIf3ELG;iuH0_$|G(D0$(- zkAfW%$TfTE@!kj}y!Wz!zVEO)PRQNTW*E`Nk9O?G(&crLxxuoiT}>OWDS5O098<@E z4o>-ovnP>#elj&NM;8BPI;I|XJb=@ak2h)Nvcjv=lkw875@68t>Vn~6C8(gjs#-6i z4!`s!mM@TXVQgu~pv%_N5L~4BgHBZ+6eRV+r{(nFDxKn>i$R52cj}o{mTj$mZhZA@^Dk3=Qktz(bcdM_TYGFJ6*5?@4!lTPo zBXMuEQ0j8~xs=VCSRAcdmT0PiiYw}7{>MqUQr`B0`}9Hl^LM9x+3IfSwwTMrN-41C zi*u^ZQU>|HE$P+!HK2u2s>b)PHl)x}^{~q5!jx`!ZukXVP%=D_kRYuK_X<8mR4!=2 zymHQ{nyMP4&lcuAP*Z}X7It5UMhQr8&MV_KVuTnxRsEbp1nZ;e{&=e?pl%D@FL{b8 zPUNRv_VCceyaAn)?5Km|Cvq#~*ItJNq0r@gSzRnV_aiOcTN|h5jD=}ZHSkvMa>xEc zm9^ulz){2Za`;T_Lr}RmAIg3_`gmK%eyH8rM$fn#KqfR>u9@u=SkAsMzC)`Ga$DUQ zoAq^}eZlwD2}?a_ytTqOJ){R$ByM;GyXirXi@~y`gf86XqFHE|(E!~irBh*&%Amf# zJ;}OE98@Xh4ubAOg%rYRSKxUGWUr449WyzJ%u=5s)V^tA%(YPgo@rg&dagUIY()=u z#mVtY29t5-;N)ZL*JKP7IyT&GNyfCh8(B@R>Z0Y`jN5Lq8tR@ODZyv`!TsSlz9oYJvR z*7oD44_G!!=)>dJZhx%W$si=MJ;DEv4!o5s?sjKZ0cVA`hq6}0A^qC-xTYOUXnG@n zX~<0;e~ctm#ZIcBb?$MM6G?iQe#Sp&(pMh?&$?{r`e=ZzdSOz^l7`5!v74ygXNaX# zIuBgG7+~^2#=l?OPUGd(D@Ql-XyLts%#SY*DWDKzm3TPoUR2axN>$1hgQl-7DY8_m zATnOZD4?YWqa07~jHw#Har1;WK_(+;4A>i*R$>GnZKQgcyp2KJYsZV=NMrD~U$*qT zZ3w$pzg6(4kRd)dUU7fjDKIhI5I>f3917<03)$arqM>uK;cbZ%NYm~|nVZwWzS-r| z^5+fkf~4N5V-Jl{{A01qMR%J%L%OK`nP|IPpaOXw~_0DcJ5a8}haL-UFec09n zBK8N_1csPF+le&!^FPdBMU$?v(&h~Od03Q@6LSWhgzfxcEM^J~t)h1}f7OSYQ?3Q> zjT-Q`{VaWIh$OIPm&-T3;K#kBSNk%mPhvpW?k6pSWPB~n%Qtq#1pVLsjU&66qljQg zTi7lO+_Trvgd%K#J3<$#nl_!qY8AD-_;VIGQxR!4GiieCgNs%+zsUIU(nalKwkNTD zzw)ssXSr~RRe6Esf*edp-#VFmLkC97#kd1Xjle^`dLU-n9I}3{M2<^aLW}Xwyi1an z5G{S-$kRAW@Zo9bzL#$aBZmTbzUNwiAe zDu7|0d1W-Qs`y1(z&Y%zJ~jmlGL%W1Bai-Xo7po~_-DlF(Z~jCe6~yCxcICU5`E{z zXv3_rk^j}x81_vN^c4KVc&Y}HY=fjz%2jkVg?KvrY=SI0*- z;6~oH$u7qV=-cZ2PP3SS(`N@i=0kcgNSz)2x)#HTjcEnv8nb zCt7l~k1V=&{%je|m#GLoA5Go{mg=v4ZhP;EV6+5-(R~#ULhaziQ8)YDQBF|*J~WJY z;0*dK!3tgc&hY3x8&lv9C-``M$5fb}13Zdx9;#%vhGefNs%meI;EbF%Wr9Zq9)&B* zP0H{g|C=#EvY9UWpTC)P<)%3*o)P|Hoo0)-_pf@e8av_ar-QB9r(Lix+9}FUz!j^T zeGbFV9{i!1#04V4uaaEEK07qa5h++a{cbLXpG zS2(!0uV0wa84fRfdYTkq3x7s-ZmJVJ12GDQCgQI&L4y1I@)?n>_%k)eE9!|3p3-wB zKUT0n34gssb8&kdp_Qmn8FoSMgu|0hp1GsN_ez@|OrF^N+;)t|%oBaCrMBE%b;G0P z8Z3=*PPnOAzT(nhTYPDu9%xUP;MB#|^y*$ktZZuLRNtill(<{twR9$6BgOl0`%YWP zPWV14knRL>!<#=1>bQdyhdNE!M^EUu<|O_r)f0AhrhmDz)f2KfPdwRn!wurZhK=Tj z9Kcf}|5W=sOQ_l`Us2nr2gfhgHjM3Jfb3s9oLOqR$Z6rS-#Ww+PgqH^ru=b0e!UaH zI2S&elDGDMU1{r!nig;YtysI$H}-f! zTgST-n>KpE0#ozvR}Nm#)G;|~7)5{|O-;c|1y7JY`YQa`9v9Gi7G=jz*Z|kaBR?ix zL%2(Bt2ZSbh1j8=_uTBXP<^{*jij$73KVRTJrL@MJ-0lgOhr79;IJrd3?i^qz|@Gn z#S4qXI+q=7ypd3T;B?g13p*XuUQzjYBKy%eS#~RDJbEbKYB<3P`4ZH#f>d?U%Yo~r z%Rzqdxw6@YQP2PmZ2GzJE|V>^Jt4h1u=alHpQYC(!D`= zc_XF$w->zB{ZJzoL%?q(C0;XEchKEg6>^u_0j>x5+zG8VgCc*m_PQGiP#T{lzfz!r z^|#6sv#BkR(bAYZ#K8fd^8WorlJ>wiC&T6iWW2B`@J&a_TW?HLzZYt*;)7gg-*#up zd83~FeU&*YVy$2~sP$!Jtq|s|a{57_E&BOb{CWS-0GVeURm;d+xc5GBpRr93y8m!q z5H7X`8kw@5+DI3m>yS89zc!Dlz|M}j9p2y|cc|=^k`J^>+yAXD@P@~G{3M>8^aA_) z-sOqXp1|r@$sN14)Zl4tkE{4uOK`rBQWBWl~Hwq>@z+W zida11?p{TXN($IFVQR6o^U@l>+t<%XIN!; z+#nve z1!ieCrwulG!`1jxdIr0_z=@i1)>O(9Ugku&*U`Da!rrHwdfaS4tJp4>#nS+G8mTig z4jhJSKB`WRm0D;Tbxc0n&k`ByIL?`EaKyE8#SN3QZuodfNigAwCzf!nOt4dXVMil_ zh1@cMN%^d80mYtp?^E!M-77czz&i5b&ie&~NZ-?%FC>G_C!< z96IcT8EM%IJXh?Hfnhbs!`cGBC#o~*cK<`&#Kl;XSmrS&ylF*3YsQtQ*paoA*Z3-WKX;kB(8j4O-Zy{ zTYyCDnO!>%OTJ>+QEg!WO>es@LQgBet~)z#r%WD4q2%-Xdajf4ITO+t(OMwQTk`h3 z2{u@!wnxI#kmezuYG~w zX~tDsQ#?G}-r@z&P&XOyk=5 zyj*IJ(SfUL{qA=zlv4);(2te!ea~7wK~#o?hMvU`eL@{5uV&40(sw`gtKHVPqomeg ztiuNPPd8_`9<{}hWnIQxaT^p9>+3zoXN42brssPPo7!I^tc*y|DfBhv zPurt*0Nox4&L3Ek2Bx4dFLf-nAvop0;*TW*ur}>7Q@dvdM^kreRkB*bbLIC{50$N; zob8*k`AIA2vShlrL);Q#R@1W!PM(217f0xmv<%@2w9*|)(1I(r>dg%&WucZ*bikyT z6%TeRbF*wbf%K+a40A_x@rkFP5m*{xI{yq$!cS9N^s&=oxo(cT#llO%_FCX_qfmWD z(;1}R*8B5nlL=mIY!G#KGeBRJax|pU#k;bH%y;ioLKfE+uiG>Hcyzd1*x;ThbXVR7 zs?8^1`UF){Hk&RCh)0Z_x}gupWq-XYOE&_kWXVlij7-2|ioL?n$posBc!Q{ij3AU_ z`hNKjeK4O@v733J3yBWVTg88g0Ezykfc%cd!h&FSPnHHLY^i(OH#n(+ z1^&(kkg0>;%amg+cb>++6NAIGGzPfCyZnb(qyY|GQc(JZ`nWOgL}}L>GV(n5&`isw zgMOn*x5GJ9QEOCfZH<#GjyLAhc5dN8rUvs~_M-hzC|wdV!zK@3>pKp5g{wg5_hNj4s|ppbPsC7@mBCOkr|$6Dew}qbN7uGP0;Y%ULOS~oKv9H@ zgM#b;DQPG@JU5dc$%QBDBy`2Gsw(dq-Ag&7UM&zm>8psVbVA4Z?UeBCrEzs7+&UV98S=nl6|}!~ZJmnWJH^cn8x){!=-p2` z5fbE{ta~ebO$HVyK5my(kAuy%eVZ>;2!kr^RvPJ@`=BGo5C5N|R-BX?(;9GS9hljJ_U# zslJ(jB>y8oclTVPh{rLw77}%2vQrGatAB?y{uP7cg2%HtT8_aZ@dOjIP*JdX|LWrP z>xW>wyP@m7^Sq$?{$vyH1{UBM>3UFosDPBgB~dXK$%G`w%)&QI+!&-2AwSe{0P8)^ zZrZg&7(K+ZgJm@@v;i)5Cvb`&)#uzRpQD-0~oH97z|SjN!wH zX11N&8GG=>#RAP;VVm&49j88O%V832#n{XAGDi5^GDn?M$O$Gyz3$axe9)6SDN}o2 z09*@f>qO5B0=48Zp_jRW@bH=M;~wh+Fn09E(bOmWKsQar`7?tX><04>t}bOc~zWIQaC-%R@OP3*gji0 zox|5`VLn0X8M*5hb!-c8{XVaH%$fzZn))8CQrZJ2coe_WU)&2+8JP()jr)MEn6C6S zJ2wn)0WDP~7mR#8Dd2sT12*p5*&i{m3r>WLHSRyk1oMiv2Qiu!!h$m`y>3jBN}3N- zZS(3YJW?K?YiC4-xyP?R+M-Q|=dXQ>D4J))v$TQ@G*&DaRQ)4xPMZ}+4Hj%wYFY4@ zr-`CB2WW?->n+ZgE(A#hY% zifV0v^R~?D(`)ll4Lz{6kpZj*V|PhBW`blhVNRh;X4sYL&bR$3Gn}XDetRK^33&38 z!M2A1=H+}}(^74PqkD?ZA1J1V-W|5N<3=;2{lN>1XPe4Mgizh9L8FfWNiUx}3HLIf{-m_b1^MlmfBKU7 zeN8&dygSpQ++Buj7Tgr1Y(y#WfUu(__;*GuTX0 zso<^=AdQe>95{kYQl?2XKd6lQqvuHJ&qQ`Q)6A0uH*OT(erc9e&;Ew%XWA4gY_2kS zB4Ugr8`TojTr@~>x>+P(YST?B{4k;Vi1L-BwEcIz+Ws0++SzagnmffL2d3SLIh%7x zveLzaqok$gy{% zGB@j3%0U(Gl#9+myEa->DHC6oQjX>~;K4zos|eR1f~D4966u5hUBAd`rk z0wG2r#L2|*b=A1`J|m;{v;rYPAtcF!)VeYXhjFbqo>d^EDTEA}kX_fg``dX{%^Vd7 zISL_9CP?eb9lZbVRNo~90w{z6nNVEU2+(Wg{^6%UoS+a&WI}mex0V+t-!F(!AWl*U z6*8f^u8&7l7#>ZBDG;Y9gc_MpU)RyX*R=G3&lLy_3ZY3RwAS@nj&qAEaZiEJrVu)0 zLU&zb@Ad8Coyk@p^e60+2qOw%OeRd$ z_2LaqUhC2i3WOAuPzm*>!c>aQ<ZR6$6@gj>!ihpSlL?n~rO$cp z_d2;zfpDb|Ze+rJU2Xl0Q@7GI4ob^Y84JnMzq3 zdxb(=B@@@y^-1eFkIBZhvDYa?0GS9}*MsLD_2)XSjSZp@!DNE6u6J_ZmpcfpX$XY~ zB@wT&iutDn6rzw!;JOZ4SQ^&-(;^B{OeQ|9D<5-k9F6{(mQaXNGEug! z!FT8miTu-Y3h|jtRIID#rGV1Q|Fn`qRFR45b=`L)sbut@)=-FAGEuj#?>F;sdm60i z7Yb2NCK~?F78G-Z7!236kwP?)iRN|nb`aS8>7Rb35Z}l|%esEMxy_3(TGQ_oqLoav zt?T$N*T_Hrw4Fk9kcrN9{d2~x%-wiRyC_6Andn(pccx6q;E1=(;pP#Cz%*t*Fp#W6UzT|ghGsxiLrGxeQ>lU;h&CEhzTwodi|KeT$i}@`c)*U^&7PI{q@BS~|^S`*6o>cqdU(Ehr{Qs13mtj?`VZ(;qb|)sH zVxxkE3Mw*7u|@vRaTVAXwuJMKKkd9*vjeUQR}&tR``k3Ij*9)@>TnI=3coj1+@`rETnnzPW{3ap z0i!Lq~P{&2YP>~O33s^4ENdjK2(4-_8K z^UD2r?<5X;5Ih(jB7CIBjh$PwrFGat;bHJ_;izRXZyOeLa@Zr_k?<&C@9QTv47r}m zVULC*;W5H#3Km^hI4r-z9t)3y#|s}C?9!?6t-=m_0z46(Bz$9dpPj3(xI660@DzBe za8T7MH=c);cG%P4>F^BUVYLqVTo_Q!Vb6qT!Lx-Iom(5QwSbqyo&!h0(Za79?*29B zOC^Ur7oG>t7cNzA<;yGKRUGyLcp&=Tz-!@k!t<^A3+yD^Jo_&12+djHY;bnrR&`JH_tKmL z&I#ud{+i*%{nVPB;oPu`a8{SocW-r6oCnSe=Mz@D!`FqHUE%z20pZi}zOTYG7laGJ zg@x13SW_lIa}l^G>?V9<_3|o3G#7)5!zF~(&h$1xpIh8v54fc8_ETXVr!|*?OT%S^ zJLT{8V6o=1uqRwjIJ1*)M^0Sj;R>*qaKfsSA(md-8?Fdf5}q?P!`|$gE5kl;6=C9TzK@LixcW-ZUG0uErqAY4o%3XIT&sQw-(-?Ec263?G(3x+rsUH=Pxa= z=78o9xINrKSX~^{nW?!W913?5-rF_RfG(Ol!(HI6!b5NRcw3sg!QJ5=!mXp~1*g#5 z6Yd4~7Tz}S{i(Zc729whxUX>a8wtm^Y3>Jy!Tp8LzB%n2p*b8L07t0#M^}XRvInZA z!S*0E|LlsvYFVT`M9n~ZsG9$s6CXr4eeU=^lVuNshr=U;UAw-V^=PQ-ijnXrc(ia# z)emEmU#wx-k?UKRL;NDQn zo&ryWrwM;s`r)EaOVt(A;TiBu;aSOo)6S@(x?&bQ8=fOPcuDy-$&YWb>?k-Io-6Dz zb5)9|sa03ZgXhBwgtH8F>UXNO>WYQ%B6zXznN=}651hYl*-PN1@G{}Q4&JW3KK_Yi zFNasaD}{4bt~aoM`?r?83SJGb5%zt2ba=(Hsw>vQ>)`dmoqPlRtOto5_6B$(yh+%- zaQQ_MkJ32o&F~g@t8l$;ReGeql*M6hgSW#wgm3x`2rsfkb;VA27ra~e=r`YIapm$m z>^<;ac%Sg#=OZfaS*5ySKYRc_D4gm3IPd5osw)n`hv6f_>nGeEo55do#ZmYed|X%^ zLLd9qR9$fbJ_*MN$7YNfKmU0JhkXh@4WALNX@70q`LODWSokb_PT1?q(e1BGsjfH= zUw|(P2mh*;%XgXTic9ci_=>Q4{jGU)h3bl{@HP0l@Ga*-Yd&wN>agSB8}Loxuce>v z4Ct)7;ud@xz9YPH-f`>APsMlPd+>eX->yG9cCl4BJb)jtI0Cvd#5 zYu=ChV?QZ=3O|FN3#)emDs*yE=cO0$OZXq*)N#ui_?1-r3Vsd05mq-rS9nxZ{1#4t z-wB6adzR}*HO24Y5Aa7}ANAAc)>Qm2{0aUnobko3{?qCx{sMo6zX_`_a9&YY@pt$K z{8KpW!Tg-<>M8yO|ArGK`QMY~p!b`LuF?FLS`=v~Rx{8}BCJ}j;CanS;bd@f;TG*H z&AF{P1)LI2CA=+d@AdaIr-swOX@wo{LVjG)oDNP8XAoAoHS(zDjBqA6vv5?stTT6O z&H_8ZS%p)s&%JfM=4^0wIEV25zB6tt)|?a01v?8*_}L-DSk1X%7dVfwbGiqO2WZX< z=Yw5^pZYGZ&{=bSxBy&GczoWF^G!7uf(ye%gpXY*-q=%fQP>SGCOo6!UmILC7l%v0 z?!x8QRcW1Gvj2zdcXy0(I;yz>>;-!Z zZ+qG2+FH#O;Yx62;pwg?To3k#>kHp(9rl|;FaT}=k{H{3{gb-8of+%-3bo4`$lPj>grkw$YfxH;THc-GZa3D0UP4uV_4!NP;umTA07 zb1S$t+(!7-Ved*KHMfP^!6Cx68vKZArMW%a0q!W=_Sltb&YDBvPH<;obpZQ)s-?IK z+!gL7eDcnYs~0tQhkL+1h4cUXly05oUT|;N7VbGb>7o*%4-{6LHt3^1&k$BO-};7Xo(a!_XA9>Gc@yE#JO_?~qlF_@ln8RxJQtn^&lh%!&G`Fc zb;S$dh43Qb6Th}r+pKvpyaZk&ZAZSe{I41W=pw)hHvgTD()Tl|23!oP&2Eq=p^lInx#FE#%= zh@>qNt3~R~aX6{4v_&#FIh;aR+9D;K3QjF7ZIK2}3#Sv7wnz_WfHMk9TV#SW!&!u- zEu7%2a5iCSi|lX?IH$0*MK0JG&Mho$;R5G@^9oB_k~kAz^8Y!f+9| zsIat!8(a)7E-YpySlvG#$EN)aH-Z}r&+0I? z58MQ93O5rj)8G9C+#GHJ2MN2qPFI-U-4YIlTm3(d&G`13-tyLP8@R1-@SPe#^jteQ z1a2>^c86gX_1Ycaj&P{3x_X#$PID)?Gu%a3U2Eo!)!Y^C26q=$JJZZFntQ-K;aK2;_lLvb0m2J1tr*RwMZg2$LBjnHE_x3S zhKImIh1I2~(^{>?pg3`c#iOkU7HHPQE)Uo zS9nEJKB3=+_`5$ zb#BNUg^$6zclRv`$!k|$~rI4_*fWLG#pT)^ala3Q#`$wlCzu$#%n;Noxzligtt zxTMLY;L>mzlgq-Ma5ZllQ{=;Qb~afDgilOg;=BfsdMe3_cE@F!>}L z1D`VaG<*h*HTf)j4nA-41^6O-$>huM75J*j*Wl}LoXI!foA51@Z^L)syC&a*@52vF zeh5E;ADjFHj)$L`{0x2$zcBeF{15!fU_^rta@H_au$sgd4@V_R1f0C^@$-)%@d9GO=2uZe7AjO-=?Uhf|oG5>5rDHaQKP7EY(8*%v+g z^`$S;D-E_YsQG7KWK_!{?M!M0+L_h-WBaqfPH;adAOPO36E(4b}*%K}Ymp8cr>;-$9 zToJAWS2o!Pt^)g-Y{3q=s>#)0Ke)QdHQ<_XEt6}*b>O-t*Mt4x`X&d!4d8|*2f}~D zjZAI~H-Vd)+zf6Gw=g*fZV3mQ+zM_Dw=uab+zt*gxjozg?r3r-+zIY%au>KO+|A_f za1Xer$-Us-ux)Z5xG&t#+gm+)3aTf_FC}~_$Yi#c+l~MCF-W`W7)^y6YxplwLV#s^*yXO20jI!7S3|A zO@_EbeJ%S891EWn4h{8BlF2j7vd_Wi;S0hCFBC}Y_c6k%@bhElK?!zN3`wn~;z9+mq zZ1#kmna5i8efR@ci_tM&0t9ZP^*%jBqC5*iy9* zhj`7g?96Z$*hzTs)Aa$-yA)@Iv%%SgbCr5+MR-S9b`CfvoJ+XqmW)TeyrM1J8O{y6 z2uIa Date: Fri, 12 Jul 2024 18:17:50 +0200 Subject: [PATCH 13/26] Update release status --- doc/Sphinx/Overview/releases.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 3f4d92fd4..b93740024 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -27,6 +27,7 @@ Changes made in the repository (not released) * Prescribed fields in AM geometry. * Particle reflective boundary conditions at Rmax in AM geometry. + * 1st order Ruyten shape function in AM geometry. * **Bug fixes**: From 41e8dfa7637293b644cf3ed7c62c1fc94b6d11c8 Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Tue, 16 Jul 2024 15:57:19 +0200 Subject: [PATCH 14/26] fix again threading with python --- src/Tools/PyTools.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Tools/PyTools.h b/src/Tools/PyTools.h index 3a165f92b..a756f3db8 100755 --- a/src/Tools/PyTools.h +++ b/src/Tools/PyTools.h @@ -44,8 +44,8 @@ // } // SMILEI_PY_RESTORE_MASTER_THREAD #if PY_MAJOR_VERSION > 2 && PY_MINOR_VERSION > 11 - #define SMILEI_PY_SAVE_MASTER_THREAD PyThreadState *_save = NULL; _Pragma("omp master") { _save = PyEval_SaveThread(); } - #define SMILEI_PY_RESTORE_MASTER_THREAD _Pragma("omp master") PyEval_RestoreThread( _save ); + #define SMILEI_PY_SAVE_MASTER_THREAD PyThreadState *_save = NULL; _Pragma("omp master") if( Py_IsInitialized() ) _save = PyEval_SaveThread(); + #define SMILEI_PY_RESTORE_MASTER_THREAD _Pragma("omp master") if( Py_IsInitialized() ) PyEval_RestoreThread( _save ); #define SMILEI_PY_ACQUIRE_GIL PyGILState_STATE _state = PyGILState_Ensure(); #define SMILEI_PY_RELEASE_GIL PyGILState_Release( _state ); #else From e633efec035521d558f6ace9a321e6582b90c0dc Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Wed, 17 Jul 2024 15:13:37 +0200 Subject: [PATCH 15/26] Correct buffer intialization for Ruyten Projector --- src/Projector/ProjectorAM1OrderRuyten.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Projector/ProjectorAM1OrderRuyten.cpp b/src/Projector/ProjectorAM1OrderRuyten.cpp index b3db66ba8..da5818595 100755 --- a/src/Projector/ProjectorAM1OrderRuyten.cpp +++ b/src/Projector/ProjectorAM1OrderRuyten.cpp @@ -82,6 +82,8 @@ void ProjectorAM1OrderRuyten::currents( ElectroMagnAM *emAM, for( unsigned int i=0; i<5; i++ ) { Sl1[i] = 0.; + } + for( unsigned int i=0; i<4; i++ ) { Sr1[i] = 0.; } Sl0[0] = 0.; From 53c761ba60808415bae3498e84688eca2e24ddcb Mon Sep 17 00:00:00 2001 From: Arnaud Beck Date: Wed, 17 Jul 2024 16:06:55 +0200 Subject: [PATCH 16/26] Typo in Ruyten interpolator --- src/Interpolator/InterpolatorAM1OrderRuyten.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Interpolator/InterpolatorAM1OrderRuyten.cpp b/src/Interpolator/InterpolatorAM1OrderRuyten.cpp index 28313aa8a..708ee5ac8 100755 --- a/src/Interpolator/InterpolatorAM1OrderRuyten.cpp +++ b/src/Interpolator/InterpolatorAM1OrderRuyten.cpp @@ -80,12 +80,12 @@ void InterpolatorAM1OrderRuyten::fields( ElectroMagn *EMfields, Particles &parti exp_mm_theta *= exp_m_theta_ ; - *( ELoc+0*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[1], El, idx_d[0], idx_p[1] )* exp_mm_theta ) ; - *( ELoc+1*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[1], Er, idx_p[0], idx_d[1] )* exp_mm_theta ) ; - *( ELoc+2*nparts ) += std::real( compute( &coeffxp[1], &coeffyp[1], Et, idx_p[0], idx_p[1] )* exp_mm_theta ) ; - *( BLoc+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[1], Bl, idx_p[0], idx_d[1] )* exp_mm_theta ) ; - *( BLoc+1*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[1], Br, idx_d[0], idx_p[1] )* exp_mm_theta ) ; - *( BLoc+2*nparts ) += std::real( compute( &coeffxd[1], &coeffyd[1], Bt, idx_d[0], idx_d[1] )* exp_mm_theta ) ; + *( ELoc+0*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[0], El, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( ELoc+1*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[0], Er, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( ELoc+2*nparts ) += std::real( compute( &coeffxp[1], &coeffyp[0], Et, idx_p[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+0*nparts ) += std::real( compute( &coeffxp[1], &coeffyd[0], Bl, idx_p[0], idx_d[1] )* exp_mm_theta ) ; + *( BLoc+1*nparts ) += std::real( compute( &coeffxd[1], &coeffyp[0], Br, idx_d[0], idx_p[1] )* exp_mm_theta ) ; + *( BLoc+2*nparts ) += std::real( compute( &coeffxd[1], &coeffyd[0], Bt, idx_d[0], idx_d[1] )* exp_mm_theta ) ; } //Translate field into the cartesian y,z coordinates From de50a84aa4b6d4dc5529202892c83b42486f0e31 Mon Sep 17 00:00:00 2001 From: Francesco Massimo Date: Fri, 19 Jul 2024 15:23:28 +0200 Subject: [PATCH 17/26] add warnings for the envelope --- src/ElectroMagn/LaserEnvelope.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/ElectroMagn/LaserEnvelope.cpp b/src/ElectroMagn/LaserEnvelope.cpp index 33c00611e..09f949978 100755 --- a/src/ElectroMagn/LaserEnvelope.cpp +++ b/src/ElectroMagn/LaserEnvelope.cpp @@ -48,7 +48,14 @@ LaserEnvelope::LaserEnvelope( Params ¶ms, Patch *patch ) : if ( (envelope_solver != "explicit") && (envelope_solver != "explicit_reduced_dispersion") ){ ERROR("Unknown envelope_solver - only 'explicit' and 'explicit_reduced_dispersion' are available. "); } - + WARNING("CFL: no general CFL condition is available for the envelope solvers. The maximum stable timestep may be lower than the CFL limit of the electromagnetic solver."); + if (envelope_solver == "explicit"){ + WARNING("Envelope solver: more accurate results (albeit at a slightly higher computational cost), expecially at longer distances, require the 'explicit_reduced_dispersion' envelope solver.") + } + if (envelope_solver == "explicit_reduced_dispersion"){ + WARNING("CFL: The maximum stable timestep of the 'explicit_reduced_dispersion' envelope solver may be lower than the one of the 'explicit' solver.") + } + PyTools::extract( "polarization_phi", polarization_phi, "LaserEnvelope" ); PyTools::extract( "ellipticity", ellipticity, "LaserEnvelope" ); if ((ellipticity != 0.) and (ellipticity != 1.)){ From af4ee6d6090657da9f55df4db9879ad43bcf4fcb Mon Sep 17 00:00:00 2001 From: Francesco Massimo Date: Fri, 19 Jul 2024 15:41:17 +0200 Subject: [PATCH 18/26] message to say that the B-TIS3 is activated --- src/Params/Params.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Params/Params.cpp b/src/Params/Params.cpp index 27e743fd1..822c4eb93 100755 --- a/src/Params/Params.cpp +++ b/src/Params/Params.cpp @@ -1361,6 +1361,9 @@ void Params::print_init() TITLE( "Geometry: " << geometry ); MESSAGE( 1, "Interpolation order : " << interpolation_order ); + if (use_BTIS3){ + MESSAGE(1, "B-TIS3 interpolation scheme activated") + } MESSAGE( 1, "Maxwell solver : " << maxwell_sol ); MESSAGE( 1, "simulation duration = " << simulation_time <<", total number of iterations = " << n_time); MESSAGE( 1, "timestep = " << timestep << " = " << timestep/dtCFL << " x CFL, time resolution = " << res_time); From eb16a048a90e0f48845327469d67ca2f5a39c100 Mon Sep 17 00:00:00 2001 From: "charles.prouveur" Date: Sun, 21 Jul 2024 03:00:12 +0200 Subject: [PATCH 19/26] Updated the doc on GPU support --- doc/Sphinx/Understand/GPU_offloading.rst | 8 +++++--- doc/Sphinx/Use/GPU_version.rst | 23 +++++++++++++++++++++++ doc/Sphinx/Use/installation.rst | 2 -- doc/Sphinx/use.rst | 1 + 4 files changed, 29 insertions(+), 5 deletions(-) create mode 100755 doc/Sphinx/Use/GPU_version.rst diff --git a/doc/Sphinx/Understand/GPU_offloading.rst b/doc/Sphinx/Understand/GPU_offloading.rst index 20fb5ac07..6b2ea599a 100644 --- a/doc/Sphinx/Understand/GPU_offloading.rst +++ b/doc/Sphinx/Understand/GPU_offloading.rst @@ -17,9 +17,11 @@ the announced exaflopic supercomputers will include GPUs. * Currently supported features: * Both AMD's GPUs and Nvidia's GPUs are supported - * Cartesian geometry in 2D and in 3D + * Cartesian geometry in 1D, 2D and in 3D , for order 2 * Diagnostics: Field, Probes, Scalar, ParticleBinning, TrackParticles - * Moving Window. + * Moving Window * A few key features remain to be implemented (AM geometry, ionization, PML, envelope, - additional physics), but the fundamentals of the code are ported. \ No newline at end of file + additional physics), but the fundamentals of the code are ported. + +* Short term roadmap We are currently working on porting on GPU the following features: AM geometry (order 2) and collisions diff --git a/doc/Sphinx/Use/GPU_version.rst b/doc/Sphinx/Use/GPU_version.rst new file mode 100755 index 000000000..2228c0ea1 --- /dev/null +++ b/doc/Sphinx/Use/GPU_version.rst @@ -0,0 +1,23 @@ +GPU version of SMILEI +------- + + +This page contains the links of this documentation to compile and run SMILEI on GPU clusters, as well as the list of the features that are supported or in development + +* :doc:`Introduction` #Improved performance using GPU offloading + +* :doc:`GPU offloading , supported features, guidelines and roadmap` + +* :doc:`Intallation and compilation` + +* :doc:`Compilation for GPU-accelerated nodes` #compilation-for-gpu-accelerated-nodes + +* :doc:`Running on GPU-equiped nodes` #running-on-gpu-equiped-nodes + +---- + +Known issues +^^^^^^^^^^^^ + +2D and 3D runs may crash with A2000 & A6000 GPUs (used in laptops and worstations respectively, +they are not 'production GPUs' which are designed for 64 bits floating point operations ) diff --git a/doc/Sphinx/Use/installation.rst b/doc/Sphinx/Use/installation.rst index 791e96dc1..6cea1737a 100755 --- a/doc/Sphinx/Use/installation.rst +++ b/doc/Sphinx/Use/installation.rst @@ -185,8 +185,6 @@ Typically ``CXXFLAGS += -ta=tesla:cc80`` for ``nvhpc`` <23.4 and .. warning:: - * We are aware of issues with CUDA >12.0, fixes are being tested but are not deployed yet. - We recommend CUDA 11.x at the moment. * The hdf5 module should be compiled with the nvidia/cray compiler; openmpi as well, but depending on the nvhpc module it might not be needed as it can be included in the nvhpc module. diff --git a/doc/Sphinx/use.rst b/doc/Sphinx/use.rst index ae2c1c9d6..fb2ec9bad 100755 --- a/doc/Sphinx/use.rst +++ b/doc/Sphinx/use.rst @@ -12,3 +12,4 @@ Use Use/post-processing Use/contribute Use/troubleshoot + Use/GPU_version From 387b3b4be54ba76b595af2f9835759d877f002dc Mon Sep 17 00:00:00 2001 From: "charles.prouveur" Date: Mon, 22 Jul 2024 14:57:45 +0200 Subject: [PATCH 20/26] Removed useless declaration in 1D GPU projector that triggered unused parameters warning --- src/Projector/Projector1D2OrderGPU.cpp | 4 +- src/Projector/Projector1D2OrderGPU.h | 70 -------------------------- 2 files changed, 2 insertions(+), 72 deletions(-) diff --git a/src/Projector/Projector1D2OrderGPU.cpp b/src/Projector/Projector1D2OrderGPU.cpp index 19493ef8d..501d7c4ce 100755 --- a/src/Projector/Projector1D2OrderGPU.cpp +++ b/src/Projector/Projector1D2OrderGPU.cpp @@ -193,8 +193,8 @@ void Projector1D2OrderGPU::currentsAndDensityWrapper( ElectroMagn *EMfields, bool diag_flag, bool is_spectral, int ispec, - int icell, - int ipart_ref ) + int icell, + int ipart_ref ) // icell and ipart_ref unused { std::vector &iold = smpi->dynamics_iold[ithread]; std::vector &delta = smpi->dynamics_deltaold[ithread]; diff --git a/src/Projector/Projector1D2OrderGPU.h b/src/Projector/Projector1D2OrderGPU.h index f35e8e4ee..5bab5d11f 100755 --- a/src/Projector/Projector1D2OrderGPU.h +++ b/src/Projector/Projector1D2OrderGPU.h @@ -47,76 +47,6 @@ class Projector1D2OrderGPU : public Projector1D LocalFields Jion ) override; - //!Wrapper for task-based implementation of Smilei - //! compiler complains otherwise even if it is completely useless - void currentsAndDensityWrapperOnBuffers( double *b_Jx, - double *b_Jy, - double *b_Jz, - double *b_rho, - int bin_width, - Particles &particles, - SmileiMPI *smpi, - int istart, - int iend, - int ithread, - bool diag_flag, - bool is_spectral, - int ispec, - int icell = 0, - int ipart_ref = 0 ) override {}; -/*#if defined( SMILEI_ACCELERATOR_GPU ) - -extern "C" void -currentDepositionKernel1DOnDevice( double *__restrict__ Jx, - double *__restrict__ Jy, - double *__restrict__ Jz, - int Jx_size, - int Jy_size, - int Jz_size, - const double *__restrict__ particle_position_x, - const double *__restrict__ particle_momentum_y, - const double *__restrict__ particle_momentum_z, - const short *__restrict__ particle_charge, - const double *__restrict__ particle_weight, - const int *__restrict__ host_bin_index, - unsigned int x_dimension_bin_count, - const double *__restrict__ invgf_, - const int *__restrict__ iold_, - const double *__restrict__ deltaold_, - double inv_cell_volume, - double dx_inv, - double dx_ov_dt, - int i_domain_begin, - int not_spectral_ ); - -extern "C" void -currentAndDensityDepositionKernel1DOnDevice( double *__restrict__ Jx, - double *__restrict__ Jy, - double *__restrict__ Jz, - double *__restrict__ rho, - int Jx_size, - int Jy_size, - int Jz_size, - int rho_size, - const double *__restrict__ particle_position_x, - const double *__restrict__ particle_momentum_y, - const double *__restrict__ particle_momentum_z, - const short *__restrict__ particle_charge, - const double *__restrict__ particle_weight, - const int *__restrict__ host_bin_index, - unsigned int x_dimension_bin_count, - const double *__restrict__ invgf_, - const int *__restrict__ iold_, - const double *__restrict__ deltaold_, - double inv_cell_volume, - double dx_inv, - double dx_ov_dt, - int i_domain_begin, - int not_spectral_ ); - -#endif*/ - - protected: double dts2_; double dts4_; From 9d445cd90b959a3de2c0bb7a56fa95e522f34ae4 Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Thu, 25 Jul 2024 11:23:53 +0200 Subject: [PATCH 21/26] fix weight units for openPMD --- src/Diagnostic/DiagnosticNewParticles.cpp | 4 ++-- src/Diagnostic/DiagnosticParticleList.cpp | 2 +- src/Params/OpenPMDparams.cpp | 4 ++++ src/Params/OpenPMDparams.h | 1 + 4 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/Diagnostic/DiagnosticNewParticles.cpp b/src/Diagnostic/DiagnosticNewParticles.cpp index 888873f22..90417eb20 100755 --- a/src/Diagnostic/DiagnosticNewParticles.cpp +++ b/src/Diagnostic/DiagnosticNewParticles.cpp @@ -93,8 +93,8 @@ void DiagnosticNewParticles::openFile( Params ¶ms, SmileiMPI *smpi ) } } if( write_weight_ ) { - loc_weight_ = newDataset( species_group, "weight", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_DENSITY ); - openPMD_->writeRecordAttributes( *loc_weight_, SMILEI_UNIT_DENSITY ); + loc_weight_ = newDataset( species_group, "weight", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_WEIGHT ); + openPMD_->writeRecordAttributes( *loc_weight_, SMILEI_UNIT_WEIGHT ); } if( write_chi_ ) { loc_chi_ = newDataset( species_group, "chi", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_NONE ); diff --git a/src/Diagnostic/DiagnosticParticleList.cpp b/src/Diagnostic/DiagnosticParticleList.cpp index bdbf45ff2..21bde8161 100755 --- a/src/Diagnostic/DiagnosticParticleList.cpp +++ b/src/Diagnostic/DiagnosticParticleList.cpp @@ -299,7 +299,7 @@ void DiagnosticParticleList::run( SmileiMPI *smpi, VectorPatch &vecPatches, int if( write_weight_ ) { fill_buffer( vecPatches, iprop, data_double ); #pragma omp master - write_scalar_double( loc_weight_, "weight", data_double[0], file_space, mem_space, SMILEI_UNIT_DENSITY ); + write_scalar_double( loc_weight_, "weight", data_double[0], file_space, mem_space, SMILEI_UNIT_WEIGHT ); } iprop++; diff --git a/src/Params/OpenPMDparams.cpp b/src/Params/OpenPMDparams.cpp index 001836e9e..d00abf457 100755 --- a/src/Params/OpenPMDparams.cpp +++ b/src/Params/OpenPMDparams.cpp @@ -74,6 +74,10 @@ OpenPMDparams::OpenPMDparams( Params &p ): } else if( unit_type == SMILEI_UNIT_ENERGY ) { unitDimension[unit_type] = { 2., 1., -2., 0., 0., 0., 0. }; unitSI[unit_type] = 8.187104382e-14; // me * c^2 + } else if( unit_type == SMILEI_UNIT_WEIGHT ) { + const auto &n = params->nDim_particle; + unitDimension[unit_type] = { -3. + n, 0., 1., 1., 0., 0., 0. }; + unitSI[unit_type] = 5.034163 * Wr*Wr * pow( 299792458. / Wr, n ); // Nr * Lr^n } } diff --git a/src/Params/OpenPMDparams.h b/src/Params/OpenPMDparams.h index 3a3216ea1..2a9500f44 100755 --- a/src/Params/OpenPMDparams.h +++ b/src/Params/OpenPMDparams.h @@ -15,6 +15,7 @@ #define SMILEI_UNIT_CHARGE 7 #define SMILEI_UNIT_TIME 8 #define SMILEI_UNIT_ENERGY 9 +#define SMILEI_UNIT_WEIGHT 10 class OpenPMDparams { From 27f316f0fd3800fc33592605eb0d966e1b993add Mon Sep 17 00:00:00 2001 From: Frederic Perez Date: Thu, 25 Jul 2024 12:05:35 +0200 Subject: [PATCH 22/26] typo in previous commit --- src/Params/OpenPMDparams.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Params/OpenPMDparams.h b/src/Params/OpenPMDparams.h index 2a9500f44..7f8b12720 100755 --- a/src/Params/OpenPMDparams.h +++ b/src/Params/OpenPMDparams.h @@ -4,7 +4,7 @@ #include "Params.h" #include "H5.h" -#define SMILEI_NUNITS 10 +#define SMILEI_NUNITS 11 #define SMILEI_UNIT_NONE 0 #define SMILEI_UNIT_EFIELD 1 #define SMILEI_UNIT_BFIELD 2 From c95a9be816abbf3cd74be564deca1b5c12774237 Mon Sep 17 00:00:00 2001 From: "charles.prouveur" Date: Fri, 26 Jul 2024 19:11:58 +0200 Subject: [PATCH 23/26] fixed necessary override in projector 1D order 2 for GPU --- src/Projector/Projector1D2OrderGPU.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Projector/Projector1D2OrderGPU.h b/src/Projector/Projector1D2OrderGPU.h index 5bab5d11f..fad2df696 100755 --- a/src/Projector/Projector1D2OrderGPU.h +++ b/src/Projector/Projector1D2OrderGPU.h @@ -29,6 +29,8 @@ class Projector1D2OrderGPU : public Projector1D int icell = 0, int ipart_ref = 0 ) override; + void currentsAndDensityWrapperOnBuffers( double *, double *, double *, double *, int , Particles &, SmileiMPI *, int, int, int, bool, bool, int, int = 0, int = 0 ) override final {}; + void susceptibility( ElectroMagn *EMfields, Particles &particles, double species_mass, @@ -54,4 +56,4 @@ class Projector1D2OrderGPU : public Projector1D unsigned int x_dimension_bin_count_; }; -#endif \ No newline at end of file +#endif From f2b7bf1114060fa3f41763bf118b47b5c57cc78a Mon Sep 17 00:00:00 2001 From: "charles.prouveur" Date: Thu, 8 Aug 2024 11:45:38 +0200 Subject: [PATCH 24/26] fix for 1D GPU with no diag every iteration --- src/Projector/Projector1D2OrderGPU.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Projector/Projector1D2OrderGPU.cpp b/src/Projector/Projector1D2OrderGPU.cpp index 501d7c4ce..c96add685 100755 --- a/src/Projector/Projector1D2OrderGPU.cpp +++ b/src/Projector/Projector1D2OrderGPU.cpp @@ -246,7 +246,7 @@ void Projector1D2OrderGPU::currentsAndDensityWrapper( ElectroMagn *EMfields, else{ #if defined( SMILEI_ACCELERATOR_GPU ) - currentDepositionKernel1DOnDevice(Jx_, Jy_, Jz_, + currentDepositionKernel1DOnDevice(EMfields->Jx_->data(), EMfields->Jz_->data(), EMfields->Jz_->data(), EMfields->Jx_->size(), EMfields->Jy_->size(), EMfields->Jz_->size(), particles.getPtrPosition( 0 ), particles.getPtrMomentum( 1 ), From 3e775f01b896fec253151298f924fdedc930e12e Mon Sep 17 00:00:00 2001 From: "charles.prouveur" Date: Thu, 8 Aug 2024 12:02:37 +0200 Subject: [PATCH 25/26] removed useless class member pointers Jx_, etc. in 1D, could be done in 2D and 3D as well --- src/Projector/Projector1D.h | 2 +- src/Projector/Projector1D2Order.cpp | 12 ++++++------ src/Projector/Projector1D4Order.cpp | 12 ++++++------ 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Projector/Projector1D.h b/src/Projector/Projector1D.h index c08c0e9a8..67eb72382 100755 --- a/src/Projector/Projector1D.h +++ b/src/Projector/Projector1D.h @@ -30,7 +30,7 @@ class Projector1D : public Projector double dx_inv_; double dx_ov_dt_; int i_domain_begin_; - double *Jx_, *Jy_, *Jz_, *rho_; + //double *Jx_, *Jy_, *Jz_, *rho_; }; #endif diff --git a/src/Projector/Projector1D2Order.cpp b/src/Projector/Projector1D2Order.cpp index 451bca539..67d050ead 100755 --- a/src/Projector/Projector1D2Order.cpp +++ b/src/Projector/Projector1D2Order.cpp @@ -315,20 +315,20 @@ void Projector1D2Order::currentsAndDensityWrapper( ElectroMagn *EMfields, Partic std::vector *delta = &( smpi->dynamics_deltaold[ithread] ); std::vector *invgf = &( smpi->dynamics_invgf[ithread] ); - Jx_ = &( *EMfields->Jx_ )( 0 ); - Jy_ = &( *EMfields->Jy_ )( 0 ); - Jz_ = &( *EMfields->Jz_ )( 0 ); - rho_ = &( *EMfields->rho_ )( 0 ); + double *Jx = &( *EMfields->Jx_ )( 0 ); + double *Jy = &( *EMfields->Jy_ )( 0 ); + double *Jz = &( *EMfields->Jz_ )( 0 ); + double *rho = &( *EMfields->rho_ )( 0 ); // If no field diagnostics this timestep, then the projection is done directly on the total arrays if( !diag_flag ) { if( !is_spectral ) { for( int ipart=istart ; ipart *delta = &( smpi->dynamics_deltaold[ithread] ); std::vector *invgf = &( smpi->dynamics_invgf[ithread] ); - Jx_ = &( *EMfields->Jx_ )( 0 ); - Jy_ = &( *EMfields->Jy_ )( 0 ); - Jz_ = &( *EMfields->Jz_ )( 0 ); - rho_ = &( *EMfields->rho_ )( 0 ); + double *Jx = &( *EMfields->Jx_ )( 0 ); + double *Jy = &( *EMfields->Jy_ )( 0 ); + double *Jz = &( *EMfields->Jz_ )( 0 ); + double *rho = &( *EMfields->rho_ )( 0 ); // If no field diagnostics this timestep, then the projection is done directly on the total arrays if( !diag_flag ) { if( !is_spectral ) { for( int ipart=istart ; ipart Date: Thu, 15 Aug 2024 14:43:12 +0200 Subject: [PATCH 26/26] update releases.rst --- doc/Sphinx/Overview/releases.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index b93740024..9e2b1445b 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -35,6 +35,7 @@ Changes made in the repository (not released) * Custom functions in ``ParticleBinning`` crashed with python 3.12. * Species-specific diagnostics in AM geometry with vectorization. * Happi's ``average`` argument would sometimes be missing the last bin. + * 1D projector on GPU without diagnostics ----