From 052a4dc05e8cb907417bd670d1107556b2d982f6 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 8 Apr 2025 16:32:59 +0200 Subject: [PATCH] Corrections to tag v2.5 Small set of corrections to the tag v2.5. Includes changes to python scripts to plot data. Includes new 'fast' setup conditions that allow to output cases in an hour or so with still a good CFL condition and grid resolution. --- .gitignore | 1 + bc_80ns_T10.csv | 6 + bc_80ns_T10Z6.csv | 5 - bc_80ns_T30.csv | 6 + bc_80ns_T30Z11.csv | 5 - bc_80ns_T6.csv | 6 + bc_80ns_T60.csv | 6 + bc_80ns_T60Z16.csv | 5 - bc_80ns_T6Z4.csv | 5 - bc_fa_T30.csv | 3 + bc_fa_T30Z11.csv | 3 - moduleOutput.f90 | 61 ++-- scripts_python/TZ_data.json | 641 ++++++++++++++++++++++++++++++++++++ scripts_python/TtoZne.py | 4 +- scripts_python/plotBC.py | 9 +- scripts_python/plotCumF.py | 21 +- scripts_python/plotMom.py | 42 ++- scripts_python/readBC.py | 3 +- scripts_python/readF.py | 5 +- vlaplex.f90 | 60 ++-- 20 files changed, 790 insertions(+), 107 deletions(-) create mode 100644 bc_80ns_T10.csv delete mode 100644 bc_80ns_T10Z6.csv create mode 100644 bc_80ns_T30.csv delete mode 100644 bc_80ns_T30Z11.csv create mode 100644 bc_80ns_T6.csv create mode 100644 bc_80ns_T60.csv delete mode 100644 bc_80ns_T60Z16.csv delete mode 100644 bc_80ns_T6Z4.csv create mode 100644 bc_fa_T30.csv delete mode 100644 bc_fa_T30Z11.csv create mode 100644 scripts_python/TZ_data.json diff --git a/.gitignore b/.gitignore index 420618c..3320179 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ vlaplex *.csv *.mod *.o +*.tar.gz diff --git a/bc_80ns_T10.csv b/bc_80ns_T10.csv new file mode 100644 index 0000000..a77c2bf --- /dev/null +++ b/bc_80ns_T10.csv @@ -0,0 +1,6 @@ +t (s),n (m^-3),u (m s^-1),T (eV) +0.000000E-000,1.000000E+024,0.000000E+000,1.000000E+001 +8.000000E-008,1.000000E+024,0.000000E+000,1.000000E+001 +8.500000E-008,1.000000E+020,0.000000E+000,5.000000E-001 +1.000000E-007,1.000000E+010,0.000000E+000,5.000000E-001 +2.000000E-006,1.000000E+010,0.000000E+000,5.000000E-001 diff --git a/bc_80ns_T10Z6.csv b/bc_80ns_T10Z6.csv deleted file mode 100644 index 04bccba..0000000 --- a/bc_80ns_T10Z6.csv +++ /dev/null @@ -1,5 +0,0 @@ -t (s),n (m^-3),u (m s^-1),T (eV),Z -0.000000E-000,1.111111E+025,0.000000E+000,1.000000E+001,6.000000E+000 -8.000000E-008,1.111111E+025,0.000000E+000,1.000000E+001,6.000000E+000 -8.500000E-008,1.000000E+020,0.000000E+000,5.000000E+000,6.000000E+000 -2.000000E-006,1.000000E+020,0.000000E+000,5.000000E+000,6.000000E+000 diff --git a/bc_80ns_T30.csv b/bc_80ns_T30.csv new file mode 100644 index 0000000..554c029 --- /dev/null +++ b/bc_80ns_T30.csv @@ -0,0 +1,6 @@ +t (s),n (m^-3),u (m s^-1),T (eV) +0.000000E-000,1.000000E+024,0.000000E+000,3.000000E+001 +8.000000E-008,1.000000E+024,0.000000E+000,3.000000E+001 +8.500000E-008,1.000000E+020,0.000000E+000,5.000000E-001 +1.000000E-007,1.000000E+010,0.000000E+000,5.000000E-001 +2.000000E-006,1.000000E+010,0.000000E+000,5.000000E-001 diff --git a/bc_80ns_T30Z11.csv b/bc_80ns_T30Z11.csv deleted file mode 100644 index a23d865..0000000 --- a/bc_80ns_T30Z11.csv +++ /dev/null @@ -1,5 +0,0 @@ -t (s),n (m^-3),u (m s^-1),T (eV),Z -0.000000E-000,1.111111E+025,0.000000E+000,3.000000E+001,1.100000E+001 -8.000000E-008,1.111111E+025,0.000000E+000,3.000000E+001,1.100000E+001 -8.500000E-008,1.000000E+020,0.000000E+000,5.000000E+000,1.100000E+001 -2.000000E-006,1.000000E+020,0.000000E+000,5.000000E+000,1.100000E+001 diff --git a/bc_80ns_T6.csv b/bc_80ns_T6.csv new file mode 100644 index 0000000..0598227 --- /dev/null +++ b/bc_80ns_T6.csv @@ -0,0 +1,6 @@ +t (s),n (m^-3),u (m s^-1),T (eV) +0.000000E-000,1.000000E+024,0.000000E+000,6.000000E+000 +8.000000E-008,1.000000E+024,0.000000E+000,6.000000E+000 +8.500000E-008,1.000000E+020,0.000000E+000,5.000000E-001 +1.000000E-007,1.000000E+010,0.000000E+000,5.000000E-001 +2.000000E-006,1.000000E+010,0.000000E+000,5.000000E-001 diff --git a/bc_80ns_T60.csv b/bc_80ns_T60.csv new file mode 100644 index 0000000..a5603d1 --- /dev/null +++ b/bc_80ns_T60.csv @@ -0,0 +1,6 @@ +t (s),n (m^-3),u (m s^-1),T (eV) +0.000000E-000,1.000000E+024,0.000000E+000,6.000000E+001 +8.000000E-008,1.000000E+024,0.000000E+000,6.000000E+001 +8.500000E-008,1.000000E+020,0.000000E+000,5.000000E-001 +1.000000E-007,1.000000E+010,0.000000E+000,5.000000E-001 +2.000000E-006,1.000000E+010,0.000000E+000,5.000000E-001 diff --git a/bc_80ns_T60Z16.csv b/bc_80ns_T60Z16.csv deleted file mode 100644 index 63c934d..0000000 --- a/bc_80ns_T60Z16.csv +++ /dev/null @@ -1,5 +0,0 @@ -t (s),n (m^-3),u (m s^-1),T (eV),Z -0.000000E-000,1.111111E+025,0.000000E+000,6.000000E+001,1.600000E+001 -8.000000E-008,1.111111E+025,0.000000E+000,6.000000E+001,1.600000E+001 -8.500000E-008,1.000000E+020,0.000000E+000,5.000000E+000,1.600000E+001 -2.000000E-006,1.000000E+020,0.000000E+000,5.000000E+000,1.600000E+001 diff --git a/bc_80ns_T6Z4.csv b/bc_80ns_T6Z4.csv deleted file mode 100644 index f0bdb43..0000000 --- a/bc_80ns_T6Z4.csv +++ /dev/null @@ -1,5 +0,0 @@ -t (s),n (m^-3),u (m s^-1),T (eV),Z -0.000000E-000,1.111111E+025,0.000000E+000,6.000000E+000,4.000000E+000 -8.000000E-008,1.111111E+025,0.000000E+000,6.000000E+000,4.000000E+000 -8.500000E-008,1.000000E+020,0.000000E+000,5.000000E+000,4.000000E+000 -2.000000E-006,1.000000E+020,0.000000E+000,5.000000E+000,4.000000E+000 diff --git a/bc_fa_T30.csv b/bc_fa_T30.csv new file mode 100644 index 0000000..2e7a922 --- /dev/null +++ b/bc_fa_T30.csv @@ -0,0 +1,3 @@ +t (s),n (m^-3),u (m s^-1),T (eV) +0.000000E-000,1.000000E+024,0.000000E+000,3.000000E+001 +2.000000E-006,1.000000E+024,0.000000E+000,3.000000E+001 diff --git a/bc_fa_T30Z11.csv b/bc_fa_T30Z11.csv deleted file mode 100644 index 2c0c946..0000000 --- a/bc_fa_T30Z11.csv +++ /dev/null @@ -1,3 +0,0 @@ -t (s),n (m^-3),u (m s^-1),T (eV),Z -0.000000E-000,1.111111E+025,0.000000E+000,3.000000E+001,1.100000E+001 -1.000000E-006,1.111111E+025,0.000000E+000,3.000000E+001,1.100000E+001 diff --git a/moduleOutput.f90 b/moduleOutput.f90 index 1602bb0..bba4e1a 100644 --- a/moduleOutput.f90 +++ b/moduleOutput.f90 @@ -163,8 +163,8 @@ module output do i = 1, nr write(dataPhi_id, '('//formatFloat//',3('//formatSep //','//formatFloat//'))') & r(i)*L_ref, & - n_i(j,i)*n_ref, & - u_i(j,i)*u_ref, & + n_i(j,i)*n_ref, & + u_i(j,i)*u_ref, & T_i(j,i)*Temp_ref/ev_to_K end do @@ -173,13 +173,13 @@ module output end do end subroutine writeOutputMom - subroutine writeOutputBoundary(t, dt, n, u, Temp, TtoZ, Zinj) + subroutine writeOutputBoundary(t, dt, n, u, Temp, Zinj) use constantParameters, only: eV_to_K use referenceValues, only: t_ref, n_ref, u_ref, Temp_ref integer, intent(in):: t real(dp), intent(in):: dt - real(dp), intent(in):: n, u, Temp, TtoZ, Zinj + real(dp), intent(in):: n, u, Temp, Zinj character(len=6), parameter:: filename = 'bc.csv' logical:: res @@ -187,42 +187,43 @@ module output if (.not. res) then write (*, '(A, A)') 'Writing: ', filename open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') - write(dataBC_id, '(A,5(' // formatSep // ',A))') 't (s)', 'n (m^-3)', 'u (m s^-1)', 'T (eV)', 'TtoZ','Zinj' + write(dataBC_id, '(A,4(' // formatSep // ',A))') 't (s)', 'n (m^-3)', 'u (m s^-1)', 'T (eV)','Zinj' close(dataBC_id) end if open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') - write(dataBC_id, '(' // formatFloat // ',5('// formatSep // ',' // formatFloat // '))') & - t*dt*t_ref, n*n_ref, u*u_ref, Temp*Temp_ref/eV_to_K, TtoZ, Zinj + write(dataBC_id, '(' // formatFloat // ',4('// formatSep // ',' // formatFloat // '))') & + t*dt*t_ref, n*n_ref, u*u_ref, Temp*Temp_ref/eV_to_K, Zinj close(dataBC_id) end subroutine writeOutputBoundary - subroutine writeOutputTime(t, time, bins) - integer, intent(in):: t - real(dp), intent(in):: time - real(dp), intent(in):: bins - character(len=8), parameter:: filename = 'time.csv' - logical:: res - - inquire(file=pathOutput // filename, exist=res) - if (.not. res) then - write (*, '(A, A)') 'Writing: ', filename - open(unit=dataTime_id, file=pathOutput // filename, action='write', position='append') - write(dataTime_id, '(A,2(' // formatSep // ',A))') 'timestep', 'duration (s)', '#bins' - close(dataTime_id) - - end if - - open(unit=dataTime_id, file=pathOutput // filename, action='write', position='append') - write(dataTime_id, '(' // formatInt // ',2('// formatSep // ',' // formatFloat // '))') & - t, time, bins - - close(dataTime_id) - - end subroutine writeOutputTime +! JG: What is this procedure? +! subroutine writeOutputTime(t, time, bins) +! integer, intent(in):: t +! real(dp), intent(in):: time +! real(dp), intent(in):: bins +! character(len=8), parameter:: filename = 'time.csv' +! logical:: res +! +! inquire(file=pathOutput // filename, exist=res) +! if (.not. res) then +! write (*, '(A, A)') 'Writing: ', filename +! open(unit=dataTime_id, file=pathOutput // filename, action='write', position='append') +! write(dataTime_id, '(A,2(' // formatSep // ',A))') 'timestep', 'duration (s)', '#bins' +! close(dataTime_id) +! +! end if +! +! open(unit=dataTime_id, file=pathOutput // filename, action='write', position='append') +! write(dataTime_id, '(' // formatInt // ',2('// formatSep // ',' // formatFloat // '))') & +! t, time, bins +! +! close(dataTime_id) +! +! end subroutine writeOutputTime subroutine writeOutputZList(nz, Z_list) integer, intent(in):: nz diff --git a/scripts_python/TZ_data.json b/scripts_python/TZ_data.json new file mode 100644 index 0000000..b72c0cc --- /dev/null +++ b/scripts_python/TZ_data.json @@ -0,0 +1,641 @@ +{ + "now": "2025-04-04T14:21:35.000Z", + "program": "ZVView, ver. 2021-09-16", + "file": "Copy of plotted data", + "f":"y(x)", + "title": "Tin Average Charge State", + "xAxis": "temperature", + "yAxis": "", + "xUnits": "eV", + "yUnits": "", + "xMin": 0.271593, + "xMax": 184099, + "xScale": 1, + "yScale": 0, + "funcs": [ + { + "ifunc":0, + "fName":"Ne = 1.0E+12", + "numCol":1, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.032965} + ,{"x":1, "y":0.92671} + ,{"x":1.5, "y":1.6823} + ,{"x":2, "y":1.9381} + ,{"x":5, "y":3.3694} + ,{"x":7, "y":3.8735} + ,{"x":10, "y":4.1766} + ,{"x":15, "y":4.9651} + ,{"x":23, "y":6.2945} + ,{"x":32, "y":7.4784} + ,{"x":52, "y":9.4762} + ,{"x":74, "y":11.197} + ,{"x":100, "y":12.921} + ,{"x":165, "y":19.431} + ,{"x":235, "y":21.062} + ,{"x":310, "y":21.488} + ,{"x":390, "y":21.742} + ,{"x":475, "y":21.969} + ,{"x":655, "y":22.552} + ,{"x":845, "y":23.504} + ,{"x":1000, "y":24.553} + ,{"x":1441, "y":28.026} + ,{"x":1925, "y":31.867} + ,{"x":2454, "y":36.133} + ,{"x":3030, "y":38.155} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.155} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.442} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.403} + ] + } + ,{ + "ifunc":1, + "fName":"Ne = 1.0E+13", + "numCol":2, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.033255} + ,{"x":1, "y":1.0488} + ,{"x":1.5, "y":1.8974} + ,{"x":2, "y":1.984} + ,{"x":5, "y":3.4598} + ,{"x":7, "y":3.8851} + ,{"x":10, "y":4.18} + ,{"x":15, "y":4.9679} + ,{"x":23, "y":6.2958} + ,{"x":32, "y":7.4791} + ,{"x":52, "y":9.4765} + ,{"x":74, "y":11.197} + ,{"x":100, "y":12.921} + ,{"x":165, "y":19.431} + ,{"x":235, "y":21.062} + ,{"x":310, "y":21.488} + ,{"x":390, "y":21.742} + ,{"x":475, "y":21.969} + ,{"x":655, "y":22.552} + ,{"x":845, "y":23.504} + ,{"x":1000, "y":24.553} + ,{"x":1441, "y":28.026} + ,{"x":1925, "y":31.867} + ,{"x":2454, "y":36.133} + ,{"x":3030, "y":38.155} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.155} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.442} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.403} + ] + } + ,{ + "ifunc":2, + "fName":"Ne = 1.0E+14", + "numCol":3, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.074263} + ,{"x":1, "y":1.7619} + ,{"x":1.5, "y":1.9889} + ,{"x":2, "y":1.9987} + ,{"x":5, "y":3.6943} + ,{"x":7, "y":3.9281} + ,{"x":10, "y":4.2053} + ,{"x":15, "y":4.9938} + ,{"x":23, "y":6.3082} + ,{"x":32, "y":7.4852} + ,{"x":52, "y":9.479} + ,{"x":74, "y":11.199} + ,{"x":100, "y":12.922} + ,{"x":165, "y":19.431} + ,{"x":235, "y":21.063} + ,{"x":310, "y":21.488} + ,{"x":390, "y":21.742} + ,{"x":475, "y":21.969} + ,{"x":655, "y":22.552} + ,{"x":845, "y":23.504} + ,{"x":1000, "y":24.553} + ,{"x":1441, "y":28.026} + ,{"x":1925, "y":31.867} + ,{"x":2454, "y":36.133} + ,{"x":3030, "y":38.155} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.155} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.442} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.403} + ] + } + ,{ + "ifunc":3, + "fName":"Ne = 1.0E+15", + "numCol":4, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.067024} + ,{"x":1, "y":1.8546} + ,{"x":1.5, "y":1.9951} + ,{"x":2, "y":2.0151} + ,{"x":5, "y":3.9202} + ,{"x":7, "y":4.027} + ,{"x":10, "y":4.3596} + ,{"x":15, "y":5.1831} + ,{"x":23, "y":6.4113} + ,{"x":32, "y":7.5405} + ,{"x":52, "y":9.5012} + ,{"x":74, "y":11.211} + ,{"x":100, "y":12.928} + ,{"x":165, "y":19.432} + ,{"x":235, "y":21.063} + ,{"x":310, "y":21.488} + ,{"x":390, "y":21.742} + ,{"x":475, "y":21.969} + ,{"x":655, "y":22.552} + ,{"x":845, "y":23.505} + ,{"x":1000, "y":24.553} + ,{"x":1441, "y":28.027} + ,{"x":1925, "y":31.867} + ,{"x":2454, "y":36.133} + ,{"x":3030, "y":38.155} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.155} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.442} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.403} + ] + } + ,{ + "ifunc":4, + "fName":"Ne = 1.0E+16", + "numCol":5, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.012541} + ,{"x":1, "y":1.4079} + ,{"x":1.5, "y":1.9792} + ,{"x":2, "y":2.0845} + ,{"x":5, "y":3.9919} + ,{"x":7, "y":4.2166} + ,{"x":10, "y":4.914} + ,{"x":15, "y":5.8148} + ,{"x":23, "y":6.8933} + ,{"x":32, "y":7.8899} + ,{"x":52, "y":9.6814} + ,{"x":74, "y":11.315} + ,{"x":100, "y":12.98} + ,{"x":165, "y":19.435} + ,{"x":235, "y":21.063} + ,{"x":310, "y":21.488} + ,{"x":390, "y":21.742} + ,{"x":475, "y":21.97} + ,{"x":655, "y":22.553} + ,{"x":845, "y":23.506} + ,{"x":1000, "y":24.555} + ,{"x":1441, "y":28.028} + ,{"x":1925, "y":31.868} + ,{"x":2454, "y":36.134} + ,{"x":3030, "y":38.156} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.155} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.442} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.404} + ] + } + ,{ + "ifunc":5, + "fName":"Ne = 1.0E+17", + "numCol":6, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.0014986} + ,{"x":1, "y":0.44787} + ,{"x":1.5, "y":1.8165} + ,{"x":2, "y":2.0197} + ,{"x":5, "y":3.9944} + ,{"x":7, "y":4.4865} + ,{"x":10, "y":5.5905} + ,{"x":15, "y":6.7704} + ,{"x":23, "y":7.9337} + ,{"x":32, "y":8.891} + ,{"x":52, "y":10.474} + ,{"x":74, "y":11.947} + ,{"x":100, "y":13.325} + ,{"x":165, "y":19.456} + ,{"x":235, "y":21.065} + ,{"x":310, "y":21.49} + ,{"x":390, "y":21.745} + ,{"x":475, "y":21.974} + ,{"x":655, "y":22.563} + ,{"x":845, "y":23.521} + ,{"x":1000, "y":24.569} + ,{"x":1441, "y":28.035} + ,{"x":1925, "y":31.872} + ,{"x":2454, "y":36.136} + ,{"x":3030, "y":38.156} + ,{"x":3655, "y":38.94} + ,{"x":4331, "y":39.397} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.156} + ,{"x":6685, "y":40.618} + ,{"x":7585, "y":41.202} + ,{"x":8546, "y":41.919} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.441} + ,{"x":50000, "y":47.852} + ,{"x":100000, "y":48.404} + ] + } + ,{ + "ifunc":6, + "fName":"Ne = 1.0E+18", + "numCol":7, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":0.00019982} + ,{"x":1, "y":0.064341} + ,{"x":1.5, "y":0.85685} + ,{"x":2, "y":1.6966} + ,{"x":5, "y":3.9218} + ,{"x":7, "y":4.4706} + ,{"x":10, "y":5.8288} + ,{"x":15, "y":7.578} + ,{"x":23, "y":9.2131} + ,{"x":32, "y":10.343} + ,{"x":52, "y":12.203} + ,{"x":74, "y":13.329} + ,{"x":100, "y":14.264} + ,{"x":165, "y":19.576} + ,{"x":235, "y":21.076} + ,{"x":310, "y":21.503} + ,{"x":390, "y":21.77} + ,{"x":475, "y":22.018} + ,{"x":655, "y":22.657} + ,{"x":845, "y":23.649} + ,{"x":1000, "y":24.693} + ,{"x":1441, "y":28.092} + ,{"x":1925, "y":31.894} + ,{"x":2454, "y":36.144} + ,{"x":3030, "y":38.158} + ,{"x":3655, "y":38.941} + ,{"x":4331, "y":39.398} + ,{"x":5060, "y":39.77} + ,{"x":5844, "y":40.156} + ,{"x":6685, "y":40.619} + ,{"x":7585, "y":41.203} + ,{"x":8546, "y":41.92} + ,{"x":10000, "y":43.042} + ,{"x":20000, "y":46.437} + ,{"x":50000, "y":47.851} + ,{"x":100000, "y":48.403} + ] + } + ,{ + "ifunc":7, + "fName":"Ne = 1.0E+19", + "numCol":8, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":4.0615e-5} + ,{"x":1, "y":0.009011} + ,{"x":1.5, "y":0.12002} + ,{"x":2, "y":0.56798} + ,{"x":5, "y":3.4289} + ,{"x":7, "y":4.0527} + ,{"x":10, "y":5.4859} + ,{"x":15, "y":7.6474} + ,{"x":23, "y":10.06} + ,{"x":32, "y":11.848} + ,{"x":52, "y":13.767} + ,{"x":74, "y":14.563} + ,{"x":100, "y":16.174} + ,{"x":165, "y":19.984} + ,{"x":235, "y":21.158} + ,{"x":310, "y":21.628} + ,{"x":390, "y":22.016} + ,{"x":475, "y":22.426} + ,{"x":655, "y":23.377} + ,{"x":845, "y":24.505} + ,{"x":1000, "y":25.5} + ,{"x":1441, "y":28.514} + ,{"x":1925, "y":32.066} + ,{"x":2454, "y":36.197} + ,{"x":3030, "y":38.169} + ,{"x":3655, "y":38.945} + ,{"x":4331, "y":39.4} + ,{"x":5060, "y":39.773} + ,{"x":5844, "y":40.159} + ,{"x":6685, "y":40.622} + ,{"x":7585, "y":41.206} + ,{"x":8546, "y":41.923} + ,{"x":10000, "y":43.04} + ,{"x":20000, "y":46.395} + ,{"x":50000, "y":47.839} + ,{"x":100000, "y":48.401} + ] + } + ,{ + "ifunc":8, + "fName":"Ne = 1.0E+20", + "numCol":9, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":2.139e-5} + ,{"x":1, "y":0.0019556} + ,{"x":1.5, "y":0.018338} + ,{"x":2, "y":0.073538} + ,{"x":5, "y":2.1985} + ,{"x":7, "y":3.1813} + ,{"x":10, "y":4.3789} + ,{"x":15, "y":6.5863} + ,{"x":23, "y":9.5938} + ,{"x":32, "y":12.154} + ,{"x":52, "y":14.35} + ,{"x":74, "y":16.437} + ,{"x":100, "y":18.34} + ,{"x":165, "y":20.65} + ,{"x":235, "y":21.779} + ,{"x":310, "y":22.783} + ,{"x":390, "y":23.65} + ,{"x":475, "y":24.425} + ,{"x":655, "y":25.796} + ,{"x":845, "y":27.004} + ,{"x":1000, "y":27.888} + ,{"x":1441, "y":30.187} + ,{"x":1925, "y":33.06} + ,{"x":2454, "y":36.547} + ,{"x":3030, "y":38.245} + ,{"x":3655, "y":38.972} + ,{"x":4331, "y":39.415} + ,{"x":5060, "y":39.785} + ,{"x":5844, "y":40.172} + ,{"x":6685, "y":40.638} + ,{"x":7585, "y":41.225} + ,{"x":8546, "y":41.941} + ,{"x":10000, "y":43.041} + ,{"x":20000, "y":46.243} + ,{"x":50000, "y":47.781} + ,{"x":100000, "y":48.391} + ] + } + ,{ + "ifunc":9, + "fName":"Ne = 1.0E+21", + "numCol":10, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":8.6796e-5} + ,{"x":1, "y":0.0012015} + ,{"x":1.5, "y":0.0061077} + ,{"x":2, "y":0.017414} + ,{"x":5, "y":1.0162} + ,{"x":7, "y":2.0528} + ,{"x":10, "y":3.0096} + ,{"x":15, "y":4.8753} + ,{"x":23, "y":7.6568} + ,{"x":32, "y":10.348} + ,{"x":52, "y":13.822} + ,{"x":74, "y":16.818} + ,{"x":100, "y":19.1} + ,{"x":165, "y":22.57} + ,{"x":235, "y":24.342} + ,{"x":310, "y":25.709} + ,{"x":390, "y":26.883} + ,{"x":475, "y":27.889} + ,{"x":655, "y":29.46} + ,{"x":845, "y":30.595} + ,{"x":1000, "y":31.34} + ,{"x":1441, "y":33.437} + ,{"x":1925, "y":35.945} + ,{"x":2454, "y":37.71} + ,{"x":3030, "y":38.614} + ,{"x":3655, "y":39.138} + ,{"x":4331, "y":39.525} + ,{"x":5060, "y":39.883} + ,{"x":5844, "y":40.28} + ,{"x":6685, "y":40.767} + ,{"x":7585, "y":41.379} + ,{"x":8546, "y":42.105} + ,{"x":10000, "y":43.18} + ,{"x":20000, "y":46.152} + ,{"x":50000, "y":47.725} + ,{"x":100000, "y":48.381} + ] + } + ,{ + "ifunc":10, + "fName":"Ne = 1.0E+22", + "numCol":11, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":1.0393} + ,{"x":1, "y":1.2424} + ,{"x":1.5, "y":1.4551} + ,{"x":2, "y":1.6071} + ,{"x":5, "y":2.1738} + ,{"x":7, "y":2.3232} + ,{"x":10, "y":2.6442} + ,{"x":15, "y":3.7366} + ,{"x":23, "y":5.8249} + ,{"x":32, "y":8.0823} + ,{"x":52, "y":11.952} + ,{"x":74, "y":14.41} + ,{"x":100, "y":17.488} + ,{"x":165, "y":24.02} + ,{"x":235, "y":27.042} + ,{"x":310, "y":29.179} + ,{"x":390, "y":30.612} + ,{"x":475, "y":31.573} + ,{"x":655, "y":33.326} + ,{"x":845, "y":34.893} + ,{"x":1000, "y":35.843} + ,{"x":1441, "y":37.5} + ,{"x":1925, "y":38.396} + ,{"x":2454, "y":38.942} + ,{"x":3030, "y":39.341} + ,{"x":3655, "y":39.695} + ,{"x":4331, "y":40.066} + ,{"x":5060, "y":40.5} + ,{"x":5844, "y":41.033} + ,{"x":6685, "y":41.673} + ,{"x":7585, "y":42.387} + ,{"x":8546, "y":43.107} + ,{"x":10000, "y":43.994} + ,{"x":20000, "y":46.237} + ,{"x":50000, "y":47.722} + ,{"x":100000, "y":48.387} + ] + } + ,{ + "ifunc":11, + "fName":"Ne = 1.0E+23", + "numCol":12, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":4} + ,{"x":1, "y":4} + ,{"x":1.5, "y":3.1272} + ,{"x":2, "y":3.1804} + ,{"x":5, "y":3.4601} + ,{"x":7, "y":3.5987} + ,{"x":10, "y":3.805} + ,{"x":15, "y":4.3015} + ,{"x":23, "y":5.5703} + ,{"x":32, "y":7.3753} + ,{"x":52, "y":10.315} + ,{"x":74, "y":12.332} + ,{"x":100, "y":14.468} + ,{"x":165, "y":19.988} + ,{"x":235, "y":26.448} + ,{"x":310, "y":29.602} + ,{"x":390, "y":31.778} + ,{"x":475, "y":33.775} + ,{"x":655, "y":36.183} + ,{"x":845, "y":37.386} + ,{"x":1000, "y":37.966} + ,{"x":1441, "y":38.913} + ,{"x":1925, "y":39.567} + ,{"x":2454, "y":40.176} + ,{"x":3030, "y":40.815} + ,{"x":3655, "y":41.501} + ,{"x":4331, "y":42.217} + ,{"x":5060, "y":42.923} + ,{"x":5844, "y":43.575} + ,{"x":6685, "y":44.14} + ,{"x":7585, "y":44.609} + ,{"x":8546, "y":44.989} + ,{"x":10000, "y":45.409} + ,{"x":20000, "y":46.694} + ,{"x":50000, "y":47.853} + ,{"x":100000, "y":48.44} + ] + } + ,{ + "ifunc":12, + "fName":"Ne = 1.0E+24", + "numCol":13, + "typSoed":1, + "typUsi":-1, + "lx":36, + "i0x":0, + "i1x":35, + "pts":[ + {"x":0.5, "y":8} + ,{"x":1, "y":8} + ,{"x":1.5, "y":8} + ,{"x":2, "y":8.0001} + ,{"x":5, "y":8.0069} + ,{"x":7, "y":8.0197} + ,{"x":10, "y":8.0501} + ,{"x":15, "y":8.1262} + ,{"x":23, "y":8.3073} + ,{"x":32, "y":7.9918} + ,{"x":52, "y":9.3773} + ,{"x":74, "y":11.177} + ,{"x":100, "y":12.915} + ,{"x":165, "y":16.245} + ,{"x":235, "y":23.165} + ,{"x":310, "y":26.987} + ,{"x":390, "y":29.485} + ,{"x":475, "y":31.555} + ,{"x":655, "y":35.008} + ,{"x":845, "y":36.799} + ,{"x":1000, "y":37.754} + ,{"x":1441, "y":39.94} + ,{"x":1925, "y":41.496} + ,{"x":2454, "y":42.632} + ,{"x":3030, "y":43.504} + ,{"x":3655, "y":44.197} + ,{"x":4331, "y":44.779} + ,{"x":5060, "y":45.238} + ,{"x":5844, "y":45.603} + ,{"x":6685, "y":45.9} + ,{"x":7585, "y":46.147} + ,{"x":8546, "y":46.365} + ,{"x":10000, "y":46.634} + ,{"x":20000, "y":47.448} + ,{"x":50000, "y":48.101} + ,{"x":100000, "y":48.558} + ] + } + ] +} \ No newline at end of file diff --git a/scripts_python/TtoZne.py b/scripts_python/TtoZne.py index cd5a329..4e7b860 100644 --- a/scripts_python/TtoZne.py +++ b/scripts_python/TtoZne.py @@ -22,9 +22,9 @@ for func in data["funcs"]: T_values[T] = [] T_values[T].append(Z) + T_sorted = sorted(T_values.keys()) -Z_avg = [np.mean(T_values[T]) for x in T_sorted] -# y_std = [np.std(x_values[x]) for x in x_sorted] +Z_avg = [np.mean(T_values[T]) for T in T_sorted] # Save to CSV diff --git a/scripts_python/plotBC.py b/scripts_python/plotBC.py index 05173b5..445f3b9 100644 --- a/scripts_python/plotBC.py +++ b/scripts_python/plotBC.py @@ -4,15 +4,14 @@ import numpy as np import readBC -fileBC = glob.glob('../2025-02-10_11.14.59/bc.csv') -time, n, u, T, TtoZ, Zinj = readBC.read(fileBC[0]) +paths = ['../polytropic_80ns_T30/'] +time, n, u, T, Zinj = readBC.read(paths[0] + 'bc.csv') fig, ax = plt.subplots() -plt.plot(time, n / n[0] , label = f"$n_i$ ($\\times {n[0] * 1e-6} \\; cm^{{-3}})$") -plt.plot(time, T / T[0], label = f"$T \\; (\\times{T[0]} \\; eV)$") -plt.plot(time, TtoZ, label = "$Z$") +plt.plot(time, n / n[0], label = f"$n_i$ ($\\times {n[0] * 1e-6:.0e} \\; cm^{{-3}})$") +plt.plot(time, T / T[0], label = f"$T \\; (\\times{T[0]:.1f} \\; eV)$") plt.plot(time, Zinj, label = "Injection species") plt.semilogy() plt.legend() diff --git a/scripts_python/plotCumF.py b/scripts_python/plotCumF.py index 8a72112..3df5eb3 100644 --- a/scripts_python/plotCumF.py +++ b/scripts_python/plotCumF.py @@ -1,5 +1,6 @@ import readPhi import readF +import readZlist import matplotlib.pyplot as plt import glob import numpy as np @@ -13,21 +14,23 @@ m_i = 1.9712e-25 # paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/', '../quasiNeutral_partialAblation/', '../Poisson_partialAblation/'] # paths = ['../2024-12-02_21.07.52/', '../Poisson_50ns_T30Z11/'] # paths = ['../2024-12-10_18.45.17/', '../Poisson_50ns_T30Z11/'] -paths = ['../2024-12-11_12.38.27/', '../Poisson_polytropic_fa_T30Z11/', '../Poisson_fa_T30Z11/'] +# paths = ['../2024-12-11_12.38.27/', '../Poisson_polytropic_fa_T30Z11/', '../Poisson_fa_T30Z11/'] # paths = ['../Poisson_partialAblation/','../Poisson_partialAblation_lowerT/','../Poisson_partialAblation_lowT/','../Poisson_partialAblation_highT/'] +paths = ['../polytropic_80ns_T60/'] labels = [path[3:-1] for path in paths] for path, label in zip(paths, labels): + Zlist = readZlist.read(path+'ZList.csv') filesCum_i = sorted(glob.glob(path+'time_*_fCum_i.csv')) - # start = 0 - # end = len(filesCum_i) - # every = 20 - # for fileCum_i in filesCum_i[start:end+1:every]: - # time, x, v, f_i = readF.read(fileCum_i) - # plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label='t = {:.1f} ns'.format(time*1e9)) + _, _, v, _ = readF.read(filesCum_i[-1]) + sumF = np.zeros(len(v)) - time, x, v, f_i = readF.read(filesCum_i[-1]) - plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label=label) + for Z in Zlist: + filesCum_i = sorted(glob.glob(path+'time_*_Z{:.0f}000_fCum_i.csv'.format(Z))) + time, x, v, f_i = readF.read(filesCum_i[-1]) + sumF += f_i[0] + plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label=Z) + plt.plot(v**2*m_i*0.5/e, sumF*e/m_i/v, label='sum', color='k', linestyle='dashed') plt.yscale('log') plt.ylim([1e16,5e27]) diff --git a/scripts_python/plotMom.py b/scripts_python/plotMom.py index f980681..23e3648 100644 --- a/scripts_python/plotMom.py +++ b/scripts_python/plotMom.py @@ -1,6 +1,7 @@ import readPhi import readMom import readF +import readZlist import matplotlib.pyplot as plt import glob import numpy as np @@ -11,29 +12,44 @@ from scipy.constants import e, k # paths = ['../quasiNeutral_partialAblation/','../Poisson_partialAblation/'] # paths = ['../2024-10-02_14.30.44/'] # paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/'] -paths = ['../2024-12-10_18.45.17/'] +paths = ['../2025-04-08_09.36.52/'] labels = [path[3:-1] for path in paths] for path, label in zip(paths, labels): + Zlist = readZlist.read(path+'ZList.csv') filesPhi = sorted(glob.glob(path+'time_*_phi.csv')) - filesMom_i = sorted(glob.glob(path+'time_*_mom_i.csv')) - start = 50 - end = len(filesMom_i) - every = 20 + start = 0 + end = 20#len(filesPhi) + every = 5 fig, ax = plt.subplots(4, sharex='all') - for fileMom_i, filePhi in zip(filesMom_i[start:end+1:every], filesPhi[start:end+1:every]): - time, r, phi, E, n_e = readPhi.read(filePhi) - time, r, n_i, u_i, T_i, Zave = readMom.read(fileMom_i) + ax[1].set_yscale('log') + ax[1].set_ylim(bottom=1e10, top=1e24) + _, r, _, _, _ = readPhi.read(filesPhi[0]) + for t in range(start,end+1,every): + sum_Zni = np.zeros(len(r)) + ave_ui = np.zeros(len(r)) + ave_Ti = np.zeros(len(r)) + for Z in Zlist: + filesMom_i = sorted(glob.glob(path+'time_*_Z{:.0f}000_mom_i.csv'.format(Z))) + fileMom_i = filesMom_i[t] + time, r, n_i, u_i, T_i, Zave = readMom.read(fileMom_i) + sum_Zni += Zave*n_i + ave_ui += Zave*n_i*u_i + ave_Ti += Zave*n_i*T_i + + ave_ui = np.divide(ave_ui, sum_Zni, out=np.zeros_like(ave_ui), where=sum_Zni!=0.0) + ave_Ti = np.divide(ave_Ti, sum_Zni, out=np.zeros_like(ave_Ti), where=sum_Zni!=0.0) + + filePhi = filesPhi[t] + time, r, phi, E, n_e = readPhi.read(filePhi) ax[0].plot(r, phi, label='t = {:.1f} ns'.format(time*1e9)) # ax[0].plot(r, E, label='t = {:.1f} ns'.format(time*1e9)) # ax[0].plot(r, (Zave*n_i - n_e), label='t = {:.1f} ns'.format(time*1e9)) - ax[1].set_yscale('log') - # ax[1].set_ylim(bottom=1e16) - ax[1].plot(r, Zave*n_i) + ax[1].plot(r, sum_Zni) ax[1].plot(r, n_e, color='k', linestyle='dashed') - ax[2].plot(r, u_i) - ax[3].plot(r, T_i) + ax[2].plot(r, ave_ui) + ax[3].plot(r, ave_Ti) ax[0].set_title(label) ax[0].legend() diff --git a/scripts_python/readBC.py b/scripts_python/readBC.py index 10be123..b3de37f 100644 --- a/scripts_python/readBC.py +++ b/scripts_python/readBC.py @@ -7,9 +7,8 @@ def read(filename): n = df['n (m^-3)'].to_numpy() u = df['u (m s^-1)'].to_numpy() T = df['T (eV)'].to_numpy() - TtoZ = df['TtoZ'].to_numpy() Zinj = df['Zinj'].to_numpy() - return time, n, u, T, TtoZ, Zinj + return time, n, u, T, Zinj diff --git a/scripts_python/readF.py b/scripts_python/readF.py index 7fe3c15..0c8cdf6 100644 --- a/scripts_python/readF.py +++ b/scripts_python/readF.py @@ -5,9 +5,12 @@ def read(filename): df = pandas.read_csv(filename,skiprows=0,nrows=1) time = df['t (s)'].to_numpy()[0] + df = pandas.read_csv(filename,skiprows=2,nrows=1) + Z = df['Z'].to_numpy()[0] + df = pandas.read_csv(filename,skiprows=2,nrows=1,header=None) x = df.to_numpy()[0][1:] - df = pandas.read_csv(filename,skiprows=3,header=None) + df = pandas.read_csv(filename,skiprows=5,header=None) f = [] for col in df: if col == 0: diff --git a/vlaplex.f90 b/vlaplex.f90 index 5031648..f22597d 100644 --- a/vlaplex.f90 +++ b/vlaplex.f90 @@ -37,6 +37,7 @@ program VlaPlEx real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons real(dp), parameter:: gamma_e_exp = 1.0_dp /(gamma_e - 1.0_dp) ! Exponent for polytropic electrons real(dp), parameter:: gamma_e_dexp = (2.0_dp - gamma_e)/(gamma_e - 1.0_dp) ! Exponent for polytropic db_dphi + real(dp), parameter:: n_epsilon = 1.0e-16_dp real(dp):: r0, rf real(dp), allocatable, dimension(:):: r @@ -46,6 +47,7 @@ program VlaPlEx real(dp):: time real(dp):: dr, dv, dt integer:: nr, nv, nt, nz + integer:: nzMin, nzMax integer:: i, iz, j, t, z_inj integer:: j0 ! First integer of positive velocity @@ -89,7 +91,7 @@ program VlaPlEx integer:: rCum_index ! Set number of threads - call omp_set_num_threads(8) + call omp_set_num_threads(16) ! Set reference numbers (in SI units) Temp_ref = 30.0_dp * eV_to_K @@ -126,9 +128,9 @@ program VlaPlEx ! Index for cumulative sum rCum_index = minloc(abs(r - rCum), 1) - v0 =-1.0e1_dp*c_s - vf = 2.0e1_dp*c_s - dv = 1.0e-1_dp + v0 =-0.5e1_dp*c_s + vf = 1.0e1_dp*c_s + dv = 2.0e-1_dp nv = nint((vf - v0) / dv) + 1 dv = (vf - v0) / float(nv-1) @@ -146,7 +148,7 @@ program VlaPlEx t0 = 0.0_dp tf = 2.0e-7_dp / t_ref ! tf = 1.0e1_dp * (rf - r0) / c_s - dt = 1.0e-2_dp*dr/c_s + dt = 5.0e-2_dp*dr/c_s nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) @@ -164,7 +166,10 @@ program VlaPlEx write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr - nz = 2 + nzMin = 3 + nzMax = 14 + nz = nzMax - nzMin + 1 + nz = nz + 1 ! Add bin for low Z plasma ! Allocate vectors allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv)) allocate(n_i(1:nz,1:nr)) @@ -178,14 +183,17 @@ program VlaPlEx f_i = 0.0_dp f_i_old = 0.0_dp n_i = 0.0_dp - sum_ni = 0.0_dp + sum_ni = 0.0_dp u_i = 0.0_dp E_i = 0.0_dp T_i = 0.0_dp n_e = 0.0_dp T_e = 0.0_dp Zave = 0.0_dp - Z_list = (/ 6.0, 12.0 /) + Z_list(1) = 1.0_dp ! Low Z bin + do iz = nzMin, nzMax + Z_list(iz-nzMin+1+1) = float(iz) + end do Zave_bc_old = 0.0_dp phi = 0.0_dp phi_old = 0.0_dp @@ -247,18 +255,24 @@ program VlaPlEx ! Main loop do t = 1, nt time = t * dt + t0 + + ! Find new \bar{Z}_i based on T and density call boundaryConditions%get(time, n_bc, u_bc, Temp_bc) + ! Reset previous value + Zave_bc_old = 0.0_dp + ! Initial guess based on average table call TtoZ%get(Temp_bc, Zave_bc) z_inj = minloc(abs(Z_list - Zave_bc),1) - Zave_bc = Z_list(z_inj) + Zave_bc = Z_list(z_inj) + ! Start iterative process based on T, n_e table do while (Zave_bc - Zave_bc_old > 0.1_dp) Zave_bc_old = Zave_bc call TtoZne%get(Temp_bc, Zave_bc * n_bc, Zave_bc) z_inj = minloc(abs(Z_list - Zave_bc),1) - Zave_bc = Z_list(z_inj) + Zave_bc = Z_list(z_inj) end do - call writeOutputBoundary(t, dt, Zave_bc*n_bc, u_bc, Temp_bc, Zave_bc, Zave_bc) u_bc = sqrt(Zave_bc * Temp_bc) + call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc) ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) @@ -279,7 +293,7 @@ program VlaPlEx sum_ni = 0.0_dp ! 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(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then cycle end if !$omp parallel do @@ -296,13 +310,15 @@ program VlaPlEx end if n_i(iz,i) = sum(f_i(iz,i,:))*dv - if (n_i(iz,i) > 1.0e-10_dp) then - u_i(iz,i) = sum(v(:) *f_i(iz,i,:))*dv / n_i(iz,i) - E_i(i) = sum(v(:)**2*f_i(iz,i,:))*dv / n_i(iz,i) - T_i(iz,i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(iz,i)**2 + if (n_i(iz,i) > n_epsilon) then + u_i(iz,i) = sum(v(:) *f_i(iz,i,:))*dv / n_i(iz,i) + E_i(i) = sum(v(:)**2*f_i(iz,i,:))*dv / n_i(iz,i) + T_i(iz,i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(iz,i)**2 else - u_i(iz,i) = 0.0_dp - T_i(iz,i) = 0.0_dp + f_i(iz,i,:) = 0.0_dp + n_i(iz,i) = 0.0_dp + u_i(iz,i) = 0.0_dp + T_i(iz,i) = 0.0_dp end if end do @@ -360,7 +376,7 @@ program VlaPlEx end if ! ! Calculate new potential to ensure 0 current at the edge - ! if (n_i(nr) > 1.0e-10_dp) then + ! if (n_i(nr) > n_epsilon) then ! phiF = phi0 + T_e * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_e/m_e))) ! ! else @@ -386,7 +402,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(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then cycle end if ! Advect in the v direction @@ -423,14 +439,14 @@ program VlaPlEx ! Reset values for next iteration 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(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then cycle end if fCum_i(iz,:) = fCum_i(iz,:) + f_i_old(iz,rCum_index,:) end do ! Write output if (mod(t,everyOutput) == 0 .or. t == nt) then - call writeOutputF(t, dt, nz, nr, r, nv, v, f_i_old, Z_list) + ! call writeOutputF(t, dt, nz, nr, r, nv, v, f_i_old, Z_list) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list) call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list)