Hi Julia community,
My name is Nikolaj, I am not an expert on Julia, but I really like the syntax and speed of Julia.
I am using Julia to determine the wavelength dependent absorption coefficient of CO2, based on line-by-line data HITEMP2010 (calculated from quantum mechanics) available from HITRANonline data: HITEMP .
I am trying to read data from a file in a multithreading loop. The file consists of many lines of 160 characters each. After reading one line, I need to index into specific entries in the string. When I use only 1 thread I get the same result every time. When I use 2 or more threads, I get a different result for each run. I suspect that this happens because the different threads get ‘mixed up’ and access each others strings?
I have attached a minimal reproducible example (you also need the file “02_2250-2500_HITEMP2010.par” to run the example which can be downloaded from Index of /hitemp/data/HITEMP-2010/CO2_line_list).
using Plots
pyplot()
using LaTeXStrings
# set the number of threads
nthreads = 1 # [-] number of threads to use in parallel
# inputs
T = 1200.0 # [K] temperature
P_bar = 1.0 # [bar] pressure
P_atm = P_bar/1.01325 # [atm] pressure
P_Pa = P_atm*101325.0 # [Pa] pressure
CO2_molFrac = 0.1 # [-] mole fraction of CO2
nu_lower = 2300.0 # [1/cm] lower limit on wavenumber
nu_upper = 2305.0 # [1/cm] upper limit on wavenumber
nu_step = 0.01 # [1/cm] stepsize for wavenumber range
# initial calculations
nu_range = nu_lower:nu_step:nu_upper # [1/cm] wavenumber range
PpartCO2_atm = P_atm*CO2_molFrac # [atm] partial pressure of CO2 in atm
PpartCO2_Pa = P_Pa*CO2_molFrac # [Pa] partial pressure of CO2 in Pa
T_ref = 296.0 # [K] reference temperature of HITEMP2010
QT_ref_CO2 = 2.8609E+02 # QT at 296 K for carbon dioxide
c2 = 1.4387769 # [cm*K] constant for calculations
V_cm3 = 1e-6 # [m^3] = 1 cm^3
R = 8.314 # [J/mol-K] universal gas constant
mol_per_cm3 = PpartCO2_Pa*V_cm3/(R*T) # (mol*cm^-3) moles of CO2 per cm^3
molecules_per_cm3 = mol_per_cm3*6.02214129e23 # (molecules*cm^-3) molecules per cm^3
QT_value_CO2 = 4738.8 # The Q(1200) value has been looked up for this case
# connect to the spectral data of carbon dioxide
f_CO2 = Vector{IOStream}() # create the variable for storage
f_CO2 = push!(f_CO2, open("02_2250-2500_HITEMP2010.par")) # HITEMP2010 data file
# read the data from the file
rangeCO2 = String[] # empty string
templines = readlines(f_CO2[1]) # get data
for n = eachindex(templines) # loop over lines
push!(rangeCO2, templines[n]) # put the lines into the CO2 range variable
end
# get the range (these values are really looked up)
index_low = 219893 # the lower index of the data we're interested in
index_upper = 234622 # the upper index of the data we're interested in
index_rangeWidth = index_upper-index_low # the range that we need
# get ready for parallel for-loop
widthNumber = 80.0 # the number of halfwidths to consider on each side of peak
kappaCO2_per_cm = zeros(length(nu_range)) # preallocate
sumKappaCO2 = zeros(length(nu_range)) # preallocate
indexCount = zeros(Int, nthreads) # preallocate
sumKappaCO2 = zeros(length(nu_range), nthreads) # preallocate
kappaCO2_per_cm = zeros(length(nu_range)) # preallocate
# enter the parallel for-loop
Threads.@threads for logicalCore in 1:nthreads
# here I try to split the range into subranges
indexCount[logicalCore] = trunc(Int, index_low+(index_rangeWidth*((logicalCore-1)/nthreads)))
@label newRun # I use goto/label instead of for-loops here
indexCount[logicalCore] += 1 # increase counter for this thread
# now see if the current index for this thread is below the maximum for this thread
if indexCount[logicalCore] <= trunc(Int, index_low+(index_rangeWidth*(logicalCore/nthreads)))
line = rangeCO2[indexCount[logicalCore]] # get a single line from the data
# indices refer to positions in the HITEMP2010 string
vacuumWavenumber = parse(Float64, line[(2+1+1):(2+1+12)]) # 12
intensity = parse(Float64, line[(2+1+12+1):(2+1+12+10)]) # 10
airBroadenedHalfWidth = parse(Float64, line[(2+1+12+10+10+1):(2+1+12+10+10+5)]) # 5
selfBroadenedHalfWidth = parse(Float64, line[(2+1+12+10+10+5+1):(2+1+12+10+10+5+5)]) # 5
lowerStateEnergy = parse(Float64, line[(2+1+12+10+10+5+5+1):(2+1+12+10+10+5+5+10)]) # 10
temperatureDependenceCoef = parse(Float64, line[(2+1+12+10+10+5+5+10+1):(2+1+12+10+10+5+5+10+4)]) # 4
airPressureInducedLineShift = parse(Float64, line[(2+1+12+10+10+5+5+10+4+1):(2+1+12+10+10+5+5+10+4+8)]) # 8
# calculate the halfwidth
gamma = (T_ref/T)^temperatureDependenceCoef * (airBroadenedHalfWidth*(P_atm-PpartCO2_atm)+selfBroadenedHalfWidth*PpartCO2_atm)
# calculate the intensity
S = (intensity*(QT_ref_CO2/QT_value_CO2)*((exp(-c2*lowerStateEnergy/T))/(exp(-c2*lowerStateEnergy/T_ref)))*
((1-exp(-c2*vacuumWavenumber/T))/(1-exp(-c2*vacuumWavenumber/T_ref))))
# calculate the shifted peak position
waveNumberShifted = vacuumWavenumber+airPressureInducedLineShift*P_atm
nuCount = 1 # set the counter for wavenumber
@label newNu
if nuCount <= length(nu_range)
nu = nu_range[nuCount]
# now we check if we are within 80 halfwidths of the 'middle' wavenumber of the current index
# we only do the calculations if we are within this range
# otherwise we keep the value at zero, which should be sufficiently accurate (based on HITRAN Online)
if (nu >= (waveNumberShifted-(widthNumber*gamma))) && (nu <= (waveNumberShifted+(widthNumber*gamma)))
# this is really the spectral absorption cross section (1/(cm^-2*molecule))
kappaCO2 = (S/pi)*((gamma)/(gamma^2 +(nu-waveNumberShifted)^2))
kappaCO2_per_cm[nuCount] = molecules_per_cm3 *kappaCO2 # (molecule*cm^-2)*(1/(cm^-2*molecule)) = (1/cm)
elseif (nu > waveNumberShifted+(widthNumber*gamma))
@goto endOfRun
end
nuCount += 1
@goto newNu
end
@label endOfRun
sumKappaCO2[:,logicalCore] += kappaCO2_per_cm[:] # add vector to sum
@goto newRun
end
end
# sum the contributions for each thread to get the total
sumSumKappaCO2 = sum(sumKappaCO2, dims=2)
plot() # now plot the sum
plot!(nu_range, sumSumKappaCO2 ,legend=false)
xlabel!(latexstring("Wavenumber \n [cm\$^{-1}\$]"))
ylabel!(latexstring("Spectral Absorption Coefficient \n [cm\$^{-1}\$]"))
