Merge branch 'DDfusion' into develop

.gitlab-ci.yml

@@ -6,7 +6,7 @@ compile:
   stage: compile
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Force workdir cleaning in case of retried
@@ -28,7 +28,7 @@ makerun1D:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir
@@ -40,7 +40,7 @@ makerun2D:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir
@@ -52,7 +52,7 @@ makerun3D:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir
@@ -64,7 +64,7 @@ makerunAM:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir
@@ -76,7 +76,7 @@ makerunV:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir
@@ -88,7 +88,7 @@ makerunCollisions:
   stage: makerun
   only:
     - develop
-    - smilei_Poisson_AM
+    - DDfusion
 
   script:
     # Move in test dir

benchmarks/tst_collisions4_DD_fusion_rate.py

@@ -0,0 +1,133 @@
+# ---------------------------------------------
+# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
+# ---------------------------------------------
+
+import math
+L0 = 2.*math.pi # conversion from normalization length to wavelength
+
+
+Main(
+	geometry = "1Dcartesian",
+
+	number_of_patches = [ 8 ],
+
+	interpolation_order = 2,
+
+	timestep = 0.05 * L0,
+	simulation_time = 10 * L0,
+
+
+	time_fields_frozen = 100000000000.,
+
+	cell_length = [5.*L0],
+	grid_length = [1600.*L0],
+
+	EM_boundary_conditions = [ ["periodic"] ],
+
+	reference_angular_frequency_SI = L0 * 3e8 /1.e-6,
+
+	random_seed = smilei_mpi_rank	
+)
+
+
+i = 0
+for Da_nppc, Db_nppc, v in [
+	[100, 100, 0.1],
+	[100,  10, 0.1],
+	[10 , 100, 0.1],
+	[100, 100, 0.00082],
+	[100, 100, 0.00201],
+	[100, 100, 0.00493],
+	[100, 100, 0.01209],
+	[100, 100, 0.02960],
+	[100, 100, 0.07236],
+	[100, 100, 0.17545],
+	[100, 100, 0.40641],
+	[100, 100, 0.76966],
+	[100, 100, 0.97376],
+	]:
+
+	Da = "Da_"+str(i)
+	Db = "Db_"+str(i)
+	He = "He_"+str(i)
+
+	Species(
+		name = Da,
+		position_initialization = "regular",
+		momentum_initialization = "maxwell-juettner",
+		particles_per_cell= Da_nppc,
+		atomic_number = 1,
+		mass = 3870.5,
+		charge = 0.,
+		number_density = 100.,
+		mean_velocity = [v, 0., 0.],
+		temperature = [0.0000001]*3,
+		time_frozen = 100000000.0,
+		boundary_conditions = [
+			["periodic", "periodic"],
+		],
+	)
+
+	Species(
+		name = Db,
+		position_initialization = "regular",
+		momentum_initialization = "maxwell-juettner",
+		particles_per_cell= Db_nppc,
+		atomic_number = 1,
+		mass = 3870.5,
+		charge = 0.,
+		number_density = 100.,
+		mean_velocity = [-v, 0., 0.],
+		temperature = [0.00000001]*3,
+		time_frozen = 100000000.0,
+		boundary_conditions = [
+			["periodic", "periodic"],
+		],
+	)
+
+	Species(
+		name = He,
+		position_initialization = "regular",
+		momentum_initialization = "maxwell-juettner",
+		particles_per_cell= 0,
+		atomic_number = 2,
+		mass = 5497.9,
+		charge = 0.,
+		number_density = 0.,
+		time_frozen = 100000000.0,
+		boundary_conditions = [
+			["periodic", "periodic"],
+		],
+	)
+
+	Collisions(
+		species1 = [Da],
+		species2 = [Db],
+		coulomb_log = 0.001,
+		nuclear_reaction = [He],
+		debug_every = 10
+	)
+
+	DiagParticleBinning(
+		deposited_quantity = "weight",
+		every = 20,
+		species = [He],
+		axes = [
+			 ["x",    0,    Main.grid_length[0],   1]
+		]
+	)
+
+	DiagParticleBinning(
+		deposited_quantity = "weight",
+		every = 20,
+		species = [Da, Db],
+		axes = [
+			 ["x",    0,    Main.grid_length[0],   1]
+		]
+	)
+
+	i+=1
+
+DiagScalar(
+	every = 50,
+)

doc/Sphinx/namelist.rst

@@ -1961,6 +1961,7 @@ Collisions
       coulomb_log = 5.,
       debug_every = 1000,
       ionizing = False,
+  #      nuclear_reaction = [],
   )
 
 
@@ -2013,7 +2014,7 @@ Collisions
 
 .. py:data:: ionizing
 
-  :default: False
+  :default: ``False``
 
   :ref:`Collisional ionization <CollIonization>` is set when this parameter is not ``False``.
   It can either be set to the name of a pre-existing electron species (where the ionized
@@ -2023,7 +2024,36 @@ Collisions
   One of the species groups must be all electrons (:py:data:`mass` = 1), and the other
   one all ions of the same :py:data:`atomic_number`.
 
+.. rst-class:: experimental
+
+.. py:data:: nuclear_reaction
+  
+  :type: a list of strings
+  :default: ``None`` (no nuclear reaction)
+
+  A list of the species names for the products of nuclear reactions
+  that may occur during collisions. You may omit product species if they are not necessary
+  for the simulation.
+
+  All members of :py:data:`species1` must be the same type of atoms, which is automatically
+  recognized by their :py:data:`mass` and :py:data:`atomic_number`. The same applies for
+  all members of :py:data:`species2`.
+
+  In the current version, only the reaction D(d,n)He³ is available.
+
+.. rst-class:: experimental
+
+.. py:data:: nuclear_reaction_multiplier
+  
+  :type: a float
+  :default: 0. (automatically adjusted)
 
+  The rate multiplier for nuclear reactions. It is a positive number that artificially
+  increases the occurence of reactions so that a good statistics is obtained. The number
+  of actual reaction products is adjusted by changing their weights in order to provide
+  a physically correct number of reactions. Leave this number to ``0.`` for an automatic
+  rate multiplier: the final number of produced macro-particles will be of the same order
+  as that of reactants.
 
 
 

src/Checkpoint/Checkpoint.cpp

@@ -26,6 +26,7 @@
 #include "DiagnosticScreen.h"
 #include "DiagnosticTrack.h"
 #include "LaserEnvelope.h"
+#include "Collisions.h"
 
 using namespace std;
 
@@ -255,7 +256,7 @@ void Checkpoint::dumpAll( VectorPatch &vecPatches, unsigned int itime,  SmileiMP
         string patchName=Tools::merge( "patch-", patch_name.str() );
         hid_t patch_gid = H5::group( fid, patchName.c_str() );
 
-        dumpPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, params, patch_gid );
+        dumpPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, vecPatches( ipatch )->vecCollisions, params, patch_gid );
 
         // Random number generator state
         H5::attr( patch_gid, "xorshift32_state", vecPatches( ipatch )->xorshift32_state );
@@ -283,7 +284,7 @@ void Checkpoint::dumpAll( VectorPatch &vecPatches, unsigned int itime,  SmileiMP
 
 }
 
-void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector<Species *> vecSpecies, Params &params, hid_t patch_gid )
+void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector<Species *> vecSpecies, std::vector<Collisions *> &vecCollisions, Params &params, hid_t patch_gid )
 {
     if (  params.geometry != "AMcylindrical" ) {
         dumpFieldsPerProc( patch_gid, EMfields->Ex_ );
@@ -442,6 +443,13 @@ void Checkpoint::dumpPatch( ElectroMagn *EMfields, std::vector<Species *> vecSpe
         H5Gclose( gid );
 
     } // End for ispec
+
+    // Manage some collisions parameters
+    std::vector<double> rate_multiplier( vecCollisions.size() );
+    for( unsigned int icoll = 0; icoll<vecCollisions.size(); icoll++ ) {
+        rate_multiplier[icoll] = vecCollisions[icoll]->NuclearReaction->rate_multiplier_;
+    }
+    H5::vect( patch_gid, "collisions_rate_multiplier", rate_multiplier );
 };
 
 
@@ -539,7 +547,7 @@ void Checkpoint::restartAll( VectorPatch &vecPatches,  SmileiMPI *smpi, SimWindo
         string patchName=Tools::merge( "patch-", patch_name.str() );
         hid_t patch_gid = H5Gopen( fid, patchName.c_str(), H5P_DEFAULT );
 
-        restartPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, params, patch_gid );
+        restartPatch( vecPatches( ipatch )->EMfields, vecPatches( ipatch )->vecSpecies, vecPatches( ipatch )->vecCollisions, params, patch_gid );
 
         // Random number generator state
         H5::getAttr( patch_gid, "xorshift32_state", vecPatches( ipatch )->xorshift32_state );
@@ -566,7 +574,7 @@ void Checkpoint::restartAll( VectorPatch &vecPatches,  SmileiMPI *smpi, SimWindo
 }
 
 
-void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector<Species *> &vecSpecies, Params &params, hid_t patch_gid )
+void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector<Species *> &vecSpecies, std::vector<Collisions *> &vecCollisions, Params &params, hid_t patch_gid )
 {
     if ( params.geometry != "AMcylindrical" ) {
         restartFieldsPerProc( patch_gid, EMfields->Ex_ );
@@ -754,6 +762,14 @@ void Checkpoint::restartPatch( ElectroMagn *EMfields, std::vector<Species *> &ve
         H5Gclose( gid );
     }
 
+    // Manage some collisions parameters
+    if( H5::getVectSize( patch_gid, "collisions_rate_multiplier" ) > 0 ) {
+        std::vector<double> rate_multiplier;
+        H5::getVect( patch_gid, "collisions_rate_multiplier", rate_multiplier, true );
+        for( unsigned int icoll = 0; icoll<rate_multiplier.size(); icoll++ ) {
+            vecCollisions[icoll]->NuclearReaction->rate_multiplier_ = rate_multiplier[icoll];
+        }
+    }
 }
 
 void Checkpoint::dumpFieldsPerProc( hid_t fid, Field *field )

src/Checkpoint/Checkpoint.h

@@ -23,6 +23,7 @@ class Field;
 class cField;
 class Species;
 class VectorPatch;
+class Collisions;
 
 #include <csignal>
 
@@ -42,7 +43,7 @@ class Checkpoint
     //! restart everything to file per processor
     void readPatchDistribution( SmileiMPI *smpi, SimWindow *simWin );
     void restartAll( VectorPatch &vecPatches,  SmileiMPI *smpi, SimWindow *simWin, Params &params, OpenPMDparams &openPMD );
-    void restartPatch( ElectroMagn *EMfields, std::vector<Species *> &vecSpecies, Params &params, hid_t patch_gid );
+    void restartPatch( ElectroMagn *EMfields, std::vector<Species *> &vecSpecies, std::vector<Collisions *> &vecCollisions, Params &params, hid_t patch_gid );
 
     //! restart field per proc
     void restartFieldsPerProc( hid_t fid, Field *field );
@@ -58,7 +59,7 @@ class Checkpoint
 
     //! dump everything to file per processor
     void dumpAll( VectorPatch &vecPatches, unsigned int itime,  SmileiMPI *smpi, SimWindow *simWin, Params &params );
-    void dumpPatch( ElectroMagn *EMfields, std::vector<Species *> vecSpecies, Params &params, hid_t patch_gid );
+    void dumpPatch( ElectroMagn *EMfields, std::vector<Species *> vecSpecies, std::vector<Collisions *> &vecCollisions, Params &params, hid_t patch_gid );
 
     //! incremental number of times we've done a dump
     unsigned int dump_number;