diff --git a/doc/Sphinx/Overview/highlights.rst b/doc/Sphinx/Overview/highlights.rst index 1e755f6e2..c49c46004 100755 --- a/doc/Sphinx/Overview/highlights.rst +++ b/doc/Sphinx/Overview/highlights.rst @@ -88,7 +88,7 @@ Following is a comparison of the accelerated electron spectra at the end of thes In the green line it is shown the result of the previously known averaged ionization model. Without the longitudinal momentum initialization, few electrons obtained through ionization are trapped and accelerated in the plasma wave. The red line shows the result with the new averaged ionization model implemented in :program:`Smilei`, which accurately reproduces -the spectrum obtained with the simulation without an envelope model (blue line). +the spectrum obtained with the *standard laser* simulation without an envelope model. The envelope simulation required an amount of computing resources orders of magnitude smaller than those required by the simulation without a @@ -145,7 +145,7 @@ an amount of resources greater by at least an order of magnitude. The laser (in red) propagates through a low density plasma and drives a nonlinear Langmuir wave (in blue) that propagates at a velocity close to that of light in its wake. In this simulation, a moving window is used -so we can follow the laser as it propagates through the plasma. We see electrons (in white) being self-injected +so we can follow the laser as it propagates through the plasma. We see electrons being self-injected in this wakefield where they see a strong electric field that accelerates them up to ultra-relativistic (GeV) energy level. An animation generated from the simulation data can be found `here `_ @@ -230,7 +230,7 @@ shows how this balancing reduces the time of the simulation. :align: center The red curve is the best situation obtained in the previous section, while -the black curve corresponds to the DLB algorithm enabled. +the curves with dynamic load balancing are indicated as *DLB*. The portion of the box belonging to each MPI process varies when the load balancing occurs. The following figure shows how each of these portions evolve with time. @@ -239,7 +239,7 @@ occurs. The following figure shows how each of these portions evolve with time. The four panels correspond to four timesteps during the simulation. The colorscale represents the log-scaled load of each patch. -The black lines show the borders of each MPI process' portion of the box. +The lines show the borders of each MPI process' portion of the box. The MPI processes that are close to the hotspot tend to handle a smaller portion of the box. @@ -266,7 +266,7 @@ with a total of :math:`\sim 1.4` billion quasi-particles in the box. The following figure (top panel) shows half of the simulation box in the y-direction, and the laser field is reported at three different times. The reflected laser pulse (at time :math:`t_2`) shows a different spectral content than -the incident pulse (at time :math:`t_0`). The plasma electron density is shown in black. +the incident pulse (at time :math:`t_0`). A close-up view of the interaction region is given in the bottom panel, illustrating the electron bunches being pulled out from the plasma surface. @@ -317,8 +317,8 @@ A 2-dimensional simulation, in conditions close to actual experiments, ran on a box size of 1024 µm x 512 µm for 10 ps with 25 billion quasi-particles. The following figure shows the evolution of the pump and seed intensities in the head-on collision at three different times. -The blue-yellow maps correspond to the plasma density while the white-red maps -correspond to the lasers intensity. +Surfaces corresponding to the laser intensity are overlaid on top of surfaces +corresponding to the plasma density. .. image:: /_static/pump_seed.jpg :align: center diff --git a/doc/Sphinx/Overview/material.rst b/doc/Sphinx/Overview/material.rst index c05384fc8..45d06b92b 100644 --- a/doc/Sphinx/Overview/material.rst +++ b/doc/Sphinx/Overview/material.rst @@ -48,7 +48,18 @@ As of November 2021, 90 papers have been published covering a broad range of top Use the python script doc/doi2publications.py to generate entries from a DOI number, and paste them here +.. [Gorlova2024] + D. A. Gorlova, I. N. Tsymbalov, I. P. Tsygvintsev and A. B. Savelev, + `THz transition radiation of electron bunches laser-accelerated in long-scale near-critical-density plasmas`, + `Laser Physics Letters 21, 035001 (2024) `_ + +.. [Seidel2024] + + A. Seidel, B. Lei, C. Zepter, M. C. Kaluza, A. Sävert, M. Zepf, and D. Seipt, + `Polarization and CEP dependence of the transverse phase space in laser driven accelerators`, + `Physical Review Research 6, 013056 (2024) `_ + .. [Gao2023b] X. Gao, diff --git a/doc/Sphinx/Overview/partners.rst b/doc/Sphinx/Overview/partners.rst index 2170321ba..69b87e746 100755 --- a/doc/Sphinx/Overview/partners.rst +++ b/doc/Sphinx/Overview/partners.rst @@ -56,7 +56,6 @@ Partners | | * `Julien Dérouillat `_ | | | * `Haithem Kallala `_ | | | * `Mathieu Lobet `_ | -| | * `Francesco Massimo `_ | | | * `Charles Prouveur `_ | | | | +------------+---------------------------------------------------------------------------------------------------------+ diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 98e30cde6..1b0bd9ba7 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -25,13 +25,21 @@ Changes made in the repository (not released) * Happi: + * In ``Scalar``, it is now possible to make an operation on scalars such as ``"Uelm+Ukin"``. + The list of available scalars can be obtained from ``getScalars()``. * New arguments ``xoffset`` and ``yoffset`` to shift plot coordinates. + * New argument ``timestep_indices`` as an alternative to ``timesteps``. * Changed coordinate reference for 2D probe in 3D or AM geometry (zero is the box origin projected orthogonally on the probe plane). +* Documentation: + + * Dark theme (click the switch on the bottom left, or set browser preferences). + * Bug fixes: * ``dump_minutes`` often failed to write some checkpoint files. + * ``"auto"`` limits in ``ParticleBinning`` could fail with only one side on ``"auto"``. ---- diff --git a/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst b/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst index 8cfdab1e6..deb171eca 100644 --- a/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst +++ b/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst @@ -320,7 +320,7 @@ Cancellation on axis """"""""""""""""""""""" The first basic principle is that a mode 0 field defined on axis can only be longitudinal otherwise it would be ill defined. -On the opposite, longitudinal fields on axis can only be of mode 0 since they do not depend on :math:`theta`. +On the opposite, longitudinal fields on axis can only be of mode 0 since they do not depend on :math:`\theta`. From this we can already state that :math:`E_r^{m=0},\ E_t^{m=0},\ B_r^{m=0},\ B_t^{m=0},\ E_l^{m>0},\ B_l^{m>0}` are zero on axis. @@ -415,7 +415,12 @@ To ensure this derivative cancels on axis we simply pick: And equation :eq:`transverse_on_axis` then gives .. math:: - \tilde{E_{\theta}}^{m=1}[2] = -i\tilde{E_r}^{m=1}[3] + \tilde{E_{\theta}}^{m=1}(r=0) = -i\tilde{E_r}^{m=1}(r=0) + +With a finite difference scheme, this is implemented as + +.. math:: + \tilde{E_{\theta}}^{m=1}[2] = -\frac{i}{8}(9\tilde{E_r}^{m=1}[3]-\tilde{E_r}^{m=1}[4]) All the equation derived here are also valid for the magnetic field. But because of a different duality, it is more convenient to use a different approach. @@ -450,7 +455,7 @@ With a similar interpolation we obtain the boundary condition on axis for :math: .. math:: B_{\theta}^{m=1}[2]=-2iB_{r}^{m=1}[2]-B_{\theta}^{m=1}[3] -Longitudinal fields on axis +Longitudinal field on axis """"""""""""""""""""""""""""" We have alreayd established that only modes :math:`m=0` of longitudinal fields are non zero on axis. @@ -468,11 +473,6 @@ Introducing this result in the standard FDTD expression of :math:`E_l` we get: Again, the :math:`n` indice indicates the time step here. -:math:`B_l^{m=0}` is independant of :math:`\theta`. If we assume it is differentiable at :math:`r=0` then its derivative along :math:`r` is zero -on axis (derivative of a pair function is zero at :math:`x=0` ). From this we get: - -.. math:: - B_{l}^{m=0}[2]=B_{l}^{m=0}[3] Below axis """""""""""""""""""""""""""" @@ -481,7 +481,8 @@ Fields "below" axis are primal fields data with indice :math:`j<2` and dual fiel These fields are not physical in the sense that they do not contribute to the reconstruction of any physical field in real space and are not obtained by solving Maxwell's equations. Nevertheless, it is numerically convenient to give them a value in order to facilitate field interpolation for macro-particles near axis. This is already what is done for dual fields in :math:`r` which cancel on axis for instance. -We extend this logic to primal fields in :math:`r`: +We extend this logic to primal fields in :math:`r` as well. +For any given field :math:`F`, the symetric of :math:`F` with respect to the axis is :math:`F` if :math:`F` is non zero on axis and :math:`-F` if :math:`F` is zero on axis: .. math:: @@ -489,19 +490,35 @@ We extend this logic to primal fields in :math:`r`: E_{l}^{m>0}[1] = -E_{l}^{m>0}[3] + E_{r}^{m\neq1}[2] = -E_{r}^{m\neq1}[3] + + E_{r}^{m=1}[2] = E_{r}^{m=1}[3] + E_{\theta}^{m\neq1}[1] = -E_{\theta}^{m\neq1}[3] E_{\theta}^{m=1}[1] = E_{\theta}^{m=1}[3] + B_{l}^{m=0}[2]=B_{l}^{m=0}[3] + + B_{l}^{m>0}[2]=-B_{l}^{m>0}[3] + B_{r}^{m\neq1}[1] = -B_{r}^{m\neq1}[3] B_{r}^{m=1}[1] = B_{r}^{m=1}[3] + B_{t}^{m\neq1}[2] = -B_{t}^{m\neq1}[3] + + B_{t}^{m=1}[2] = B_{t}^{m=1}[3] + Currents near axis """"""""""""""""""""" A specific treatment must be applied to charge and current densities near axis because the projector deposits charge and current "below" axis. -Quantities below axis must be brought back in the "physical" terms on and above axis. +Quantities below axis must be folded back onto their symetric "above" axis. +Mode 0 contributions below axis are added to their "above axis" counterparts. +On the opposite, strictly positive modes contributions are deduced. +This ensures that a particle sitting on axis will have a net contribution to the domain only in mode 0 as expected because in that case :math:`\theta` is not defined and therefore +the deposition can not be a function of :math:`\theta`. Using the continuity equation instead of Gauss law for transverse current of mode :math:`m=1` on axis, we can derive the exact same boundary conditions on axis for current density as for the electric field. diff --git a/doc/Sphinx/Use/installation.rst b/doc/Sphinx/Use/installation.rst index bf83d58a0..a2f872e36 100755 --- a/doc/Sphinx/Use/installation.rst +++ b/doc/Sphinx/Use/installation.rst @@ -105,9 +105,6 @@ Advanced compilation options make config=vtune # For Intel Vtune make config=inspector # For Intel Inspector make config=detailed_timers # More detailed timers, but somewhat slower execution - make config=omptasks # use OpenMP task parallelization, not supported by old compilers - make config=part_event_tracing_tasks_off # trace the use particle operators, without task parallelization - make config=part_event_tracing_tasks_on # trace the use particle operators, with OpenMP task parallelization make config="gpu_nvidia noopenmp" # For Nvidia GPU acceleration make config="gpu_amd" # For AMD GPU acceleration diff --git a/doc/Sphinx/Use/post-processing.rst b/doc/Sphinx/Use/post-processing.rst index 78f5083c7..8f90f56f6 100755 --- a/doc/Sphinx/Use/post-processing.rst +++ b/doc/Sphinx/Use/post-processing.rst @@ -101,12 +101,16 @@ about the corresponding diagnostics in the simulation. Returns a list of available diagnostics of the given type * ``diagType``: The diagnostic type (``"Field"``, ``"Probe"``, etc.) + +.. rubric:: Information on specific diagnostics + +.. py:method:: getScalars() + + Returns a list of available scalars. .. py:method:: getTrackSpecies() Returns a list of available tracked species. - -.. rubric:: Information on specific diagnostics .. py:method:: fieldInfo(diag) @@ -145,12 +149,16 @@ Open a Scalar diagnostic .. py:method:: Scalar(scalar=None, timesteps=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``scalar``: The name of the scalar. - | If not given, then a list of available scalars is printed. - * ``timesteps``: The requested timestep(s). - | If omitted, all timesteps are used. - | If one number given, the nearest timestep available is used. - | If two numbers given, all the timesteps in between are used. + * ``scalar``: The name of the scalar, or an operation on scalars, such as ``"Uelm+Ukin"``. + * ``timesteps`` or ``timestep_indices``: The requested range of timesteps. + + * If omitted, all timesteps are used. + * If one number given, the nearest timestep available is used. + * If two numbers given, all the timesteps in between are used. + + When using ``timesteps``, provide the timesteps themselves, but + when using ``timestep_indices``, provide their indices in the list + of the available timesteps. * ``units``: A unit specification (see :ref:`units`) * ``data_log``: | If ``True``, then :math:`\log_{10}` is applied to the output. @@ -170,7 +178,7 @@ Open a Field diagnostic .. py:method:: Field(diagNumber=None, field=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, moving=False, export_dir=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``: same as before. * ``diagNumber``: number or ``name`` of the fields diagnostic | If not given, then a list of available diagnostic numbers is printed. * ``field``: The name of a field (``"Ex"``, ``"Ey"``, etc.) @@ -228,7 +236,7 @@ Open a Probe diagnostic .. py:method:: Probe(probeNumber=None, field=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``probeNumber``: number or ``name`` of the probe (the first one has number 0). | If not given, a list of available probes is printed. * ``field``: name of the field (``"Bx"``, ``"By"``, ``"Bz"``, ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Jx"``, ``"Jy"``, ``"Jz"`` or ``"Rho"``). @@ -256,7 +264,7 @@ Open a ParticleBinning diagnostic .. py:method:: ParticleBinning(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``: number or ``name`` of the particle binning diagnostic (starts at 0). | If not given, a list of available diagnostics is printed. | It can also be an operation between several diagnostics. @@ -306,7 +314,7 @@ Open a Screen diagnostic .. py:method:: Screen(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``, ``subset`` and ``average``: identical to that of ParticleBinning diagnostics. * See also :ref:`otherkwargs` @@ -323,7 +331,7 @@ Open a RadiationSpectrum diagnostic .. py:method:: ParticleBinning(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``, ``subset`` and ``average``: identical to that of ParticleBinning diagnostics. * See also :ref:`otherkwargs` @@ -344,7 +352,7 @@ Open a TrackParticles diagnostic .. py:method:: TrackParticles(species=None, select="", axes=[], timesteps=None, sort=True, length=None, units=[""], **kwargs) - * ``timesteps``, ``units``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``export_dir``: same as before. * ``species``: the name of a tracked-particle species. If omitted, a list of available tracked-particle species is printed. * ``select``: Instructions for selecting particles among those available. diff --git a/doc/Sphinx/_static/LWFA_Plas@Par.jpg b/doc/Sphinx/_static/LWFA_Plas@Par.jpg index 002204fcd..3360fbdf1 100644 Binary files a/doc/Sphinx/_static/LWFA_Plas@Par.jpg and b/doc/Sphinx/_static/LWFA_Plas@Par.jpg differ diff --git a/doc/Sphinx/_static/PWFA.jpg b/doc/Sphinx/_static/PWFA.jpg index 88f0c5130..58de8a641 100644 Binary files a/doc/Sphinx/_static/PWFA.jpg and b/doc/Sphinx/_static/PWFA.jpg differ diff --git a/doc/Sphinx/_static/smileiLogo.svg b/doc/Sphinx/_static/smileiLogo.svg index 1b70b4a5c..2da3028e4 100755 --- a/doc/Sphinx/_static/smileiLogo.svg +++ b/doc/Sphinx/_static/smileiLogo.svg @@ -1,14 +1,7 @@ - + viewBox="50 380 700 240"> {{ toc }} -
    +
      {%- else %}
    • {{ entry }} @@ -162,16 +162,15 @@ + style="opacity:1;fill:none;stroke:#ffffff;stroke:var(--header_text);stroke-width:10;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" /> + style="opacity:1;fill:#ffffff;fill:var(--header_text);fill-opacity:1;stroke:none;" /> + style="opacity:1;fill:none;stroke:#ffffff;stroke:var(--header_text);stroke-width:3;stroke-linecap:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1" /> @@ -180,21 +179,18 @@ xmlns="http://www.w3.org/2000/svg" viewBox="0 0 80 120"> - - + transform="translate(0,-932.36216)" + style="fill:none;stroke:#ffffff;stroke:var(--header_text);stroke-width:10;stroke-linecap:round;stroke-linejoin:miter;stroke-opacity:1;stroke-miterlimit:4;stroke-dasharray:none"> + + -
      +
      - + @@ -226,6 +222,9 @@ {%- block footer %}
      +
      + +
      @@ -261,6 +260,30 @@ for( var i=0; i div { display: inline-block; } +#themebutton { + height: 1.6em; width: 1.6em; border-radius: 0.8em; cursor: pointer; margin: 0 auto; + background: linear-gradient(90deg, var(--bg) 50%, var(--alt-bg) 50%) no-repeat; + border: 0.2em solid; border-color: var(--alt-bg); + margin: 0.5em; +} +#themebutton:hover { filter: opacity(0.8); } +#themebutton:active { filter: opacity(1.2); } + div.logo { - padding-right: 2em; - margin: 0; + height: 10mm; text-align: left; - height:{{ theme_head_height }}mm; + margin: 1mm 2em 1mm 1mm; + padding: 0 2mm; + border-radius: 2mm; + background-color: {{ theme_header_text }}; background-color: var(--header_text); + height:{{ theme_head_height| float -2 }}mm; } div.logo a, @@ -91,7 +142,12 @@ div.logo a:hover { } .version { - position: relative; + position: relative; + margin: 0; +} + +.version img { + margin: 0 9%; } .version::before { @@ -101,18 +157,18 @@ div.logo a:hover { bottom: 5%; right: 5%; font: bold 1.6em BPreplay; - color: {{ theme_main_color_bold }};; + color:{{ theme_main_color_bold }}; color: var(--main_bold); } /* -- body styles ----------------------------------------------------------- */ a { - color: {{ theme_main_color_bold }}; + color:{{ theme_bold_text }}; color: var(--bold_text); text-decoration: none; } a:hover { - color: {{ theme_main_color_med }}; + color:{{ theme_main_color_med }}; color: var(--main_med); text-decoration: underline; } @@ -126,13 +182,13 @@ div.body h6 { font-weight: bold; margin: 10px 0px 10px 0px; padding: 0; - color: {{ theme_main_color_bold }}; + color:{{ theme_bold_text }}; color: var(--bold_text); } div.body h1 { font-size: 240%; padding: 0.1em 0.3em 0.1em 1em; - background-color:{{ theme_main_color_pale }}; + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); display: inline-block; border-radius:0.1em } @@ -150,6 +206,7 @@ div.body h6 { font-size: 100%; } h1:target, h2:target, h3:target, dt:target { background: {{ theme_main_color_pale }} content-box !important; + background-color: var(--main_pale); } a.headerlink { @@ -216,11 +273,11 @@ div.admonition p.last { } div.highlight { - background-color: white; + background-color: #fff; background-color: var(--bg); } .highlight { - background: {{ theme_main_color_pale }}; + background:{{ theme_main_color_pale }}; background: var(--main_pale); } div.note { @@ -283,7 +340,7 @@ table.fancy { position: relative; z-index:0; overflow-x: auto; - border:1px solid {{ theme_main_color_med }}; + border:1px solid {{ theme_main_color_med }}; border:1px solid var(--main_med); white-space: nowrap; } @@ -292,7 +349,7 @@ table.fancy * { } table.fancy > tbody > tr.row-odd { - background-color: {{ theme_main_color_pale }}; + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); } table.fancy > tbody > tr > td { @@ -391,7 +448,7 @@ li.toctree-l1 > ul { } pre { - background: {{ theme_pre_bg }}; + background: rgba(120,120,120,0.15); padding: 7px 1em; margin: 15px 0px; line-height: 1.3em; @@ -402,13 +459,13 @@ div.highlight > pre { } code { - background-color: {{ theme_pre_bg }}; - color: #222; + background-color: rgba(160,160,160,0.2); + color: {{ theme_body_text }}; color: var(--text) word-break: break-word; } code.xref { - background-color: {{ theme_pre_bg }}; + background-color: rgba(160,160,160,0.2); border-bottom: 1px solid transparent; } @@ -417,18 +474,18 @@ a.reference { } a.reference:hover { - border-bottom: 1px solid {{ theme_main_color_med }}; + border-bottom: 1px solid {{ theme_main_color_med }}; border-bottom: 1px solid var(--main_med); } a.footnote-reference { text-decoration: none; font-size: 0.7em; vertical-align: top; - border-bottom: 1px dotted {{ theme_main_color_bold }}; + border-bottom: 1px dotted {{ theme_main_color_bold }}; border-bottom: 1px dotted var(--main_bold); } a.footnote-reference:hover { - border-bottom: 1px solid {{ theme_main_color_med }}; + border-bottom: 1px solid {{ theme_main_color_med }}; border-bottom: 1px solid var(--main_med); } a:hover code { @@ -441,6 +498,10 @@ a:hover code { display: none!important; } +dl.py > dt { + overflow-x: scroll; +} + /* Make nested-list/multi-paragraph items look better in Releases changelog * pages. Without this, docutils' magical list fuckery causes inconsistent * formatting between different release sub-lists. @@ -472,7 +533,7 @@ div.headercolor { top:0; height: {{ theme_head_height }}mm; width: 100%; - background-color:{{ theme_main_color_bold }}; + background-color:{{ theme_header_color }}; background-color: var(--header_color); } div.hpositioner { @@ -486,7 +547,6 @@ div.header { position: relative; height: 100%; font-size: 1em; - color: white; display: inline-block; text-align: left; } @@ -521,52 +581,48 @@ div.header > * { right:0; cursor:pointer; text-align: center; - background-color: #eee; - color: #888; - border-right: 2px solid #aaa; - border-left: 2px solid #aaa; - border-bottom: 2px solid #aaa; - border-bottom-left-radius: 4px; - border-bottom-right-radius: 4px; + background-color: {{ theme_nav_bg }}; background-color: var(--nav_bg); + color: {{ theme_body_text }}; color: var(--text); margin-right:12px; padding:0.2em 1em; font-size:80%; font-weight:bold; z-index:10; + border-left: 2px solid rgba(128,128,128,0.6); + border-right: 2px solid rgba(128,128,128,0.6); + border-bottom: 2px solid rgba(128,128,128,0.6); + border-bottom-left-radius: 4px; + border-bottom-right-radius: 4px; } #nav_button.pushed { display:none; } #nav_list { display:none; - background-color: transparent; + background-color: {{ theme_nav_bg }}; background-color: var(--nav_bg); position: relative; + border-left: 2px solid rgba(128,128,128,0.6); + border-right: 2px solid rgba(128,128,128,0.6); + border-bottom: 2px solid rgba(128,128,128,0.6); + border-bottom-left-radius: 4px; + border-bottom-right-radius: 4px; } #nav_list > ul a { padding:0.3em 1em; } #nav_title { line-height: {{ theme_head_height |float *0.8 }}mm; - border-left: 2px solid #aaa; - border-right: 2px solid #aaa; - background-color: #eee; padding: 0 1em; font-weight: bold; font-size:120%; position: relative; } #nav_title > a { - color: #444; + color: {{ theme_body_text }}; color: var(--text); } #nav_list > ul { position: relative; right:0; - background-color: #eee; - border-left: 2px solid #aaa; - border-right: 2px solid #aaa; - border-bottom: 2px solid #aaa; - border-bottom-left-radius: 4px; - border-bottom-right-radius: 4px; } .menu { @@ -578,7 +634,7 @@ div.header > * { .menu > div > a, .menu > div > a:hover { - color: white; + color: {{ theme_header_text }}; color: var(--header_text); text-decoration: none; } @@ -611,13 +667,13 @@ div.header > * { } .menuButton.pushed { - background-color: {{ theme_main_color_med }}; + background-color:{{ theme_header_pushed }}; background-color: var(--header_pushed); } .menu div.toctree-smilei { - border-left: 1px solid {{ theme_main_color_bold }}; - border-right: 1px solid {{ theme_main_color_bold }}; - border-bottom: 1px solid {{ theme_main_color_bold }}; + border-left: 1px solid {{ theme_main_color_bold }}; border-left: 1px solid var(--main_bold); + border-right: 1px solid {{ theme_main_color_bold }}; border-right: 1px solid var(--main_bold); + border-bottom: 1px solid {{ theme_main_color_bold }}; border-bottom: 1px solid var(--main_bold); border-bottom-left-radius: 4px; border-bottom-right-radius: 4px; } @@ -647,11 +703,11 @@ div.header > * { } #smallScreenMenuButton.pushed { - background-color: {{ theme_main_color_pale }}; + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); } #smallScreenMenuButton.pushed svg g { - fill:{{ theme_main_color_bold }} !important; + fill:{{ theme_main_color_bold }} !important; fill:var(--main_bold) !important; } #smallScreenMenu { @@ -662,7 +718,7 @@ div.header > * { left:0; padding:0.3em 0 0 0; overflow-y: auto; - background-color:{{ theme_main_color_pale }}; + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); border: 0; z-index:500; } @@ -730,15 +786,18 @@ ul.experimental p:first-of-type::after /* --------------------------------------------------------------------*/ /* Logos and icons */ img.logo { - height:{{ theme_head_height }}mm; + height:{{ theme_head_height | float -2 }}mm; } a.icon { display:inline-block; text-align:center; margin:0.3em; + padding: 0.5em 0; width:100px; text-decoration:none; - color:{{ theme_body_text }}; + color:{{ theme_body_text }}; color: var(--text); + background-color: rgba(180,180,180,0.2); + border-radius: 1em; } a.icon * { margin:0 auto; @@ -748,11 +807,11 @@ a.icon img { height:3em; } a.icon:hover { - color:{{ theme_body_text }}; - background-color:{{ theme_main_color_pale }}; + color:{{ theme_body_text }}; color: var(--text); + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); } a.icon:active { - background-color:{{ theme_main_color_med }}; + background-color:{{ theme_main_color_med }}; background-color: var(--main_med); } /* --------------------------------------------------------------------*/ @@ -771,14 +830,14 @@ div.lablogos > div { /* --------------------------------------------------------------------*/ /* table of contents */ div.toctree-smilei { - background-color: {{ theme_main_color_pale }}; + background-color:{{ theme_main_color_pale }}; background-color: var(--main_pale); font-weight:normal; } div.toctree-smilei a { text-decoration: none; border:0 !important; - color: black; + color: {{ theme_body_text }}; color: var(--text); display:block; } div.toctree-smilei > ul { @@ -790,17 +849,17 @@ div.toctree-smilei > ul > li { padding: 0; } div.toctree-smilei > ul > li > a { - background-color: {{ theme_main_color_bold }}; - color: {{ theme_main_color_pale }}; + background-color:{{ theme_bold_text }}; background-color: var(--bold_text); + color:{{ theme_main_color_pale }}; color: var(--main_pale); padding: 0.4em 0 0.4em 1em; } div.toctree-smilei li.outer > a { background-color: transparent; - color: {{ theme_main_color_bold }}; + color:{{ theme_bold_text }}; color: var(--bold_text); } div.toctree-smilei li.outer > a:hover { - background-color: {{ theme_main_color_med }}; - color: {{ theme_main_color_pale }}; + background-color:{{ theme_main_color_med }}; background-color: var(--main_med); + color:{{ theme_pale_text }}; color: var(--pale_text); } div.toctree-smilei > ul > li > ul { margin: 0; @@ -810,7 +869,7 @@ div.toctree-smilei > ul > li > ul > li > a { padding:0.3em 0 0.3em 2.5em; } div.toctree-smilei > ul > li > ul > li > a:hover { - background-color: #ddd; + background-color: rgba(128,128,128,0.3); } div.toctree-smilei > ul > li > ul > li > ul { @@ -873,11 +932,11 @@ dl.citation p { hr { margin-top:1em; border: 0; - border-top: 2px solid {{ theme_main_color_bold }}; + border-top: 2px solid {{ theme_bold_text }}; border-top: 2px solid var(--bold_text); } div.toctree-smilei hr { margin: 0.2em 0; - background: {{ theme_main_color_bold }}; + background:{{ theme_main_color_bold }}; background: var(--main_bold); width:100%; } @@ -940,8 +999,8 @@ div.toctree-smilei hr { margin:{{ theme_head_height|float * 0.1 }}mm 0; border-radius:0.5em; border:0; - color: #555; - background-color: white; + color: #555; color: var(--text); + background-color: white; background-color: var(--bg); } #search-results { padding-top: 5px; @@ -1017,7 +1076,7 @@ div.toctree-smilei hr { } #searchbox { display:block !important; - left:{{ theme_head_height|float * 2.8 }}mm; + left:{{ theme_head_height|float * 3 }}mm; right:{{ theme_head_height|float * 1.9 }}mm; } #searchinput { diff --git a/doc/Sphinx/smilei_theme/theme.conf b/doc/Sphinx/smilei_theme/theme.conf index 048723647..648f66bc2 100755 --- a/doc/Sphinx/smilei_theme/theme.conf +++ b/doc/Sphinx/smilei_theme/theme.conf @@ -10,30 +10,46 @@ main_color_bold = #0672BA main_color_med = #409ddb main_color_pale = #EAF4FD keyword_color = #05baa8 +body_text = #3E4349 +bold_text = #0672BA +pale_text = #EAF4FD +header_color = #0672BA +header_pushed = #409ddb +header_text = white +nav_bg = #eee + +dark_main_color_bold = #0672BA +dark_main_color_med = #409ddb +dark_main_color_pale = #26313A +dark_keyword_color = #05baa8 +dark_body_text = #ccc +dark_bold_text = #7ECDFF +dark_pale_text = #26313A +dark_header_color = #003765 +dark_header_pushed = #0672BA +dark_header_text = #C5DCEA +dark_nav_bg = #252525 + +footer_text = #777 page_width = 210 head_height = 12 -body_text = #3E4349 -footer_text = #888 - anchor = #DDD anchor_hover_fg = #444 anchor_hover_bg = transparent -note_bg = #EEE -note_border = #CCC - -seealso_bg = #EEE -seealso_border = #CCC +note_bg = rgba(123,123,123,0.2) +note_border = rgba(123,123,123,0.6) -warn_bg = #EAF4FD -warn_border = #0672BA +seealso_bg = rgba(123,123,123,0.2) +seealso_border = rgba(123,123,123,0.6) -footnote_bg = #FDFDFD -footnote_border = #EEE +warn_bg = rgba(123,123,123,0.2) +warn_border = rgba(123,123,123,0.2) -pre_bg = #EEE +footnote_bg = rgba(123,123,123,0.2) +footnote_border = rgba(123,123,123,0.2) code_font_size = 0.9em code_font_family = 'Consolas', "BPmonoRound", 'Menlo', 'Deja Vu Sans Mono', 'Bitstream Vera Sans Mono', monospace diff --git a/happi/_Diagnostics/Diagnostic.py b/happi/_Diagnostics/Diagnostic.py index b156736bf..b17e4c449 100755 --- a/happi/_Diagnostics/Diagnostic.py +++ b/happi/_Diagnostics/Diagnostic.py @@ -26,7 +26,6 @@ def __init__(self, simulation, *args, **kwargs): self._log = [] self._data_log = False self._data_transform = None - self._error = [] self._xoffset = 0. # The 'simulation' is a SmileiSimulation object. It is passed as an instance attribute @@ -43,9 +42,7 @@ def __init__(self, simulation, *args, **kwargs): # Reload the simulation, in case it has been updated self.simulation.reload() - if not self.simulation.valid: - self._error += ["Invalid Smilei simulation"] - return + assert self.simulation.valid, "Invalid Smilei simulation" # Copy some parameters from the simulation self._results_path = self.simulation._results_path @@ -65,21 +62,20 @@ def __init__(self, simulation, *args, **kwargs): self.units = kwargs.pop("units", [""]) if type(self.units) in [list, tuple]: self.units = Units(*self.units , verbose = self._verbose) if type(self.units) is dict : self.units = Units(verbose = self._verbose, **self.units) - if type(self.units) is not Units: - self._error += ["Could not understand the 'units' argument"] - return + assert type(self.units) is Units, "Could not understand the 'units' argument" self.units._initRegistry(self._ureg) + # Manage the 'timesteps' or 'timestep_indices' + assert not ( "timesteps" in kwargs and "timestep_indices" in kwargs ),\ + "Cannot define both `timesteps` and `timestep_indices`" + # Call the '_init' function of the child class remaining_kwargs = self._init(*args, **kwargs) - if remaining_kwargs is not None and len(remaining_kwargs) > 0: - self.valid = False - self._error += ["The following keyword-arguments are unknown: "+", ".join(remaining_kwargs.keys())] - return + assert remaining_kwargs is None or len(remaining_kwargs) == 0,\ + "The following keyword-arguments are unknown: "+", ".join(remaining_kwargs.keys()) self.dim = len(self._shape) - if self.valid: - self._prepareUnits() + self._prepareUnits() # Prepare data_log output self._dataAtTime = self._dataLogAtTime if self._data_log else self._dataLinAtTime @@ -88,19 +84,6 @@ def __init__(self, simulation, *args, **kwargs): def __repr__(self): self.info() return "" - - # Method to verify everything was ok during initialization - def _validate(self): - try: - self.simulation.valid - except Exception as e: - print("No valid Smilei simulation selected") - return False - if not self.simulation.valid or not self.valid: - print("***ERROR*** - Diagnostic is invalid") - print("\n".join(self._error)) - return False - return True # Prepare units for axes def _prepareUnits(self): @@ -143,7 +126,7 @@ def limits(self): # Method to print info on this diag def info(self): - if self._validate() and self._verbose: + if self._verbose: print(self._info()) # Method to get only the arrays of data @@ -159,7 +142,6 @@ def getData(self, timestep=None): A list of arrays: each array corresponding to the diagnostic data at a given timestep. """ - if not self._validate(): return self._prepare1() # prepare the vfactor data = [] @@ -175,7 +157,6 @@ def getData(self, timestep=None): def getTimesteps(self): """Obtains the list of timesteps selected in this diagnostic""" - if not self._validate(): return [] return self._timesteps def getTimes(self): @@ -184,7 +165,6 @@ def getTimes(self): By default, times are in the code's units, but are converted to the diagnostic's units defined by the `units` argument, if provided. """ - if not self._validate(): return [] return self.units.tcoeff * self.timestep * self._np.array(self._timesteps) def _getCenters(self, axis_index, timestep): @@ -234,7 +214,6 @@ def get(self, timestep=None): !!! Deprecated !!! Use functions `getData`, `getTimesteps`, `getTimes` and `getAxis` instead. """ - if not self._validate(): return # obtain the data arrays data = self.getData(timestep=timestep) # format the results into a dictionary @@ -286,7 +265,6 @@ def plot(self, timestep=None, saveAs=None, axes=None, dpi=200, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).plot(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -334,7 +312,6 @@ def streak(self, saveAs=None, axes=None, dpi=200, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).streak(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -422,7 +399,6 @@ def animate(self, movie="", fps=15, dpi=200, saveAs=None, axes=None, **kwargs): This takes the particle binning diagnostic #1 and plots the resulting array in figure 1 from 0 to 3e14. """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -482,7 +458,6 @@ def slide(self, axes=None, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).slide(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return ax = self._make_axes(axes) @@ -516,20 +491,26 @@ def update(t): self.info() + # Method to select specific timesteps (or indices) among those available + def _selectTimesteps(self, timestep_selection, timestep_index_selection, timesteps): + if timestep_selection is not None: + ts = self._np.array(timestep_selection, ndmin=1, dtype=float) + if ts.size==2: # get all times in between bounds + return timesteps[ (timesteps>=ts[0]) * (timesteps<=ts[1]) ] + elif ts.size==1: # get nearest time + return self._np.array([timesteps[(self._np.abs(timesteps-ts)).argmin()]]) + else: + raise Exception("Argument `timesteps` must be one or two non-negative integers") + elif timestep_index_selection is not None: + ts = self._np.array(timestep_index_selection, ndmin=1, dtype=int) + if ts.size==2: + return timesteps[ ts[0]:ts[1] ] + elif ts.size==1: + return self._np.array(timesteps[ts[0]]) + else: + raise Exception("Argument `timestep_indices` must be one or two non-negative integers") + return timesteps - # Method to select specific timesteps among those available in times - def _selectTimesteps(self, timesteps, times): - ts = self._np.array(self._np.double(timesteps),ndmin=1) - if ts.size==2: - # get all times in between bounds - times = times[ (times>=ts[0]) * (times<=ts[1]) ] - elif ts.size==1: - # get nearest time - times = self._np.array([times[(self._np.abs(times-ts)).argmin()]]) - else: - raise - return times - # Method to select portion of a mesh based on a slice def _selectSubset(self, portion, meshpoints, axisname, axisunits, operation): try: @@ -877,7 +858,6 @@ def _dataLogAtTime(self, t): # Convert data to VTK format def toVTK(self, numberOfPieces=1): - if not self._validate(): return # prepare vfactor self._prepare1() diff --git a/happi/_Diagnostics/Field.py b/happi/_Diagnostics/Field.py index cf65a5faa..4b9a678c8 100755 --- a/happi/_Diagnostics/Field.py +++ b/happi/_Diagnostics/Field.py @@ -176,11 +176,8 @@ def fieldTranslator(f): for i,t in enumerate(self._timesteps): self._data.update({ t : i }) # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) # Need at least one timestep if self._timesteps.size < 1: @@ -316,7 +313,6 @@ def _getCenters(self, axis_index, timestep, h5item = None): # get the value of x_moved for a requested timestep def getXmoved(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -330,7 +326,6 @@ def getXmoved(self, t): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -372,7 +367,6 @@ def _getDataAtTime(self, t): # Method to obtain the data only # Specific to cylindrical geometry, to reconstruct the plane at some theta def _theta_getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -429,7 +423,6 @@ def _theta_getDataAtTime(self, t): # Specific to cylindrical geometry, to reconstruct a 3d box def _build3d_getDataAtTime(self, t): from scipy.interpolate import RegularGridInterpolator - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/NewParticles.py b/happi/_Diagnostics/NewParticles.py index 8bc6285bc..6b45abc32 100755 --- a/happi/_Diagnostics/NewParticles.py +++ b/happi/_Diagnostics/NewParticles.py @@ -133,13 +133,14 @@ def _init(self, species=None, select="", axes=[], chunksize=20000000, **kwargs): print("Kept "+str(self.nselectedParticles)+" particles") # Finish constructor + assert "timesteps" not in kwargs and "timestep_indices" not in kwargs,\ + "The NewParticles diagnostic does not accept 'timesteps' or 'timestep_indices' as arguments." 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 diff --git a/happi/_Diagnostics/ParticleBinning.py b/happi/_Diagnostics/ParticleBinning.py index 851442ada..bb4cfbf2d 100755 --- a/happi/_Diagnostics/ParticleBinning.py +++ b/happi/_Diagnostics/ParticleBinning.py @@ -95,13 +95,16 @@ def _init(self, diagNumber=None, timesteps=None, subset=None, average=None, sum= self._timesteps = {} self._alltimesteps = {} self._indexOfTime = {} + self._h5files = [] self._h5items = {} + timestep_indices = kwargs.pop("timestep_indices", None) for d in self._diags: # Gather data from all timesteps, and the list of timesteps items = {} for path in self._results_path: try: f = self._h5py.File(path+self._os.sep+self._diagType+str(d)+'.h5', 'r') + self._h5files += [f] except: continue items.update( dict(f) ) @@ -113,12 +116,8 @@ def _init(self, diagNumber=None, timesteps=None, subset=None, average=None, sum= self._indexOfTime.update({ d:{} }) for i,t in enumerate(self._timesteps[d]): self._indexOfTime[d].update({ t : i }) - # If timesteps is None, then keep all timesteps, otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps[d] = self._selectTimesteps(timesteps, self._timesteps[d]) - except Exception as e: - raise Exception("Argument 'timesteps' must be one or two non-negative integers") + # Select timesteps if requested + self._timesteps[d] = self._selectTimesteps(timesteps, timestep_indices, self._timesteps[d]) # Verify that timesteps are the same for all diagnostics if (self._timesteps[d] != self._timesteps[self._diags[0]]).any() : raise Exception( @@ -493,7 +492,6 @@ def getAvailableTimesteps(self, diagNumber=None): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Auto axes require recalculation of bin size and centers if self.auto_axes: self._updateAxes(t) diff --git a/happi/_Diagnostics/Performances.py b/happi/_Diagnostics/Performances.py index 460fe44e9..c7a8df254 100755 --- a/happi/_Diagnostics/Performances.py +++ b/happi/_Diagnostics/Performances.py @@ -191,16 +191,10 @@ def _init(self, raw=None, map=None, histogram=None, timesteps=None, data_log=Fal self._data = {} for i,t in enumerate(self._timesteps): self._data.update({ t : i }) - # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") - - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") + # Select timesteps if requested + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" # 3 - Manage axes # ------------------------------------------------------------------- @@ -270,7 +264,6 @@ def getAvailableQuantities(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/Probe.py b/happi/_Diagnostics/Probe.py index 1fd4958eb..3c6b2c73d 100755 --- a/happi/_Diagnostics/Probe.py +++ b/happi/_Diagnostics/Probe.py @@ -92,16 +92,10 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av # ------------------------------------------------------------------- # If timesteps is None, then keep all timesteps otherwise, select timesteps self._timesteps = self._alltimesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") - - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") - + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" + # 3 - Manage axes # ------------------------------------------------------------------- # Fabricate all axes values @@ -362,7 +356,6 @@ def getFields(self): # get the value of x_moved for a requested timestep def getXmoved(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -379,7 +372,6 @@ def getAvailableTimesteps(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+t+" not found in this diagnostic") diff --git a/happi/_Diagnostics/Scalar.py b/happi/_Diagnostics/Scalar.py index 550420e9d..bbacb6484 100755 --- a/happi/_Diagnostics/Scalar.py +++ b/happi/_Diagnostics/Scalar.py @@ -6,39 +6,55 @@ class Scalar(Diagnostic): def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None, **kwargs): # Get available scalars - scalars = self.getScalars() + scalars = self.simulation.getScalars() # If no scalar chosen, only print the available scalars if scalar is None: + error = [] if len(scalars)>0: - self._error += ["Error: no scalar chosen"] - self._error += ["Printing available scalars:"] - self._error += ["---------------------------"] + error += ["Error: no scalar chosen"] + error += ["Printing available scalars:"] + error += ["---------------------------"] l = [""] for s in scalars: if len(s)>4 and s[:2]!=l[-1][:2] and s[-2:]!=l[-1][-2:]: - if l!=[""]: self._error += ["\t".join(l)] + if l!=[""]: error += ["\t".join(l)] l = [] l.append(s) - if l!=[""]: self._error += ["\t".join(l)] + if l!=[""]: error += ["\t".join(l)] else: - self._error += ["No scalars found"] - return + error += ["No scalars found"] + raise Exception("\n".join(error)) # 1 - verifications, initialization # ------------------------------------------------------------------- - # Find which scalar is requested - if scalar not in scalars: - fs = list(filter(lambda x:scalar in x, scalars)) - if len(fs)==0: - self._error += ["No scalar `"+scalar+"` found"] - return - if len(fs)>1: - self._error += ["Several scalars match: "+(' '.join(fs))] - self._error += ["Please be more specific and retry."] - return - scalar = fs[0] - self._scalarname = scalar + # Find which scalar is requested & associated units + self.operation = scalar + def scalarTranslator(s): + units = "??" + if s == "time": + units = "T_r" + elif s == "Ubal_norm": + units = "" + else: + units = { + "U":"K_r * N_r * L_r^%i" % self._ndim_particles, + "P":"K_r * N_r * L_r^%i" % self._ndim_particles, + "D":"N_r * L_r^%i" % self._ndim_particles, + "E":"E_r", + "B":"B_r", + "J":"J_r", + "R":"Q_r * N_r", + "Z":"Q_r", + "N":"", + }[s[0]] + return units, "S['%s']"%s, s + self._operation = Operation(scalar, scalarTranslator, self._ureg) + self._scalarname = self._operation.variables + self._vunits = self._operation.translated_units + self._title = self._operation.title + if not self._scalarname: + raise Exception("String "+self.operation+" does not seem to include any scalar") # Put data_log as object's variable self._data_log = data_log @@ -47,26 +63,31 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None # Already get the data from the file # Loop file line by line self._alltimesteps = [] - self._values = [] - times_values = {} + S = { s:[] for s in self._scalarname } for path in self._results_path: try: - with open(path+'/scalars.txt') as f: - for line in f: - line = line.strip() - if line[0]!="#": break - prevline = line - scalars = prevline[1:].strip().split() # list of scalars - scalarindex = scalars.index(scalar) # index of the requested scalar - line = str(line.strip()).split() - times_values[ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] = float(line[scalarindex]) - for line in f: - line = str(line.strip()).split() - times_values[ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] = float(line[scalarindex]) + f = open(path+'/scalars.txt') except: continue - self._alltimesteps = self._np.array(sorted(times_values.keys())) - self._values = self._np.array([times_values[k] for k in self._alltimesteps]) + for line in f: + if line[0]!="#": break + prevline = line + scalars = prevline[1:].strip().split() # list of scalars + scalarindexes = {s:scalars.index(s) for s in self._scalarname} # indexes of the requested scalars + line = str(line.strip()).split() + self._alltimesteps += [ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] + for s,i in scalarindexes.items(): + S[s] += [float(line[i])] + for line in f: + line = str(line.strip()).split() + self._alltimesteps += [ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] + for s,i in scalarindexes.items(): + S[s] += [float(line[i])] + f.close() + self._alltimesteps = self._np.array(self._alltimesteps) + for s in S: + S[s] = self._np.array(S[s]) + self._values = self._operation.eval(locals()) self._timesteps = self._np.copy(self._alltimesteps) # 2 - Manage timesteps @@ -76,39 +97,11 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None for i,t in enumerate(self._timesteps): self._data.update({ t : i }) # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - self._error += ["Argument `timesteps` must be one or two non-negative integers"] - return + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) # Need at least one timestep - if self._timesteps.size < 1: - self._error += ["Timesteps not found"] - return - - - # 3 - Build units - # ------------------------------------------------------------------- - self._vunits = "??" - if self._scalarname == "time": - self._vunits = "T_r" - elif self._scalarname == "Ubal_norm": - self._vunits = "" - else: - self._vunits = { - "U":"K_r * N_r * L_r^%i" % self._ndim_particles, - "P":"K_r * N_r * L_r^%i" % self._ndim_particles, - "D":"N_r * L_r^%i" % self._ndim_particles, - "E":"E_r", - "B":"B_r", - "J":"J_r", - "R":"Q_r * N_r", - "Z":"Q_r", - "N":"", - }[self._scalarname[0]] - self._title =self._scalarname + assert self._timesteps.size > 0, "Timesteps not found" # Finish constructor self.valid = True @@ -116,37 +109,7 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None # Method to print info on included scalars def _info(self): - return "Scalar "+self._scalarname - - # get all available scalars - def getScalars(self): - allScalars = None - for path in self._results_path: - try: - file = path+'/scalars.txt' - f = open(file, 'r') - except Exception as e: - continue - try: - # Find last commented line - prevline = "" - for line in f: - line = line.strip() - if line[0]!="#": break - prevline = line[1:].strip() - scalars = str(prevline).split() # list of scalars - scalars = scalars[1:] # remove first, which is "time" - except Exception as e: - scalars = [] - f.close() - if allScalars is None: - allScalars = scalars - else: - allScalars = self._np.intersect1d(allScalars, scalars) - if allScalars is None: - self._error += ["Cannot open 'scalars.txt'"] - return [] - return allScalars + return "Scalar "+self.operation # get all available timesteps def getAvailableTimesteps(self): @@ -154,7 +117,6 @@ def getAvailableTimesteps(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/TrackParticles.py b/happi/_Diagnostics/TrackParticles.py index 2d4e8bfb4..de5a63c55 100755 --- a/happi/_Diagnostics/TrackParticles.py +++ b/happi/_Diagnostics/TrackParticles.py @@ -27,6 +27,8 @@ class TrackParticles(ParticleList): def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sorted_as="", length=None, chunksize=20000000, **kwargs): + timestep_indices = kwargs.pop("timestep_indices", None) + # If argument 'species' not provided, then print available species and leave if species is None: species = self.simulation.getTrackSpecies() @@ -89,13 +91,13 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor self._timesteps = self._np.array(sorted(self._locationForTime)) # If specific timesteps requested, narrow the selection - if sorted_as and timesteps is not None: - self._timesteps = self._filterTimesteps( self._timesteps, timesteps ) - + if sorted_as: + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" + self._alltimesteps = self._np.copy(self._timesteps) - if not self._locationForTime: - raise Exception("No data found for this diagnostic") + assert self._locationForTime, "No data found for this diagnostic" # List available properties try: # python 2 @@ -132,11 +134,8 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor if self._timesteps.size == 0: raise Exception("No tracked particles found") # If specific timesteps requested, narrow the selection - if timesteps is not None: - self._timesteps = self._filterTimesteps( self._timesteps, timesteps ) - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" # Select particles @@ -522,7 +521,6 @@ def _orderFiles( self, fileOrdered, chunksize, sort ): # Method to generate the raw data (only done once) def _generateRawData(self, times=None): - if not self._validate(): return self._prepare1() # prepare the vfactor if self._sort: @@ -594,7 +592,6 @@ def _generateRawData(self, times=None): # We override the get and getData methods def getData(self, timestep=None): - if not self._validate(): return self._prepare1() # prepare the vfactor if timestep is None: @@ -629,7 +626,6 @@ def getData(self, timestep=None): # Iterator on UNSORTED particles for a given timestep def iterParticles(self, timestep, chunksize=1): - if not self._validate(): return self._prepare1() # prepare the vfactor if timestep not in self._timesteps: @@ -765,7 +761,6 @@ def _animateOnAxes_2D(self, ax, t, cax_id=0): # Convert data to VTK format def toVTK(self, rendering="trajectory", data_format="xml"): - if not self._validate(): return if self._ndim_particles != 3: print ("Cannot export tracked particles of a "+str(self._ndim_particles)+"D simulation to VTK") diff --git a/happi/_Factories.py b/happi/_Factories.py index 8f47325ec..b862c48b8 100755 --- a/happi/_Factories.py +++ b/happi/_Factories.py @@ -39,10 +39,8 @@ def __init__(self, simulation, scalar=None): # If not a specific scalar (root level), build a list of scalar shortcuts if scalar is None: if simulation._verbose: print("Scanning for Scalar diagnostics") - # Create a temporary, empty scalar diagnostic - tmpDiag = Scalar(simulation) # Get a list of scalars - scalars = tmpDiag.getScalars() + scalars = simulation.getScalars() # Create scalars shortcuts for scalar in scalars: setattr(self, scalar, ScalarFactory(simulation, scalar)) @@ -53,7 +51,25 @@ def __init__(self, simulation, scalar=None): def __call__(self, *args, **kwargs): return Scalar(self._simulation, *(self._additionalArgs+args), **kwargs) - + + def __repr__(self): + msg = object.__repr__(self) + if len(self._additionalArgs) == 0 and self._simulation._scan: + scalars = self._simulation.getScalars() + if scalars: + msg += "\nAvailable Scalar diagnostics:\n" + l = [""] + for s in scalars: + if len(s)>4 and s[:2]!=l[-1][:2] and s[-2:]!=l[-1][-2:]: + if l!=[""]: + msg += "\t".join(l) + "\n" + l = [] + l.append(s) + if l!=[""]: + msg += "\t".join(l) + "\n" + else: + msg += "\nNo Scalar diagnostics available" + return msg class FieldFactory(object): """Import and analyze a Field diagnostic from a Smilei simulation diff --git a/happi/_SmileiSimulation.py b/happi/_SmileiSimulation.py index b8c3a3fbe..279f9f512 100644 --- a/happi/_SmileiSimulation.py +++ b/happi/_SmileiSimulation.py @@ -281,6 +281,35 @@ def _getParticleListSpecies(self, filePrefix): if s: species += [ s.groups()[0] ] return list(set(species)) # unique species + # get all available scalars + def getScalars(self): + allScalars = None + for path in self._results_path: + try: + file = path+'/scalars.txt' + f = open(file, 'r') + except Exception as e: + continue + try: + # Find last commented line + prevline = "" + for line in f: + line = line.strip() + if line[0]!="#": break + prevline = line[1:].strip() + scalars = str(prevline).split() # list of scalars + scalars = scalars[1:] # remove first, which is "time" + except Exception as e: + scalars = [] + f.close() + if allScalars is None: + allScalars = scalars + else: + allScalars = self._np.intersect1d(allScalars, scalars) + if allScalars is None: + return [] + return allScalars + def getTrackSpecies(self): """ List the available tracked species """ return self._getParticleListSpecies("TrackParticlesDisordered") diff --git a/happi/_Utils.py b/happi/_Utils.py index 0b7221545..9fd35a757 100755 --- a/happi/_Utils.py +++ b/happi/_Utils.py @@ -331,8 +331,7 @@ class Operation(object): Parameters: ----------- - operation: the user's requested operation - pattern: the regexp pattern to find the variables in the operation + operation: the user's requested operation (string containing variables and operators) QuantityTranslator: function that takes a string as argument (a quantity name) and outputs its units + its replacement string + its displayed name ureg: Pint's unit registry or None for no unit awareness diff --git a/scripts/compile_tools/machine/joliot_curie_rome b/scripts/compile_tools/machine/joliot_curie_rome index d77047de0..c2974e380 100755 --- a/scripts/compile_tools/machine/joliot_curie_rome +++ b/scripts/compile_tools/machine/joliot_curie_rome @@ -4,12 +4,12 @@ # # Load the correct modules: # - +#module load mpi #module load hdf5 #module switch flavor/hdf5/serial flavor/hdf5/parallel #export HDF5_ROOT_DIR=${HDF5_ROOT} #export OMPI_MCA_btl_portals4_use_rdma=0 -#module load python3/3.7.5 +#module load python3/3.8.10 #export PYTHONEXE=python3 #export PYTHONHOME=$PYTHON3_ROOT #export PATH=$PYTHONHOME/bin:$PATH diff --git a/src/Checkpoint/Checkpoint.cpp b/src/Checkpoint/Checkpoint.cpp index 1ab5aa605..13c3d28a5 100755 --- a/src/Checkpoint/Checkpoint.cpp +++ b/src/Checkpoint/Checkpoint.cpp @@ -35,7 +35,6 @@ #include "LaserEnvelope.h" #include "BinaryProcesses.h" #include "CollisionalNuclearReaction.h" -#include using namespace std; diff --git a/src/Diagnostic/DiagnosticFields.cpp b/src/Diagnostic/DiagnosticFields.cpp index 94ec4edaa..2f6b20444 100755 --- a/src/Diagnostic/DiagnosticFields.cpp +++ b/src/Diagnostic/DiagnosticFields.cpp @@ -154,10 +154,11 @@ DiagnosticFields::DiagnosticFields( Params ¶ms, SmileiMPI *smpi, VectorPatch for( unsigned int ipatch=0 ; ipatchEMfields->allFields_avg.resize( diag_n+1 ); if( time_average > 1 ) { - for( unsigned int ifield=0; ifieldEMfields->allFields_avg[diag_n].push_back( vecPatches( ipatch )->EMfields->createField( fields_names[ifield], params ) ); + } } } } @@ -333,9 +334,16 @@ void DiagnosticFields::run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, // Write H5Write dset = writeField( iteration_group_, fields_names[ifield] ); // Attributes for openPMD + Field *f = vecPatches( 0 )->EMfields->allFields[fields_indexes[ifield]]; + vector stagger( f->dims().size() ); + for( unsigned int i = 0; i < stagger.size(); i++ ) { + stagger[i] = 0.5 * (double) f->isDual(i); + } + bool ends_with_m = 0 == fields_names[ifield].compare( fields_names[ifield].length()-2, 2, "_m" ); + double stagger_t = ends_with_m ? vecPatches( 0 )->EMfields->timestep*0.5 : 0.; openPMD_->writeFieldAttributes( dset, subgrid_start_, subgrid_step_ ); - openPMD_->writeRecordAttributes( dset, field_type[ifield] ); - openPMD_->writeFieldRecordAttributes( dset ); + openPMD_->writeRecordAttributes( dset, field_type[ifield], stagger_t ); + openPMD_->writeFieldRecordAttributes( dset, stagger ); openPMD_->writeComponentAttributes( dset, field_type[ifield] ); } #pragma omp barrier diff --git a/src/Diagnostic/DiagnosticParticleBinningBase.cpp b/src/Diagnostic/DiagnosticParticleBinningBase.cpp index 33288631d..421956ccf 100755 --- a/src/Diagnostic/DiagnosticParticleBinningBase.cpp +++ b/src/Diagnostic/DiagnosticParticleBinningBase.cpp @@ -246,8 +246,8 @@ void DiagnosticParticleBinningBase::calculate_auto_limits( Patch *patch, SimWind if( !std::isnan( axis->min ) && !std::isnan( axis->max ) ) { continue; } - double axis_min = numeric_limits::max(); - double axis_max = numeric_limits::lowest(); + double axis_min = std::isnan( axis->min ) ? numeric_limits::max() : numeric_limits::lowest(); + double axis_max = std::isnan( axis->max ) ? numeric_limits::lowest() : numeric_limits::max(); for( unsigned int i_s=0; i_svecSpecies[species_indices[i_s]]; unsigned int n = s->getNbrOfParticles(); diff --git a/src/Diagnostic/Histogram.cpp b/src/Diagnostic/Histogram.cpp index 0da6251b4..6090e171d 100755 --- a/src/Diagnostic/Histogram.cpp +++ b/src/Diagnostic/Histogram.cpp @@ -71,7 +71,7 @@ void Histogram::digitize( vector species, } } else { // if the particles out of the "box" must be included - + for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { // skip already discarded particles if( int_buffer[ipart] < 0 ) { @@ -82,13 +82,12 @@ void Histogram::digitize( vector species, // move out-of-range indexes back into range if( ind < 0 ) { ind = 0; - } - if( ind >= axis->nbins ) { + } else if( ind >= axis->nbins ) { ind = axis->nbins-1; } int_buffer[ipart] += ind; } - + } } // loop axes diff --git a/src/ElectroMagn/ElectroMagn.h b/src/ElectroMagn/ElectroMagn.h index f21bf53d6..b5bd3e5bc 100755 --- a/src/ElectroMagn/ElectroMagn.h +++ b/src/ElectroMagn/ElectroMagn.h @@ -305,7 +305,6 @@ class ElectroMagn //! Number of elements in arrays without duplicated borders unsigned int bufsize[3][2]; - //!\todo should this be just an integer??? //! Oversize domain to exchange less particles (from params) std::vector oversize; @@ -380,9 +379,6 @@ class ElectroMagn cField *p_AM_; cField *Ap_AM_; - //! \todo check time_dual or time_prim (MG) -// //! method used to solve Maxwell's equation (takes current time and time-step as input parameter) -// virtual void solveMaxwellAmpere() = 0; //! Maxwell Ampere Solver Solver *MaxwellAmpereSolver_; //! Maxwell Faraday Solver diff --git a/src/ElectroMagnSolver/MA_SolverAM_norm.cpp b/src/ElectroMagnSolver/MA_SolverAM_norm.cpp index 11d239a53..6c5633fcd 100755 --- a/src/ElectroMagnSolver/MA_SolverAM_norm.cpp +++ b/src/ElectroMagnSolver/MA_SolverAM_norm.cpp @@ -79,7 +79,7 @@ void MA_SolverAM_norm::operator()( ElectroMagn *fields ) ( *El )( i, j-1 )=-( *El )( i, j+1 ); } for( unsigned int i=0 ; i 1 diff --git a/src/Field/Field.h b/src/Field/Field.h index aecd63a15..669106245 100755 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -86,7 +86,6 @@ class Field virtual void shift_x( unsigned int delta ) = 0; //! vector containing the dimensions of the Field - //! \todo private/friend/modify (JD) std::vector dims_; //! keep track ofwich direction of the Field is dual diff --git a/src/Params/OpenPMDparams.cpp b/src/Params/OpenPMDparams.cpp index b7f4405e8..001836e9e 100755 --- a/src/Params/OpenPMDparams.cpp +++ b/src/Params/OpenPMDparams.cpp @@ -18,16 +18,12 @@ OpenPMDparams::OpenPMDparams( Params &p ): // Grid parameters string xyz = "xyz"; - gridGlobalOffset.resize( params->nDim_field ); - gridOffset .resize( params->nDim_field ); - position .resize( params->nDim_field ); + gridGlobalOffset.resize( params->nDim_field, 0. ); + gridOffset .resize( params->nDim_field, 0. ); gridSpacing .resize( params->nDim_field ); for( unsigned int idim=0; idimnDim_field; idim++ ) { axisLabels.addString( xyz.substr( idim, 1 ) ); - gridGlobalOffset[idim] = 0.; - gridOffset [idim] = 0.; - position [idim] = 0.; - gridSpacing [idim] = params->cell_length[idim]; + gridSpacing[idim] = params->cell_length[idim]; } fieldSolverParameters = ""; if( params->maxwell_sol == "Yee" ) { @@ -192,15 +188,15 @@ void OpenPMDparams::writeSpeciesAttributes( H5Write & ) { } -void OpenPMDparams::writeRecordAttributes( H5Write &location, unsigned int unit_type ) +void OpenPMDparams::writeRecordAttributes( H5Write &location, unsigned int unit_type, double timeOffset ) { location.attr( "unitDimension", unitDimension[unit_type] ); - location.attr( "timeOffset", 0. ); + location.attr( "timeOffset", timeOffset ); } -void OpenPMDparams::writeFieldRecordAttributes( H5Write &location ) +void OpenPMDparams::writeFieldRecordAttributes( H5Write &location, vector &stagger ) { - location.attr( "position", position ); + location.attr( "position", stagger ); } void OpenPMDparams::writeComponentAttributes( H5Write &location, unsigned int unit_type ) diff --git a/src/Params/OpenPMDparams.h b/src/Params/OpenPMDparams.h index 8971360a6..3a3216ea1 100755 --- a/src/Params/OpenPMDparams.h +++ b/src/Params/OpenPMDparams.h @@ -49,8 +49,6 @@ class OpenPMDparams DividedString fieldBoundary, fieldBoundaryParameters, particleBoundary, particleBoundaryParameters; //! current smoothing description std::string currentSmoothing, currentSmoothingParameters; - //! position of the fields within a cell - std::vector position; //! Returns a time string in the openPMD format std::string getLocalTime(); @@ -74,10 +72,10 @@ class OpenPMDparams void writeSpeciesAttributes( H5Write& ); //! Write the attributes for a record - void writeRecordAttributes( H5Write&, unsigned int ); + void writeRecordAttributes( H5Write&, unsigned int, double timeOffset = 0 ); //! Write the attributes for a field record - void writeFieldRecordAttributes( H5Write& ); + void writeFieldRecordAttributes( H5Write&, std::vector &stagger ); //! Write the attributes for a component void writeComponentAttributes( H5Write&, unsigned int ); diff --git a/src/PartCompTime/PartCompTimeFactory.h b/src/PartCompTime/PartCompTimeFactory.h index fccb45038..e8fd5478f 100644 --- a/src/PartCompTime/PartCompTimeFactory.h +++ b/src/PartCompTime/PartCompTimeFactory.h @@ -70,7 +70,7 @@ class PartCompTimeFactory } } else { - ERROR_NAMELIST( "Unknwon parameters setting the particle computing time evaluation operator: " + ERROR_NAMELIST( "Unknown parameters setting the particle computing time evaluation operator: " << params.geometry << ", Order : " << params.interpolation_order, LINK_NAMELIST + std::string("#main-variables")); } diff --git a/src/Particles/ParticleCreator.cpp b/src/Particles/ParticleCreator.cpp index 05a26fb57..270d0fee1 100644 --- a/src/Particles/ParticleCreator.cpp +++ b/src/Particles/ParticleCreator.cpp @@ -285,13 +285,13 @@ int ParticleCreator::create( struct SubSpace sub_space, } // In AM, normalization of weights might be required - bool renormalize = true; + bool regular_weight = false; if( position_initialization_ == "regular" ) { - renormalize = false; + regular_weight = true; } else if( position_initialization_on_species_ ) { unsigned int ispec = species_->position_initialization_on_species_index_; if( patch->vecSpecies[ispec]->position_initialization_ == "regular" ) { - renormalize = false; + regular_weight = true; } } @@ -328,7 +328,7 @@ int ParticleCreator::create( struct SubSpace sub_space, ParticleCreator::createPosition( position_initialization_, regular_number_array_, particles_, species_, nPart, iPart, indexes, params, patch->rand_ ); } ParticleCreator::createMomentum( momentum_initialization_, particles_, species_, nPart, iPart, &temp[0], &vel[0], patch->rand_ ); - ParticleCreator::createWeight( particles_, nPart, iPart, density( i, j, k ), params, renormalize ); + ParticleCreator::createWeight( particles_, nPart, iPart, density( i, j, k ), params, regular_weight ); ParticleCreator::createCharge( particles_, species_, nPart, iPart, charge( i, j, k ) ); iPart += nPart; @@ -931,33 +931,28 @@ void ParticleCreator::createMomentum( std::string momentum_initialization, //! For all (nPart) particles in a mesh initialize its numerical weight (equivalent to a number density) // --------------------------------------------------------------------------------------------------------------------- void ParticleCreator::createWeight( Particles * particles, unsigned int nPart, unsigned int iPart, double n_real_particles, - Params ¶ms, bool renormalize ) + Params ¶ms, bool regular_weight ) { double w = n_real_particles / nPart; for( unsigned int p=iPart; pweight( p ) = w ; } - // In AM, we have a correction to make : multiply by radius + // In AM, we have a correction to make : multiply by cell radius // because n_real_particles was computed with the cell section, not cell volume // See above : density( i, j, k ) *= params.cell_volume; - // where params.cell_volume is 2*pi*dR*dL (in AM only) + // where params.cell_volume is 2*pi*dR*dL (in AM only). Multiplication by R is missing. if( params.geometry == "AMcylindrical" ) { - double radius; - for( unsigned int p=iPart; pposition(1,p)*particles->position(1,p) + particles->position(2,p)*particles->position(2,p)); - particles->weight( p ) *= radius; - } - // We also need to renormalize in case total weight is not exaclty the same anymore - if( renormalize ) { - double total_weight = 0.; + if (regular_weight){ // Multiply each weight by the radius of the particle. Particles towards the axis are lighter. for( unsigned int p=iPart; pweight( p ); + double radius = sqrt(particles->position(1,p)*particles->position(1,p) + particles->position(2,p)*particles->position(2,p)); + particles->weight( p ) *= radius; } + } else { // Multiply each weight by the radius of the cell. All particles of the cell have identical weights. + double radius = sqrt(particles->position(1,iPart)*particles->position(1,iPart) + particles->position(2,iPart)*particles->position(2,iPart)); double cell_radius = params.cell_length[1] * (floor(radius/params.cell_length[1]) + 0.5); - double coeff = n_real_particles * cell_radius / total_weight; for( unsigned int p=iPart; pweight( p ) *= coeff; + particles->weight( p ) *= cell_radius; } } } diff --git a/src/Particles/ParticleCreator.h b/src/Particles/ParticleCreator.h index 336e98ad6..544acda65 100644 --- a/src/Particles/ParticleCreator.h +++ b/src/Particles/ParticleCreator.h @@ -81,7 +81,7 @@ class ParticleCreator unsigned int iPart, double n_real_particles, Params ¶ms, - bool renormalize ); + bool regular_weight ); // For all particles in a mesh initialize its charge state static void createCharge( Particles * particles, Species * species, diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 232a2a474..0c2fbb036 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -1402,6 +1402,10 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int for( unsigned int i=0; ipatches_mins[ipatch][i] ); + } + } + for( unsigned int i=0; ipatches_maxs[ipatch][i] ); } } diff --git a/src/Projector/ProjectorAM2Order.cpp b/src/Projector/ProjectorAM2Order.cpp index caaeff63a..5487995df 100755 --- a/src/Projector/ProjectorAM2Order.cpp +++ b/src/Projector/ProjectorAM2Order.cpp @@ -77,7 +77,7 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, // arrays used for the Esirkepov projection method double Sl0[5], Sl1[5], Sr0[5], Sr1[5], DSl[5], DSr[5]; complex Jl_p[5], Jr_p[5]; - complex e_delta, e_delta_m1, e_delta_inv, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; + complex e_delta, e_delta_m1, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; complex *Jl, *Jr, *Jt, *rho; for( unsigned int i=0; i<5; i++ ) { @@ -107,7 +107,7 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, //calculate exponential coefficients double rp = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); - std::complex theta_old = array_eitheta_old[0]; + std::complex eitheta_old = array_eitheta_old[0]; std::complex eitheta = ( particles.position( 1, ipart ) + Icpx * particles.position( 2, ipart ) ) / rp ; //exp(i theta) e_bar = 1.; // locate the particle on the primal grid at current time-step & calculate coeff. S1 @@ -138,8 +138,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, double r_bar = ((jpo + j_domain_begin_ + deltaold[1*nparts])*dr + rp) * 0.5; // r at t = t0 - dt/2 - e_delta_m1 = std::sqrt(eitheta * (2.*std::real(theta_old) - theta_old)); // std::sqrt keeps the root with positive real part which is what we need here. - e_bar_m1 = theta_old * e_delta_m1; + //e_delta_m1 = std::sqrt(eitheta * (2.*std::real(theta_old) - theta_old)); // std::sqrt keeps the root with positive real part which is what we need here. + e_delta_m1 = std::sqrt(eitheta * std::conj(eitheta_old)); //ei^(theta-theta_old)/2. std::sqrt keeps the root with positive real part which is what we need here. + e_bar_m1 = eitheta_old * e_delta_m1; //ei^(theta+theta_old)/2. ipo -= 2; //This minus 2 come from the order 2 scheme, based on a 5 points stencil from -2 to +2. // i/j/kpo stored with - i/j/k_domain_begin in Interpolator @@ -183,9 +184,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, } e_delta = 0.5; - e_delta_inv = 0.5; + //e_delta_inv = 0.5; - //Compute division by R in advance for Jt and rho evaluation. + //Compute division by R in advance for Jt and rho evaluation. for( unsigned int j=0 ; j<5 ; j++ ) { Sr0[j] *= invR__local[j]; Sr1[j] *= invR__local[j]; @@ -197,8 +198,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, e_delta *= e_delta_m1; e_bar *= e_bar_m1; C_m = 2. * e_bar ; //multiply modes > 0 by 2 and C_m = 1 otherwise. - e_delta_inv =1./e_delta - 1.; + //e_delta_inv =1./e_delta - 1.; crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.*r_bar; + //crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.; } // Add contribution J_p to global array @@ -247,7 +249,8 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, iloc = ( i+ipo )*nprimr_ + jpo; for( unsigned int j=0 ; j<5 ; j++ ) { linindex = iloc+j; - Jt [linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. )); + //Jt [linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. )); + Jt [linindex] -= crt_p*(Sr1[j]*Sl1[i]*( e_delta-1. ) - Sr0[j]*Sl0[i]*(std::conj(e_delta) - 1.)); } } @@ -354,8 +357,6 @@ void ProjectorAM2Order::basicForComplexOnBuffer( complex *rhoj, Particle // This function also assumes that particles position is evaluated at the same time as currents which is usually not true (half time-step difference). // It will therefore fail to evaluate the current accurately at t=0 if a plasma is already in the box. - - // ------------------------------------- // Variable declaration & initialization // ------------------------------------- @@ -366,17 +367,6 @@ void ProjectorAM2Order::basicForComplexOnBuffer( complex *rhoj, Particle if( type > 0 ) { //if current density ERROR("This projector can be used only for charge density at the moment."); - // charge_weight *= 1./sqrt( 1.0 + particles.momentum( 0, ipart )*particles.momentum( 0, ipart ) - // + particles.momentum( 1, ipart )*particles.momentum( 1, ipart ) - // + particles.momentum( 2, ipart )*particles.momentum( 2, ipart ) ); - // if( type == 1 ) { //if Jl - // charge_weight *= particles.momentum( 0, ipart ); - // } else if( type == 2 ) { //if Jr - // charge_weight *= ( particles.momentum( 1, ipart )*particles.position( 1, ipart ) + particles.momentum( 2, ipart )*particles.position( 2, ipart ) )/ r ; - // nr++; - // } else { //if Jt - // charge_weight *= ( -particles.momentum( 1, ipart )*particles.position( 2, ipart ) + particles.momentum( 2, ipart )*particles.position( 1, ipart ) ) / r ; - // } } complex e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; @@ -460,8 +450,9 @@ void ProjectorAM2Order::axisBC(ElectroMagnAM *emAM, bool diag_flag ) void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ) { - double sign = 1.; - for (unsigned int i=0; i< imode; i++) sign *= -1; + // Mode 0 contribution "below axis" is added. + // Non zero modes are substracted because a particle sitting exactly on axis has a non defined theta and can not contribute to a theta dependent mode. + double sign = (imode == 0) ? 1 : -1 ; if (diag_flag && rhoj) { for( unsigned int i=2 ; i *rhoj,std::complex 0){ rhoj[i] = 0.; - rhoj[i-1] = - rhoj[i+1]; + rhoj[i-1] = - rhoj[i+1]; // Zero Jl mode > 0 on axis. } else { - //This is just for cosmetics on the picture, rho has no influence on the results - rhoj[i] = (4.*rhoj[i+1] - rhoj[i+2])/3.; - rhoj[i-1] = rhoj[i+1]; + rhoj[i] = rhoj[i+1]; //This smoothing is just for cosmetics on the picture, rho has no influence on the results. + rhoj[i-1] = rhoj[i+1]; // Non zero Jl mode > 0 on axis. } } } @@ -489,10 +479,11 @@ void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex 0){ Jl [i] = 0. ; - Jl[i-1] = -Jl[i+1]; + Jl[i-1] = -Jl[i+1]; // Zero Jl mode > 0 on axis. } else { - //Jl mode 1 on axis is left as is. It looks over estimated but it is necessary to conserve a correct divergence and a proper evaluation on the field on axis. - Jl [i-1] = Jl [i+1] ; + //Jl mode 0 on axis should be left as is. It looks over estimated but it might be necessary to conserve a correct divergence and a proper evaluation on the field on axis. + //Jl [i] = (4.*Jl [i+1] - Jl [i+2])/3. ; //Here we use smoothig on axis for results appearing more physical. + Jl [i-1] = Jl [i+1] ; // Non zero Jl mode 0 on axis. } } } @@ -503,15 +494,15 @@ void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex