From a309ccdca82c7595813d3a708b1994661c667923 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Tue, 21 Sep 2021 12:48:06 -0400 Subject: [PATCH] Update with optional effects due to area expansion (#4) * initial changes from PC * clean up calc_c1 * add new keywords; remove redundante statements * clean up some comments * update the readme * add a tests folder for quick checks * update readme * cleanup comments * last docstring cleanup --- README.md | 22 ++++++- calc_c1.pro | 36 +++++------ ebtel2.pro | 98 +++++++++++++++++++----------- tests/.gitignore | 138 +++++++++++++++++++++++++++++++++++++++++++ tests/README.md | 26 ++++++++ tests/conftest.py | 20 +++++++ tests/test_cases.py | 141 ++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 428 insertions(+), 53 deletions(-) create mode 100644 tests/.gitignore create mode 100644 tests/README.md create mode 100644 tests/conftest.py create mode 100644 tests/test_cases.py diff --git a/README.md b/README.md index 09b6ea7..1985d84 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # EBTEL + [![ascl:1203.007](https://img.shields.io/badge/ascl-1203.007-blue.svg?colorB=262255)](http://ascl.net/1203.007) +[![DOI](https://zenodo.org/badge/31337216.svg)](https://zenodo.org/badge/latestdoi/31337216) Enthalpy-Based Thermal Evolution of Loops (EBTEL) model coded in the IDL language by J. A. Klimchuk, S. Patsourakos, and P.J. Cargill. @@ -12,31 +14,41 @@ For more information regarding the EBTEL model see: See also [ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus), a C++ implementation of the EBTEL model that treats the electron and ion populations separately. Note that while the IDL version of EBTEL only solves the single fluid equations, there is little or no difference between the two-fluid and single-fluid models below roughly 5 MK. ## Word of Caution + In its integration scheme, EBTEL IDL uses a fixed timestep with a default value of 1 second. In order to accurately model a strong, impulsive heating event, it is necessary to resolve the timescales of the relevant physical processes, most notably thermal conduction, which can be on the order of milliseconds for a typical coronal plasma (e.g. [Bradshaw and Cargill, 2013](http://adsabs.harvard.edu/abs/2013ApJ...770...12B)). Thus the EBTEL IDL integration routine with the default timestep may under-resolve portions of the loop's evolution, particularly in the early heating and conductive cooling phases. -[ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus) addresses this issue by providing an adaptive timestep routine that ensures the timestep is always sufficiently small compared to the timescales of the relevant physical processes. ebtel++ also provides an additional correction to the _c_1 parameter during the conduction phase and treats the electrons and the ions separately (see appendices of [Barnes et al., 2016](http://adsabs.harvard.edu/abs/2016ApJ...829...31B) for more details). Both of these corrections are likely to be important in the case of strong, impulsive heating. Thus, for sufficiently short heating timescales and large heating rates, use of ebtel++ is recommended. +[ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus) addresses this issue by providing an adaptive timestep routine that ensures the timestep is always sufficiently small compared to the timescales of the relevant physical processes. This correction is likely to be important in the case of strong, impulsive heating. Thus, for sufficiently short heating timescales and large heating rates, use of ebtel++ is recommended. ## Terms of Use + This code is the authorized version of EBTEL dated January 2013. Updates will be made as and when necessary.This code is distributed as is. Modifications by the user are unsupported. If you believe you have encountered an error in the original (i.e. unaltered) code, create an issue or submit a pull request with the requested changes. Use of EBTEL should be acknowledged by referencing Papers 1 & 2 as listed above. ## Purpose + Compute the evolution of spatially-averaged (along the field) loop quantities using simplified equations. The instantaneous differential emission measure of the transition region is also computed. This version incorporates all the modifications from Paper 2 and is written in modular form for clarity. DEM parts unchanged except for TR pressure correction (see Paper 2). ## Inputs + + ttime = time array (s) + heat = heating rate array (erg cm^-3 s^-1) (direct heating only; the first element, heat(0), determines the initial static equilibrium) + length = loop half length (top of chromosphere to apex) (cm) ## Optional Keyword Inputs + + classical = set to use the UNsaturated classical heat flux + dem_old = set to use old technique of computing DEM(T) in the trans. reg. (weighted average of demev, demcon, and demeq) + flux_nt = energy flux array for nonthermal electrons impinging the chromosphere (input as a positive quantity; erg cm^-2 s^-1) + energy_nt = mean energy of the nonthermal electrons in keV (default is 50 keV) + rtv = set to use Rosner, Tucker, Vaiana radiative loss function (Winebarger form) ++ a_c = normalised area factor for the corona (average) as described in literature (see below). (default is 1) ++ a_0 = normalised area factor for the top of the TR as described in literature (see below). (default is 1) ++ a_tr = normalised area factor for the TR (average) as described in literature (see below). (default is 1) ++ l_fact = l_cor/length, coronal fraction of loop length recommended value is 0.85 (default is 1) ## Outputs + + t (t_a) = temperature array corresponding to time (avg. over coronal section of loop / apex) + n (n_a) = electron number density array (cm^-3) (coronal avg. / apex) + p (p_a) = pressure array (dyn cm^-2) (coronal avg. /apex) @@ -51,6 +63,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us + rad_cor = coronal radiative loss ## Correspondence with Variables in ApJ Articles + + Paper 1 + r1 = c_3 + r2 = c_2 @@ -61,6 +74,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us + dem_eq = DEM_se ## Usage + + To include the transition region DEM: + `IDL> ebtel2, ttime, heat, t, n, p, v, dem_tr, dem_cor, logtdem` + To exclude the transition region DEM (faster): @@ -71,6 +85,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us + `IDL> ebtel2, ttime, heat, t, n, p, v, dem_tr, dem_cor, logtdem, f_ratio, rad_ratio` ### Example Run + + Define the necessary parameters + `IDL> time = findgen(10000)` Define time array + `IDL> heat = fltarr(10000)` Define corresponding heating array @@ -109,6 +124,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us + `IDL> plot, logtdem, alog10(dem60), xtit='log T (K)',ytit='log DEM (cm!U-5!N K!U-1!N)', tit='1000-1059 s Integration', xran=[5.,7.], /ynoz` ## Intensities + For observations in which temperature response function, G(T), has units of DN s^-1 pix^-1 cm^5 and the loop diameter, d, is larger than the pixel dimension, l_pix: + I_cor_perp = d/(2L) * Int{G(T) * dem_cor(T) * dT} @@ -118,6 +134,7 @@ For observations in which temperature response function, G(T), has units of DN s for lines-of-sight perpendicular and parallel to the loop axis. I_tr_perp assumes that the transition region is thinner than l_pix. See `intensity_ebtel.pro` for more information. ## Miscellaneous Comments + + A 1 sec time step is generally adequate, but in applications where exceptionally strong conductive cooling is expected (e.g., intense short duration heating events, especially in short loops), a shorter time step may be necessary. If there is a question, users should compare runs with different time steps and verify that there are no significant differences in the results. + Runs much more quickly if the transition region DEM is not computed. + Speed can be increased by increasing the minimum DEM temperature from 10^4 to, say, 10^5 K or by decreasing the maximum DEM temperature from 10^8.5 to, say, 10^7.5 (search on 450 and 451). @@ -125,8 +142,11 @@ for lines-of-sight perpendicular and parallel to the loop axis. I_tr_perp assume + To have equal amounts of thermal and nonthermal heating: flux_nt = heat*length. + It is desirable to have a low-level background heating during the cooling phase so that the coronal temperature does not drop below values at which the corona DEM is invalid. + v = (c_2/c_3)(t_tr/t)v_0 = (r2/r1)(t_tr/t)(v/r4) at temperature t_tr in the transition region, where t is the average coronal temperature. ++ The non-uniform area version assumes A_c >= A_0 >= A_TR. There are limits on how large A_c can be. Too large area factors lead to violation of EBTEL assumptions, in particular that TR thickness is small compared to the loop length. An initial recommendation is A_c < 5. See Cargill et al, MNRAS, 2021 for more information. ## CHANGELOG + ++ September 2021, Modified to include non-uniform loop area (see above for caveats) + May 2012. PC version. Modular form. + 2013 Jan 15, JAK, Fixed a bug in the compution of the radiation function at 10^4 K; important for computing radiation losses based on dem_tr; ge vs. gt in computing rad; lt vs. le in computing rad_dem diff --git a/calc_c1.pro b/calc_c1.pro index 8c286b6..cf99a2a 100644 --- a/calc_c1.pro +++ b/calc_c1.pro @@ -1,13 +1,12 @@ pro calc_c1, temp, den, length, rad, c1 common params, k_b, mp, kappa_0 + common area, a0_ac, atr_ac, l_c, l_tr, l_star ; Calculate scale height calc_lambda, temp, sc - ; r3_rad_0: Radiative phase value, no gravity - ; r3_eqm_0: Value in eqm with no gravity, -2/3 loss power law. - ; - r3_rad_0 = 0.6 - r3_eqm_0 = 2. + + r3_rad_0 = 0.6 ; Radiative phase value, no gravity + r3_eqm_0 = 2. ; Value in eqm with no gravity, -2/3 loss power law. l_fact_eq = 5. l_fact_rad = 5. @@ -17,8 +16,8 @@ pro calc_c1, temp, den, length, rad, c1 ; Adjust values for gravity ; - r3_eqm_g = r3_eqm_0*exp(2*2*sin(3.14159/l_fact_eq)*length/3.14159/sc) - r3_radn_g = r3_rad_0*exp(2*2*sin(3.14159/l_fact_rad)*length/3.14159/sc) + r3_eqm_g = r3_eqm_0*exp(2*2*sin(3.14159/l_fact_eq)*l_c/3.14159/sc) + r3_radn_g = r3_rad_0*exp(2*2*sin(3.14159/l_fact_rad)*l_c/3.14159/sc) ; ; Adjust for loss function ; @@ -27,20 +26,21 @@ pro calc_c1, temp, den, length, rad, c1 ; These lines turn off loss function fix ; r3_eqm = r3_eqm_g ; r3_radn = r3_radn_g + ; ;Calculate over/under density - - n_eq_2 = kappa_0/3.5/r3_eqm/rad/length/length*(temp/r2)^(7./2.) + ; + n_eq_2 = (a0_ac/atr_ac)*kappa_0/3.5/r3_eqm/rad*(l_star/length/l_c^2)*(temp/r2)^(7./2.)/(1.-l_tr/l_c/r3_eqm) + noneq2=den^2./n_eq_2 - - ; Hot loops: use eqm C1 + + r3_cond = 6. + ; r3_cond = 2. ; uncomment this to turn off the conductive phase fix if (noneq2 lt 1.) then begin - r3=r3_eqm - endif - - ; Radiative loops: transition from eqm. - if (noneq2 ge 1.) then begin - r3=(2.*r3_eqm+r3_radn*(noneq2-1.))/(1+noneq2) - endif + r3=(2.*r3_eqm+r3_cond*(1./noneq2-1.))/(1+1./noneq2) + endif else begin + ; Radiative loops: transition from eqm. + r3=(2.*r3_eqm+r3_radn*(noneq2-1.))/(1+noneq2) + endelse ; Final value c1=r3 diff --git a/ebtel2.pro b/ebtel2.pro index 0bf14cd..55fb5c1 100644 --- a/ebtel2.pro +++ b/ebtel2.pro @@ -1,7 +1,8 @@ pro ebtel2, ttime, heat, length, t, n, p, v, ta, na, pa, c11, dem_tr, dem_cor, $ logtdem, f_ratio, rad_ratio, cond, rad_cor,$ classical=classical, dynamic=dynamic, dem_old=dem_old, $ - flux_nt=flux_nt, energy_nt=energy_nt, rtv=rtv + flux_nt=flux_nt, energy_nt=energy_nt, rtv=rtv, $ + a_c=a_c, a_0=a_0, a_tr=a_tr, l_fact=l_fact ; ; NAME: Enthalpy-Based Thermal Evolution of Loops (EBTEL) ; @@ -22,31 +23,39 @@ ; OPTIONAL KEYWORD INPUTS: ; classical = set to use the UNsaturated classical heat flux ; dynamic = set to use dynamical r1 and r2 (NOT recommended, especially when T > 10 MK). Now redundant. - ; dem_old = set to use old technique of computing DEM(T) in the trans. reg. + ; dem_old = set to use old technique of computing DEM(T) in the transition region ; (weighted average of demev, demcon, and demeq) ; flux_nt = energy flux array for nonthermal electrons impinging the chromosphere ; (input as a positive quantity; erg cm^-2 s^-1) ; energy_nt = mean energy of the nonthermal electrons in keV (default is 50 keV) ; rtv = set to use Rosner, Tucker, Vaiana radiative loss function (Winebarger form) + ; a_c = normalised area factor for the corona (average) as + ; described in literature (see below). (default is 1) + ; a_0 = normalised area factor for the top of the TR as + ; described in literature (see below). (default is 1) + ; a_tr = normalised area factor for the TR (average) as + ; described in literature (see below). (default is 1) + ; l_fact = l_cor/length, coronal fraction of loop length recommended value is 0.85 + ; (default is 1) ; ; OUTPUTS: - ; t (t_a) = temperature array corresponding to time (avg. over coronal section of loop / apex) - ; n (n_a) = electron number density array (cm^-3) (coronal avg. / apex) - ; p (p_a) = pressure array (dyn cm^-2) (coronal avg. /apex) - ; v = velocity array (cm s^-1) (r4 * velocity at base of corona) - ; c11 = C1 (or r3 in this code) - ; dem_tr = differential emission measure of transition region, dem(time,T), both legs - ; (dem = n^2 * ds/dT cm^-5 K^-1) - ; (Note: dem_tr is not reliable when a nonthermal electron flux is used.) - ; dem_cor = differential emission measure of corona, dem(time,T), both legs - ; (Int{dem_cor+dem_tr dT} gives the total emission measure of a - ; loop strand having a cross sectional area of 1 cm^2) - ; logtdem = logT array corresponding to dem_tr and dem_cor - ; f_ratio = ratio of heat flux to equilibrium heat flux - ; (ratio of heat flux to tr. reg. radiative loss rate) + ; t (t_a) = temperature array corresponding to time (avg. over coronal section of loop / apex) + ; n (n_a) = electron number density array (cm^-3) (coronal avg. / apex) + ; p (p_a) = pressure array (dyn cm^-2) (coronal avg. /apex) + ; v = velocity array (cm s^-1) (r4 * velocity at base of corona) + ; c11 = C1 (or r3 in this code) + ; dem_tr = differential emission measure of transition region, dem(time,T), both legs + ; (dem = n^2 * ds/dT cm^-5 K^-1) + ; (Note: dem_tr is not reliable when a nonthermal electron flux is used.) + ; dem_cor = differential emission measure of corona, dem(time,T), both legs + ; (Int{dem_cor+dem_tr dT} gives the total emission measure of a + ; loop strand having a cross sectional area of 1 cm^2) + ; logtdem = logT array corresponding to dem_tr and dem_cor + ; f_ratio = ratio of heat flux to equilibrium heat flux + ; (ratio of heat flux to tr. reg. radiative loss rate) ; rad_ratio = ratio of tr. reg. radiative loss rate from dem_tr and from r3*(coronal rate) - ; cond = conductive loss from corona - ; rad_cor = coronal radiative loss + ; cond = conductive loss from corona + ; rad_cor = coronal radiative loss ; ; CORRESPONDENCE WITH VARIABLES IN ASTROPHYSICAL JOURNAL ARTICLES: ; (Klimchuk et al., 2008, ApJ, 682, 1351; Cargill et al., ApJ 2012) @@ -100,9 +109,14 @@ ; the same time that r1 = 0.7 provides a more accurate radiative cooling. ; v = (c_2/c_3)*(t_tr/t)*v_0 = (r2/r1)*(t_tr/t)*(v/r4) at temperature t_tr in the transition ; region, where t is the average coronal temperature. + ; The non-uniform area version assumes A_c >= A_0 >= A_TR. There are limits on how large A_c + ; can be. Too large area factors lead to violation of EBTEL assumptions, in particular that TR + ; thickness is small compared to the loop length. An initial recommendation is A_c < 5. See + ; Cargill et al, MNRAS, 2021 for more information. ; ; HISTORY: + ; September 2021, Modified to include non-uniform loop area (see above for caveats) ; May 2012. PC version. Modular form. ; See original ebtel.pro for many additional comments. ; 2013 Jan 15, JAK, Fixed a bug in the compution of the radiation function at 10^4 K; @@ -110,7 +124,22 @@ ; ge vs. gt in computing rad; lt vs. le in computing rad_dem ; --------------------------------------------------------------------- common params, k_b, mp, kappa_0 - + common area, a0_ac, atr_ac, l_c, l_tr, l_star + + if not keyword_set(a_0) then a_0 = 1. + if not keyword_set(a_c) then a_c = 1. + if not keyword_set(a_tr) then a_tr = 1. + if not keyword_set(l_fact) then l_fact = 1. + + ; Area ratios + a0_ac=a_0/a_c + atr_ac=a_tr/a_c + + ; Modified lengths + l_c = l_fact*length + l_tr=length - l_c + l_star = l_c + l_tr*atr_ac + ntot = n_elements(ttime) ; Physical constants Can comment out Hydrad lines if needed. @@ -242,11 +271,13 @@ ; Iterate on TT and r3 + atr_a0 = a_tr/a_0 + for i=1,100 do begin calc_c1, tt_old, nn, length, rad, r3 - tt_new = r2*(3.5*r3/(1. + r3)*length*length*q(0)/kappa_0)^(2./7.) + tt_new = r2*(3.5*(r3-l_tr/l_c)/(1. + r3*atr_ac)*l_c^2.*q(0)*atr_a0/kappa_0)^(2./7.) radloss, rad, tt_new, 1, rtv=rtv - nn = (q(0)/((1. + r3)*rad))^0.5 + nn = (l_star/l_c*q(0)/((1. + atr_ac*r3)*rad))^0.5 err = tt_new - tt_old err_n = nn_old - nn if abs(err) lt 1e3 then begin @@ -258,7 +289,7 @@ endfor tt = tt_old - nn = (q(0)/((1. + r3)*rad))^0.5 + nn = (l_star/l_c*q(0)/((1. + atr_ac*r3)*rad))^0.5 ; If want to fix out of eqm start, e.g. cooling flare. Section 4.2, Paper 3. ; tt=1.e7*r2 @@ -277,7 +308,7 @@ v(0) = 0. ta(0) = t(0)/r2 calc_lambda, t(0), sc - na(0) = n(0)*r2*exp(-2*length*(1.-sin(3.14159/5.))/3.14159/sc); + na(0) = n(0)*r2*exp(-2*l_c*(1.-sin(3.14159/5.))/3.14159/sc) pa(0) = 2*k_b*na(0)*ta(0) print, 'Initial static equilibrium' @@ -318,7 +349,7 @@ ; Thermal conduction flux at base - f_cl = c1*(t(i)/r2)^3.5/length + f_cl = c1*(t(i)/r2)^3.5/l_c if keyword_set(classical) then begin f = f_cl @@ -343,16 +374,15 @@ ; Equilibrium thermal conduction flux at base (-R_tr in ApJ paper) f_eq = -r3*n(i)*n(i)*rad*length - ; pv = 0.4*(f_eq - f) - pv = 0.4*(f_eq - f - flux_nt(i)) - ; dn = pv*0.5/(r12*k_b*t(i)*length)*dt - dn = (pv*0.5/(r12*k_b*t(i)*length) + j_nt(i)/length)*dt + pv = 0.4*(atr_ac*f_eq/r3/length*(l_c^2/l_star)*(r3-l_tr/l_c) - a0_ac*f - flux_nt(i)) + + dn = (pv*0.5/(r12*k_b*t(i)*l_c) + j_nt(i)/length)*dt n(i+1) = n(i) + dn - - ; dp = 2./3.*(q(i) + (1. + 1./r3)*f_eq/length)*dt - dp = 2./3.*(q(i) + (1. + 1./r3)*f_eq/length $ - - (1. - 1.5*k_b*t(i)/energy_nt)*flux_nt(i)/length)*dt + + dp = 2./3.*(q(i) + (l_c/l_star)*(atr_ac + 1./r3)*f_eq/length $ + - (1. - 1.5*k_b*t(i)/energy_nt)*flux_nt(i)/length)*dt + p(i+1) = p(i) + dp t(i+1) = p(i+1)/(n(i+1)*2.*k_b) @@ -367,7 +397,7 @@ ; calculate apex quantities ; ta(i+1) = t(i+1)/r2; - na(i+1) = n(i+1)*r2*exp(-2.*length*(1.-sin(3.14159/5.))/3.14159/sc); + na(i+1) = n(i+1)*r2*exp(-2.*l_c*(1.-sin(3.14159/5.))/3.14159/sc); pa(i+1) = 2*k_b*na(i+1)*ta(i+1) ; Differential emission measure @@ -480,7 +510,7 @@ endif - cond(i)=f + cond(i)=f*a0_ac rad_cor(i)=f_eq/r3 endfor diff --git a/tests/.gitignore b/tests/.gitignore new file mode 100644 index 0000000..5391d87 --- /dev/null +++ b/tests/.gitignore @@ -0,0 +1,138 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ \ No newline at end of file diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000..63c0bef --- /dev/null +++ b/tests/README.md @@ -0,0 +1,26 @@ +# Testing with Python + +This directory contains a number of tests that can be run with the +[`pytest`](https://docs.pytest.org/en/stable/) testing framework in Python. +To install the needed dependencies, + +```shell +$ pip install pytest numpy hissw matplotlib +``` + +You will also need to properly configure [`hissw`](https://wtbarnes.github.io/hissw/) for your particular IDL installation. Then, to run the tests, + +```shell +$ pytest --ebtel_idl_path=/path/to/EBTEL +``` + +where `/path/to/EBTEL` is the path to the root of the EBTEL repository on your machine. +Alternatively, you can also plot the results from the tests to visually confirm +that everything is working as expected, + +```shell +$ pytest --ebtel_idl_path=/path/to/EBTEL --show_plots +``` + +Note that these tests do not necessarily test the correctness of the code, but simply +test that the code runs without crashing for several different input configurations. diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..9ca2e6c --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,20 @@ +import pytest + + +def pytest_addoption(parser): + parser.addoption( + "--ebtel_idl_path", action="store", default=None, help="Path to EBTEL IDL code" + ) + parser.addoption( + "--show_plots", action="store_true", default=False, + ) + + +@pytest.fixture +def ebtel_idl_path(request): + return request.config.getoption("--ebtel_idl_path") + + +@pytest.fixture +def show_plots(request): + return request.config.getoption("--show_plots") diff --git a/tests/test_cases.py b/tests/test_cases.py new file mode 100644 index 0000000..5eab03c --- /dev/null +++ b/tests/test_cases.py @@ -0,0 +1,141 @@ +import copy +from collections import OrderedDict +from subprocess import run + +import pytest + +import numpy as np +import matplotlib.pyplot as plt +import hissw + + +def run_idl(ebtel_idl_path, config): + # Flags and keywords + flags = [] + if 'dem' not in config or not config['dem']['use_new_method']: + flags += ['dem_old'] + if not config['use_flux_limiting']: + flags += ['classical'] + keys = [ + 'flux_nt', + 'energy_nt', + 'a_tr', + 'a_c', + 'a_0', + 'l_fact', + ] + keywords = [f'{k}={config[k]}' for k in keys if k in config] + + # Heating + time = np.arange(0, config['total_time']-config['tau'], config['tau']) + heat = np.ones(time.shape) * config['heating']['background'] + for _e in config['heating']['events']: + e = _e['event'] + # Rise + i = np.where(np.logical_and(time >= e['rise_start'], time < e['rise_end'])) + heat[i] += e['magnitude'] * (time[i] - e['rise_start']) / (e['rise_end'] - e['rise_start']) + # Plateau + i = np.where(np.logical_and(time >= e['rise_end'], time < e['decay_start'])) + heat[i] += e['magnitude'] + # Decay + i = np.where(np.logical_and(time >= e['decay_start'], time <= e['decay_end'])) + heat[i] += e['magnitude'] * (e['decay_end'] - time[i])/(e['decay_end'] - e['decay_start']) + + args = { + 'time': time.tolist(), + 'loop_length': config['loop_length'], + 'heat': heat.tolist(), + 'flags': flags, + 'keywords': keywords, + 'return_vars': config['return_vars'], + } + idl = hissw.Environment(extra_paths=[ebtel_idl_path]) + script = """time = {{ time }} +heat = {{ heat }} +loop_length = {{ loop_length }} +ebtel2,time,heat,loop_length,{{ return_vars | join(',') }}{% if flags %}, /{{ flags | join(', /') }}{% endif %}{% if keywords %}, {{ keywords | join(',') }}{% endif %} + """ + + return idl.run(script, args=args, save_vars=config['return_vars']) + + +def plot_results(result): + fig = plt.figure(figsize=(10,10)) + for i,v in enumerate(['temperature', 'density', 'pressure']): + ax = fig.add_subplot(2,2,i+1) + ax.plot(result['time'], result[v]) + ax.set_xlim(result['time'][[0,-1]]) + ax = fig.add_subplot(224) + ax.plot(result['temperature'], result['density'],) + ax.set_yscale('log') + ax.set_xscale('log') + plt.show() + + +# Set up configuration dictionaries for all the cases we want to consider +# When adding a new test case, put it here and then add it to the list of +# test cases below +base_config = { + 'total_time': 5e3, + 'tau': 1.0, + 'loop_length': 4e9, + 'use_flux_limiting': False, + 'heating': OrderedDict({ + 'background': 3e-5, + 'events': [] + }), + 'return_vars': [ + 'temperature', + 'density', + 'pressure', + 'velocity', + ] +} +# Case 1: No events +no_events = copy.deepcopy(base_config) +# Case 2: 1 event +one_event = copy.deepcopy(base_config) +one_event['heating']['events'] = [ + {'event': { + 'rise_start':0, + 'rise_end':100, + 'decay_start':100, + 'decay_end':200, + 'magnitude':0.1, + }} +] +# Case 3: Non-uniform cross-sections +varying_xs = copy.deepcopy(one_event) +varying_xs['a_c'] = 3.0 +varying_xs['a_0'] = 2.0 +varying_xs['a_tr'] = 1.0 +varying_xs['l_fact'] = 0.85 +# Case 4: +mulitple_returns = copy.deepcopy(base_config) +mulitple_returns['return_vars'] = mulitple_returns['return_vars'] + [ + 'ta', + 'na', + 'pa', + 'c11', + 'dem_tr', + 'dem_cor', + 'logtdem', + 'f_ratio', + 'rad_ratio', + 'cond', + 'rad_cor', +] + + +@pytest.mark.parametrize('config', [ + no_events, + one_event, + varying_xs, + mulitple_returns, +]) +def test_all_cases(ebtel_idl_path, show_plots, config): + r = run_idl(ebtel_idl_path, config) + for k in config['return_vars']: + assert k in r + if show_plots: + plot_results(r)