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/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index e271b32c5..9e2b1445b 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -16,13 +16,27 @@ Get Smilei You can find older, `unsupported versions here `_ -.. -.. ---- +---- + +.. _latestVersion: + +Changes made in the repository (not released) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +* **Features**: + + * Prescribed fields in AM geometry. + * Particle reflective boundary conditions at Rmax in AM geometry. + * 1st order Ruyten shape function in AM geometry. + +* **Bug fixes**: -.. .. _latestVersion: + * Tunnel ionization was wrong in some cases for high atomic numbers. + * 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 -.. Changes made in the repository (not released) -.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ---- 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/namelist.rst b/doc/Sphinx/Use/namelist.rst index a07f19005..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"`` @@ -2009,8 +2013,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 +2059,25 @@ This feature is accessible using the ``PrescribedField`` block:: .. py:data:: field - Field name: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_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. + +.. 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 ` @@ -3641,4 +3657,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/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 diff --git a/happi/_Diagnostics/Diagnostic.py b/happi/_Diagnostics/Diagnostic.py index 39604f0ec..0a79de56d 100755 --- a/happi/_Diagnostics/Diagnostic.py +++ b/happi/_Diagnostics/Diagnostic.py @@ -578,8 +578,8 @@ def _selectRange(self, portion, meshpoints, axisname, axisunits, operation, edge info = operation+" for "+axisname+" from "+axismin+" to "+axismax else: info = operation+" for "+axisname+" from "+str(meshpoints[indices[0]])+" to "+str(meshpoints[indices[-1]])+" "+axisunits - selection = slice(indices[0],indices[-1]) - finalShape = indices[-1] - indices[0] + selection = slice(indices[0], indices[-1] + 1) + finalShape = indices[-1] - indices[0] + 1 return info, selection, finalShape # Method to prepare some data before plotting 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/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/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 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.)){ diff --git a/src/Interpolator/InterpolatorAM1OrderRuyten.cpp b/src/Interpolator/InterpolatorAM1OrderRuyten.cpp new file mode 100755 index 000000000..708ee5ac8 --- /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[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 + 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/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; 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..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 @@ -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 { diff --git a/src/Params/Params.cpp b/src/Params/Params.cpp index 69973d104..822c4eb93 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 @@ -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); 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/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 ); diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 42f4dd3d8..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. // --------------------------------------------------------------------------------------------------------------------- @@ -1443,10 +1442,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 ); @@ -4586,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/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 &iold = smpi->dynamics_iold[ithread]; std::vector &delta = smpi->dynamics_deltaold[ithread]; @@ -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 ), diff --git a/src/Projector/Projector1D2OrderGPU.h b/src/Projector/Projector1D2OrderGPU.h index f35e8e4ee..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, @@ -47,76 +49,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_; @@ -124,4 +56,4 @@ currentAndDensityDepositionKernel1DOnDevice( double *__restrict__ Jx, unsigned int x_dimension_bin_count_; }; -#endif \ No newline at end of file +#endif diff --git a/src/Projector/Projector1D4Order.cpp b/src/Projector/Projector1D4Order.cpp index ea4eafa4a..a10664c24 100755 --- a/src/Projector/Projector1D4Order.cpp +++ b/src/Projector/Projector1D4Order.cpp @@ -361,20 +361,20 @@ void Projector1D4Order::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 +#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.; + } + for( unsigned int i=0; i<4; i++ ) { + 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/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; 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/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 ); 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") diff --git a/src/Tools/PyTools.h b/src/Tools/PyTools.h index 055875718..a756f3db8 100755 --- a/src/Tools/PyTools.h +++ b/src/Tools/PyTools.h @@ -32,6 +32,29 @@ #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 +#if PY_MAJOR_VERSION > 2 && PY_MINOR_VERSION > 11 + #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 + #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 class PyTools 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 000000000..4b2dfccf0 Binary files /dev/null and b/validation/references/tstAM_04d_laser_propagation_Order1.py.txt differ