diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5a509eaac..c4262426f 100755 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -16,6 +16,7 @@ install: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Force workdir cleaning in case of retried @@ -35,6 +36,7 @@ compile_default: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -47,6 +49,7 @@ runQuick: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -59,6 +62,7 @@ run1D: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -72,6 +76,7 @@ run2D: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -87,6 +92,7 @@ run3D: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -103,6 +109,7 @@ runAM: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -115,6 +122,7 @@ runCollisions: only: - develop - BTIS3_for_develop + - diag_new_particles script: # Move in test dir @@ -127,6 +135,7 @@ compile_picsar: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei @@ -138,6 +147,7 @@ run_picsar: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation @@ -148,6 +158,7 @@ compile_debug: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei @@ -159,6 +170,7 @@ compile_no_mpi_threadmultiple: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei @@ -170,6 +182,7 @@ compile_no_openmp: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei @@ -181,6 +194,7 @@ compile_omptasks: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei @@ -192,6 +206,7 @@ run_omptasks: only: - develop - BTIS3_for_develop + - diag_new_particles script: - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation diff --git a/benchmarks/tst1d_05_tunnel_ionisation.py b/benchmarks/tst1d_05_tunnel_ionisation.py index 7ce923b46..c0bc4e711 100755 --- a/benchmarks/tst1d_05_tunnel_ionisation.py +++ b/benchmarks/tst1d_05_tunnel_ionisation.py @@ -116,4 +116,10 @@ def Bz(t): species = "electron", every = [1,1000,30], attributes = ["x","px","py","pz","w","Wy"] -) \ No newline at end of file +) + +DiagNewParticles( + species = "electron", + every = 100, + attributes = ["x","py","w","q"], +) diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 8f32c6ec0..9f12d7ea7 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -23,30 +23,34 @@ You can find older, `unsupported versions here 0``. + It is possible to make logical operations: ``+`` is *OR*; ``*`` is *AND*; ``~`` is *NOT*. + + | **Example:** ``select="(x>1)*(x<2)"`` + + It is also possible to select directly a list of IDs. + + | **Example:** ``select=[ID1, ID2, ...]`` ---- diff --git a/doc/Sphinx/smilei_theme/static/smilei_theme.css_t b/doc/Sphinx/smilei_theme/static/smilei_theme.css_t index 0668da71b..1a8d7bf78 100755 --- a/doc/Sphinx/smilei_theme/static/smilei_theme.css_t +++ b/doc/Sphinx/smilei_theme/static/smilei_theme.css_t @@ -709,8 +709,8 @@ div.header > * { font-size:1.3em; } -div.experimental > h1:first-of-type::after, -div.experimental > h2:first-of-type::after, +.experimental > h1:first-of-type::after, +.experimental > h2:first-of-type::after, dl.experimental > dt:first-of-type::after, ul.experimental p:first-of-type::after { diff --git a/happi/_Diagnostics/NewParticles.py b/happi/_Diagnostics/NewParticles.py new file mode 100755 index 000000000..b258c8594 --- /dev/null +++ b/happi/_Diagnostics/NewParticles.py @@ -0,0 +1,293 @@ +from .Diagnostic import Diagnostic +from .ParticleList import ParticleList +from .._Utils import * + +class NewParticles(ParticleList): + """Class for loading a NewParticles diagnostic""" + + _diagType = "NewParticles" + + def _init(self, species=None, select="", axes=[], chunksize=20000000, **kwargs): + + # If argument 'species' not provided, then print available species and leave + if species is None: + species = self.simulation.getNewParticlesSpecies() + error = ["Argument `species` not provided"] + if len(species)>0: + error += ["Printing available species for NewParticles:"] + error += ["--------------------------------------------"] + error += ["\n".join(species)] + else: + error += ["No files found for NewParticles"] + raise Exception("\n".join(error)) + + # Get info from the hdf5 files + verifications + # ------------------------------------------------------------------- + self.species = species + files = self._findFiles("NewParticles") + self._h5files = [self._h5py.File(file, "r") for file in files] + + # List available properties + self._short_properties_from_raw["birth_time"] = "t" + self._raw_properties_from_short = {v:k for k,v in self._short_properties_from_raw.items()} + T0 = self._h5files[0]["data/0/particles/"+self.species] + self.available_properties = {v for k,v in self._short_properties_from_raw.items() if k in T0} + + # Find the overlap between files + if len(self._h5files) == 1: + # If only 1 file, there is 1 range containing the whole file + self.nParticles = self._h5files[0]["iteration_npart"][-1,1] + self.ranges_in_file = {0:[self._np.s_[0:self.nParticles]]} + else: + # If several files, this is more complicated + iterations_in_file = [] + npart_in_file = [] + all_iterations = self._np.array([], dtype=int) + # for each file, get the list of iterations dumped and the corresponding number of particles for each iteration + for i,f in enumerate(self._h5files): + iterations_in_file += [f["iteration_npart"][:,0]] + npart_in_file += [f["iteration_npart"][:,1]] + # all_iteration is the union of all iterations in all files + all_iterations = self._np.union1d( all_iterations, iterations_in_file[-1] ) + # Now, assign one file to each existing iteration. If 2 files contain the same iteration, the last file prevails + file_vs_iteration = self._np.empty_like( all_iterations, dtype = int ) + for i,f in enumerate(self._h5files): + locations = self._np.isin( all_iterations, iterations_in_file[i], assume_unique=True ) + file_vs_iteration[locations] = i + # find for which iterations there is a change in the file number + change_locs = 1 + self._np.flatnonzero(self._np.diff(file_vs_iteration)) + change_locs = self._np.concatenate(([0], change_locs, [len(all_iterations)])) + # Now the goal is to calculate several ranges of particles that each file should handle + # Instead of looping every iteration, we loop on the iterations where a file change occurs + self.nParticles = 0 + self.ranges_in_file = {i:[] for i in range(len(self._h5files))} + for i in range(len(change_locs)-1): + # Get the locations of relevant iterations + istart, istop = change_locs[i], change_locs[i+1] + # Get the corresponding file + ifile = file_vs_iteration[istart] + # Get the relevant iterations + iterstart = all_iterations[istart] + iterstop = all_iterations[-1]+1 if istop >= len(all_iterations) else all_iterations[istop] + # Get the iteration locations in the file + iterindexstart, iterindexstop = self._np.searchsorted(iterations_in_file[ifile], [iterstart, iterstop]) + # Get the corresponding particle indices in the file + ipartstart = 0 if iterindexstart == 0 else npart_in_file[ifile][iterindexstart-1] + ipartstop = npart_in_file[ifile][iterindexstop -1] + # Add this range of particles + self.nParticles += ipartstop - ipartstart + self.ranges_in_file[ifile] += [self._np.s_[ipartstart:ipartstop]] + + # Manage axes + self._initAxes(axes) + + # Select particles + self.selectedParticles = None + self.nselectedParticles = self.nParticles + if type(select) is str: + # Transform the operation + try: + particleSelector = select + doubleProps = [] + int16Props = [] + for prop in self.available_properties: + particleSelector, nsubs = self._re.subn(r"\b"+prop+r"\b", "particle_data['"+prop+"']", particleSelector) + if nsubs > 0: + if prop == "q" : int16Props += [prop] + else : doubleProps += [prop] + except: + raise Exception("Error: `select` not understood") + # Nothing to select if empty `select` does not contain any particle property + if len(int16Props) + len(doubleProps) > 0.: + # Get particle selection chunk by chunk + self.selectedParticles = [] + for particle_data in self._iterAllParticles(chunksize): + self.selectedParticles += [self._np.flatnonzero( eval(particleSelector) )] + # Concatenate all chunks + self.selectedParticles = self._np.concatenate( self.selectedParticles ) + self.nselectedParticles = len(self.selectedParticles) + + # Otherwise, the selection can be a list of particle IDs + elif select: + # Loop only IDs + self.selectedParticles = [] + for particle_data in self._iterAllParticles(chunksize, axes=["Id"]): + self.selectedParticles += [self._np.flatnonzero(self._np.in1d(particle_data["Id"], select))] + self.selectedParticles = self._np.sort( self._np.concatenate( self.selectedParticles ) ) + self.nselectedParticles = len(self.selectedParticles) + + # Now split the selected particles per the different files where they belong + if self.nselectedParticles < self.nParticles: + sel = {} + offset = 0 + for ifile, ranges in self.ranges_in_file.items(): + sel[ifile] = [] + for r in ranges: + istart, istop = self._np.searchsorted( self.selectedParticles, [offset + r.start, offset + r.stop] ) + sel[ifile] += [self.selectedParticles[istart:istop]-offset] + offset += r.stop - r.start + sel[ifile] = self._np.concatenate( sel[ifile] ) + self.selectedParticles = sel + + if select and self._verbose: + print("Kept "+str(self.nselectedParticles)+" particles") + + # Finish constructor + self._timesteps = [None] + self.valid = True + return kwargs + + # We override the get and getData methods + def getData(self): + if not self._validate(): return + for data in self.iterParticles(chunksize=self.nParticles): + return data + + # Iterator on particles + def iterParticles(self, chunksize=1): + if self.nselectedParticles < self.nParticles: + return self._iterSelectedParticles(chunksize, self.axes) + else: + return self._iterAllParticles(chunksize, self.axes) + + # Iterator for all particles + def _iterAllParticles(self, chunksize=1, axes=None): + if axes is None: + axes = self.axes + # Make buffers + properties = self._raw_properties_from_short + data = {} + for axis in axes: + if axis == "Id": + data[axis] = self._np.empty((chunksize,), dtype=self._np.uint64) + elif axis == "q": + data[axis] = self._np.empty((chunksize,), dtype=self._np.int16 ) + else: + data[axis] = self._np.empty((chunksize,), dtype=self._np.double) + + # Loop files + npart_in_chunk = 0 + for ifile,f in enumerate(self._h5files): + # Get the main group containing all datasets + try: + group = f["data/0/particles/"+self.species] + except: + continue + # Loop on the ranges of particles + for r in self.ranges_in_file[ifile]: + nremaining_in_chunk = chunksize - npart_in_chunk + start = r.start + # Yield chunks until the current range of particles is complete + while r.stop - start > nremaining_in_chunk: + for axis in axes: + group[properties[axis]].read_direct(data[axis], source_sel=self._np.s_[start:start+nremaining_in_chunk], dest_sel=self._np.s_[npart_in_chunk:]) + yield data + # Start new chunk + start += nremaining_in_chunk + npart_in_chunk = 0 + nremaining_in_chunk = chunksize + # If the range still has a few particles, start filling the next chunk + if start < r.stop: + for axis in axes: + group[properties[axis]].read_direct(data[axis], source_sel=self._np.s_[start:], dest_sel=self._np.s_[npart_in_chunk:npart_in_chunk+r.stop-start]) + npart_in_chunk += r.stop-start + # If the last chunk was not complete, yield that + if npart_in_chunk > 0: + # Resize buffers + for axis in axes: + data[axis].resize((npart_in_chunk,)) + yield data + + # Iterator for selected particles + def _iterSelectedParticles(self, chunksize=1, axes=None): + # Make buffers + properties = self._raw_properties_from_short + data = {} + for axis in axes: + if axis == "Id": + data[axis] = self._np.empty((chunksize,), dtype=self._np.uint64) + elif axis == "q": + data[axis] = self._np.empty((chunksize,), dtype=self._np.int16 ) + else: + data[axis] = self._np.empty((chunksize,), dtype=self._np.double) + + # Loop files + npart_in_chunk = 0 + for ifile,f in enumerate(self._h5files): + # Get the main group containing all datasets + try: + group = f["data/0/particles/"+self.species] + except: + continue + #selection for this file + sel = self.selectedParticles[ifile] + # Yield chunks until the current selection of particles is complete + nremaining_in_chunk = chunksize - npart_in_chunk + start = 0; stop = sel.size + while stop - start > nremaining_in_chunk: + print("chunk %d"%chunkstart) + for axis in axes: + group[properties[axis]].read_direct(data[axis], source_sel=sel[start:start+nremaining_in_chunk], dest_sel=self._np.s_[npart_in_chunk:]) + yield data + # Start new chunk + start += nremaining_in_chunk + npart_in_chunk = 0 + nremaining_in_chunk = chunksize + # If the selection still has a few particles, start filling the next chunk + if start < stop: + for axis in axes: + group[properties[axis]].read_direct(data[axis], source_sel=sel[start:], dest_sel=self._np.s_[npart_in_chunk:npart_in_chunk+stop-start]) + npart_in_chunk += stop-start + # If the last chunk was not complete, yield that + if npart_in_chunk > 0: + # Resize buffers + for axis in axes: + data[axis].resize((npart_in_chunk,)) + yield data + + # Method to get info + def _info(self): + info = "New particles diagnostic: species '"+self.species+"' containing "+str(self.nParticles)+" particles" + if self.selectedParticles: + info += "\n\twith selection of "+str(self.nselectedParticles)+" particles" + info += "\n\tAxes: " + ", ".join(self.axes) + return info + + # Prevent using some plotting functions + def streak(self, *args, **kwargs): + raise("Error: Cannot use streak in NewParticles") + def animate(self, *args, **kwargs): + raise("Error: Cannot use animate in NewParticles") + def slide(self, *args, **kwargs): + raise("Error: Cannot use slide in NewParticles") + def toVTK(self, *args, **kwargs): + raise("Error: Cannot use toVTK in NewParticles") + + # Override plotting functions + def _prepare3(self): + if self._tmpdata is None: + A = self.getData() + self._tmpdata = [] + for axis in self.axes: self._tmpdata.append( A[axis] ) + return True + + def _plotOnAxes_0D(self, ax, t, cax_id=0): + pass + + def _plotOnAxes_1D(self, ax, t, cax_id=0): + pass + + def _plotOnAxes_2D(self, ax, t, cax_id=0): + xfactor = self.options.xfactor or 1. + yfactor = self.options.yfactor or 1. + try : ax.set_prop_cycle (None) + except: ax.set_color_cycle(None) + ax.plot(xfactor*self._tmpdata[0], yfactor*self._tmpdata[1], '.', **self.options.plot) + # Add labels and options + ax.set_xlabel(self._xlabel) + ax.set_ylabel(self._ylabel) + self._setLimits(ax, xmin=self.options.xmin, xmax=self.options.xmax, ymin=self.options.ymin, ymax=self.options.ymax) + self._setTitle(ax) + self._setAxesOptions(ax) + return 1 + diff --git a/happi/_Diagnostics/ParticleList.py b/happi/_Diagnostics/ParticleList.py new file mode 100755 index 000000000..8f6ef6cb5 --- /dev/null +++ b/happi/_Diagnostics/ParticleList.py @@ -0,0 +1,91 @@ +from .Diagnostic import Diagnostic +from .._Utils import * + +class ParticleList(Diagnostic): + """Generic class for a diagnostic with listed particles (trackParticles or newParticles)""" + + _short_properties_from_raw = { + "id":"Id", "position/x":"x", "position/y":"y", "position/z":"z", + "momentum/x":"px", "momentum/y":"py", "momentum/z":"pz", + "charge":"q", "weight":"w", "chi":"chi", + "E/x":"Ex", "E/y":"Ey", "E/z":"Ez", "B/x":"Bx", "B/y":"By", "B/z":"Bz", + "W/x":"Wx", "W/y":"Wy", "W/z":"Wz" + } + + def _initAxes(self, axes): + + if type(axes) is not list: + raise Exception("Error: Argument 'axes' must be a list") + + # if axes provided, verify them + if len(axes)>0: + self.axes = axes + unknown_axes = set(axes) - self.available_properties + if unknown_axes: + raise Exception( + "Error: Argument 'axes' has unknown item(s): "+", ".join(unknown_axes)+".\n" + + " Available axes are: "+", ".join(sorted(self.available_properties)) + ) + # otherwise use default + else: + self.axes = list(self.available_properties) + + # Then figure out axis units + self._factors = [] + for axis in self.axes: + axisunits = "" + if axis == "Id": + self._centers.append( [0, 281474976710655] ) + elif axis in ["x" , "y" , "z", "moving_x"]: + axisunits = "L_r" + self._centers.append( [0., self.namelist.Main.grid_length[{"x":0,"y":1,"z":-1}[axis[-1]]]] ) + elif axis in ["px", "py", "pz"]: + axisunits = "P_r" + self._centers.append( [-1., 1.] ) + elif axis == "w": + axisunits = "N_r * L_r^%i" % self._ndim_particles + self._centers.append( [0., 1.] ) + elif axis == "q": + axisunits = "Q_r" + self._centers.append( [-10., 10.] ) + elif axis == "chi": + axisunits = "1" + self._centers.append( [0., 2.] ) + elif axis[0] == "E": + axisunits = "E_r" + self._centers.append( [-1., 1.] ) + elif axis[0] == "B": + axisunits = "B_r" + self._centers.append( [-1., 1.] ) + elif axis[0] == "W": + axisunits = "K_r" + self._centers.append( [0., 1.] ) + self._log += [False] + self._label += [axis] + self._units += [axisunits] + if axis == "Id": + self._factors += [1] + else: + factor, _ = self.units._convert(axisunits, None) + self._factors += [factor] + self._title = self._diagType+" '"+self.species+"'" + self._shape = [0]*len(self.axes) + self._centers = [self._np.array(c) for c in self._centers] + self._type = self.axes + + # Hack to work with 1 axis + self._vunits = self._units[0] if len(axes)==1 else "" + + def get(self): + return self.getData() + + # Get a list of files + def _findFiles(self, filePrefix): + files = [] + for path in self._results_path: + file = path+self._os.sep+filePrefix+"_"+self.species+".h5" + if self._os.path.isfile(file): + files += [file] + if not files: + raise Exception("No files "+filePrefix+" found") + return files diff --git a/happi/_Diagnostics/TrackParticles.py b/happi/_Diagnostics/TrackParticles.py index a6afaa54d..2d4e8bfb4 100755 --- a/happi/_Diagnostics/TrackParticles.py +++ b/happi/_Diagnostics/TrackParticles.py @@ -1,4 +1,5 @@ from .Diagnostic import Diagnostic +from .ParticleList import ParticleList from .._Utils import * # Define a function that finds the next closing character in a string @@ -19,11 +20,13 @@ def findClosingCharacter(string, character, start=0): raise Exception("Error in selector syntax: missing `"+character+"`") -class TrackParticles(Diagnostic): +class TrackParticles(ParticleList): """Class for loading a TrackParticles diagnostic""" - + + _diagType = "TrackParticles" + def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sorted_as="", length=None, chunksize=20000000, **kwargs): - + # If argument 'species' not provided, then print available species and leave if species is None: species = self.simulation.getTrackSpecies() @@ -47,14 +50,7 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor # ------------------------------------------------------------------- self.species = species self._h5items = {} - disorderedfiles = self._findDisorderedFiles() - self._short_properties_from_raw = { - "id":"Id", "position/x":"x", "position/y":"y", "position/z":"z", - "momentum/x":"px", "momentum/y":"py", "momentum/z":"pz", - "charge":"q", "weight":"w", "chi":"chi", - "E/x":"Ex", "E/y":"Ey", "E/z":"Ez", "B/x":"Bx", "B/y":"By", "B/z":"Bz", - "W/x":"Wx", "W/y":"Wy", "W/z":"Wz" - } + disorderedfiles = self._findFiles("TrackParticlesDisordered") # Get x_moved and add moving_x in the list of properties self._XmovedForTime = {} @@ -64,7 +60,7 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor x_moved = val.attrs.get("x_moved") if x_moved is not None: self._XmovedForTime[int(t)] = x_moved - extra_properties = ["moving_x"] if self._XmovedForTime else [] + extra_properties = {"moving_x"} if self._XmovedForTime else set() # If sorting allowed, find out if ordering needed needsOrdering = False @@ -110,7 +106,7 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor self._raw_properties_from_short = {v:k for k,v in self._short_properties_from_raw.items()} f, _ = next(iter(self._locationForTime.values())) T0 = next(iter(f["data"].values()))["particles/"+self.species] - self.available_properties = [v for k,v in self._short_properties_from_raw.items() if k in T0] + extra_properties + self.available_properties = {v for k,v in self._short_properties_from_raw.items() if k in T0} | extra_properties # If sorting allowed, then do the sorting if sort: @@ -125,7 +121,7 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor "Ex", "Ey", "Ez", "Bx", "By", "Bz", "Wx", "Wy", "Wz"]: if prop in self._lastfile: self._h5items[prop] = self._lastfile[prop] - self.available_properties = list(self._h5items.keys()) + extra_properties + self.available_properties = set(self._h5items.keys()) | extra_properties # Memorize the locations of timesteps in the files self._locationForTime = {t:it for it, t in enumerate(self._lastfile["Times"])} self._timesteps = self._np.array(sorted(self._lastfile["Times"])) @@ -169,69 +165,7 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor if self._verbose: print("Kept "+str(self.nselectedParticles)+" particles") # Manage axes - # ------------------------------------------------------------------- - if type(axes) is not list: - raise Exception("Error: Argument 'axes' must be a list") - # if axes provided, verify them - if len(axes)>0: - self.axes = axes - for axis in axes: - if axis not in self.available_properties: - raise Exception( - "Error: Argument 'axes' has item '"+str(axis)+"' unknown.\n" - + " Available axes are: "+(", ".join(sorted(self.available_properties))) - ) - # otherwise use default - else: - self.axes = self.available_properties - - - # Then figure out axis units - self._type = self.axes - self._factors = [] - for axis in self.axes: - axisunits = "" - if axis == "Id": - self._centers.append( [0, 281474976710655] ) - elif axis in ["x" , "y" , "z", "moving_x"]: - axisunits = "L_r" - self._centers.append( [0., self.namelist.Main.grid_length[{"x":0,"y":1,"z":-1}[axis[-1]]]] ) - elif axis in ["px", "py", "pz"]: - axisunits = "P_r" - self._centers.append( [-1., 1.] ) - elif axis == "w": - axisunits = "N_r * L_r^%i" % self._ndim_particles - self._centers.append( [0., 1.] ) - elif axis == "q": - axisunits = "Q_r" - self._centers.append( [-10., 10.] ) - elif axis == "chi": - axisunits = "1" - self._centers.append( [0., 2.] ) - elif axis[0] == "E": - axisunits = "E_r" - self._centers.append( [-1., 1.] ) - elif axis[0] == "B": - axisunits = "B_r" - self._centers.append( [-1., 1.] ) - elif axis[0] == "W": - axisunits = "K_r" - self._centers.append( [0., 1.] ) - self._log += [False] - self._label += [axis] - self._units += [axisunits] - if axis == "Id": - self._factors += [1] - else: - factor, _ = self.units._convert(axisunits, None) - self._factors += [factor] - self._title = "Track particles '"+species+"'" - self._shape = [0]*len(self.axes) - self._centers = [self._np.array(c) for c in self._centers] - - # Hack to work with 1 axis - if len(axes)==1: self._vunits = self._units[0] - else: self._vunits = "" + self._initAxes(axes) # Set the directory in case of exporting self._exportPrefix = "TrackParticles_"+self.species+"_"+"".join(self.axes) @@ -452,17 +386,6 @@ def _readUnstructuredH5(self, dataset, indices, first_time, last_time=None): def getAvailableTimesteps(self): return self._alltimesteps - # Get a list of disordered files - def _findDisorderedFiles(self): - disorderedfiles = [] - for path in self._results_path: - file = path+self._os.sep+"TrackParticlesDisordered_"+self.species+".h5" - if self._os.path.isfile(file): - disorderedfiles += [file] - if not disorderedfiles: - raise Exception("No TrackParticles files") - return disorderedfiles - # Make the particles ordered by Id in the file, in case they are not def _orderFiles( self, fileOrdered, chunksize, sort ): if self._verbose: @@ -704,9 +627,6 @@ def getData(self, timestep=None): data[t][axis] = self._rawData[t][axis] * factor return data - def get(self): - return self.getData() - # Iterator on UNSORTED particles for a given timestep def iterParticles(self, timestep, chunksize=1): if not self._validate(): return @@ -719,7 +639,7 @@ def iterParticles(self, timestep, chunksize=1): properties = {"moving_x":"x"} properties.update( self._raw_properties_from_short ) - disorderedfiles = self._findDisorderedFiles() + disorderedfiles = self._findFiles("TrackParticlesDisordered") for file in disorderedfiles: f = self._h5py.File(file, "r") # This is the timestep for which we want to produce an iterator diff --git a/happi/_Factories.py b/happi/_Factories.py index d3897f547..8f47325ec 100755 --- a/happi/_Factories.py +++ b/happi/_Factories.py @@ -5,6 +5,7 @@ from ._Diagnostics.Screen import Screen from ._Diagnostics.RadiationSpectrum import RadiationSpectrum from ._Diagnostics.TrackParticles import TrackParticles +from ._Diagnostics.NewParticles import NewParticles from ._Diagnostics.Performances import Performances class ScalarFactory(object): @@ -511,15 +512,14 @@ class TrackParticlesFactory(object): To get a list of available tracked species, simply omit this argument. select: string (optional) Instructions for selecting particles among those available. - Syntax 1: select="any(times, condition)" - Syntax 2: select="all(times, condition)" + Syntax 1: select="any(times, condition)" -> particles satisfying `condition` for at least one of the `times` + Syntax 2: select="all(times, condition)" -> particles satisfying `condition` at all `times` + Syntax 3: select=[ID1, ID2, ...] -> particles selected by their ID `times` is a selection of timesteps t, for instance `t>50`. `condition` is a condition on particles properties (x, y, z, px, py, pz), for instance `px>0`. - Syntax 1 selects particles satisfying `condition` for at least one of the `times`. - Syntax 2 selects particles satisfying `condition` at all `times`. Example: select="all(t<40, px<0.1)" selects particles that kept px<0.1 until timestep 40. Example: select="any(t>0, px>1.)" selects particles that reached px>1 at some point. - It is possible to make logical operations: + is OR; * is AND; - is NOT. + It is possible to make logical operations: + is OR; * is AND; ~ is NOT. Example: select="any((t>30)*(t<60), px>1) + all(t>0, (x>1)*(x<2))" timesteps : int or [int, int] (optional) If omitted, all timesteps are used. @@ -578,3 +578,66 @@ def __repr__(self): else: msg += "\nNo TrackParticles diagnostics available" return msg + + + +class NewParticlesFactory(object): + """Import and analyze new particles from a Smilei simulation + + Parameters: + ----------- + species : string (optional) + Name of a species in a NewParticles diagnostic. + To get a list of available species, simply omit this argument. + select: string (optional) + Instructions for selecting particles among those available. + It must be a condition on particles properties (x, y, z, px, py, pz), for instance `px>0`. + It is possible to make logical operations: + is OR; * is AND; ~ is NOT. + Example: select="(x>1)*(x<2)" + It is also possible to select directly a list of IDs. + Example: select=[ID1, ID2, ...] + units : A units specification such as ["m","second"] + axes: A list of required particle properties. + Each axis is "x", "y", "z", "px", "py" or "pz", "chi". + + Usage: + ------ + S = happi.Open("path/to/simulation") # Load the simulation + np = S.NewParticles(...) # Load the tracked-particle diagnostic + np.get() # Obtain the data + """ + + def __init__(self, simulation, species=None): + self._simulation = simulation + self._additionalKwargs = dict() + self._children = [] + if not simulation._scan: return + + # If not a specific species (root level), build a list of species shortcuts + if species is None: + if simulation._verbose: print("Scanning for new particle diagnostics") + # Get a list of species + specs = self._simulation.getNewParticlesSpecies() + # Create species shortcuts + for spec in specs: + child = NewParticlesFactory(simulation, spec) + setattr(self, spec, child) + self._children += [child] + + else: + # the species is saved for generating the object in __call__ + self._additionalKwargs.update( {"species":species} ) + + def __call__(self, *args, **kwargs): + kwargs.update(self._additionalKwargs) + return NewParticles(self._simulation, *args, **kwargs) + + def __repr__(self): + msg = object.__repr__(self) + if len(self._additionalKwargs) == 0 and self._simulation._scan: + specs = self._simulation.getNewParticlesSpecies() + if specs: + msg += "\nAvailable NewParticles diagnostics:\n" + ", ".join(specs) + else: + msg += "\nNo NewParticles diagnostics available" + return msg diff --git a/happi/_SmileiSimulation.py b/happi/_SmileiSimulation.py index 340978968..b8c3a3fbe 100644 --- a/happi/_SmileiSimulation.py +++ b/happi/_SmileiSimulation.py @@ -1,4 +1,4 @@ -from ._Factories import ScalarFactory, FieldFactory, ProbeFactory, ParticleBinningFactory, RadiationSpectrumFactory, PerformancesFactory, ScreenFactory, TrackParticlesFactory +from ._Factories import ScalarFactory, FieldFactory, ProbeFactory, ParticleBinningFactory, RadiationSpectrumFactory, PerformancesFactory, ScreenFactory, TrackParticlesFactory, NewParticlesFactory from ._Utils import * @@ -109,6 +109,7 @@ def __init__(self, results_path=".", reference_angular_frequency_SI=None, show=T self.Performances = PerformancesFactory(self) self.Screen = ScreenFactory(self) self.TrackParticles = TrackParticlesFactory(self) + self.NewParticles = NewParticlesFactory(self) def _openNamelist(self, path): # empty class to store the namelist variables @@ -270,16 +271,24 @@ def __repr__(self): files = "\n\t".join(files) return "Smilei simulation with input file(s) located at:\n\t"+files - def getTrackSpecies(self): - """ List the available tracked species """ + def _getParticleListSpecies(self, filePrefix): + """ List the available species in diagnostics of type ParticleList """ species = [] for path in self._results_path: - files = self._glob(path+self._os.sep+"TrackParticles*.h5") + files = self._glob(path+self._os.sep+filePrefix+"*.h5") for file in files: - s = self._re.search("^TrackParticlesDisordered_(.+).h5",self._os.path.basename(file)) + s = self._re.search("^"+filePrefix+"_(.+).h5",self._os.path.basename(file)) if s: species += [ s.groups()[0] ] return list(set(species)) # unique species + def getTrackSpecies(self): + """ List the available tracked species """ + return self._getParticleListSpecies("TrackParticlesDisordered") + + def getNewParticlesSpecies(self): + """ List the available NewParticles species """ + return self._getParticleListSpecies("NewParticles") + def fieldInfo(self, diag): """ Information on a specific Field diagnostic diff --git a/src/Checkpoint/Checkpoint.cpp b/src/Checkpoint/Checkpoint.cpp index cfd9ef776..ef5059ba7 100755 --- a/src/Checkpoint/Checkpoint.cpp +++ b/src/Checkpoint/Checkpoint.cpp @@ -323,7 +323,7 @@ void Checkpoint::dumpAll( VectorPatch &vecPatches, Region ®ion, unsigned int for( unsigned int idiag=0; idiag( vecPatches.localDiags[idiag] ) ) { ostringstream n( "" ); - n<< "latest_ID_" << vecPatches( 0 )->vecSpecies[track->speciesId_]->name_; + n<< "latest_ID_" << track->species_name_; f.attr( n.str(), track->latest_Id, H5T_NATIVE_UINT64 ); } } @@ -533,49 +533,17 @@ void Checkpoint::dumpPatch( Patch *patch, Params ¶ms, H5Write &g ) s.attr( "radiatedEnergy", spec->nrj_radiated_ ); if( spec->getNbrOfParticles()>0 ) { - - for( unsigned int i=0; iparticles->Position.size(); i++ ) { - ostringstream my_name( "" ); - my_name << "Position-" << i; - s.vect( my_name.str(), spec->particles->Position[i] );//, dump_deflate ); - } - - for( unsigned int i=0; iparticles->Momentum.size(); i++ ) { - ostringstream my_name( "" ); - my_name << "Momentum-" << i; - s.vect( my_name.str(),spec->particles->Momentum[i] );//, dump_deflate ); - } - - s.vect( "Weight", spec->particles->Weight );//, dump_deflate ); - s.vect( "Charge", spec->particles->Charge );//, dump_deflate ); - - if( spec->particles->tracked ) { - s.vect( "Id", spec->particles->Id, H5T_NATIVE_UINT64 );//, dump_deflate ); - } - - // Monte-Carlo process - if (spec->particles->has_Monte_Carlo_process) { - s.vect( "Tau", spec->particles->Tau );//, dump_deflate ); - } - - // Copy interpolated fields that must be accumulated over time - if( spec->particles->interpolated_fields_ ) { - if( spec->particles->interpolated_fields_->mode_[6] == 2 ) { - s.vect( "Wx", spec->particles->interpolated_fields_->F_[6] ); - } - if( spec->particles->interpolated_fields_->mode_[7] == 2 ) { - s.vect( "Wy", spec->particles->interpolated_fields_->F_[7] ); - } - if( spec->particles->interpolated_fields_->mode_[8] == 2 ) { - s.vect( "Wz", spec->particles->interpolated_fields_->F_[8] ); - } - } - + dumpParticles( s, *spec->particles ); s.vect( "first_index", spec->particles->first_index ); s.vect( "last_index", spec->particles->last_index ); - - } // End if partSize - + } + + // Copy birth records that haven't been written yet + if( spec->birth_records_ ) { + H5Write b = s.group( "birth_records" ); + b.vect( "birth_time", spec->birth_records_->birth_time_ ); + dumpParticles( b, spec->birth_records_->p_ ); + } } // End for ispec // Save some scalars @@ -722,7 +690,7 @@ void Checkpoint::restartAll( VectorPatch &vecPatches, Region ®ion, SmileiMPI for( unsigned int idiag=0; idiag( vecPatches.localDiags[idiag] ) ) { ostringstream n( "" ); - n<< "latest_ID_" << vecPatches( 0 )->vecSpecies[track->speciesId_]->name_; + n<< "latest_ID_" << track->species_name_; if( f.hasAttr( n.str() ) ) { f.attr( n.str(), track->latest_Id, H5T_NATIVE_UINT64 ); } else { @@ -987,50 +955,23 @@ void Checkpoint::restartPatch( Patch *patch, Params ¶ms, H5Read &g ) s.attr( "radiatedEnergy", spec->nrj_radiated_ ); if( partSize>0 ) { - for( unsigned int i=0; iparticles->Position.size(); i++ ) { - ostringstream namePos( "" ); - namePos << "Position-" << i; - s.vect( namePos.str(), spec->particles->Position[i] ); - } - - for( unsigned int i=0; iparticles->Momentum.size(); i++ ) { - ostringstream namePos( "" ); - namePos << "Momentum-" << i; - s.vect( namePos.str(), spec->particles->Momentum[i] ); - } - - s.vect( "Weight", spec->particles->Weight ); - - s.vect( "Charge", spec->particles->Charge ); - - if( spec->particles->tracked ) { - s.vect( "Id", spec->particles->Id, H5T_NATIVE_UINT64 ); - } - - if (spec->particles->has_Monte_Carlo_process) { - s.vect( "Tau", spec->particles->Tau ); - } - - // Retrieve interpolated fields that must be accumulated over time - if( spec->particles->interpolated_fields_ ) { - if( spec->particles->interpolated_fields_->mode_[6] == 2 ) { - s.vect( "Wx", spec->particles->interpolated_fields_->F_[6] ); - } - if( spec->particles->interpolated_fields_->mode_[7] == 2 ) { - s.vect( "Wy", spec->particles->interpolated_fields_->F_[7] ); - } - if( spec->particles->interpolated_fields_->mode_[8] == 2 ) { - s.vect( "Wz", spec->particles->interpolated_fields_->F_[8] ); - } - } + restartParticles( s, *spec->particles ); if( ! params.cell_sorting_ ) { s.vect( "first_index", spec->particles->first_index, true ); s.vect( "last_index", spec->particles->last_index, true ); } // When cell sorting is activated, indexes are recomputed directly after the restart. - } + + // Read birth records that haven't been written yet + if( spec->birth_records_ && s.has( "birth_records" ) ) { + H5Read b = s.group( "birth_records" ); + b.vect( "birth_time", spec->birth_records_->birth_time_, true ); + spec->birth_records_->p_.initialize( spec->birth_records_->birth_time_.size(), *spec->particles ); + restartParticles( b, spec->birth_records_->p_ ); + } + } // Load some scalars @@ -1103,6 +1044,85 @@ void Checkpoint::restart_cFieldsPerProc( H5Read &g, Field *field ) g.vect( field->name, *cfield->cdata_, H5T_NATIVE_DOUBLE ); } +void Checkpoint::dumpParticles( H5Write& s, Particles &p ) +{ + for( unsigned int i=0; imode_[6] == 2 ) { + s.vect( "Wx", p.interpolated_fields_->F_[6] ); + } + if( p.interpolated_fields_->mode_[7] == 2 ) { + s.vect( "Wy", p.interpolated_fields_->F_[7] ); + } + if( p.interpolated_fields_->mode_[8] == 2 ) { + s.vect( "Wz", p.interpolated_fields_->F_[8] ); + } + } +} + +void Checkpoint::restartParticles( H5Read& s, Particles &p ) +{ + for( unsigned int i=0; imode_[6] == 2 ) { + s.vect( "Wx", p.interpolated_fields_->F_[6] ); + } + if( p.interpolated_fields_->mode_[7] == 2 ) { + s.vect( "Wy", p.interpolated_fields_->F_[7] ); + } + if( p.interpolated_fields_->mode_[8] == 2 ) { + s.vect( "Wz", p.interpolated_fields_->F_[8] ); + } + } +} + void Checkpoint::dumpMovingWindow( H5Write &f, SimWindow *simWin ) { f.attr( "x_moved", simWin->getXmoved() ); diff --git a/src/Checkpoint/Checkpoint.h b/src/Checkpoint/Checkpoint.h index be0647d41..7dec58d26 100755 --- a/src/Checkpoint/Checkpoint.h +++ b/src/Checkpoint/Checkpoint.h @@ -51,13 +51,6 @@ class Checkpoint void restartAll( VectorPatch &vecPatches, Region ®ion, SmileiMPI *smpi, Params ¶ms ); void restartPatch( Patch *patch, Params ¶ms, H5Read &g ); - //! restart field per proc - void restartFieldsPerProc( H5Read &g, Field *field ); - void restart_cFieldsPerProc( H5Read &g, Field *field ); - - //! load moving window parameters - void restartMovingWindow( H5Read &f, SimWindow *simWindow ); - //! test before writing everything to file per processor //bool dump(unsigned int itime, double time, Params ¶ms); void dump( VectorPatch &vecPatches, Region ®ion, unsigned int itime, SmileiMPI *smpi, SimWindow *simWindow, Params ¶ms ); @@ -108,12 +101,17 @@ class Checkpoint //! initialize the time zero of the simulation void initDumpCases(); - //! dump field per proc + //! dump/restart field per proc void dumpFieldsPerProc( H5Write &g, Field *field ); void dump_cFieldsPerProc( H5Write &g, Field *field ); - - //! dump moving window parameters + void restartFieldsPerProc( H5Read &g, Field *field ); + void restart_cFieldsPerProc( H5Read &g, Field *field ); + //! dump/restart a particles object + void dumpParticles( H5Write& s, Particles &p ); + void restartParticles( H5Read& s, Particles &p ); + //! dump/restart moving window parameters void dumpMovingWindow( H5Write &f, SimWindow *simWindow ); + void restartMovingWindow( H5Read &f, SimWindow *simWindow ); //! function that returns elapsed time from creator (uses private var time_reference) //double time_seconds(); diff --git a/src/Collisions/CollisionalIonization.cpp b/src/Collisions/CollisionalIonization.cpp index 61e896c26..96d0489b1 100755 --- a/src/Collisions/CollisionalIonization.cpp +++ b/src/Collisions/CollisionalIonization.cpp @@ -271,7 +271,7 @@ void CollisionalIonization::calculate( double gamma_s, double gammae, double gam // Finish the ionization (moves new electrons in place) -void CollisionalIonization::finish( Params ¶ms, Patch *patch, std::vector &localDiags, bool, std::vector, std::vector, int ) +void CollisionalIonization::finish( Params ¶ms, Patch *patch, std::vector &localDiags, bool, std::vector, std::vector, int itime ) { - patch->vecSpecies[ionization_electrons_]->importParticles( params, patch, new_electrons, localDiags ); + patch->vecSpecies[ionization_electrons_]->importParticles( params, patch, new_electrons, localDiags, ( itime + 0.5 ) * params.timestep ); } diff --git a/src/Collisions/CollisionalNuclearReaction.cpp b/src/Collisions/CollisionalNuclearReaction.cpp index 2e5fb6ac3..249d82f9b 100644 --- a/src/Collisions/CollisionalNuclearReaction.cpp +++ b/src/Collisions/CollisionalNuclearReaction.cpp @@ -148,11 +148,11 @@ void CollisionalNuclearReaction::apply( Random *random, BinaryProcessData &D ) // Finish the reaction void CollisionalNuclearReaction::finish( Params ¶ms, Patch *patch, std::vector &localDiags, - bool intra_collisions, vector sg1, vector sg2, int + bool intra_collisions, vector sg1, vector sg2, int itime ) { // Move new particles in place for( unsigned int i=0; ivecSpecies[product_ispecies_[i]]->importParticles( params, patch, *product_particles_[i], localDiags ); + patch->vecSpecies[product_ispecies_[i]]->importParticles( params, patch, *product_particles_[i], localDiags, ( itime + 0.5 ) * params.timestep ); } // Remove reactants that have fully reacted (very rare) diff --git a/src/Diagnostic/Diagnostic.h b/src/Diagnostic/Diagnostic.h index 2c64a260f..140edee4a 100755 --- a/src/Diagnostic/Diagnostic.h +++ b/src/Diagnostic/Diagnostic.h @@ -5,8 +5,8 @@ #include "Patch.h" #include "Timers.h" -class Params; -class OpenPMDparams; +#include "OpenPMDparams.h" + class SmileiMPI; class VectorPatch; class TimeSelection; diff --git a/src/Diagnostic/DiagnosticFactory.h b/src/Diagnostic/DiagnosticFactory.h index 1f6299347..9379b9e47 100755 --- a/src/Diagnostic/DiagnosticFactory.h +++ b/src/Diagnostic/DiagnosticFactory.h @@ -9,6 +9,7 @@ #include "DiagnosticProbes.h" #include "DiagnosticScalar.h" #include "DiagnosticTrack.h" +#include "DiagnosticNewParticles.h" #include "DiagnosticPerformances.h" #include "DiagnosticFields1D.h" @@ -52,16 +53,16 @@ class DiagnosticFactory std::vector vecDiagnostics; vecDiagnostics.push_back( new DiagnosticScalar( params, smpi, vecPatches( 0 ) ) ); - for( unsigned int n_diag_particles = 0; n_diag_particles < PyTools::nComponents( "DiagParticleBinning" ); n_diag_particles++ ) { - vecDiagnostics.push_back( new DiagnosticParticleBinning( params, smpi, vecPatches( 0 ), n_diag_particles ) ); + for( unsigned int i = 0, n = PyTools::nComponents( "DiagParticleBinning" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticParticleBinning( params, smpi, vecPatches( 0 ), i ) ); } - for( unsigned int n_diag_screen = 0; n_diag_screen < PyTools::nComponents( "DiagScreen" ); n_diag_screen++ ) { - vecDiagnostics.push_back( new DiagnosticScreen( params, smpi, vecPatches( 0 ), n_diag_screen ) ); + for( unsigned int i = 0, n = PyTools::nComponents( "DiagScreen" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticScreen( params, smpi, vecPatches( 0 ), i ) ); } - - for (unsigned int n_diag_rad_spectrum = 0; n_diag_rad_spectrum < PyTools::nComponents("DiagRadiationSpectrum"); n_diag_rad_spectrum++) { - vecDiagnostics.push_back( new DiagnosticRadiationSpectrum(params, smpi, vecPatches(0), radiation_tables_ , n_diag_rad_spectrum) ); + + for( unsigned int i = 0, n = PyTools::nComponents( "DiagRadiationSpectrum" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticRadiationSpectrum(params, smpi, vecPatches(0), radiation_tables_ , i) ); } return vecDiagnostics; @@ -73,19 +74,21 @@ class DiagnosticFactory static std::vector createLocalDiagnostics( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, OpenPMDparams &openPMD ) { std::vector vecDiagnostics; - //MESSAGE("in create local diags: global dims after declaring vecdiag " << vecPatches(0)->EMfields->Jx_s[1]->number_of_points_); - for( unsigned int n_diag_fields = 0; n_diag_fields < PyTools::nComponents( "DiagFields" ); n_diag_fields++ ) { - vecDiagnostics.push_back( DiagnosticFieldsFactory::create( params, smpi, vecPatches, n_diag_fields, openPMD ) ); - // MESSAGE("in create local diags: global dims after creating and pushing back field diag " << vecPatches(0)->EMfields->Jx_s[1]->number_of_points_); + for( unsigned int i = 0, n = PyTools::nComponents( "DiagFields" ); i < n; i++ ) { + vecDiagnostics.push_back( DiagnosticFieldsFactory::create( params, smpi, vecPatches, i, openPMD ) ); + } + + for( unsigned int i = 0, n = PyTools::nComponents( "DiagProbe" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticProbes( params, smpi, vecPatches, i ) ); } - for( unsigned int n_diag_probe = 0; n_diag_probe < PyTools::nComponents( "DiagProbe" ); n_diag_probe++ ) { - vecDiagnostics.push_back( new DiagnosticProbes( params, smpi, vecPatches, n_diag_probe ) ); + for( unsigned int i = 0, n = PyTools::nComponents( "DiagTrackParticles" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticTrack( params, smpi, vecPatches, i, vecDiagnostics.size(), openPMD ) ); } - for( unsigned int n_diag_track = 0; n_diag_track < PyTools::nComponents( "DiagTrackParticles" ); n_diag_track++ ) { - vecDiagnostics.push_back( new DiagnosticTrack( params, smpi, vecPatches, n_diag_track, vecDiagnostics.size(), openPMD ) ); + for( unsigned int i = 0, n = PyTools::nComponents( "DiagNewParticles" ); i < n; i++ ) { + vecDiagnostics.push_back( new DiagnosticNewParticles( params, smpi, vecPatches, i, vecDiagnostics.size(), openPMD ) ); } if( PyTools::nComponents( "DiagPerformances" ) > 0 ) { @@ -102,7 +105,7 @@ class DiagnosticFactory { std::vector probes( 0 ); - for( unsigned int n_probe = 0; n_probe < PyTools::nComponents( "DiagProbe" ); n_probe++ ) { + for( unsigned int i = 0, n = PyTools::nComponents( "DiagProbe" ); i < n; i++ ) { probes.push_back( new ProbeParticles() ); } diff --git a/src/Diagnostic/DiagnosticNewParticles.cpp b/src/Diagnostic/DiagnosticNewParticles.cpp new file mode 100755 index 000000000..888873f22 --- /dev/null +++ b/src/Diagnostic/DiagnosticNewParticles.cpp @@ -0,0 +1,285 @@ +#include "PyTools.h" + +#include +#include + +#include "ParticleData.h" +#include "PeekAtSpecies.h" +#include "DiagnosticNewParticles.h" +#include "VectorPatch.h" +#include "Params.h" + +using namespace std; + +DiagnosticNewParticles::DiagnosticNewParticles( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, unsigned int iDiagnosticNewParticles, unsigned int, OpenPMDparams &oPMD ) : + DiagnosticParticleList( params, smpi, vecPatches, "DiagNewParticles", "NewParticles_", iDiagnosticNewParticles, oPMD ) +{ + write_id_ = vecPatches.species( 0, species_index_ )->particles->tracked; + + // Inform each patch about this diag + for( unsigned int ipatch=0; ipatchbirth_records_ = new BirthRecords( *s->particles ); + // Find out the other species that may create this one by ionization + for( auto s: vecPatches( ipatch )->vecSpecies ) { + if( s->Ionize && s->electron_species_index == species_index_ ) { + s->Ionize->save_ion_charge_ = true; + } + } + } +} + +DiagnosticNewParticles::~DiagnosticNewParticles() +{ + closeFile(); +} + +void DiagnosticNewParticles::openFile( Params ¶ms, SmileiMPI *smpi ) +{ + // Create HDF5 file + file_ = new H5Write( filename, &smpi->world() ); + file_->attr( "name", diag_name_ ); + + // Groups for openPMD + H5Write data_group( file_, "data" ); + H5Write iteration_group( &data_group, "0" ); + H5Write particles_group( &iteration_group, "particles" ); + H5Write species_group( &particles_group, species_name_ ); + + // Attributes for openPMD + openPMD_->writeRootAttributes( *file_, "no_meshes", "particles/" ); + openPMD_->writeBasePathAttributes( iteration_group, 0 ); + openPMD_->writeParticlesAttributes( particles_group ); + openPMD_->writeSpeciesAttributes( species_group ); + + // PositionOffset (for OpenPMD) + string xyz = "xyz"; + H5Write positionoffset_group = species_group.group( "positionOffset" ); + openPMD_->writeRecordAttributes( positionoffset_group, SMILEI_UNIT_POSITION ); + vector np = {0}; + for( unsigned int idim=0; idimwriteComponentAttributes( xyz_group, SMILEI_UNIT_POSITION ); + xyz_group.attr( "value", 0. ); + xyz_group.attr( "shape", np, H5T_NATIVE_UINT64 ); + } + + // Make empty datasets + H5Space file_space( 0, 0, 0, 1000, true ); + if( write_id_ ) { + loc_id_ = newDataset( species_group, "id", H5T_NATIVE_UINT64, file_space, SMILEI_UNIT_NONE ); + openPMD_->writeRecordAttributes( *loc_id_, SMILEI_UNIT_NONE ); + } + if( write_charge_ ) { + loc_charge_ = newDataset( species_group, "charge", H5T_NATIVE_SHORT, file_space, SMILEI_UNIT_CHARGE ); + openPMD_->writeRecordAttributes( *loc_charge_, SMILEI_UNIT_CHARGE ); + } + if( write_any_position_ ) { + H5Write position_group( &species_group, "position" ); + openPMD_->writeRecordAttributes( position_group, SMILEI_UNIT_POSITION ); + for( unsigned int idim=0; idimwriteRecordAttributes( momentum_group, SMILEI_UNIT_MOMENTUM ); + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_momentum_[idim] ) { + loc_momentum_[idim] = newDataset( momentum_group, xyz.substr( idim, 1 ).c_str(), H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_MOMENTUM ); + } + } + } + if( write_weight_ ) { + loc_weight_ = newDataset( species_group, "weight", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_DENSITY ); + openPMD_->writeRecordAttributes( *loc_weight_, SMILEI_UNIT_DENSITY ); + } + if( write_chi_ ) { + loc_chi_ = newDataset( species_group, "chi", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_NONE ); + openPMD_->writeRecordAttributes( *loc_chi_, SMILEI_UNIT_NONE ); + } + if( write_any_E_ ) { + H5Write E_group( &species_group, "E" ); + openPMD_->writeRecordAttributes( E_group, SMILEI_UNIT_EFIELD ); + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_E_[idim] ) { + loc_E_[idim] = newDataset( E_group, xyz.substr( idim, 1 ).c_str(), H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_EFIELD ); + } + } + } + if( write_any_B_ ) { + H5Write B_group( &species_group, "B" ); + openPMD_->writeRecordAttributes( B_group, SMILEI_UNIT_BFIELD ); + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_B_[idim] ) { + loc_B_[idim] = newDataset( B_group, xyz.substr( idim, 1 ).c_str(), H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_BFIELD ); + } + } + } + if( write_any_W_ ) { + H5Write W_group( &species_group, "W" ); + openPMD_->writeRecordAttributes( W_group, SMILEI_UNIT_ENERGY ); + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_W_[idim] ) { + loc_W_[idim] = newDataset( W_group, xyz.substr( idim, 1 ).c_str(), H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_ENERGY ); + } + } + } + loc_birth_time_ = newDataset( species_group, "birth_time", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_TIME ); + openPMD_->writeRecordAttributes( *loc_birth_time_, SMILEI_UNIT_TIME ); + + // Make a dataset supposed to contain the number of particles written at every output iteration + hsize_t ntimes = timeSelection->howManyTimesBefore( params.n_time ) + 1; + H5Space fs( {0, 2}, {0, 0}, {0, 2}, {ntimes, 2}, {true, false} ); + iteration_npart_ = new H5Write( file_, "iteration_npart", H5T_NATIVE_INT64, &fs ); + + file_->flush(); +} + +void DiagnosticNewParticles::closeFile() +{ + if( file_ ) { + for( auto d : loc_position_ ) delete d; + for( auto d : loc_momentum_ ) delete d; + delete loc_id_; + delete loc_charge_; + delete loc_weight_; + delete loc_chi_; + for( auto d : loc_E_ ) delete d; + for( auto d : loc_B_ ) delete d; + for( auto d : loc_W_ ) delete d; + delete loc_birth_time_; + delete iteration_npart_; + delete file_; + file_ = nullptr; + } +} + +void DiagnosticNewParticles::init( Params ¶ms, SmileiMPI *smpi, VectorPatch & ) +{ + // create the file + openFile( params, smpi ); +} + + +H5Space * DiagnosticNewParticles::prepareH5( SimWindow *, SmileiMPI *smpi, int itime, uint32_t nParticles_local, uint64_t nParticles_global, uint64_t offset ) +{ + // Resize datasets + hsize_t new_size = nParticles_written + nParticles_global; + loc_birth_time_->extend( new_size ); + if( write_id_ ) { + loc_id_->extend( new_size ); + } + if( write_charge_ ) { + loc_charge_->extend( new_size ); + } + if( write_any_position_ ) { + for( unsigned int idim=0; idimextend( new_size ); + } + } + } + if( write_any_momentum_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_momentum_[idim] ) { + loc_momentum_[idim]->extend( new_size ); + } + } + } + if( write_weight_ ) { + loc_weight_->extend( new_size ); + } + if( write_chi_ ) { + loc_chi_->extend( new_size ); + } + if( write_any_E_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_E_[idim] ) { + loc_E_[idim]->extend( new_size ); + } + } + } + if( write_any_B_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_B_[idim] ) { + loc_B_[idim]->extend( new_size ); + } + } + } + if( write_any_W_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_W_[idim] ) { + loc_W_[idim]->extend( new_size ); + } + } + } + // update iteration_npart_ + iteration_npart_->extend( {nTimes_written+1, 2} ); + if( smpi->isMaster() ) { + H5Space filespace( {nTimes_written+1, 2}, {nTimes_written, 0}, {1, 2} ); + H5Space memspace( 2 ); + uint64_t i_n[2] = { (uint64_t) itime, new_size }; + iteration_npart_->write( i_n[0], H5T_NATIVE_UINT64, &filespace, &memspace, true ); + } + nTimes_written += 1; + + // Filespace + hsize_t full_offset = nParticles_written + offset; + nParticles_written = new_size; + return new H5Space( new_size, full_offset, nParticles_local ); +} + +void DiagnosticNewParticles::writeOther( VectorPatch &vecPatches, size_t iprop, H5Space *file_space, H5Space *mem_space ) +{ + fill_buffer( vecPatches, iprop, data_double ); + #pragma omp master + write_scalar_double( loc_birth_time_, "birth_time", data_double[0], file_space, mem_space, SMILEI_UNIT_NONE ); + + // Delete the records + #pragma omp for schedule(static) + for( unsigned int ipatch=0 ; ipatchbirth_records_->clear(); + } +} + + +void DiagnosticNewParticles::write_scalar_uint64( H5Write * location, string /*name*/, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_UINT64, file_space, mem_space ); +} +void DiagnosticNewParticles::write_scalar_short( H5Write * location, string /*name*/, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_SHORT, file_space, mem_space ); +} +void DiagnosticNewParticles::write_scalar_double( H5Write * location, string /*name*/, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_DOUBLE, file_space, mem_space ); +} + +void DiagnosticNewParticles::write_component_uint64( H5Write * location, string /*name*/, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_UINT64, file_space, mem_space ); +} +void DiagnosticNewParticles::write_component_short( H5Write * location, string /*name*/, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_SHORT, file_space, mem_space ); +} +void DiagnosticNewParticles::write_component_double( H5Write * location, string /*name*/, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int /*unit_type*/ ) +{ + location->write( buffer, H5T_NATIVE_DOUBLE, file_space, mem_space ); +} + + +// SUPPOSED TO BE EXECUTED ONLY BY MASTER MPI +uint64_t DiagnosticNewParticles::getDiskFootPrint( int, int, Patch * ) +{ + return 0; +} + +Particles * DiagnosticNewParticles::getParticles( Patch * patch ) +{ + return &patch->vecSpecies[species_index_]->birth_records_->p_; +} diff --git a/src/Diagnostic/DiagnosticNewParticles.h b/src/Diagnostic/DiagnosticNewParticles.h new file mode 100755 index 000000000..d79f20542 --- /dev/null +++ b/src/Diagnostic/DiagnosticNewParticles.h @@ -0,0 +1,65 @@ +#ifndef DIAGNOSTICNEWPARTICLES_H +#define DIAGNOSTICNEWPARTICLES_H + +#include "DiagnosticParticleList.h" + +class Patch; +class Params; +class SmileiMPI; + +class DiagnosticNewParticles : public DiagnosticParticleList +{ + +public : + //! Default constructor + DiagnosticNewParticles( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, unsigned int, unsigned int, OpenPMDparams & ); + //! Default destructor + ~DiagnosticNewParticles() override; + + void openFile( Params ¶ms, SmileiMPI *smpi ) override; + + void closeFile() override; + + void init( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches ) override; + + //! Get disk footprint of current diagnostic + uint64_t getDiskFootPrint( int istart, int istop, Patch *patch ) override; + + //! Returns the Particles object of interest in a given patch + Particles * getParticles( Patch * patch ) override; + + //! Prepare all HDF5 groups, datasets and spaces + H5Space * prepareH5( SimWindow *simWindow, SmileiMPI *smpi, int itime, uint32_t nParticles_local, uint64_t nParticles_global, uint64_t offset ) override; + + //! Write extra particles properties + void writeOther( VectorPatch &, size_t, H5Space *, H5Space * ) override; + + //! Write a dataset + void write_scalar_uint64( H5Write * location, std::string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_scalar_short ( H5Write * location, std::string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_scalar_double( H5Write * location, std::string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_uint64( H5Write * location, std::string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_short ( H5Write * location, std::string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_double( H5Write * location, std::string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + + H5Write * newDataset( H5Write &group, std::string name, hid_t dtype, H5Space &file_space, unsigned int unit_type ) { + H5Write * d = new H5Write( &group, name, dtype, &file_space ); + openPMD_->writeComponentAttributes( *d, unit_type ); + return d; + }; + + +private : + + //! Number of particles previously written in the file + hsize_t nParticles_written = 0; + + //! Number of times the file was previously written + hsize_t nTimes_written = 0; + + //! Storage for the number of particles at each output iteration + H5Write * iteration_npart_ = nullptr; +}; + +#endif + diff --git a/src/Diagnostic/DiagnosticParticleList.cpp b/src/Diagnostic/DiagnosticParticleList.cpp new file mode 100755 index 000000000..31daeed72 --- /dev/null +++ b/src/Diagnostic/DiagnosticParticleList.cpp @@ -0,0 +1,481 @@ +#include "PyTools.h" + +#include +#include + +#include "ParticleData.h" +#include "PeekAtSpecies.h" +#include "DiagnosticParticleList.h" +#include "VectorPatch.h" +#include "Params.h" + +using namespace std; + +DiagnosticParticleList::DiagnosticParticleList( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, string diag_type, string file_prefix, unsigned int idiag_of_this_type, OpenPMDparams &oPMD ) : + Diagnostic( &oPMD, diag_type, idiag_of_this_type ), + nDim_particle( params.nDim_particle ) +{ + + // Extract the species + PyTools::extract( "species", species_name_, diag_type, idiag_of_this_type ); + vector species_names = {species_name_}; + vector species_ids = Params::FindSpecies( vecPatches( 0 )->vecSpecies, species_names ); + if( species_ids.size() > 1 ) { + ERROR( diag_type << " #" << idiag_of_this_type << " corresponds to more than 1 species" ); + } + if( species_ids.size() < 1 ) { + ERROR( diag_type << " #" << idiag_of_this_type << " does not correspond to any existing species" ); + } + species_index_ = species_ids[0]; + + ostringstream name( "" ); + name << diag_type << " with species '" << species_name_ << "'"; + + // Get parameter "every" which describes an iteration selection + timeSelection = new TimeSelection( PyTools::extract_py( "every", diag_type, idiag_of_this_type ), name.str() ); + + // Get parameter "flush_every" which decides the file flushing time selection + flush_timeSelection = new TimeSelection( PyTools::extract_py( "flush_every", diag_type, idiag_of_this_type ), name.str() ); + + // Get parameter "filter" which gives a python function to select particles + filter = PyTools::extract_py( "filter", diag_type, idiag_of_this_type ); + has_filter = ( filter != Py_None ); + if( has_filter ) { +#ifdef SMILEI_USE_NUMPY + // Test the filter with temporary, "fake" particles + name << " filter:"; + bool *dummy = NULL; + ParticleData test( nDim_particle, filter, name.str(), dummy ); +#else + ERROR( name.str() << " with a filter requires the numpy package" ); +#endif + } + + // Get the parameter "attributes": a list of attribute name that must be written + vector attributes( 0 ); + if( !PyTools::extractV( "attributes", attributes, diag_type, idiag_of_this_type ) ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": argument `attribute` must be a list of strings" ); + } + if( attributes.size() == 0 ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": argument `attribute` must have at least one element" ); + } + ostringstream attr_list( "" ); + InterpolatedFields * interpolated_fields = vecPatches( 0 )->vecSpecies[species_index_]->particles->interpolated_fields_; + for( unsigned int i=0; i1 ) { + write_position_[1] = true; + } else { + continue; + } + } else if( attributes[i] == "z" ) { + if( nDim_particle>2 ) { + write_position_[2] = true; + } else { + continue; + } + } else if( attributes[i] == "px" ) { + write_momentum_[0] = true; + } else if( attributes[i] == "py" ) { + write_momentum_[1] = true; + } else if( attributes[i] == "pz" ) { + write_momentum_[2] = true; + } else if( attributes[i] == "charge" || attributes[i] == "q" ) { + write_charge_ = true; + } else if( attributes[i] == "weight" || attributes[i] == "w" ) { + write_weight_ = true; + } else if( attributes[i] == "chi" ) { + write_chi_ = true; + } else if( attributes[i] == "Ex" ) { + write_E_[0] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[0] > 0 ); + } else if( attributes[i] == "Ey" ) { + write_E_[1] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[1] > 0 ); + } else if( attributes[i] == "Ez" ) { + write_E_[2] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[2] > 0 ); + } else if( attributes[i] == "Bx" ) { + write_B_[0] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[3] > 0 ); + } else if( attributes[i] == "By" ) { + write_B_[1] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[4] > 0 ); + } else if( attributes[i] == "Bz" ) { + write_B_[2] = true; + interpolate_ = interpolate_ || !( interpolated_fields && interpolated_fields->mode_[5] > 0 ); + } else if( attributes[i] == "Wx" ) { + write_W_[0] = true; + if( ! interpolated_fields || interpolated_fields->mode_[6] == 0 ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": attribute Wx requires `keep_interpolated_fields` to include 'Wx'" ); + } + } else if( attributes[i] == "Wy" ) { + write_W_[1] = true; + if( ! interpolated_fields || interpolated_fields->mode_[7] == 0 ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": attribute Wy requires `keep_interpolated_fields` to include 'Wy'" ); + } + } else if( attributes[i] == "Wz" ) { + write_W_[2] = true; + if( ! interpolated_fields || interpolated_fields->mode_[8] == 0 ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": attribute Wz requires `keep_interpolated_fields` to include 'Wz'" ); + } + } else { + ERROR( diag_type << " #" << idiag_of_this_type << ": attribute `" << attributes[i] << "` unknown" ); + } + attr_list << "," << attributes[i]; + } + write_any_position_ = write_position_[0] || write_position_[1] || write_position_[2]; + write_any_momentum_ = write_momentum_[0] || write_momentum_[1] || write_momentum_[2]; + write_any_E_ = write_E_[0] || write_E_[1] || write_E_[2]; + write_any_B_ = write_B_[0] || write_B_[1] || write_B_[2]; + write_any_W_ = write_W_[0] || write_W_[1] || write_W_[2]; + if( write_chi_ && ! vecPatches( 0 )->vecSpecies[species_index_]->particles->has_quantum_parameter ) { + ERROR( diag_type << " #" << idiag_of_this_type << ": attribute `chi` not available for this species" ); + } + + // Create the filename + ostringstream hdf_filename( "" ); + hdf_filename << file_prefix << species_name_ << ".h5" ; + filename = hdf_filename.str(); + + // Print some info + if( smpi->isMaster() ) { + MESSAGE( 1, "Created " << diag_type << " #" << idiag_of_this_type << ": species " << species_name_ ); + MESSAGE( 2, attr_list.str() ); + } +} + +DiagnosticParticleList::~DiagnosticParticleList() +{ + delete timeSelection; + delete flush_timeSelection; + Py_DECREF( filter ); +} + +bool DiagnosticParticleList::prepare( int itime ) +{ + return timeSelection->theTimeIsNow( itime ); +} + + +void DiagnosticParticleList::run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, SimWindow *simWindow, Timers & ) +{ + uint64_t nParticles_global = 0; + string xyz = "xyz"; + + H5Space *file_space=NULL, *mem_space=NULL; + #pragma omp master + { + // Obtain the particle partition of all the patches in this MPI + nParticles_local = 0; + patch_start.resize( vecPatches.size() ); + + if( has_filter ) { +#ifdef SMILEI_USE_NUMPY + patch_selection.resize( vecPatches.size() ); + PyArrayObject *ret; + ParticleData particleData( 0 ); + for( unsigned int ipatch=0 ; ipatchnumberOfParticles(); + if( npart > 0 ) { + // Expose particle data as numpy arrays + particleData.resize( npart ); + particleData.set( p ); + // run the filter function + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( filter, particleData.get(), NULL ); + PyTools::checkPyError(); + particleData.clear(); + if( ret == NULL ) { + ERROR( "A DiagTrackParticles filter has not provided a correct result" ); + } + // Loop the return value and store the selected particle index + bool *arr = ( bool * ) PyArray_GETPTR1( ret, 0 ); + for( unsigned int i=0; inumberOfParticles(); + } + } + + // Specify the memory dataspace (the size of the local buffer) + mem_space = new H5Space( (hsize_t)nParticles_local ); + + // Get the number of offset for this MPI rank + uint64_t np_local = nParticles_local, offset; + MPI_Scan( &np_local, &offset, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD ); + nParticles_global = offset; + offset -= np_local; + MPI_Bcast( &nParticles_global, 1, MPI_UNSIGNED_LONG_LONG, smpi->getSize()-1, MPI_COMM_WORLD ); + + // Prepare all HDF5 groups and datasets + file_space = prepareH5( simWindow, smpi, itime, nParticles_local, nParticles_global, offset ); + } + + // Id + if( write_id_ ) { + #pragma omp master + data_uint64.resize( nParticles_local ); + fill_buffer( vecPatches, 0, data_uint64 ); + #pragma omp master + { + write_scalar_uint64( loc_id_, "id", data_uint64[0], file_space, mem_space, SMILEI_UNIT_NONE ); + data_uint64.resize( 0 ); + } + } + + // Charge + if( write_charge_ ) { + #pragma omp master + data_short.resize( nParticles_local ); + fill_buffer( vecPatches, 0, data_short ); + #pragma omp master + { + write_scalar_short( loc_charge_, "charge", data_short[0], file_space, mem_space, SMILEI_UNIT_CHARGE ); + data_short.resize( 0 ); + } + } + + #pragma omp master + data_double.resize( nParticles_local ); + + // Position + if( write_any_position_ ) { + for( unsigned int idim=0; idimvecSpecies[species_index_]->mass_ != 1. && + vecPatches( 0 )->vecSpecies[species_index_]->mass_ > 0) { + for( unsigned int ip=0; ipvecSpecies[species_index_]->mass_; + } + } + write_component_double( loc_momentum_[idim], xyz.substr( idim, 1 ).c_str(), data_double[0], file_space, mem_space, SMILEI_UNIT_MOMENTUM ); + } + } + } + } + + // Properties included in Particles::double_prop are ordered a certain way that must be followed with care. + // iprop is the index of the property currently being treated + // positions (x, y, z) are indexed from 0 to nDim_particles + // momenta (px, py, pz) are indexed from nDim_particles to nDim_particles+3 + // etc + size_t iprop = nDim_particle+3; + + // Weight + 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 ); + } + + iprop++; + // If position_old exist, skip its components + if( ! getParticles( vecPatches( 0 ) )->Position_old.empty() ) { + iprop += nDim_particle; + } + + // Chi - quantum parameter + if( write_chi_ ) { + fill_buffer( vecPatches, iprop, data_double ); + #pragma omp master + write_scalar_double( loc_chi_, "chi", data_double[0], file_space, mem_space, SMILEI_UNIT_NONE ); + } + + #pragma omp barrier + if( getParticles( vecPatches( 0 ) )->has_quantum_parameter ) { + iprop++; + } + if( getParticles( vecPatches( 0 ) )->has_Monte_Carlo_process ) { + iprop++; + } + + InterpolatedFields * interpolated_fields = getParticles( vecPatches( 0 ) )->interpolated_fields_; + + // If field interpolation necessary + if( interpolate_ ) { + + #pragma omp master + data_double.resize( nParticles_local*6 ); + + // Do the interpolation + unsigned int nPatches = vecPatches.size(); + #pragma omp barrier + + if( has_filter ) { + #pragma omp for schedule(static) + for( unsigned int ipatch=0 ; ipatchInterp->fieldsSelection( + vecPatches.emfields( ipatch ), + *getParticles( vecPatches( ipatch ) ), + &data_double[patch_start[ipatch]], + ( int ) nParticles_local, + &patch_selection[ipatch] + ); + } + } else { + #pragma omp for schedule(static) + for( unsigned int ipatch=0 ; ipatchInterp->fieldsSelection( + vecPatches.emfields( ipatch ), + *getParticles( vecPatches( ipatch ) ), + &data_double[patch_start[ipatch]], + ( int ) nParticles_local, + NULL + ); + } + } + #pragma omp barrier + + // Write out the fields + #pragma omp master + { + if( write_any_E_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_E_[idim] ) { + write_component_double( loc_E_[idim], xyz.substr( idim, 1 ).c_str(), data_double[idim*nParticles_local], file_space, mem_space, SMILEI_UNIT_EFIELD ); + } + } + } + + if( write_any_B_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_B_[idim] ) { + write_component_double( loc_B_[idim], xyz.substr( idim, 1 ).c_str(), data_double[( 3+idim )*nParticles_local], file_space, mem_space, SMILEI_UNIT_BFIELD ); + } + } + } + } + + if( interpolated_fields ) { + iprop += count( interpolated_fields->mode_.begin(), interpolated_fields->mode_.end(), 1 ); + } + + // If field interpolation not necessary, maybe fields have been stored using `keep_interpolated_fields` + } else if( write_any_E_ || write_any_B_ ) { + + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_E_[idim] ) { + fill_buffer( vecPatches, iprop, data_double ); + #pragma omp master + write_component_double( loc_E_[idim], xyz.substr( idim, 1 ).c_str(), data_double[0], file_space, mem_space, SMILEI_UNIT_EFIELD ); + } + if( interpolated_fields->mode_[idim] > 0 ) { + iprop++; + } + } + + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_B_[idim] ) { + fill_buffer( vecPatches, iprop, data_double ); + #pragma omp master + write_component_double( loc_B_[idim], xyz.substr( idim, 1 ).c_str(), data_double[0], file_space, mem_space, SMILEI_UNIT_BFIELD ); + } + if( interpolated_fields->mode_[3+idim] > 0 ) { + iprop++; + } + } + + // If no fields required, just count the offset they create + } else if( interpolated_fields ) { + iprop += count( interpolated_fields->mode_.begin(), interpolated_fields->mode_.end(), 1 ); + } + + if( write_any_W_ ) { + for( unsigned int idim=0; idim<3; idim++ ) { + if( write_W_[idim] ) { + fill_buffer( vecPatches, iprop, data_double ); + #pragma omp master + write_component_double( loc_W_[idim], xyz.substr( idim, 1 ).c_str(), data_double[0], file_space, mem_space, SMILEI_UNIT_ENERGY ); + } + if( interpolated_fields->mode_[6+idim] > 0 ) { + iprop++; + } + } + } else if( interpolated_fields ) { + iprop += count( interpolated_fields->mode_.begin()+6, interpolated_fields->mode_.end(), 2 ); + } + + writeOther( vecPatches, iprop, file_space, mem_space ); + + #pragma omp master + { + data_double.resize( 0 ); + + // Close and flush + patch_selection.resize( 0 ); + + delete file_space; + delete mem_space; + deleteH5(); + + if( flush_timeSelection->theTimeIsNow( itime ) ) { + file_->flush(); + } + } + #pragma omp barrier +} + +template +void DiagnosticParticleList::fill_buffer( VectorPatch &vecPatches, size_t iprop, vector &buffer ) +{ + unsigned int patch_nParticles, i, j, nPatches=vecPatches.size(); + vector *property = NULL; + + #pragma omp barrier + if( has_filter ) { + #pragma omp for schedule(runtime) + for( unsigned int ipatch=0 ; ipatchgetProperty( iprop, property ); + i=0; + j=patch_start[ipatch]; + while( inumberOfParticles(); + p->getProperty( iprop, property ); + copy( property->begin(), property->begin() + patch_nParticles, buffer.begin() + patch_start[ipatch] ); + } + } +} diff --git a/src/Diagnostic/DiagnosticParticleList.h b/src/Diagnostic/DiagnosticParticleList.h new file mode 100755 index 000000000..acae40184 --- /dev/null +++ b/src/Diagnostic/DiagnosticParticleList.h @@ -0,0 +1,122 @@ +#ifndef DIAGNOSTICPARTICLELIST_H +#define DIAGNOSTICPARTICLELIST_H + +#include "Diagnostic.h" + +class Patch; +class Params; +class SmileiMPI; + + +class DiagnosticParticleList : public Diagnostic +{ + +public : + //! Default constructor + DiagnosticParticleList( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, std::string, std::string, unsigned int, OpenPMDparams & ); + //! Default destructor + ~DiagnosticParticleList() override; + + virtual void openFile( Params ¶ms, SmileiMPI *smpi ) = 0; + + virtual void closeFile() = 0; + + virtual void init( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches ) = 0; + + bool prepare( int itime ) override; + + void run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, SimWindow *simWindow, Timers &timers ) override; + + //! Get memory footprint of current diagnostic + int getMemFootPrint() override + { + return 0; + } + + //! Get disk footprint of current diagnostic + virtual uint64_t getDiskFootPrint( int istart, int istop, Patch *patch ) = 0; + + //! Returns the Particles object of interest in a given patch + virtual Particles * getParticles( Patch * patch ) = 0; + + //! Prepare all HDF5 groups, datasets and spaces + virtual H5Space * prepareH5( SimWindow *simWindow, SmileiMPI *smpi, int itime, uint32_t nParticles_local, uint64_t nParticles_global, uint64_t offset ) = 0; + + //! Close HDF5 groups, datasets and spaces + virtual void deleteH5() {}; + + //! Modify the filtered particles + virtual void modifyFiltered( VectorPatch &, unsigned int ) {}; + + //! Write extra particles properties + virtual void writeOther( VectorPatch &, size_t, H5Space *, H5Space * ) {}; + + //! Fills a buffer with the required particle property + template void fill_buffer( VectorPatch &vecPatches, size_t iprop, std::vector &buffer ); + + //! Write a dataset + virtual void write_scalar_uint64( H5Write * location, std::string name, uint64_t &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + virtual void write_scalar_short ( H5Write * location, std::string name, short &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + virtual void write_scalar_double( H5Write * location, std::string name, double &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + virtual void write_component_uint64( H5Write * location, std::string name, uint64_t &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + virtual void write_component_short ( H5Write * location, std::string name, short &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + virtual void write_component_double( H5Write * location, std::string name, double &, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) = 0; + + //! Index of the species used + unsigned int species_index_; + + //! Name of the species used + std::string species_name_; + +protected: + + //! Number of spatial dimensions + unsigned int nDim_particle; + + //! Current particle partition among the patches own by current MPI + std::vector patch_start; + + //! Tells whether this diag includes a particle filter + bool has_filter; + + //! Tells whether this diag includes a particle filter + PyObject *filter; + + //! Selection of the filtered particles in each patch + std::vector > patch_selection; + + //! Buffer for the output of double array + std::vector data_double; + //! Buffer for the output of short array + std::vector data_short; + //! Buffer for the output of uint64 array + std::vector data_uint64; + + //! Approximate total number of particles + double npart_total; + + //! Number of particles shared among patches in this proc + uint32_t nParticles_local; + + //! HDF5 locations where attributes must be written + bool write_id_ = false; H5Write* loc_id_ = nullptr; + bool write_charge_ = false; H5Write* loc_charge_ = nullptr; + bool write_weight_ = false; H5Write* loc_weight_ = nullptr; + bool write_chi_ = false; H5Write* loc_chi_ = nullptr; + std::vector write_position_ = {false, false, false}; std::vector loc_position_ = {nullptr, nullptr, nullptr}; + std::vector write_momentum_ = {false, false, false}; std::vector loc_momentum_ = {nullptr, nullptr, nullptr}; + std::vector write_E_ = {false, false, false}; std::vector loc_E_ = {nullptr, nullptr, nullptr}; + std::vector write_B_ = {false, false, false}; std::vector loc_B_ = {nullptr, nullptr, nullptr}; + std::vector write_W_ = {false, false, false}; std::vector loc_W_ = {nullptr, nullptr, nullptr}; + bool interpolate_ = false; + bool write_any_position_; + bool write_any_momentum_; + bool write_any_E_; + bool write_any_B_; + bool write_any_W_; + + H5Write* loc_birth_time_ = nullptr; +}; + +#endif + diff --git a/src/Diagnostic/DiagnosticTrack.cpp b/src/Diagnostic/DiagnosticTrack.cpp index 4f2446304..4d311f891 100755 --- a/src/Diagnostic/DiagnosticTrack.cpp +++ b/src/Diagnostic/DiagnosticTrack.cpp @@ -12,160 +12,19 @@ using namespace std; DiagnosticTrack::DiagnosticTrack( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches, unsigned int iDiagTrackParticles, unsigned int idiag, OpenPMDparams &oPMD ) : - Diagnostic( &oPMD, "DiagTrackParticles", iDiagTrackParticles ), - IDs_done( params.restart ), - nDim_particle( params.nDim_particle ) + DiagnosticParticleList( params, smpi, vecPatches, "DiagTrackParticles", "TrackParticlesDisordered_", iDiagTrackParticles, oPMD ), + IDs_done( params.restart ) { - - // Extract the species - string species_name; - PyTools::extract( "species", species_name, "DiagTrackParticles", iDiagTrackParticles ); - vector species_names = {species_name}; - vector species_ids = Params::FindSpecies( vecPatches( 0 )->vecSpecies, species_names ); - if( species_ids.size() > 1 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << " corresponds to more than 1 species" ); - } - if( species_ids.size() < 1 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << " does not correspond to any existing species" ); - } - speciesId_ = species_ids[0]; - - ostringstream name( "" ); - name << "Tracking species '" << species_name << "'"; - - // Get parameter "every" which describes an iteration selection - timeSelection = new TimeSelection( PyTools::extract_py( "every", "DiagTrackParticles", iDiagTrackParticles ), name.str() ); - - // Get parameter "flush_every" which decides the file flushing time selection - flush_timeSelection = new TimeSelection( PyTools::extract_py( "flush_every", "DiagTrackParticles", iDiagTrackParticles ), name.str() ); + write_id_ = true; // Inform each patch about this diag for( unsigned int ipatch=0; ipatchvecSpecies[speciesId_]->tracking_diagnostic = idiag; - } - - // Get parameter "filter" which gives a python function to select particles - filter = PyTools::extract_py( "filter", "DiagTrackParticles", iDiagTrackParticles ); - has_filter = ( filter != Py_None ); - if( has_filter ) { -#ifdef SMILEI_USE_NUMPY - // Test the filter with temporary, "fake" particles - name << " filter:"; - bool *dummy = NULL; - ParticleData test( nDim_particle, filter, name.str(), dummy ); -#else - ERROR( name.str() << " with a filter requires the numpy package" ); -#endif - } - - // Get the parameter "attributes": a list of attribute name that must be written - vector attributes( 0 ); - if( !PyTools::extractV( "attributes", attributes, "DiagTrackParticles", iDiagTrackParticles ) ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": argument `attribute` must be a list of strings" ); - } - if( attributes.size() == 0 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": argument `attribute` must have at least one element" ); - } - ostringstream attr_list( "" ); - attr_list << "id"; - write_position.resize( 3, false ); - write_momentum.resize( 3, false ); - write_charge = false; - write_weight = false; - write_chi = false; - write_E.resize( 3, false ); - write_B.resize( 3, false ); - write_W.resize( 3, false ); - interpolate = false; - InterpolatedFields * interpolated_fields = vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_; - for( unsigned int i=0; i1 ) { - write_position[1] = true; - } else { - continue; - } - } else if( attributes[i] == "z" ) { - if( nDim_particle>2 ) { - write_position[2] = true; - } else { - continue; - } - } else if( attributes[i] == "px" ) { - write_momentum[0] = true; - } else if( attributes[i] == "py" ) { - write_momentum[1] = true; - } else if( attributes[i] == "pz" ) { - write_momentum[2] = true; - } else if( attributes[i] == "charge" || attributes[i] == "q" ) { - write_charge = true; - } else if( attributes[i] == "weight" || attributes[i] == "w" ) { - write_weight = true; - } else if( attributes[i] == "chi" ) { - write_chi = true; - } else if( attributes[i] == "Ex" ) { - write_E[0] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[0] > 0 ); - } else if( attributes[i] == "Ey" ) { - write_E[1] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[1] > 0 ); - } else if( attributes[i] == "Ez" ) { - write_E[2] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[2] > 0 ); - } else if( attributes[i] == "Bx" ) { - write_B[0] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[3] > 0 ); - } else if( attributes[i] == "By" ) { - write_B[1] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[4] > 0 ); - } else if( attributes[i] == "Bz" ) { - write_B[2] = true; - interpolate = interpolate || !( interpolated_fields && interpolated_fields->mode_[5] > 0 ); - } else if( attributes[i] == "Wx" ) { - write_W[0] = true; - if( ! interpolated_fields || interpolated_fields->mode_[6] == 0 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": attribute Wx requires `keep_interpolated_fields` to include 'Wx'" ); - } - } else if( attributes[i] == "Wy" ) { - write_W[1] = true; - if( ! interpolated_fields || interpolated_fields->mode_[7] == 0 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": attribute Wy requires `keep_interpolated_fields` to include 'Wy'" ); - } - } else if( attributes[i] == "Wz" ) { - write_W[2] = true; - if( ! interpolated_fields || interpolated_fields->mode_[8] == 0 ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": attribute Wz requires `keep_interpolated_fields` to include 'Wz'" ); - } - } else { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": attribute `" << attributes[i] << "` unknown" ); - } - attr_list << "," << attributes[i]; - } - write_any_position = write_position[0] || write_position[1] || write_position[2]; - write_any_momentum = write_momentum[0] || write_momentum[1] || write_momentum[2]; - write_any_E = write_E[0] || write_E[1] || write_E[2]; - write_any_B = write_B[0] || write_B[1] || write_B[2]; - write_any_W = write_W[0] || write_W[1] || write_W[2]; - if( write_chi && ! vecPatches( 0 )->vecSpecies[speciesId_]->particles->has_quantum_parameter ) { - ERROR( "DiagTrackParticles #" << iDiagTrackParticles << ": attribute `chi` not available for this species" ); - } - - // Create the filename - ostringstream hdf_filename( "" ); - hdf_filename << "TrackParticlesDisordered_" << species_name << ".h5" ; - filename = hdf_filename.str(); - - // Print some info - if( smpi->isMaster() ) { - MESSAGE( 1, "Created TrackParticles #" << iDiagTrackParticles << ": species " << species_name ); - MESSAGE( 2, attr_list.str() ); + vecPatches( ipatch )->vecSpecies[species_index_]->tracking_diagnostic = idiag; } // Obtain the approximate number of particles in the species if( params.print_expected_disk_usage ) { - PeekAtSpecies peek( params, speciesId_ ); + PeekAtSpecies peek( params, species_index_ ); npart_total = peek.totalNumberofParticles(); } else { npart_total = 0; @@ -174,39 +33,32 @@ DiagnosticTrack::DiagnosticTrack( Params ¶ms, SmileiMPI *smpi, VectorPatch & DiagnosticTrack::~DiagnosticTrack() { - delete timeSelection; - delete flush_timeSelection; - Py_DECREF( filter ); closeFile(); } - void DiagnosticTrack::openFile( Params &, SmileiMPI *smpi ) { // Create HDF5 file file_ = new H5Write( filename, &smpi->world() ); - file_->attr( "name", diag_name_ ); // Attributes for openPMD openPMD_->writeRootAttributes( *file_, "no_meshes", "particles/" ); - data_group = new H5Write( file_, "data" ); + data_group_ = new H5Write( file_, "data" ); file_->flush(); } - void DiagnosticTrack::closeFile() { if( file_ ) { - delete data_group; + delete data_group_; delete file_; file_ = NULL; } } - void DiagnosticTrack::init( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches ) { // Set the IDs of the particles @@ -222,418 +74,119 @@ void DiagnosticTrack::init( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPat // create the file openFile( params, smpi ); - } -bool DiagnosticTrack::prepare( int itime ) +H5Space * DiagnosticTrack::prepareH5( SimWindow *simWindow, SmileiMPI *smpi, int itime, uint32_t nParticles_local, uint64_t nParticles_global, uint64_t offset ) { - return timeSelection->theTimeIsNow( itime ); -} - - -void DiagnosticTrack::run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, SimWindow *simWindow, Timers & ) -{ - uint64_t nParticles_global = 0; - string xyz = "xyz"; - - H5Write *momentum_group=NULL, *position_group=NULL, *species_group=NULL, *E_group=NULL, *B_group=NULL, *W_group=NULL; - H5Space *file_space=NULL, *mem_space=NULL; - #pragma omp master - { - // Obtain the particle partition of all the patches in this MPI - nParticles_local = 0; - patch_start.resize( vecPatches.size() ); - - if( has_filter ) { - -#ifdef SMILEI_USE_NUMPY - patch_selection.resize( vecPatches.size() ); - PyArrayObject *ret; - ParticleData particleData( 0 ); - for( unsigned int ipatch=0 ; ipatchvecSpecies[speciesId_]->particles; - unsigned int npart = p->numberOfParticles(); - if( npart > 0 ) { - // Expose particle data as numpy arrays - particleData.resize( npart ); - particleData.set( p ); - // run the filter function - ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( filter, particleData.get(), NULL ); - PyTools::checkPyError(); - particleData.clear(); - if( ret == NULL ) { - ERROR( "A DiagTrackParticles filter has not provided a correct result" ); - } - // Loop the return value and store the particle IDs - bool *arr = ( bool * ) PyArray_GETPTR1( ret, 0 ); - for( unsigned int i=0; iid( i ) & 72057594037927935) == 0 ) { - p->id( i ) += ++latest_Id; - } - } - } - Py_DECREF( ret ); - } - patch_start[ipatch] = nParticles_local; - nParticles_local += patch_selection[ipatch].size(); - } -#endif - - } else { - for( unsigned int ipatch=0 ; ipatchvecSpecies[speciesId_]->getNbrOfParticles(); - } - } - - // Specify the memory dataspace (the size of the local buffer) - mem_space = new H5Space( (hsize_t)nParticles_local ); - - // Get the number of offset for this MPI rank - uint64_t np_local = nParticles_local, offset; - MPI_Scan( &np_local, &offset, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD ); - nParticles_global = offset; - offset -= np_local; - MPI_Bcast( &nParticles_global, 1, MPI_UNSIGNED_LONG_LONG, smpi->getSize()-1, MPI_COMM_WORLD ); - - // Make a new group for this iteration - ostringstream t( "" ); - t << setfill( '0' ) << setw( 10 ) << itime; - // Create "data" group for openPMD compatibility - H5Write iteration_group = data_group->group( t.str() ); - H5Write particles_group = iteration_group.group( "particles" ); - species_group = new H5Write( &particles_group, vecPatches( 0 )->vecSpecies[speciesId_]->name_ ); - - // Add openPMD attributes ( "basePath" ) - openPMD_->writeBasePathAttributes( iteration_group, itime ); - // Add openPMD attributes ( "particles" ) - openPMD_->writeParticlesAttributes( particles_group ); - // Add openPMD attributes ( path of a given species ) - openPMD_->writeSpeciesAttributes( *species_group ); - - // Write x_moved - iteration_group.attr( "x_moved", simWindow ? simWindow->getXmoved() : 0. ); - - // Filespace and chunks - hsize_t chunk = 0; - if( nParticles_global>0 ) { - // Set the chunk size - unsigned int maximum_chunk_size = 100000000; - unsigned int number_of_chunks = nParticles_global/maximum_chunk_size; - if( nParticles_global%maximum_chunk_size != 0 ) { - number_of_chunks++; - } - if( number_of_chunks <= 1 ) { - chunk = 0; - } else { - unsigned int chunk_size = nParticles_global/number_of_chunks; - if( nParticles_global%number_of_chunks != 0 ) { - chunk_size++; - } - chunk = chunk_size; - } - } - file_space = new H5Space( nParticles_global, offset, nParticles_local, chunk ); - - // Create the "latest_IDs" dataset - // Create file space and select one element for each proc - iteration_group.vect( "latest_IDs", latest_Id, smpi->getSize(), H5T_NATIVE_UINT64, smpi->getRank(), 1 ); - + // Make a new group for this iteration + ostringstream t( "" ); + t << setfill( '0' ) << setw( 10 ) << itime; + // Create groups for openPMD compatibility + H5Write iteration_group = data_group_->group( t.str() ); + H5Write particles_group = iteration_group.group( "particles" ); + H5Write * species_group = new H5Write( &particles_group, species_name_ ); + loc_id_ = species_group; + loc_charge_ = species_group; + loc_weight_ = species_group; + loc_chi_ = species_group; + if( write_any_position_ ) { + fill( loc_position_.begin(), loc_position_.end(), new H5Write( species_group, "position" ) ); + openPMD_->writeRecordAttributes( *loc_position_[0], SMILEI_UNIT_POSITION ); } - - // Id - #pragma omp master - data_uint64.resize( nParticles_local, 1 ); - #pragma omp barrier - fill_buffer( vecPatches, 0, data_uint64 ); - #pragma omp master - { - write_scalar( species_group, "id", data_uint64[0], H5T_NATIVE_UINT64, file_space, mem_space, SMILEI_UNIT_NONE ); - data_uint64.resize( 0 ); + if( write_any_momentum_ ) { + fill( loc_momentum_.begin(), loc_momentum_.end(), new H5Write( species_group, "momentum" ) ); + openPMD_->writeRecordAttributes( *loc_momentum_[0], SMILEI_UNIT_MOMENTUM ); } - - // Charge - if( write_charge ) { - #pragma omp master - data_short.resize( nParticles_local, 0 ); - #pragma omp barrier - fill_buffer( vecPatches, 0, data_short ); - #pragma omp master - { - write_scalar( species_group, "charge", data_short[0], H5T_NATIVE_SHORT, file_space, mem_space, SMILEI_UNIT_CHARGE ); - data_short.resize( 0 ); - } + if( write_any_E_ ) { + fill( loc_E_.begin(), loc_E_.end(), new H5Write( species_group, "E" ) ); + openPMD_->writeRecordAttributes( *loc_E_[0], SMILEI_UNIT_EFIELD ); } - - #pragma omp master - data_double.resize( nParticles_local, 0 ); - - // Position - if( write_any_position ) { - #pragma omp master - { - position_group = new H5Write( species_group, "position" ); - openPMD_->writeRecordAttributes( *position_group, SMILEI_UNIT_POSITION ); - } - for( unsigned int idim=0; idimwriteRecordAttributes( *momentum_group, SMILEI_UNIT_MOMENTUM ); - } - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_momentum[idim] ) { - #pragma omp barrier - fill_buffer( vecPatches, nDim_particle+idim, data_double ); - #pragma omp master - { - // Multiply by the mass to obtain an actual momentum (except for photons (mass = 0)) - if( vecPatches( 0 )->vecSpecies[speciesId_]->mass_ != 1. && - vecPatches( 0 )->vecSpecies[speciesId_]->mass_ > 0) { - for( unsigned int ip=0; ipvecSpecies[speciesId_]->mass_; - } - } - write_component( momentum_group, xyz.substr( idim, 1 ).c_str(), data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_MOMENTUM ); - } - } - } - #pragma omp master - delete momentum_group; + if( write_any_B_ ) { + fill( loc_B_.begin(), loc_B_.end(), new H5Write( species_group, "B" ) ); + openPMD_->writeRecordAttributes( *loc_B_[0], SMILEI_UNIT_BFIELD ); } - - // Properties included in Particles::double_prop are ordered a certain way that must be followed with care. - // iprop is the index of the property currently being treated - // positions (x, y, z) are indexed from 0 to nDim_particles - // momenta (px, py, pz) are indexed from nDim_particles to nDim_particles+3 - // etc - size_t iprop = nDim_particle+3; - - // Weight - if( write_weight ) { - #pragma omp barrier - fill_buffer( vecPatches, iprop, data_double ); - #pragma omp master - write_scalar( species_group, "weight", data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_DENSITY ); - } - - iprop++; - // If position_old exist, skip its components - if( ! vecPatches( 0 )->vecSpecies[speciesId_]->particles->Position_old.empty() ) { - iprop += nDim_particle; + if( write_any_W_ ) { + fill( loc_W_.begin(), loc_W_.end(), new H5Write( species_group, "W" ) ); + openPMD_->writeRecordAttributes( *loc_W_[0], SMILEI_UNIT_ENERGY ); } - // Chi - quantum parameter - if( write_chi ) { - #pragma omp barrier - fill_buffer( vecPatches, iprop, data_double ); - #pragma omp master - write_scalar( species_group, "chi", data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_NONE ); - } - - #pragma omp barrier - if( vecPatches( 0 )->vecSpecies[speciesId_]->particles->has_quantum_parameter ) { - iprop++; - } - if( vecPatches( 0 )->vecSpecies[speciesId_]->particles->has_Monte_Carlo_process ) { - iprop++; + // PositionOffset (for OpenPMD) + string xyz = "xyz"; + H5Write positionoffset_group = species_group->group( "positionOffset" ); + openPMD_->writeRecordAttributes( positionoffset_group, SMILEI_UNIT_POSITION ); + vector np = {nParticles_global}; + for( unsigned int idim=0; idimwriteComponentAttributes( xyz_group, SMILEI_UNIT_POSITION ); + xyz_group.attr( "value", 0. ); + xyz_group.attr( "shape", np, H5T_NATIVE_UINT64 ); } - // If field interpolation necessary - if( interpolate ) { - - #pragma omp master - data_double.resize( nParticles_local*6 ); - - // Do the interpolation - unsigned int nPatches=vecPatches.size(); - #pragma omp barrier - - if( has_filter ) { - #pragma omp for schedule(static) - for( unsigned int ipatch=0 ; ipatchInterp->fieldsSelection( - vecPatches.emfields( ipatch ), - *( vecPatches.species( ipatch, speciesId_ )->particles ), - &data_double[patch_start[ipatch]], - ( int ) nParticles_local, - &patch_selection[ipatch] - ); - } - } else { - #pragma omp for schedule(static) - for( unsigned int ipatch=0 ; ipatchInterp->fieldsSelection( - vecPatches.emfields( ipatch ), - *( vecPatches.species( ipatch, speciesId_ )->particles ), - &data_double[patch_start[ipatch]], - ( int ) nParticles_local, - NULL - ); - } - } - #pragma omp barrier - - // Write out the fields - #pragma omp master - { - if( write_any_E ) { - H5Write Efield_group = species_group->group( "E" ); - openPMD_->writeRecordAttributes( Efield_group, SMILEI_UNIT_EFIELD ); - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_E[idim] ) { - write_component( &Efield_group, xyz.substr( idim, 1 ).c_str(), data_double[idim*nParticles_local], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_EFIELD ); - } - } - } - - if( write_any_B ) { - H5Write Bfield_group = species_group->group( "B" ); - openPMD_->writeRecordAttributes( Bfield_group, SMILEI_UNIT_BFIELD ); - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_B[idim] ) { - write_component( &Bfield_group, xyz.substr( idim, 1 ).c_str(), data_double[( 3+idim )*nParticles_local], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_BFIELD ); - } - } - } - } - - InterpolatedFields * interpolated_fields = vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_; - if( interpolated_fields ) { - iprop += count( interpolated_fields->mode_.begin(), interpolated_fields->mode_.end(), 1 ); - } - - // If field interpolation not necessary, it means fields have been stored using `keep_interpolated_fields` - } else if( write_any_E || write_any_B ) { - - #pragma omp master - if( write_any_E ) { - E_group = new H5Write( species_group, "E" ); - openPMD_->writeRecordAttributes( *E_group, SMILEI_UNIT_EFIELD ); - } - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_E[idim] ) { - #pragma omp barrier - fill_buffer( vecPatches, iprop, data_double ); - #pragma omp master - write_component( E_group, xyz.substr( idim, 1 ).c_str(), data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_EFIELD ); - } - if( vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_->mode_[idim] > 0 ) { - iprop++; - } - } - #pragma omp master - if( write_any_E ) { - delete E_group; - } - - #pragma omp master - if( write_any_B ) { - E_group = new H5Write( species_group, "B" ); - openPMD_->writeRecordAttributes( *B_group, SMILEI_UNIT_BFIELD ); - } - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_B[idim] ) { - #pragma omp barrier - fill_buffer( vecPatches, iprop, data_double ); - #pragma omp master - write_component( B_group, xyz.substr( idim, 1 ).c_str(), data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_BFIELD ); - } - if( vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_->mode_[3+idim] > 0 ) { - iprop++; - } - } - #pragma omp master - if( write_any_B ) { - delete B_group; - } - - } else { - - InterpolatedFields * interpolated_fields = vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_; - if( interpolated_fields ) { - iprop += count( interpolated_fields->mode_.begin(), interpolated_fields->mode_.end(), 1 ); - } - } + // Attributes for openPMD + openPMD_->writeBasePathAttributes( iteration_group, itime ); + openPMD_->writeParticlesAttributes( particles_group ); + openPMD_->writeSpeciesAttributes( *species_group ); - if( write_any_W ) { - #pragma omp master - { - W_group = new H5Write( species_group, "W" ); - openPMD_->writeRecordAttributes( *W_group, SMILEI_UNIT_ENERGY ); + // Write x_moved + iteration_group.attr( "x_moved", simWindow ? simWindow->getXmoved() : 0. ); + + // Create the "latest_IDs" dataset + // Create file space and select one element for each proc + iteration_group.vect( "latest_IDs", latest_Id, smpi->getSize(), H5T_NATIVE_UINT64, smpi->getRank(), 1 ); + + // Filespace and chunks + hsize_t chunk = 0; + if( nParticles_global>0 ) { + // Set the chunk size + unsigned int maximum_chunk_size = 100000000; + unsigned int number_of_chunks = nParticles_global/maximum_chunk_size; + if( nParticles_global%maximum_chunk_size != 0 ) { + number_of_chunks++; } - for( unsigned int idim=0; idim<3; idim++ ) { - if( write_W[idim] ) { - #pragma omp barrier - fill_buffer( vecPatches, iprop, data_double ); - #pragma omp master - write_component( W_group, xyz.substr( idim, 1 ).c_str(), data_double[0], H5T_NATIVE_DOUBLE, file_space, mem_space, SMILEI_UNIT_ENERGY ); - } - if( vecPatches( 0 )->vecSpecies[speciesId_]->particles->interpolated_fields_->mode_[6+idim] > 0 ) { - iprop++; + if( number_of_chunks <= 1 ) { + chunk = 0; + } else { + unsigned int chunk_size = nParticles_global/number_of_chunks; + if( nParticles_global%number_of_chunks != 0 ) { + chunk_size++; } + chunk = chunk_size; } - #pragma omp master - delete W_group; } - - #pragma omp master - { - data_double.resize( 0 ); - - // PositionOffset (for OpenPMD) - H5Write positionoffset_group = species_group->group( "positionOffset" ); - openPMD_->writeRecordAttributes( positionoffset_group, SMILEI_UNIT_POSITION ); - vector np = {nParticles_global}; - for( unsigned int idim=0; idimwriteComponentAttributes( xyz_group, SMILEI_UNIT_POSITION ); - xyz_group.attr( "value", 0. ); - xyz_group.attr( "shape", np, H5T_NATIVE_UINT64 ); - } - - // Close and flush - patch_selection.resize( 0 ); - - delete file_space; - delete mem_space; - delete species_group; - - if( flush_timeSelection->theTimeIsNow( itime ) ) { - file_->flush(); + return new H5Space( nParticles_global, offset, nParticles_local, chunk ); +} + +void DiagnosticTrack::deleteH5() +{ + delete loc_position_[0]; + delete loc_momentum_[0]; + delete loc_E_[0]; + delete loc_B_[0]; + delete loc_W_[0]; + delete loc_id_; +} + +void DiagnosticTrack::modifyFiltered( VectorPatch &vecPatches, unsigned int ipatch ) +{ + Particles *p = getParticles( vecPatches( ipatch ) ); + for( auto ipart: patch_selection[ipatch] ) { + // If particle not tracked before ( the 7 first bytes (ID<2^56) == 0 ), then set its ID + if(( p->id( ipart ) & 72057594037927935) == 0 ) { + p->id( ipart ) += ++latest_Id; } } - #pragma omp barrier } - void DiagnosticTrack::setIDs( Patch *patch ) { // If filter, IDs are set on-the-fly if( has_filter ) { return; } - unsigned int s = patch->vecSpecies[speciesId_]->getNbrOfParticles(); + unsigned int s = patch->vecSpecies[species_index_]->getNbrOfParticles(); for( unsigned int iPart=0; iPartvecSpecies[speciesId_]->particles->id( iPart ) = ++latest_Id; + patch->vecSpecies[species_index_]->particles->id( iPart ) = ++latest_Id; } } @@ -654,54 +207,38 @@ void DiagnosticTrack::setIDs( Particles &particles ) } -template -void DiagnosticTrack::fill_buffer( VectorPatch &vecPatches, unsigned int iprop, vector &buffer ) +void DiagnosticTrack::write_scalar_uint64( H5Write * location, string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) { - unsigned int patch_nParticles, i, j, nPatches=vecPatches.size(); - vector *property = NULL; - - if( has_filter ) { - #pragma omp for schedule(runtime) - for( unsigned int ipatch=0 ; ipatchvecSpecies[speciesId_]->particles->getProperty( iprop, property ); - i=0; - j=patch_start[ipatch]; - while( ivecSpecies[speciesId_]->getNbrOfParticles(); - vecPatches( ipatch )->vecSpecies[speciesId_]->particles->getProperty( iprop, property ); - i=0; - j=patch_start[ipatch]; - while( iarray( name, buffer, H5T_NATIVE_UINT64, file_space, mem_space ); + openPMD_->writeRecordAttributes( a, unit_type ); + openPMD_->writeComponentAttributes( a, unit_type ); } - - -template -void DiagnosticTrack::write_scalar( H5Write * location, string name, T &buffer, hid_t dtype, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) +void DiagnosticTrack::write_scalar_short( H5Write * location, string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) +{ + H5Write a = location->array( name, buffer, H5T_NATIVE_SHORT, file_space, mem_space ); + openPMD_->writeRecordAttributes( a, unit_type ); + openPMD_->writeComponentAttributes( a, unit_type ); +} +void DiagnosticTrack::write_scalar_double( H5Write * location, string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) { - H5Write a = location->array( name, buffer, dtype, file_space, mem_space ); + H5Write a = location->array( name, buffer, H5T_NATIVE_DOUBLE, file_space, mem_space ); openPMD_->writeRecordAttributes( a, unit_type ); openPMD_->writeComponentAttributes( a, unit_type ); } -template -void DiagnosticTrack::write_component( H5Write * location, string name, T &buffer, hid_t dtype, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) +void DiagnosticTrack::write_component_uint64( H5Write * location, string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) { - H5Write a = location->array( name, buffer, dtype, file_space, mem_space ); + H5Write a = location->array( name, buffer, H5T_NATIVE_UINT64, file_space, mem_space ); + openPMD_->writeComponentAttributes( a, unit_type ); +} +void DiagnosticTrack::write_component_short( H5Write * location, string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) +{ + H5Write a = location->array( name, buffer, H5T_NATIVE_SHORT, file_space, mem_space ); + openPMD_->writeComponentAttributes( a, unit_type ); +} +void DiagnosticTrack::write_component_double( H5Write * location, string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) +{ + H5Write a = location->array( name, buffer, H5T_NATIVE_DOUBLE, file_space, mem_space ); openPMD_->writeComponentAttributes( a, unit_type ); } @@ -729,3 +266,8 @@ uint64_t DiagnosticTrack::getDiskFootPrint( int istart, int istop, Patch * ) return footprint; } + +Particles * DiagnosticTrack::getParticles( Patch * patch ) +{ + return patch->vecSpecies[species_index_]->particles; +} diff --git a/src/Diagnostic/DiagnosticTrack.h b/src/Diagnostic/DiagnosticTrack.h index 02dd08df6..a6bcc65a3 100755 --- a/src/Diagnostic/DiagnosticTrack.h +++ b/src/Diagnostic/DiagnosticTrack.h @@ -1,14 +1,13 @@ #ifndef DIAGNOSTICTRACK_H #define DIAGNOSTICTRACK_H -#include "Diagnostic.h" +#include "DiagnosticParticleList.h" class Patch; class Params; class SmileiMPI; - -class DiagnosticTrack : public Diagnostic +class DiagnosticTrack : public DiagnosticParticleList { public : @@ -23,37 +22,35 @@ public : void init( Params ¶ms, SmileiMPI *smpi, VectorPatch &vecPatches ) override; - bool prepare( int itime ) override; - - void run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, SimWindow *simWindow, Timers &timers ) override; - - //! Get memory footprint of current diagnostic - int getMemFootPrint() override - { - return 0; - } - //! Get disk footprint of current diagnostic uint64_t getDiskFootPrint( int istart, int istop, Patch *patch ) override; - //! Fills a buffer with the required particle property - template void fill_buffer( VectorPatch &vecPatches, unsigned int iprop, std::vector &buffer ); + //! Returns the Particles object of interest in a given patch + Particles * getParticles( Patch * patch ) override; + + //! Prepare all HDF5 groups, datasets and spaces + H5Space * prepareH5( SimWindow *simWindow, SmileiMPI *smpi, int itime, uint32_t nParticles_local, uint64_t nParticles_global, uint64_t offset ) override; - //! Write a scalar dataset with the given buffer - template void write_scalar( H5Write*, std::string, T &, hid_t, H5Space*, H5Space*, unsigned int ); + //! Close HDF5 groups, datasets and spaces + void deleteH5() override; - //! Write a vector component dataset with the given buffer - template void write_component( H5Write*, std::string, T &, hid_t, H5Space*, H5Space*, unsigned int ); + //! Modify the filtered particles (apply new ID) + void modifyFiltered( VectorPatch &, unsigned int ) override; - //! Set a given patch's particles with the required IDs + //! Write a dataset + void write_scalar_uint64( H5Write * location, std::string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_scalar_short ( H5Write * location, std::string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_scalar_double( H5Write * location, std::string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_uint64( H5Write * location, std::string name, uint64_t &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_short ( H5Write * location, std::string name, short &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + void write_component_double( H5Write * location, std::string name, double &buffer, H5Space *file_space, H5Space *mem_space, unsigned int unit_type ) override; + + //! Set a given patch's particles with the required IDs (used at initialization & simWindow) void setIDs( Patch * ); - //! Set a given particles with the required IDs + //! Set a given particles with the required IDs (used by importParticles) void setIDs( Particles & ); - //! Index of the species used - unsigned int speciesId_; - //! Last ID assigned to a particle by this MPI domain uint64_t latest_Id; @@ -62,52 +59,7 @@ public : private : - H5Write *data_group; - - //! Number of spatial dimensions - unsigned int nDim_particle; - - //! Current particle partition among the patches own by current MPI - std::vector patch_start; - - //! Tells whether this diag includes a particle filter - bool has_filter; - - //! Tells whether this diag includes a particle filter - PyObject *filter; - - //! Selection of the filtered particles in each patch - std::vector > patch_selection; - - //! Buffer for the output of double array - std::vector data_double; - //! Buffer for the output of short array - std::vector data_short; - //! Buffer for the output of uint64 array - std::vector data_uint64; - - //! Approximate total number of particles - double npart_total; - - //! Number of particles shared among patches in this proc - uint32_t nParticles_local; - - //! Booleans to determine which attributes to write out - std::vector write_position; - std::vector write_momentum; - bool write_charge; - bool write_weight; - bool write_chi ; - std::vector write_E; - std::vector write_B; - std::vector write_W; - bool interpolate; - bool write_any_position; - bool write_any_momentum; - bool write_any_E; - bool write_any_B; - bool write_any_W; - + H5Write * data_group_; }; #endif diff --git a/src/Ionization/Ionization.cpp b/src/Ionization/Ionization.cpp index 7f4ec7735..debc48b71 100755 --- a/src/Ionization/Ionization.cpp +++ b/src/Ionization/Ionization.cpp @@ -1,6 +1,10 @@ +#include + #include "Ionization.h" #include "Species.h" +using namespace std; + Ionization::Ionization( Params ¶ms, Species *species ) { @@ -20,6 +24,7 @@ Ionization::Ionization( Params ¶ms, Species *species ) #ifdef _OMPTASKS new_electrons_per_bin = new Particles[species->Nbins]; + ion_charge_per_bin_.resize( species->Nbins ); #endif } @@ -29,33 +34,40 @@ Ionization::~Ionization() } -void Ionization::joinNewElectrons(unsigned int Nbins) +void Ionization::joinNewElectrons( unsigned int Nbins ) { - // if tasks on bins are used for ionization, join the lists of new electrons // created in each bin, to have the list of new electrons for this species and patch - for( unsigned int ibin = 0 ; ibin < Nbins ; ibin++ ) { - // number of particles to add from the bin - unsigned int n_particles_to_create = new_electrons_per_bin[ibin].size(); - new_electrons.createParticles(n_particles_to_create); - - for (unsigned int ipart = 0; ipart < n_particles_to_create ; ipart++){ - - int idNew = (new_electrons.size() - n_particles_to_create) + ipart; - - for( unsigned int i=0; i ion_charge_; + std::vector > ion_charge_per_bin_; + protected: double eV_to_au; diff --git a/src/Ionization/IonizationFromRate.cpp b/src/Ionization/IonizationFromRate.cpp index 3314302e2..57358f425 100755 --- a/src/Ionization/IonizationFromRate.cpp +++ b/src/Ionization/IonizationFromRate.cpp @@ -94,6 +94,10 @@ void IonizationFromRate::operator()( Particles *particles, unsigned int ipart_mi new_electrons.weight( idNew )=double( k_times )*particles->weight( ipart ); new_electrons.charge( idNew )=-1; + if( save_ion_charge_ ) { + ion_charge_.push_back( particles->charge( ipart ) ); + } + // Increase the charge of the particle particles->charge( ipart ) += k_times; } diff --git a/src/Ionization/IonizationTunnel.cpp b/src/Ionization/IonizationTunnel.cpp index 4c8ee1aec..6d396a980 100755 --- a/src/Ionization/IonizationTunnel.cpp +++ b/src/Ionization/IonizationTunnel.cpp @@ -170,6 +170,10 @@ void IonizationTunnel::operator()( Particles *particles, unsigned int ipart_min, new_electrons.weight( idNew )=double( k_times )*particles->weight( ipart ); new_electrons.charge( idNew )=-1; + if( save_ion_charge_ ) { + ion_charge_.push_back( particles->charge( ipart ) ); + } + // Increase the charge of the particle particles->charge( ipart ) += k_times; } @@ -300,6 +304,10 @@ void IonizationTunnel::ionizationTunnelWithTasks( Particles *particles, unsigned new_electrons_per_bin[ibin].weight( idNew )=double( k_times )*particles->weight( ipart ); new_electrons_per_bin[ibin].charge( idNew )=-1; + if( save_ion_charge_ ) { + ion_charge_per_bin_[ibin].push_back( particles->charge( ipart ) ); + } + // // Increase the charge of the particle particles->charge( ipart ) += k_times; } diff --git a/src/Ionization/IonizationTunnelEnvelopeAveraged.cpp b/src/Ionization/IonizationTunnelEnvelopeAveraged.cpp index 83ebf3e6a..f53ec02f2 100644 --- a/src/Ionization/IonizationTunnelEnvelopeAveraged.cpp +++ b/src/Ionization/IonizationTunnelEnvelopeAveraged.cpp @@ -242,6 +242,10 @@ void IonizationTunnelEnvelopeAveraged::envelopeIonization( Particles *particles, } + if( save_ion_charge_ ) { + ion_charge_.push_back( particles->charge( ipart ) ); + } + // weight and charge of the new electron new_electrons.weight( idNew )= particles->weight( ipart ); new_electrons.charge( idNew )= -1; @@ -302,6 +306,10 @@ void IonizationTunnelEnvelopeAveraged::envelopeIonization( Particles *particles, } + if( save_ion_charge_ ) { + ion_charge_per_bin_[ibin].push_back( particles->charge( ipart ) ); + } + // weight and charge of the new electron new_electrons_per_bin[ibin].weight( idNew )=particles->weight( ipart ); new_electrons_per_bin[ibin].charge( idNew )=-1; diff --git a/src/Params/OpenPMDparams.cpp b/src/Params/OpenPMDparams.cpp index bb1565b38..b7f4405e8 100755 --- a/src/Params/OpenPMDparams.cpp +++ b/src/Params/OpenPMDparams.cpp @@ -48,48 +48,35 @@ OpenPMDparams::OpenPMDparams( Params &p ): Wr=1.; } for( unsigned int unit_type=0; unit_type> &pold, size_t start, size_t n, size_t buffer_size, double mass_ ); //! Methods to obtain any property, given its index in the arrays double_prop_, uint64_prop_, or short_prop_ - void getProperty( unsigned int iprop, std::vector *&prop ) + void getProperty( size_t iprop, std::vector *&prop ) { prop = uint64_prop_[iprop]; } - void getProperty( unsigned int iprop, std::vector *&prop ) + void getProperty( size_t iprop, std::vector *&prop ) { prop = short_prop_[iprop]; } - void getProperty( unsigned int iprop, std::vector *&prop ) + void getProperty( size_t iprop, std::vector *&prop ) { prop = double_prop_[iprop]; } diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 5d5997317..554e41989 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -771,7 +771,7 @@ void VectorPatch::injectParticlesFromBoundaries(Params ¶ms, Timers &timers, for (unsigned int i_injector=0 ; i_injectorparticle_injector_vector_.size() ; i_injector++) { ParticleInjector * particle_injector = patch->particle_injector_vector_[i_injector]; Species * injector_species = species( ipatch, particle_injector->getSpeciesNumber() ); - injector_species->importParticles( params, patches_[ipatch], local_particles_vector[i_injector], localDiags ); + injector_species->importParticles( params, patches_[ipatch], local_particles_vector[i_injector], localDiags, ( itime + 0.5 ) * params.timestep ); } } // end for ipatch diff --git a/src/Profiles/Function.h b/src/Profiles/Function.h index 0b8154d8b..ce7f934ab 100755 --- a/src/Profiles/Function.h +++ b/src/Profiles/Function.h @@ -168,7 +168,7 @@ class Function_File : public Function Function_File( std::string path, std::string dataset_name, H5Read *file, std::vector cell_length ) : path_( path ), dataset_name_( dataset_name ), file_( file ), cell_length_( cell_length ) { - opened_file_count_ = new int( 0 ); + opened_file_count_ = new int( 1 ); }; Function_File( Function_File *f ) :path_( f->path_ ), dataset_name_( f->dataset_name_ ), file_( f->file_ ), cell_length_( f->cell_length_ ) diff --git a/src/Python/pycontrol.py b/src/Python/pycontrol.py index f255736ff..7ce08c20b 100755 --- a/src/Python/pycontrol.py +++ b/src/Python/pycontrol.py @@ -37,7 +37,8 @@ def _smilei_check(): # Verify classes were not overriden for CheckClassName in ["SmileiComponent","Species", "Laser","Collisions", "DiagProbe","DiagParticleBinning", "DiagScalar","DiagFields", - "DiagTrackParticles","DiagPerformances","ExternalField","PrescribedField", + "DiagTrackParticles","DiagNewParticles","DiagPerformances", + "ExternalField","PrescribedField", "SmileiSingleton","Main","Checkpoints","LoadBalancing","MovingWindow", "RadiationReaction", "ParticleData", "MultiphotonBreitWheeler", "Vectorization", "MultipleDecomposition"]: @@ -159,9 +160,8 @@ def _keep_python_running(): if len(LoadBalancing)>0 and len(MultipleDecomposition)>0: return True # Verify the tracked species that require a particle selection - for d in DiagTrackParticles: - if d.filter is not None: - return True + if any([d.filter for d in DiagTrackParticles]) or any([d.filter for d in DiagNewParticles]): + return True # Verify the particle binning having a function for deposited_quantity or axis type for d in DiagParticleBinning._list + DiagScreen._list: if type(d.deposited_quantity) is not str: diff --git a/src/Python/pyinit.py b/src/Python/pyinit.py index b5cfca1ba..84b9fe79a 100755 --- a/src/Python/pyinit.py +++ b/src/Python/pyinit.py @@ -561,6 +561,15 @@ class DiagTrackParticles(SmileiComponent): filter = None attributes = ["x", "y", "z", "px", "py", "pz", "w"] +class DiagNewParticles(SmileiComponent): + """Track diagnostic""" + name = "" + species = None + every = 0 + flush_every = 1 + filter = None + attributes = ["x", "y", "z", "px", "py", "pz", "w"] + class DiagPerformances(SmileiSingleton): """Performances diagnostic""" every = 0 diff --git a/src/Species/BirthRecords.h b/src/Species/BirthRecords.h new file mode 100755 index 000000000..98a28051b --- /dev/null +++ b/src/Species/BirthRecords.h @@ -0,0 +1,36 @@ +#ifndef BIRTHRECORDS_H +#define BIRTHRECORDS_H + +#include + +#include "Particles.h" +#include "Ionization.h" + +//! Records some information when the particles are created by ionization and other mechanisms +//! Required for DiagNewParticles +struct BirthRecords +{ + BirthRecords( Particles &source_particles ) { + p_.initialize( 0, source_particles ); + p_.double_prop_.push_back( &birth_time_ ); // the birth time is the last property + }; + void clear() { + birth_time_.resize( 0 ); + p_.resize( 0 ); + }; + void update( Particles &source_particles, size_t npart, double time_dual, Ionization *I ) { + size_t prev_size = birth_time_.size(); + birth_time_.resize( prev_size + npart, time_dual ); + source_particles.copyParticles( 0, npart, p_, prev_size ); + // If electrons come from ionization, then the charge is replaced by that of the ionized ions + if( I && I->save_ion_charge_ ) { + copy( I->ion_charge_.begin(), I->ion_charge_.end(), &p_.Charge[prev_size] ); + I->ion_charge_.clear(); + } + } + + std::vector birth_time_; + Particles p_; +}; + +#endif diff --git a/src/Species/Species.cpp b/src/Species/Species.cpp index 129484305..a4617f57c 100755 --- a/src/Species/Species.cpp +++ b/src/Species/Species.cpp @@ -81,7 +81,6 @@ Species::Species( Params ¶ms, Patch *patch ) : position_initialization_on_species_( false ), position_initialization_on_species_index_( -1 ), electron_species( NULL ), - electron_species_index( -1 ), photon_species_( nullptr ), //photon_species_index(-1), radiation_photon_species( "" ), @@ -407,42 +406,21 @@ void Species::initOperators( Params ¶ms, Patch *patch ) // --------------------------------------------------------------------------------------------------------------------- Species::~Species() { - delete particles; delete particles_to_move; - + delete Push; delete Interp; delete Proj; - - if( Merge ) { - delete Merge; - } - - if( Ionize ) { - delete Ionize; - } - if( Radiate ) { - delete Radiate; - } - if( part_comp_time_ ) { - delete part_comp_time_; - } - if( Multiphoton_Breit_Wheeler_process ) { - delete Multiphoton_Breit_Wheeler_process; - } - if( partBoundCond ) { - delete partBoundCond; - } - if( particles_per_cell_profile_ ) { - delete particles_per_cell_profile_; - } - if( charge_profile_ ) { - delete charge_profile_; - } - if( density_profile_ ) { - delete density_profile_; - } + delete Merge; + delete Ionize; + delete Radiate; + delete part_comp_time_; + delete Multiphoton_Breit_Wheeler_process; + delete partBoundCond; + delete particles_per_cell_profile_; + delete charge_profile_; + delete density_profile_; for( unsigned int i=0; iimportParticles( params, patch, Ionize->new_electrons, localDiags ); + electron_species->importParticles( params, patch, Ionize->new_electrons, localDiags, time_dual, Ionize ); } // if moving particle - if( time_dual>time_frozen_ ) { // moving particle + if( time_dual>time_frozen_ ) { - // Radiation losses - if( Radiate ) { - // If creation of macro-photon, we add them to photon_species - if( photon_species_ ) { - photon_species_->importParticles( params, - patch, - *radiated_photons_, - localDiags ); - } + // Radiation losses: If creation of macro-photon, we add them to photon_species + if( Radiate && photon_species_ ) { + photon_species_->importParticles( params, patch, *radiated_photons_, localDiags, time_dual ); } // Multiphoton Breit-Wheeler if( Multiphoton_Breit_Wheeler_process ) { // Addition of the electron-positron particles for( int k=0; k<2; k++ ) { - mBW_pair_species_[k]->importParticles( params, - patch, - *mBW_pair_particles_[k], - localDiags ); + mBW_pair_species_[k]->importParticles( params, patch, *mBW_pair_particles_[k], localDiags, time_dual ); } } - }//END if time vs. time_frozen_ + } } @@ -1904,7 +1869,7 @@ void Species::countSortParticles( Params ¶ms ) } // Move all particles from another species to this one -void Species::importParticles( Params ¶ms, Patch *patch, Particles &source_particles, vector &localDiags ) +void Species::importParticles( Params ¶ms, Patch *patch, Particles &source_particles, vector &localDiags, double time_dual, Ionization *I ) { unsigned int npart = source_particles.size(); unsigned int nbin = particles->numberOfBins(); @@ -1914,7 +1879,12 @@ void Species::importParticles( Params ¶ms, Patch *patch, Particles &source_p if( particles->tracked ) { dynamic_cast( localDiags[tracking_diagnostic] )->setIDs( source_particles ); } - + + // If there is a diagnostic for recording particle birth, then copy new particles to the buffer + if( birth_records_ ) { + birth_records_->update( source_particles, npart, time_dual, I ); + } + // Move particles vector src_bin_keys( npart, 0 ); for( unsigned int i=0; i & ); + virtual void importParticles( Params &, Patch *, Particles &, std::vector &, double time_dual, Ionization *I = nullptr ); //! This method eliminates the space gap between the bins //! (presence of empty particles between the bins) diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index d8c6be7f1..4c970d008 100755 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -1101,7 +1101,11 @@ class SpeciesFactory new_species->particles->interpolated_fields_ = new InterpolatedFields(); new_species->particles->interpolated_fields_->mode_ = species->particles->interpolated_fields_->mode_; } - + + if( species->birth_records_ ) { + new_species->birth_records_ = new BirthRecords( *species->particles ); + } + return new_species; } // End Species* clone() @@ -1136,166 +1140,131 @@ class SpeciesFactory // Loop species to find related species for( unsigned int ispec1 = 0; ispec1 < tot_species_number; ispec1++ ) { - + Species &s1 = *patch->vecSpecies[ispec1]; + // Ionizable species - if( patch->vecSpecies[ispec1]->Ionize ) { - // Loop all other species - for( unsigned int ispec2 = 0; ispec2vecSpecies.size(); ispec2++ ) { - if( patch->vecSpecies[ispec1]->ionization_electrons == patch->vecSpecies[ispec2]->name_ ) { - if( ispec1==ispec2 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_<<"' ionization_electrons must be a distinct species", - LINK_NAMELIST + std::string("#species") ); - } - if( patch->vecSpecies[ispec2]->mass_!=1 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_<<"' ionization_electrons must be a species with mass==1", - LINK_NAMELIST + std::string("#species") ); - } - patch->vecSpecies[ispec1]->electron_species_index = ispec2; - patch->vecSpecies[ispec1]->electron_species = patch->vecSpecies[ispec2]; - - // int max_eon_number = - // patch->vecSpecies[ispec1]->getNbrOfParticles() - // * ( patch->vecSpecies[ispec1]->atomic_number_ || patch->vecSpecies[ispec1]->maximum_charge_state_ ); - // patch->vecSpecies[ispec1]->Ionize->new_electrons.initializeReserve( - // max_eon_number, *patch->vecSpecies[ispec1]->electron_species->particles - patch->vecSpecies[ispec1]->Ionize->new_electrons.initialize( - 0, *patch->vecSpecies[ispec1]->electron_species->particles - ); -#ifdef _OMPTASKS - unsigned int Nbins = patch->vecSpecies[ispec1]->Nbins; - for (unsigned int ibin = 0 ; ibin < Nbins ; ibin++){ - patch->vecSpecies[ispec1]->Ionize->new_electrons_per_bin[ibin].initializeReserve(0, *patch->vecSpecies[ispec1]->electron_species->particles); - } -#endif - break; - } + if( s1.Ionize ) { + // Find electron species + auto s2 = std::find_if( patch->vecSpecies.begin(), patch->vecSpecies.end(), [&s1](Species *s) { return s->name_ == s1.ionization_electrons; } ); + if( s2 == patch->vecSpecies.end() ) { + ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->electron_species_index==-1 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_ - <<"' ionization_electrons named " << patch->vecSpecies[ispec1]->ionization_electrons - << " could not be found", - LINK_NAMELIST + std::string("#species") ); + s1.electron_species = *s2; + s1.electron_species_index = std::distance( patch->vecSpecies.begin(), s2 ); + if( s1.name_ == s1.electron_species->name_ ) { + ERROR_NAMELIST( "For species '"<mass_!=1 ) { + ERROR_NAMELIST( "For species '"<new_electrons.initializeReserve( max_eon_number, *s1.electron_species->particles + s1.Ionize->new_electrons.initialize( 0, *s1.electron_species->particles ); +#ifdef _OMPTASKS + for( unsigned int ibin = 0 ; ibin < s1.Nbins ; ibin++ ){ + s1.Ionize->new_electrons_per_bin[ibin].initializeReserve( 0, *s1.electron_species->particles ); + } +#endif } // Radiating species - if( patch->vecSpecies[ispec1]->Radiate ) { + if( s1.Radiate ) { // No emission of discrete photon, only scalar diagnostics are updated - if( patch->vecSpecies[ispec1]->radiation_photon_species.empty() ) { - patch->vecSpecies[ispec1]->photon_species_index = -1; - patch->vecSpecies[ispec1]->photon_species_ = nullptr; - patch->vecSpecies[ispec1]->radiated_photons_ = nullptr; + if( s1.radiation_photon_species.empty() ) { + s1.photon_species_index = -1; + s1.photon_species_ = nullptr; + s1.radiated_photons_ = nullptr; } // Else, there will be emission of macro-photons. else { unsigned int ispec2 = 0; for( ispec2 = 0; ispec2vecSpecies.size(); ispec2++ ) { - if( patch->vecSpecies[ispec1]->radiation_photon_species == patch->vecSpecies[ispec2]->name_ ) { + if( s1.radiation_photon_species == patch->vecSpecies[ispec2]->name_ ) { if( ispec1==ispec2 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_ - <<"' radiation_photon_species must be a distinct photon species", + ERROR_NAMELIST( "For species '"<vecSpecies[ispec2]->mass_!=0 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_ - <<"' radiation_photon_species must be a photon species with mass==0", + ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->photon_species_index = ispec2; - patch->vecSpecies[ispec1]->photon_species_ = patch->vecSpecies[ispec2]; - patch->vecSpecies[ispec1]->radiated_photons_ = ParticlesFactory::create( params ); - // patch->vecSpecies[ispec1]->radiated_photons_->initializeReserve( - // patch->vecSpecies[ispec1]->getNbrOfParticles(), - // *patch->vecSpecies[ispec1]->photon_species_->particles - // ); - patch->vecSpecies[ispec1]->radiated_photons_->initialize( - 0, - *patch->vecSpecies[ispec1]->photon_species_->particles - ); + s1.photon_species_index = ispec2; + s1.photon_species_ = patch->vecSpecies[ispec2]; + s1.radiated_photons_ = ParticlesFactory::create( params ); + // s1.radiated_photons_->initializeReserve( s1.getNbrOfParticles(), *s1.photon_species_->particles ); + s1.radiated_photons_->initialize( 0, *s1.photon_species_->particles ); #ifdef _OMPTASKS - unsigned int Nbins = patch->vecSpecies[ispec1]->Nbins; - for (unsigned int ibin = 0 ; ibin < Nbins ; ibin++){ - patch->vecSpecies[ispec1]->Radiate->new_photons_per_bin_[ibin].initializeReserve( - patch->vecSpecies[ispec1]->getNbrOfParticles(), - *patch->vecSpecies[ispec1]->photon_species_->particles - ); + for (unsigned int ibin = 0 ; ibin < s1.Nbins ; ibin++){ + s1.Radiate->new_photons_per_bin_[ibin].initializeReserve( s1.getNbrOfParticles(), *s1.photon_species_->particles ); } #endif break; } } if( ispec2 == patch->vecSpecies.size() ) { - ERROR_NAMELIST( "Species '" << patch->vecSpecies[ispec1]->radiation_photon_species - << "' does not exist.", + ERROR_NAMELIST( "Species '" << s1.radiation_photon_species << "' does not exist.", LINK_NAMELIST + std::string("#species") ) } } } // Breit-Wheeler species - if( patch->vecSpecies[ispec1]->Multiphoton_Breit_Wheeler_process ) { + if( s1.Multiphoton_Breit_Wheeler_process ) { unsigned int ispec2; for( int k=0; k<2; k++ ) { ispec2 = 0; while( ispec2vecSpecies.size()) { // We look for the pair species mBW_pair_species_names_[k] - if( patch->vecSpecies[ispec1]->mBW_pair_species_names_[k] == patch->vecSpecies[ispec2]->name_ ) { + if( s1.mBW_pair_species_names_[k] == patch->vecSpecies[ispec2]->name_ ) { if( ispec1==ispec2 ) { - ERROR_NAMELIST( "For species '" << patch->vecSpecies[ispec1]->name_ + ERROR_NAMELIST( "For species '" << s1.name_ << "' pair species must be a distinct particle species", LINK_NAMELIST + std::string("#species") ); } if( patch->vecSpecies[ispec2]->mass_ != 1 ) { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_ + ERROR_NAMELIST( "For species '"<vecSpecies[ispec2]->charge_profile_->getProfileName() != "constant") { - ERROR_NAMELIST( "For species '"<vecSpecies[ispec1]->name_ + ERROR_NAMELIST( "For species '"<vecSpecies[ispec2]->max_charge_) != 1) { - ERROR_NAMELIST( "For species ``"<vecSpecies[ispec1]->name_ + ERROR_NAMELIST( "For species ``"<vecSpecies[ispec2]->max_charge_ << ") is not correct.", LINK_NAMELIST + std::string("#species") ); } - patch->vecSpecies[ispec1]->mBW_pair_species_index_[k] = ispec2; - patch->vecSpecies[ispec1]->mBW_pair_species_[k] = patch->vecSpecies[ispec2]; + s1.mBW_pair_species_index_[k] = ispec2; + s1.mBW_pair_species_[k] = patch->vecSpecies[ispec2]; - patch->vecSpecies[ispec1]->mBW_pair_particles_[k] = ParticlesFactory::create( params ); + s1.mBW_pair_particles_[k] = ParticlesFactory::create( params ); - // patch->vecSpecies[ispec1]->mBW_pair_particles_[k]->initializeReserve( - // patch->vecSpecies[ispec1]->getNbrOfParticles(), - // *patch->vecSpecies[ispec1]->mBW_pair_species_[k]->particles - // ); - patch->vecSpecies[ispec1]->mBW_pair_particles_[k]->initialize( - 0, - *patch->vecSpecies[ispec1]->mBW_pair_species_[k]->particles - ); + // s1.mBW_pair_particles_[k]->initializeReserve( s1.getNbrOfParticles(), *s1.mBW_pair_species_[k]->particles ); + s1.mBW_pair_particles_[k]->initialize( 0, *s1.mBW_pair_species_[k]->particles ); #ifdef _OMPTASKS - unsigned int Nbins = patch->vecSpecies[ispec1]->Nbins; - for (unsigned int ibin = 0; ibin < Nbins; ibin++){ - patch->vecSpecies[ispec1]->Multiphoton_Breit_Wheeler_process->new_pair_per_bin[ibin][k].initializeReserve( - patch->vecSpecies[ispec1]->getNbrOfParticles(), - *patch->vecSpecies[ispec1]->mBW_pair_species_[k]->particles - ); + for (unsigned int ibin = 0; ibin < s1.Nbins; ibin++){ + s1.Multiphoton_Breit_Wheeler_process->new_pair_per_bin[ibin][k].initializeReserve( s1.getNbrOfParticles(), *s1.mBW_pair_species_[k]->particles ); } #endif ispec2 = patch->vecSpecies.size() + 1; - } ispec2++ ; } // This means that one of the pair species has not been found if( ispec2 == patch->vecSpecies.size() ) { - ERROR_NAMELIST( "In Species `" << patch->vecSpecies[ispec1]->name_ << "`," - << " the pair species `" << patch->vecSpecies[ispec1]->mBW_pair_species_names_[k] + ERROR_NAMELIST( "In Species `" << s1.name_ << "`," + << " the pair species `" << s1.mBW_pair_species_names_[k] << "` does not exist.", LINK_NAMELIST + std::string("#species") ) } @@ -1322,19 +1291,20 @@ class SpeciesFactory patch->vecSpecies[i]->initParticles( params, patch, with_particles, vector_species[i]->particles ); patch->vecSpecies[i]->initOperators( params, patch ); } - + // Ionization for( unsigned int i=0; ivecSpecies.size(); i++ ) { - if( patch->vecSpecies[i]->Ionize ) { - patch->vecSpecies[i]->electron_species_index = vector_species[i]->electron_species_index; - patch->vecSpecies[i]->electron_species = patch->vecSpecies[patch->vecSpecies[i]->electron_species_index]; - patch->vecSpecies[i]->Ionize->new_electrons.initialize( 0, *patch->vecSpecies[i]->electron_species->particles ); + Species &s = *patch->vecSpecies[i]; + if( s.Ionize ) { + s.electron_species_index = vector_species[i]->electron_species_index; + s.electron_species = patch->vecSpecies[s.electron_species_index]; + s.Ionize->new_electrons.initialize( 0, *s.electron_species->particles ); #ifdef _OMPTASKS - unsigned int Nbins = patch->vecSpecies[i]->Nbins; - for (unsigned int ibin = 0 ; ibin < Nbins ; ibin++){ - patch->vecSpecies[i]->Ionize->new_electrons_per_bin[ibin].initialize( 0, *patch->vecSpecies[i]->electron_species->particles ); + for (unsigned int ibin = 0 ; ibin < s.Nbins ; ibin++){ + s.Ionize->new_electrons_per_bin[ibin].initialize( 0, *s.electron_species->particles ); } #endif + s.Ionize->save_ion_charge_ = vector_species[i]->Ionize->save_ion_charge_; } } diff --git a/src/Species/SpeciesV.cpp b/src/Species/SpeciesV.cpp index f4286b56f..a5df2ca15 100755 --- a/src/Species/SpeciesV.cpp +++ b/src/Species/SpeciesV.cpp @@ -1640,7 +1640,7 @@ void SpeciesV::computeParticleCellKeys( Params ¶ms ) } -void SpeciesV::importParticles( Params ¶ms, Patch *, Particles &source_particles, vector &localDiags ) +void SpeciesV::importParticles( Params ¶ms, Patch *, Particles &source_particles, vector &localDiags, double time_dual, Ionization *I ) { unsigned int npart = source_particles.size(), ncells=particles->first_index.size(); @@ -1649,7 +1649,12 @@ void SpeciesV::importParticles( Params ¶ms, Patch *, Particles &source_parti if( particles->tracked ) { dynamic_cast( localDiags[tracking_diagnostic] )->setIDs( source_particles ); } - + + // If there is a diagnostic for recording particle birth, then copy new particles to the buffer + if( birth_records_ ) { + birth_records_->update( source_particles, npart, time_dual, I ); + } + // compute cell keys of new parts vector src_cell_keys( npart, 0 ); vector src_count( ncells, 0 ); diff --git a/src/Species/SpeciesV.h b/src/Species/SpeciesV.h index 8a39bd88c..39dc45089 100755 --- a/src/Species/SpeciesV.h +++ b/src/Species/SpeciesV.h @@ -81,7 +81,7 @@ class SpeciesV : public Species } //! Method to import particles in this species while conserving the sorting among bins - void importParticles( Params &, Patch *, Particles &, std::vector & )override; + void importParticles( Params &, Patch *, Particles &, std::vector &, double, Ionization *I = nullptr )override; //! Method performing the merging of particles virtual void mergeParticles( double time_dual )override; diff --git a/src/Species/SpeciesVAdaptiveMixedSort.cpp b/src/Species/SpeciesVAdaptiveMixedSort.cpp index c12340369..cc809d8c3 100755 --- a/src/Species/SpeciesVAdaptiveMixedSort.cpp +++ b/src/Species/SpeciesVAdaptiveMixedSort.cpp @@ -135,13 +135,13 @@ void SpeciesVAdaptiveMixedSort::resizeCluster( Params ¶ms ) // // } -void SpeciesVAdaptiveMixedSort::importParticles( Params ¶ms, Patch *patch, Particles &source_particles, vector &localDiags ) +void SpeciesVAdaptiveMixedSort::importParticles( Params ¶ms, Patch *patch, Particles &source_particles, vector &localDiags, double time_dual, Ionization *I ) { if( vectorized_operators ) { - importParticles( params, patch, source_particles, localDiags ); + importParticles( params, patch, source_particles, localDiags, time_dual, I ); } else { - Species::importParticles( params, patch, source_particles, localDiags ); + Species::importParticles( params, patch, source_particles, localDiags, time_dual, I ); } } diff --git a/src/Species/SpeciesVAdaptiveMixedSort.h b/src/Species/SpeciesVAdaptiveMixedSort.h index bf52ffaa0..4a3658e75 100755 --- a/src/Species/SpeciesVAdaptiveMixedSort.h +++ b/src/Species/SpeciesVAdaptiveMixedSort.h @@ -50,7 +50,7 @@ class SpeciesVAdaptiveMixedSort : public SpeciesV // void computeParticleCellKeys( Params ¶ms ) override; //! Method to import particles in this species while conserving the sorting among bins - void importParticles( Params &, Patch *, Particles &, std::vector & )override; + void importParticles( Params &, Patch *, Particles &, std::vector &, double, Ionization *I = nullptr )override; private: diff --git a/src/Tools/H5.cpp b/src/Tools/H5.cpp index f426aacc3..6c76fc444 100644 --- a/src/Tools/H5.cpp +++ b/src/Tools/H5.cpp @@ -66,6 +66,7 @@ void H5::init( std::string file, unsigned access, MPI_Comm * comm, bool _raise ) H5Pset_dxpl_mpio( dxpl_, H5FD_MPIO_COLLECTIVE ); H5Pset_alloc_time( dcr_, H5D_ALLOC_TIME_EARLY ); } + H5Pset_fill_time( dcr_, H5D_FILL_TIME_NEVER ); // Check error if( H5Eget_num( H5E_DEFAULT ) > 0 ) { id_ = -1; @@ -107,7 +108,7 @@ H5::~H5() } -//! 1D +//! 1D without selection or chunks H5Space::H5Space( hsize_t size ) { dims_ = { size }; global_ = size; @@ -119,10 +120,11 @@ H5Space::H5Space( hsize_t size ) { } //! 1D -H5Space::H5Space( hsize_t size, hsize_t offset, hsize_t npoints, hsize_t chunk ) { +H5Space::H5Space( hsize_t size, hsize_t offset, hsize_t npoints, hsize_t chunk, bool extendable ) { dims_ = { size }; global_ = size; - sid_ = H5Screate_simple( 1, &size, NULL ); + const hsize_t unlimited[] = { H5S_UNLIMITED }; + sid_ = H5Screate_simple( 1, &size, extendable ? unlimited : NULL ); if( size <= 0 || npoints <= 0 ) { H5Sselect_none( sid_ ); } else { @@ -137,13 +139,21 @@ H5Space::H5Space( hsize_t size, hsize_t offset, hsize_t npoints, hsize_t chunk ) } //! ND -H5Space::H5Space( std::vector size, std::vector offset, std::vector npoints, std::vector chunk ) { +H5Space::H5Space( std::vector size, std::vector offset, std::vector npoints, std::vector chunk, std::vector extendable ) { dims_ = size; global_ = 1; for( unsigned int i=0; i maxsize( size.size() ); + for( size_t i = 0; i < size.size(); i++ ) { + maxsize[i] = extendable[i] ? H5S_UNLIMITED : size[i]; + } + sid_ = H5Screate_simple( size.size(), &size[0], &maxsize[0] ); + } if( global_ <= 0 ) { H5Sselect_none( sid_ ); } else if( ! offset.empty() || ! npoints.empty() ) { diff --git a/src/Tools/H5.h b/src/Tools/H5.h index 3c517e1c5..e844a96e2 100644 --- a/src/Tools/H5.h +++ b/src/Tools/H5.h @@ -43,10 +43,10 @@ class H5Space //! 1D H5Space( hsize_t size ); - H5Space( hsize_t size, hsize_t offset, hsize_t npoints, hsize_t chunk = 0 ); + H5Space( hsize_t size, hsize_t offset, hsize_t npoints, hsize_t chunk = 0, bool extendable = false ); //! ND - H5Space( std::vector size, std::vector offset = {}, std::vector npoints = {}, std::vector chunk = {} ); + H5Space( std::vector size, std::vector offset = {}, std::vector npoints = {}, std::vector chunk = {}, std::vector extendable = {} ); ~H5Space() { H5Sclose( sid_ ); @@ -378,6 +378,18 @@ class H5Write : public H5 } } + // extend an open dataset (1D only) if it is extendable + void extend( hsize_t size ) + { + H5Dset_extent( id_, &size ); + } + + // extend an open dataset (ND) if it is extendable + void extend( std::vector size ) + { + H5Dset_extent( id_, &size[0] ); + } + //! Write a multi-dimensional array of uints H5Write array( std::string name, unsigned int &v, H5Space *filespace, H5Space *memspace, bool independent = false, hid_t file_type = -1 ) { diff --git a/src/Tools/Tools.cpp b/src/Tools/Tools.cpp index afb0e9ea1..12a7b8286 100755 --- a/src/Tools/Tools.cpp +++ b/src/Tools/Tools.cpp @@ -10,12 +10,8 @@ void Tools::printMemFootPrint( std::string tag ) { - - int val[4]; - char filename[80]; char sbuf[1024]; - char *S; //long lmem; pid_t numpro; @@ -31,26 +27,19 @@ void Tools::printMemFootPrint( std::string tag ) } // Peak resident set size - S=strstr( sbuf, "VmRSS:" )+6; - val[1] = ( int )atoi( S ); + int VmRSS = atoi( strstr( sbuf, "VmRSS:" )+6 ); // Peak virtual memory usage - S=strstr( sbuf, "VmSize:" )+6; - val[2] = atoi( S ); + // int VmSize = atoi( strstr( sbuf, "VmSize:" )+7 ); std::cout << "=== Mem usage === " << std::setw( 20 ) << tag << "\t=== " << std::setw( 6 ) - << "\t VmRSS << " << ( int )( ( double )val[1]/1024. ) << " Mb" << std::endl; + << "\t VmRSS << " << ( int )( ( double )VmRSS/1024. ) << " Mb" << std::endl; } double Tools::getMemFootPrint(int type_of_memory) { - - int val[4]; - char filename[80]; char sbuf[1024]; - char *S; - //long lmem; pid_t numpro; numpro = getpid(); @@ -70,17 +59,15 @@ double Tools::getMemFootPrint(int type_of_memory) } // Peak resident set size + int Vm = 0; if (type_of_memory == 0){ - S=strstr( sbuf, "VmRSS:" )+6; + Vm = atoi( strstr( sbuf, "VmRSS:" )+6 ); } else if (type_of_memory == 1){ - S=strstr( sbuf, "VmHWM:" )+6; + Vm = atoi( strstr( sbuf, "VmHWM:" )+6 ); } - - val[1] = ( int )atoi( S ); // Return RSS in GB - return ( double )val[1]/1024./1024.; - + return ( double )Vm/1024./1024.; } diff --git a/validation/analyses/validate_tst1d_05_tunnel_ionisation.py b/validation/analyses/validate_tst1d_05_tunnel_ionisation.py index 9e6f5143a..19148bd74 100755 --- a/validation/analyses/validate_tst1d_05_tunnel_ionisation.py +++ b/validation/analyses/validate_tst1d_05_tunnel_ionisation.py @@ -1,6 +1,8 @@ import os, re, numpy as np, math -from scipy.interpolate import interp1d as interp +from scipy.interpolate import interp1d from scipy.ndimage import gaussian_filter, maximum_filter1d +from h5py import File +from glob import glob import happi S = happi.Open(["./restart*"], verbose=False) @@ -57,42 +59,47 @@ def calculate_ionization(Ip, l): return times, Zstar # hydrogen -charge = S.ParticleBinning.Diag0().get() -charge_distribution = np.array( charge["data"] ) -charge_distribution /= charge_distribution[0,:].sum() -n1, n2 = charge_distribution.shape -mean_charge = (charge_distribution * np.outer(np.ones((n1,)), np.arange(n2))).sum(axis=1) +charge_distribution = S.ParticleBinning.Diag0().getData() +charge_distribution /= charge_distribution[0].sum() +n = charge_distribution[0].size +mean_charge = [np.sum(d*np.arange(n)) for d in charge_distribution] # # theory # Ip = np.array([13.5984])/27.2114 # l = np.array([0]) # t, Zs = calculate_ionization(Ip, l) # times = charge["times"]*S.namelist.Main.timestep -# Zs_theory = interp(t, Zs) (times) -Validate("Hydrogen mean charge vs time", mean_charge, 0.1) +# Zs_theory = interp1d(t, Zs) (times) +Validate("Hydrogen mean charge vs time", mean_charge, 0.03) -# carbon (does not work yet) -charge = S.ParticleBinning.Diag1().get() -charge_distribution = np.array( charge["data"] ) -charge_distribution /= charge_distribution[0,:].sum() -n1, n2 = charge_distribution.shape -mean_charge = (charge_distribution * np.outer(np.ones((n1,)), np.arange(n2))).sum(axis=1) +# carbon +charge_distribution = S.ParticleBinning.Diag1().getData() +charge_distribution /= charge_distribution[0].sum() +n = charge_distribution[0].size +mean_charge = [np.sum(d*np.arange(n)) for d in charge_distribution] # # theory # Ip = np.array([11.2602,24.3845,47.8877,64.4935,392.0905,489.9931]) /27.2114 # l = np.array([1,1,0,0,0,0]) # t, Zs = calculate_ionization(Ip, l) # times = charge["times"]*S.namelist.Main.timestep -# Zs_theory = interp(t, Zs) (times) -Validate("Carbon mean charge vs time", mean_charge, 0.1) +# Zs_theory = interp1d(t, Zs) (times) +Validate("Carbon mean charge vs time", mean_charge, 0.03) # SCALARS RELATED TO SPECIES Validate("Scalar Dens_electron", S.Scalar.Dens_electron().getData(), 0.003) Validate("Scalar Ntot_electron", S.Scalar.Ntot_electron().getData(), 100.) Validate("Scalar Zavg_carbon" , S.Scalar.Zavg_carbon ().getData(), 0.2) -# # TRACKING DIAGNOSTIC +# TRACKING DIAGNOSTIC d = S.TrackParticles("electron", axes=["Id","x","Wy"], timesteps=1000).getData() keep = d["Id"] > 0 order = np.argsort(d["x"][keep]) Validate("Track electron x", d["x"][keep][order][::200], 1e-4) Validate("Track electron Wy", gaussian_filter(maximum_filter1d(d["Wy"][keep][order],20),200)[::200], 1e-5) +# NEW PARTICLES DIAGNOSTIC +d = S.NewParticles.electron().get() +t = d["t"] +q = d["q"] +Validate("DiagNewParticles: number of particles", t.size, 5. ) +tavg = [np.mean(t[q==i]) for i in [0,1,2,3]] +Validate("DiagNewParticles: time vs ionization state", tavg, 0.01 ) diff --git a/validation/references/tst1d_05_tunnel_ionisation.py.txt b/validation/references/tst1d_05_tunnel_ionisation.py.txt index bc08a2342..f42d7e3ea 100755 Binary files a/validation/references/tst1d_05_tunnel_ionisation.py.txt and b/validation/references/tst1d_05_tunnel_ionisation.py.txt differ