From f00a5c15c2c316e15363f6da88648f4fca4e3e11 Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Mon, 7 Apr 2025 13:35:06 +0200 Subject: [PATCH 1/4] Change Z string and reduce output files --- moduleOutput.f90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/moduleOutput.f90 b/moduleOutput.f90 index 1602bb0..b2d5fc7 100644 --- a/moduleOutput.f90 +++ b/moduleOutput.f90 @@ -145,8 +145,11 @@ module output character(len=10) :: ZString do j = 1, nz + if (all(n_i(j,i)*n_ref < 1.0e-16_dp)) then + cycle + end if write (timeString, formatTime) t - write(ZString, '(F6.3)') Z_list(j) + write(ZString, '(F6.1)') Z_list(j) ZString = adjustl(trim(ZString)) ZString = adjustl(ZString) @@ -278,8 +281,11 @@ module output do j = 1, nz + if (all(f(j,:)*n_ref/u_ref < 1.0e+20_dp)) then + cycle + end if write (timeString, formatTime) t - write(ZString, '(F6.3)') Z_list(j) + write(ZString, '(F6.1)') Z_list(j) ZString = adjustl(trim(ZString)) ZString = adjustl(ZString) From ff9047cee59acf87c41f35490d56ba3bab2db39a Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Mon, 7 Apr 2025 13:45:09 +0200 Subject: [PATCH 2/4] Revert Cycle changes, and use F4.1 string for Z --- moduleOutput.f90 | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/moduleOutput.f90 b/moduleOutput.f90 index b2d5fc7..3fbaa85 100644 --- a/moduleOutput.f90 +++ b/moduleOutput.f90 @@ -145,11 +145,8 @@ module output character(len=10) :: ZString do j = 1, nz - if (all(n_i(j,i)*n_ref < 1.0e-16_dp)) then - cycle - end if write (timeString, formatTime) t - write(ZString, '(F6.1)') Z_list(j) + write(ZString, '(F4.1)') Z_list(j) ZString = adjustl(trim(ZString)) ZString = adjustl(ZString) @@ -281,11 +278,8 @@ module output do j = 1, nz - if (all(f(j,:)*n_ref/u_ref < 1.0e+20_dp)) then - cycle - end if write (timeString, formatTime) t - write(ZString, '(F6.1)') Z_list(j) + write(ZString, '(F4.1)') Z_list(j) ZString = adjustl(trim(ZString)) ZString = adjustl(ZString) From 03b0619991c37a4789c6cbc3d22e34244ceb3369 Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Mon, 26 May 2025 10:31:54 +0200 Subject: [PATCH 3/4] Add recombination and ionization to model --- vlaplex.f90 | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 2 deletions(-) diff --git a/vlaplex.f90 b/vlaplex.f90 index 5031648..04d3163 100644 --- a/vlaplex.f90 +++ b/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 From 4d89e0e94289e31f92695978d649f158d1c7b06c Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Mon, 26 May 2025 10:44:36 +0200 Subject: [PATCH 4/4] Add python script for rate coefficients conversion --- scripts_python/rate_coef.py | 66 +++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 scripts_python/rate_coef.py diff --git a/scripts_python/rate_coef.py b/scripts_python/rate_coef.py new file mode 100644 index 0000000..4c908e6 --- /dev/null +++ b/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 +