Add PRL LULI simu. example

input-deck

@@ -91,17 +91,17 @@
 #########################################################################
 #########################################################################
 #
-#simu Vlasov
+#simu NJP-Academic-case
 #
-#N_threads   1
+#N_threads   0
 #
 #hll_order           2
 #
-#implicit_coll   .false.
+#implicit_coll   .true.
 #
 #cfl               0.9
 #
-#bi_temp        .true.
+#bi_temp        .false.
 #
 #magnetic_diff .true.
 #
@@ -111,23 +111,23 @@
 #
 #Kalpha .false.
 #
-#L_t    150.
+#L_t    3500.
 #
-#Delta_t_diag 30.
+#Delta_t_diag 35.
 #
-#d_z      3.
+#d_z      1.
 #
-#d_x      6.
+#d_x      1.
 #
-#L_z      9.
+#L_z    100.
 #
-#L_x     60.
+#L_x    100.
 #
-#d_eps     10.
+#d_eps     5.
 #
-#eps_min   10.
+#eps_min   20.
 #
-#L_eps     9.e3
+#L_eps     1.5e3
 #
 #########################################################################
 #########################################################################
@@ -192,15 +192,15 @@
 #########################################################################
 #########################################################################
 #
-#E_tot      2.
+#E_tot      10.
 #
-#Delta_t    30.
+#Delta_t   1177.41
 #
-#Delta_r    6.
+#Delta_r    23.5482
 #
-#spectrum   5
+#spectrum   1
 #
-#eps0       700.
+#eps0      1000.
 #
 #sigma_eps   50.
 #
@@ -214,9 +214,9 @@
 #
 #eps1        3.e2
 #
-#angle0       20.
+#angle0       0.
 #
-#Delta_theta  40.
+#Delta_theta  0.
 #
 #########################################################################
 #########################################################################
@@ -280,19 +280,19 @@
 #########################################################################
 #########################################################################
 #
-#Material  2
+#Material  1
 #
 #Mat Al
 #
-#Z0  13.
+#Z0  1.
 #
-#A0  26.98
+#A0  1.008
 #
 #tabulated_plasma .false.
 #
-#rho   2.7
+#rho  50.
 #
-#T_ini 2.6e-2
+#T_ini 1.e0
 #
 #tabulated_resistivity .false.
 #

readme.md

@@ -26,7 +26,21 @@
 To achieve this, it computes the two first momentum angular moments of the beam electrons Vlasov-Fokker-Planck-Beliaev-Budker (VFPBB) kinetic equation (the two first order equations of the kinetic equation Cartesian tensor scalar product expansion), completed with the Minerbo maximum angular entropy closure to express the distribution function second order momentum angular moment needed in the first order equation. The resulting kinetic equations are coupled with the target electrons and ions Magneto-Hydrodynamic (MHD) equations considering time scales on which the target electron return current has already set up. [AMoRE](https://github.com/michaeltouati/AMoRE) thus takes into account both collective effects, with the self-generated electromagnetic fields, and collisional effects, with the slowing down of beam electrons in collisions with plasmons, bound and free electrons and their angular scattering on both ions and electrons. The kinetic energy cut-off separating target electrons and beam electrons is assumed to be around 10 keV for the beam electrons collisional friction and diffusion VFPBB terms and the solid/plasma species MHD equations to be valid. More pieces of information can be found in my [PhD manuscript](https://tel.archives-ouvertes.fr/tel-01238782/document) and in this [peer-reviewed article](https://iopscience.iop.org/article/10.1088/1367-2630/16/7/073014/pdf). In the latter, there is no mention of updates in relation with the possible beam electrons specular reflection boundary conditions (refluxing), thermal capacities valid at ambient temperatures according to the Einstein and Debye models as well as the distinction between collisions of free electrons with s, p or d-band electrons in the electron-electron collisions contribution on transport coefficients. Citations to these references are recommended and appreciated for publications of scientific results using [AMoRE](https://github.com/michaeltouati/AMoRE) in peer-reviewed journals. 
 
 Python scripts, using the Matplotlib and Numpy packages, are provided to automatically extract and plot the simulation results. The simulation parameters are described in the [input-deck](https://github.com/michaeltouati/AMoRE/blob/main/input-deck) and they can be modified without having to recompile the code. Compilation rules can be modified in the [makefile](https://github.com/michaeltouati/AMoRE/blob/main/makefile) depending on the user compiler preferences. Tools for testing the compilation of the code and tools for checking the simulation parameters are provided. X-ray Kα and Kβ photon emission from the target can be computed using a provided [X-ray Kα and Kβ table](https://github.com/michaeltouati/AMoRE/blob/master/sources/data/Kalpha_tab.dat). The user can use its own [tabulated electrical resistivity](https://github.com/michaeltouati/AMoRE/blob/master/sources/user/resistivity_tab.dat) for the solid/plasma and/or its own [tabulated plasma material-density-temperature profile](https://github.com/michaeltouati/AMoRE/blob/master/sources/user/plasma_tab.dat) and/or its own [tabulated electron beam kinetic energy spectrum](https://github.com/michaeltouati/AMoRE/blob/master/sources/user/spectrum_tab.dat).
-
+
+# Simulation plot examples
+
+<p align="center">
+  <img width="400" height="325" src="test-cases/NJP-Academic-case/nb_51.png">
+  <img width="400" height="325" src="test-cases/NJP-Academic-case/fb_x51_51.png">
+</p>
+
+<p align="center">
+  <img width="200" height="325" src="test-cases/PRL-LULI-C-vitreous/nb_51.png">
+  <img width="200" height="325" src="test-cases/PRL-LULI-C-vitreous/Te_51.png">
+  <img width="200" height="325" src="test-cases/PRL-LULI-C-vitreous/By_51.png">
+  <img width="200" height="325" src="test-cases/PRL-LULI-C-vitreous/nKa_51.png">
+</p>
+
 # Compiling the code
 
 Modify the [makefile](https://github.com/michaeltouati/AMoRE/blob/main/makefile) as a function of the Fortran compiler installed on your computer and then type

sources/Initialization.f90

@@ -373,13 +373,13 @@ subroutine initialize(x_a, z_a, eps_a, ni, A, Z, Ex, Ez, By_n, By_np1,&
           end do
           close(600)
           N_plasma_tab = i
-          allocate(plasma_tab(1:N_plasma_tab,1:4))
+          allocate(plasma_tab(1:N_plasma_tab,1:6))
           !  Read the file 'sources/user/plasma_tab.dat' and define the table plasma_tab with the 
           !  corresponding values :
           open (unit=700,file='sources/user/plasma_tab.dat',form='formatted',status='unknown')      
           read(700,*)
           do i =1,N_plasma_tab,1
-            read(700,*) plasma_tab(i,1),plasma_tab(i,2),plasma_tab(i,3),plasma_tab(i,4)
+            read(700,*) plasma_tab(i,1),plasma_tab(i,2),plasma_tab(i,3),plasma_tab(i,4),plasma_tab(i,5),plasma_tab(i,6)
           end do
           close(700)
           ! find the number of cells on x-axis and deduce the
@@ -440,39 +440,39 @@ subroutine initialize(x_a, z_a, eps_a, ni, A, Z, Ex, Ez, By_n, By_np1,&
                         / (plasma_tab((im-1)*Nz_plasma_tab+km+1,2) - plasma_tab((im-1)*Nz_plasma_tab+km,2))
                 beta_z  = (z_a(k) - plasma_tab((im-1)*Nz_plasma_tab+km,2))&
                         / (plasma_tab((im-1)*Nz_plasma_tab+km+1,2) - plasma_tab((im-1)*Nz_plasma_tab+km,2))
-                r = (alpha_x * ((alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,3))&
-                               +(beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,3))))&
-                  + (beta_x  * ((alpha_z * plasma_tab( im   *Nz_plasma_tab+km  ,3))&
-                               +(beta_z  * plasma_tab( im   *Nz_plasma_tab+km+1,3))))
-                t = (alpha_x * ((alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,4))&
-                               +(beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,4))))&
-                  + (beta_x  * ((alpha_z * plasma_tab( im   *Nz_plasma_tab+km  ,4))&
-                               +(beta_z  * plasma_tab( im   *Nz_plasma_tab+km+1,4))))
+                r = (alpha_x * ((alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,5))&
+                               +(beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,5))))&
+                  + (beta_x  * ((alpha_z * plasma_tab( im   *Nz_plasma_tab+km  ,5))&
+                               +(beta_z  * plasma_tab( im   *Nz_plasma_tab+km+1,5))))
+                t = (alpha_x * ((alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,6))&
+                               +(beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,6))))&
+                  + (beta_x  * ((alpha_z * plasma_tab( im   *Nz_plasma_tab+km  ,6))&
+                               +(beta_z  * plasma_tab( im   *Nz_plasma_tab+km+1,6))))
               else if ((im == Nx_plasma_tab).and.(km.ne.Nz_plasma_tab)) then
                 alpha_z = (plasma_tab((im-1)*Nz_plasma_tab+km+1,2) - z_a(k))&
                         / (plasma_tab((im-1)*Nz_plasma_tab+km+1,2) - plasma_tab((im-1)*Nz_plasma_tab+km,2))
                 beta_z  = (z_a(k) - plasma_tab((im-1)*Nz_plasma_tab+km,2))&
                         / (plasma_tab((im-1)*Nz_plasma_tab+km+1,2) - plasma_tab((im-1)*Nz_plasma_tab+km,2))
-                r = (alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,3))&
-                  + (beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,3))
-                t = (alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,4))&
-                  + (beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,4))
+                r = (alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,5))&
+                  + (beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,5))
+                t = (alpha_z * plasma_tab((im-1)*Nz_plasma_tab+km  ,6))&
+                  + (beta_z  * plasma_tab((im-1)*Nz_plasma_tab+km+1,6))
               else if ((im.ne.Nx_plasma_tab).and.(km == Nz_plasma_tab)) then
                 alpha_x = (plasma_tab(im * Nz_plasma_tab+km,1) - x_a(i))&
                         / (plasma_tab(im * Nz_plasma_tab+km,1) - plasma_tab((im-1)*Nz_plasma_tab+km,1))
                 beta_x  = (x_a(i) - plasma_tab((im-1)*Nz_plasma_tab+km,1))&
                         / (plasma_tab(im*Nz_plasma_tab+km,1) - plasma_tab((im-1)*Nz_plasma_tab+km,1))
-                r = (alpha_x * plasma_tab((im-1)*Nz_plasma_tab+km,3))&
-                  + (beta_x  * plasma_tab( im   *Nz_plasma_tab+km,3))
-                t = (alpha_x * plasma_tab((im-1)*Nz_plasma_tab+km,4))&
-                  + (beta_x  * plasma_tab( im   *Nz_plasma_tab+km,4))
+                r = (alpha_x * plasma_tab((im-1)*Nz_plasma_tab+km,5))&
+                  + (beta_x  * plasma_tab( im   *Nz_plasma_tab+km,5))
+                t = (alpha_x * plasma_tab((im-1)*Nz_plasma_tab+km,6))&
+                  + (beta_x  * plasma_tab( im   *Nz_plasma_tab+km,6))
               else
-                r = plasma_tab((im-1)*Nz_plasma_tab+km,3)
-                t = plasma_tab((im-1)*Nz_plasma_tab+km,4)
+                r = plasma_tab((im-1)*Nz_plasma_tab+km,5)
+                t = plasma_tab((im-1)*Nz_plasma_tab+km,6)
               end if
-              A(i,k)  = A1
-              Z(i,k)  = Z1
-              ni(i,k) = r / (A1 * mu)
+              A(i,k)  = plasma_tab((im-1)*Nz_plasma_tab+km,3)
+              Z(i,k)  = plasma_tab((im-1)*Nz_plasma_tab+km,4)
+              ni(i,k) = r / (A(i,k) * mu)
               Te(i,k) = t * eV
               Ti(i,k) = t * eV
             end do

sources/Transport_coef.f90

@@ -639,17 +639,17 @@ function resis(eta_tab, Z, ni, Te, Ti)
     real(PR), dimension(1:N_eta_tab,1:2), intent(in) :: eta_tab
     real(PR), intent(in)                             :: Z, ni, Te, Ti
     real(PR)                                         :: resis
-    if (tabulated_resistivity) then
-      resis = resistivity_tab(eta_tab, Z, ni, Te, Ti)
-    else
-      ! Metals
-      if ((Z == Z_Al).or.(Z == Z_Cu).or.(Z == Z_Ta)) then
-        resis = resistivity_metals(Z, ni, Te, Ti)
-      ! Hydrogen plasmas
-      else if (Z == 1) then
-        resis = resistivity_H(Z, ni, Te, Ti)
-      ! Other plasmas
-      else 
+    ! Metals
+    if ((Z == Z_Al).or.(Z == Z_Cu).or.(Z == Z_Ta)) then
+      resis = resistivity_metals(Z, ni, Te, Ti)
+    ! Hydrogen plasmas
+    else if (Z == 1) then
+      resis = resistivity_H(Z, ni, Te, Ti)
+    ! Other plasmas
+    else 
+      if (tabulated_resistivity) then
+        resis = resistivity_tab(eta_tab, Z, ni, Te, Ti)
+      else
         resis = resistivity_Spitzer(Z, ni, Te, Ti)
       end if
     end if

sources/plot/library.py

@@ -37,6 +37,18 @@
         'size': FONT_SIZE,
         }
 
+def get_string_parameter(param):
+    """
+    Return AMoRE simulation input deck param value
+    """
+    with open('input-deck', 'r', encoding='utf-8') as file :
+        for line in file:
+            array     = line.strip().split()
+            if array[0] == '#'+param :
+                param_value = array[1]
+                break
+    return param_value
+
 def get_results_dir():
     """
     Return AMoRE simulation name
@@ -48,6 +60,7 @@ def get_results_dir():
             if array[0] == '#simu' :
                 name = array[1]
                 break
+    name = get_string_parameter('simu')
     to_print      = ' ' + name + ' AMoRE SIMULATION PLOTS :'
     line_to_print = ' '
     for _ in range(0,len(to_print)-1):
@@ -562,7 +575,7 @@ def construct_distribution_and_plot(**kwargs):
                                     xmap_max = np.matrix(kwargs['pz_map']).max(),
                                     ymap_min = np.matrix(kwargs['px_map']).min(),
                                     ymap_max = np.matrix(kwargs['px_map']).max(),
-                                    zmap_min = np.matrix(fb_map).min(),
+                                    zmap_min = -15.,#np.matrix(fb_map).min(),
                                     zmap_max = np.matrix(fb_map).max(),
                                     xlabel   = r'$p_z\,(m_e c)$',
                                     ylabel   = r'$p_x\,(m_e c)$',