Add python script for rate coefficients conversion

moduleOutput.f90

@@ -146,7 +146,7 @@ module output
 
       do j = 1, nz
         write (timeString, formatTime) t
-        write(ZString, '(F6.3)') Z_list(j)
+        write(ZString, '(F4.1)') Z_list(j)
         ZString = adjustl(trim(ZString))
 
         ZString = adjustl(ZString)
@@ -279,7 +279,7 @@ module output
 
       do j = 1, nz
         write (timeString, formatTime) t
-        write(ZString, '(F6.3)') Z_list(j)
+        write(ZString, '(F4.1)') Z_list(j)
         ZString = adjustl(trim(ZString))
 
         ZString = adjustl(ZString)

scripts_python/rate_coef.py

@@ -0,0 +1,66 @@
+import re
+import pandas as pd
+
+file_path = "rate_coef_Sn.dat"
+
+with open(file_path, "r") as f:
+    content = f.read()
+
+lines = content.splitlines()
+
+# Dictionary to store tables with their respective dataframes
+tables = {}
+current_table = None
+headers = []
+data_rows = []
+
+# Mapping for the short names to full descriptions
+full_name_map = {
+    "ci": "electron_impact_direct_ionization",
+    "ea": "electron_impact_excitation_autoionization",
+    "rr": "radiative_recombination",
+    "dr": "dielectronic_recombination",
+    "bb": "bound_bound_radiative_loss_rates",
+    "bf": "recombination_radiative_loss_rates",
+    "csd": "charge_state_distribution"
+}
+
+# Regular expression to detect a table header (e.g., "ci temp Sn0+ Sn1+ ...")
+header_pattern = re.compile(r'^(ci|ea|rr|dr|bb|bf|csd)\s+temp')
+
+# Iterate through each line in the file
+for line in lines:
+    if header_pattern.match(line):
+        if current_table and headers and data_rows:
+            headers = ['temp' if h == 'temp' else re.sub(r'\D', '', h) for h in headers]
+            df = pd.DataFrame(data_rows, columns=headers)
+            df.iloc[:, 0] = pd.to_numeric(df.iloc[:, 0], errors='coerce')  # Convert temp column to numeric
+            df.set_index(headers[0], inplace=True)
+            tables[current_table] = df
+
+        parts = line.strip().split()
+        current_table = parts[0]
+        headers = parts[1:]
+        data_rows = []
+    elif current_table and line.strip() and not line.strip().startswith("#"):
+        parts = line.strip().split()
+        if len(parts) == len(headers) + 1:
+            data_rows.append(parts)
+        elif len(parts) == len(headers): 
+            data_rows.append(parts)
+
+# Save the last table
+if current_table and headers and data_rows:
+    df = pd.DataFrame(data_rows, columns=headers)
+    df.iloc[:, 0] = pd.to_numeric(df.iloc[:, 0], errors='coerce')
+    df.set_index(headers[0], inplace=True)
+    tables[current_table] = df
+
+# Save all tables to CSV files with full descriptive names
+csv_paths = {}
+for name, df in tables.items():
+    full_name = full_name_map.get(name, name)  # Default to name if no mapping found
+    path = f"{full_name}.csv"
+    df.to_csv(path)
+    csv_paths[name] = path
+

vlaplex.f90

@@ -20,14 +20,36 @@
 
 end module eos
 
+ module Temperature
+   use constantParameters, only: dp
+   implicit none
+
+   private
+   public:: T_function
+
+   contains
+     pure function T_function(T0,n0, ne, gamma_e) result(T)
+       use constantParameters, only: eV_to_K
+       implicit none
+
+       real(dp), intent(in):: T0, n0, ne, gamma_e
+       real(dp):: T
+
+      T = T0 * (ne/n0)**(gamma_e -1.0_dp)
+
+     end function T_function
+
+end module Temperature
 program VlaPlEx
   use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI
   use output
   use referenceValues
   use eos,                only: T_to_Z
+  use Temperature,         only: T_function
   use moduleTableBC
   use moduleTableTtoZ
   use moduleTableTtoZne
+  use moduleTableRateCoeff
   use omp_lib
   implicit none
 
@@ -60,6 +82,10 @@ program VlaPlEx
   character(:), allocatable:: TtoZ_file
   type(tableTtoZne):: TtoZne
   character(:), allocatable:: TtoZne_file
+  type(tableRateCoeff):: recombination_rateCoeff
+  character(:), allocatable:: recombination_rate_file
+  type(tableRateCoeff):: ionization_rateCoeff
+  character(:), allocatable:: ionization_rate_file
 
   real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
   real(dp), allocatable, dimension(:):: f0   ! Boundary at r = x_0
@@ -81,6 +107,8 @@ program VlaPlEx
   real(dp):: phiConv
   real(dp):: phi0
   real(dp):: T_e
+  real(dp):: recombination_rate
+  real(dp):: ionization_rate
   ! real(dp):: phiF
   integer:: k
 
@@ -107,6 +135,10 @@ program VlaPlEx
   call TtoZ%init(TtoZ_file)
   TtoZne_file = 'Zave_values_per_Ne.csv'
   call TtoZne%init(TtoZne_file)
+  recombination_rate_file = 'dr.csv'
+  call recombination_rateCoeff%init(recombination_rate_file)
+  ionization_rate_file = 'ci.csv'
+  call ionization_rateCoeff%init(ionization_rate_file)
 
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
@@ -184,6 +216,8 @@ program VlaPlEx
   T_i     = 0.0_dp
   n_e     = 0.0_dp
   T_e     = 0.0_dp
+  recombination_rate = 0.0_dp
+  ionization_rate = 0.0_dp
   Zave    = 0.0_dp
   Z_list  = (/ 6.0, 12.0 /)
   Zave_bc_old = 0.0_dp
@@ -277,9 +311,37 @@ program VlaPlEx
     f_i_old(:,:,1)  = 0.0_dp
     f_i_old(:,:,nv) = 0.0_dp
     sum_ni = 0.0_dp
+
+    ! --- Recombination Step --- !
+    do iz = 2, nz
+      do i = 2, nr
+        call recombination_rateCoeff%get(T_function(Temp_bc, n_e(1), n_e(i), gamma_e), Z_list(iz), recombination_rate)
+        ! Recombination of current charge state
+        f_i(iz,i,:) = f_i_old(iz,i,:) - f_i_old(iz,i,:) * dt * recombination_rate * n_e(i)
+
+        ! Transfer recombination density to lower charge state
+        f_i(iz - 1,i,:) = f_i_old(iz - 1,i,:) + f_i_old(iz,i,:) * dt * recombination_rate * n_e(i)
+      end do
+      f_i_old = f_i
+    end do
+
+
+    ! --- Ionization Step --- !
+    do iz = 1, nz - 1
+      do i = 2, nr
+        call ionization_rateCoeff%get(T_function(Temp_bc, n_e(1), n_e(i), gamma_e), Z_list(iz), ionization_rate)
+        ! Ionization of current charge state
+        f_i(iz,i,:) = f_i_old(iz,i,:) - f_i_old(iz,i,:) * dt * ionization_rate * n_e(i)
+
+        ! Transfer ionization density to higher charge state
+        f_i(iz + 1,i,:) = f_i_old(iz + 1,i,:) + f_i_old(iz,i,:) * dt * ionization_rate * n_e(i)
+      end do
+      f_i_old = f_i
+    end do
+
     ! Advect in the r direction
     do iz = 1, nz
-      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+      if (all(f_i(iz,:,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
         cycle
       end if
       !$omp parallel do
@@ -386,7 +448,7 @@ program VlaPlEx
     f_i_old = f_i
 
     do iz = 1, nz
-      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+      if (all(f_i(iz,:,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
         cycle
       end if
       ! Advect in the v direction