Merge pull request #137 from n-claes/feature/cooling_curves

Addition of new cooling curves
n-claes · Apr 27, 2023 · 9fd0bdc · 9fd0bdc
2 parents 6a2f6da + 77b5c08
commit 9fd0bdc
Show file tree

Hide file tree

Showing 8 changed files with 284 additions and 3 deletions.
diff --git a/docs/general/parameter_file.md b/docs/general/parameter_file.md
@@ -58,7 +58,7 @@ This namelist includes all physics-related variables.
 | flow  | logical | inclusion of background flow effects | `.false.` |
 | radiative_cooling | logical | whether to include optically thin radiative losses | `.false.` |
 | ncool | int | number of points used when interpolating cooling curves | 4000 |
-| cooling_curve | string | which cooling curve to use, can be `{"jc_corona", "dalgarno", "ml_solar", "spex", "spex_dalgarno", "rosner"}` | `"jc_corona"` |
+| cooling_curve | string | which cooling curve to use, can be `{"jc_corona", "dalgarno", "dalgarno2", "ml_solar", "spex", "spex_dalgarno", "rosner", "colgan", "colgan_dm"}` | `"jc_corona"` |
 | heating | logical | whether to include background heating | `.false.` |
 | force_thermal_balance | logical | whether to set the heating in such a way to enforce thermal equilibrium | `.true.` |
 | external_gravity | logical | whether to include external gravity | `.false.` |

diff --git a/docs/physics/equations.md b/docs/physics/equations.md
@@ -84,6 +84,26 @@ The function $\HH(\rho, T)$ specifies the heating function. If thermal balance i
 constant in time but possibly varying in space and as such that it balances out the cooling contribution to ensure thermal equilibrium. $\HH(\rho, T)$ can be user-specified, in which case its density and
 temperature derivatives should be provided as well.
 
+#### Cooling curves
+At the moment 8 tabulated cooling curves are implemented, along with 1 analytical cooling curve. The list below gives an overview of the cooling curves that are available, along with their references.
+
+| Name    		 		| Notes   | Approximate $\log_{10}(T)$ range | Reference     |
+| :---         		| :---    | :---: 		|   :---  		 |
+| `rosner`     		| Analytical cooling curve | $3.891 - 7.605$ | [Rosner et al. (1978)](https://ui.adsabs.harvard.edu/abs/1978ApJ...220..643R/abstract) |
+| `jc_corona` 	  | Cooling curve for coronal conditions | $4.000 - 7.954$ | [Colgan et al. (2008)](https://ui.adsabs.harvard.edu/abs/2008ApJ...689..585C/abstract) |
+| `dalgarno`   	  | Cooling curve for low temperatures | $2.000 - 9.000$ | [Dalgarno & McCray (1972)](https://ui.adsabs.harvard.edu/abs/1972ARA%26A..10..375D/abstract) |
+| `dalgarno2` 		| Cooling curve for very low temperatures | $1.000 - 4.000$ | [Dalgarno & McCray (1972)](https://ui.adsabs.harvard.edu/abs/1972ARA%26A..10..375D/abstract) |
+| `ml_solar`  		| Cooling curve for solar abundances | $2.000 - 9.000$ | [Mellema & Lundqvist (2002)](https://ui.adsabs.harvard.edu/abs/2002A%26A...394..901M/abstract) |
+| `spex`			 		| Cooling curve for solar abundances | $3.800 - 8.160$ | [Schure et al. (2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A...508..751S/abstract) |
+| `spex_dalgarno` | `spex` supplemented with `dalgarno` for lower temperatures | $2.000 - 8.160$ | - |
+| `colgan`     		| The original cooling curve | $4.065 - 9.065$ | [Colgan & Feldman (2008)](https://ui.adsabs.harvard.edu/abs/2008ApJ...689..585C/abstract) |
+| `colgan_dm`  		| `colgan` supplemented with `dalgarno2` for lower temperatures | $1.000 - 9.065$ | - |
+
+All cooling curves are interpolated at high resolution using a local second order polynomial and then sampled on the grid, with the exception of the analytical curves.
+Note that all tabulated cooling curves are extended to higher temperatures (i.e. outside of their tabulated range) by assuming pure Brehmsstrahlung ($\Lambda(T) \sim \sqrt{T}$). For temperatures
+below the tabulated range the cooling is set to zero.
+
+
 ### Thermal conduction
 The thermal conduction prescription relies on a tensor representation to model the anisotropy in MHD and is given by
 

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -45,11 +45,13 @@ set (sources
         dataIO/mod_input.f08
         physics/cooling_curves/mod_cooling_curve_names.f08
         physics/cooling_curves/mod_data_dalgarno.f08
+        physics/cooling_curves/mod_data_dalgarno2.f08
         physics/cooling_curves/mod_data_jccorona.f08
         physics/cooling_curves/mod_data_mlsolar.f08
         physics/cooling_curves/mod_data_rosner.f08
         physics/cooling_curves/mod_data_spex.f08
         physics/cooling_curves/mod_data_spex_enh.f08
+        physics/cooling_curves/mod_data_colgan.f08
         physics/cooling_curves/mod_cooling_curves.f08
         physics/mod_physics_utils.f08
         physics/mod_physics.f08

diff --git a/src/physics/cooling_curves/mod_cooling_curve_names.f08 b/src/physics/cooling_curves/mod_cooling_curve_names.f08
@@ -6,12 +6,16 @@ module mod_cooling_curve_names
 
   character(len=str_len), parameter :: JC_CORONA = "jc_corona"
   character(len=str_len), parameter :: DALGARNO = "dalgarno"
+  character(len=str_len), parameter :: DALGARNO2 = "dalgarno2"
   character(len=str_len), parameter :: ML_SOLAR = "ml_solar"
   character(len=str_len), parameter :: SPEX = "spex"
   character(len=str_len), parameter :: SPEX_DALGARNO = "spex_dalgarno"
   character(len=str_len), parameter :: ROSNER = "rosner"
+  character(len=str_len), parameter :: COLGAN = "colgan"
+  character(len=str_len), parameter :: COLGAN_DM = "colgan_dm"
 
-  character(len=str_len), parameter :: KNOWN_CURVES(6) = [ &
-    JC_CORONA, DALGARNO, ML_SOLAR, SPEX, SPEX_DALGARNO, ROSNER &
+  character(len=str_len), parameter :: KNOWN_CURVES(9) = [ &
+    JC_CORONA, DALGARNO, DALGARNO2, ML_SOLAR, SPEX, &
+    SPEX_DALGARNO, ROSNER, COLGAN, COLGAN_DM &
   ]
 end module mod_cooling_curve_names
diff --git a/src/physics/cooling_curves/mod_cooling_curves.f08 b/src/physics/cooling_curves/mod_cooling_curves.f08
@@ -215,11 +215,13 @@ end function is_valid_cooling_curve
   subroutine get_cooling_table(name, table_T, table_lambda)
     use mod_cooling_curve_names
     use mod_data_dalgarno
+    use mod_data_dalgarno2
     use mod_data_jccorona
     use mod_data_mlsolar
     use mod_data_rosner
     use mod_data_spex
     use mod_data_spex_enh
+    use mod_data_colgan
 
     character(len=*), intent(in) :: name
     real(dp), intent(out), allocatable :: table_T(:)
@@ -235,6 +237,10 @@ subroutine get_cooling_table(name, table_T, table_lambda)
       table_n = n_dalgarno
       table_T = logT_dalgarno
       table_lambda = logL_dalgarno
+    case(DALGARNO2)
+      table_n = n_dalgarno2
+      table_T = logT_dalgarno2
+      table_lambda = logL_dalgarno2
     case(ML_SOLAR)
       table_n = n_mlsolar
       table_T = logT_mlsolar
@@ -257,6 +263,18 @@ subroutine get_cooling_table(name, table_T, table_lambda)
       table_lambda(n_spex_enh_dalgarno:) = ( &
         logL_spex(6:n_spex) + log10(L_spex_enh(6:n_spex)) &
       )
+    case(COLGAN)
+      table_n = n_colgan
+      table_T = logT_colgan
+      table_lambda = logL_colgan
+    case(COLGAN_DM)
+      table_n = n_colgan + n_dalgarno2
+      allocate(table_T(table_n))
+      allocate(table_lambda(table_n))
+      table_T(1:n_dalgarno2) = logT_dalgarno2(1:n_dalgarno2)
+      table_T(n_dalgarno2+1:) = logT_colgan(1:n_colgan)
+      table_lambda(1:n_dalgarno2) = logL_dalgarno2(1:n_dalgarno2)
+      table_lambda(n_dalgarno2+1:) = logL_colgan(1:n_colgan)
     end select
   end subroutine get_cooling_table
 

diff --git a/src/physics/cooling_curves/mod_data_colgan.f08 b/src/physics/cooling_curves/mod_data_colgan.f08
@@ -0,0 +1,39 @@
+module mod_data_colgan
+    use mod_global_variables, only: dp
+    implicit none
+
+    public
+
+    integer, parameter :: n_colgan = 55
+    !> log10 temperature values from Colgan (2008)
+    real(dp), protected :: logT_colgan(n_colgan)
+    !> log10 luminosity values from Colgan (2008)
+    real(dp), protected :: logL_colgan(n_colgan)
+
+    data logT_colgan / &
+      4.06460772_dp, 4.14229559_dp, 4.21995109_dp, 4.29760733_dp, 4.37527944_dp, 4.45293587_dp, &
+      4.53060946_dp, 4.60826923_dp, 4.68592974_dp, 4.76359269_dp, 4.79704583_dp, 4.83049243_dp, &
+      4.86394114_dp, 4.89738514_dp, 4.93083701_dp, 4.96428321_dp, 4.99773141_dp, 5.03116600_dp, &
+      5.06460772_dp, 5.17574368_dp, 5.28683805_dp, 5.39795738_dp, 5.50906805_dp, 5.62017771_dp, &
+      5.73129054_dp, 5.84240328_dp, 5.95351325_dp, 6.06460772_dp, 6.17574368_dp, 6.28683805_dp, &
+      6.39795738_dp, 6.50906805_dp, 6.62017771_dp, 6.73129054_dp, 6.84240328_dp, 6.95351325_dp, &
+      7.06460772_dp, 7.17574368_dp, 7.28683805_dp, 7.39795738_dp, 7.50906805_dp, 7.62017771_dp, &
+      7.73129054_dp, 7.84240328_dp, 7.95351325_dp, 8.06460772_dp, 8.17574368_dp, 8.28683805_dp, &
+      8.39795738_dp, 8.50906805_dp, 8.62017771_dp, 8.73129054_dp, 8.84240328_dp, 8.95351325_dp, & 
+      9.06460772_dp &
+    /
+    data logL_colgan / &
+      -22.18883401_dp, -21.78629635_dp, -21.60383554_dp, -21.68480662_dp, -21.76444630_dp, &
+      -21.67935529_dp, -21.54217864_dp, -21.37958284_dp, -21.25171892_dp, -21.17584161_dp, &
+      -21.15783402_dp, -21.14491111_dp, -21.13526945_dp, -21.12837453_dp, -21.12485189_dp, &
+      -21.12438898_dp, -21.12641785_dp, -21.12802448_dp, -21.12547760_dp, -21.08964778_dp, &
+      -21.08812360_dp, -21.19542445_dp, -21.34582346_dp, -21.34839251_dp, -21.31700703_dp, &
+      -21.29072156_dp, -21.28900309_dp, -21.34104468_dp, -21.43122351_dp, -21.62448270_dp, &
+      -21.86694036_dp, -22.02897478_dp, -22.08050874_dp, -22.06057061_dp, -22.01973295_dp, &
+      -22.00000434_dp, -22.05161149_dp, -22.22175466_dp, -22.41451671_dp, -22.52581288_dp, &
+      -22.56913516_dp, -22.57485721_dp, -22.56150512_dp, -22.53968863_dp, -22.51490350_dp, &
+      -22.48895932_dp, -22.46071057_dp, -22.42908363_dp, -22.39358639_dp, -22.35456791_dp, &
+      -22.31261375_dp, -22.26827428_dp, -22.22203698_dp, -22.17422996_dp, -22.12514145_dp &
+    /
+
+  end module mod_data_colgan
diff --git a/src/physics/cooling_curves/mod_data_dalgarno2.f08 b/src/physics/cooling_curves/mod_data_dalgarno2.f08
@@ -0,0 +1,48 @@
+module mod_data_dalgarno2
+  use mod_global_variables, only: dp
+  implicit none
+
+  public
+
+  integer, parameter :: n_dalgarno2 = 76
+  !> log10 temperature values from Dalgarno and McCray (1978) for low temperatures (DM2)
+  real(dp), protected :: logT_dalgarno2(n_dalgarno2)
+  !> log10 luminosity values from Dalgarno and McCray (1978) for low temperatures (DM2)
+  real(dp), protected :: logL_dalgarno2(n_dalgarno2)
+
+  data logT_dalgarno2 / &
+    1.00_dp, 1.04_dp, 1.08_dp, 1.12_dp, 1.16_dp, 1.20_dp, &
+    1.24_dp, 1.28_dp, 1.32_dp, 1.36_dp, 1.40_dp, &
+    1.44_dp, 1.48_dp, 1.52_dp, 1.56_dp, 1.60_dp, &
+    1.64_dp, 1.68_dp, 1.72_dp, 1.76_dp, 1.80_dp, &
+    1.84_dp, 1.88_dp, 1.92_dp, 1.96_dp, 2.00_dp, &
+    2.04_dp, 2.08_dp, 2.12_dp, 2.16_dp, 2.20_dp, &
+    2.24_dp, 2.28_dp, 2.32_dp, 2.36_dp, 2.40_dp, &
+    2.44_dp, 2.48_dp, 2.52_dp, 2.56_dp, 2.60_dp, & 
+    2.64_dp, 2.68_dp, 2.72_dp, 2.76_dp, 2.80_dp, &
+    2.84_dp, 2.88_dp, 2.92_dp, 2.96_dp, 3.00_dp, &
+    3.04_dp, 3.08_dp, 3.12_dp, 3.16_dp, 3.20_dp, &
+    3.24_dp, 3.28_dp, 3.32_dp, 3.36_dp, 3.40_dp, &
+    3.44_dp, 3.48_dp, 3.52_dp, 3.56_dp, 3.60_dp, & 
+    3.64_dp, 3.68_dp, 3.72_dp, 3.76_dp, 3.80_dp, &
+    3.84_dp, 3.88_dp, 3.92_dp, 3.96_dp, 4.00_dp &
+  /
+  data logL_dalgarno2 / &
+    -30.0377_dp, -29.7062_dp, -29.4055_dp, -29.1331_dp, -28.8864_dp, -28.6631_dp, &
+    -28.4614_dp, -28.2791_dp, -28.1146_dp, -27.9662_dp, -27.8330_dp, &
+    -27.7129_dp, -27.6052_dp, -27.5088_dp, -27.4225_dp, -27.3454_dp, &
+    -27.2767_dp, -27.2153_dp, -27.1605_dp, -27.1111_dp, -27.0664_dp, &
+    -27.0251_dp, -26.9863_dp, -26.9488_dp, -26.9119_dp, -26.8742_dp, &
+    -26.8353_dp, -26.7948_dp, -26.7523_dp, -26.7080_dp, -26.6619_dp, &
+    -26.6146_dp, -26.5666_dp, -26.5183_dp, -26.4702_dp, -26.4229_dp, &
+    -26.3765_dp, -26.3317_dp, -26.2886_dp, -26.2473_dp, -26.2078_dp, &
+    -26.1704_dp, -26.1348_dp, -26.1012_dp, -26.0692_dp, -26.0389_dp, &
+    -26.0101_dp, -25.9825_dp, -25.9566_dp, -25.9318_dp, -25.9083_dp, &
+    -25.8857_dp, -25.8645_dp, -25.8447_dp, -25.8259_dp, -25.8085_dp, &
+    -25.7926_dp, -25.7778_dp, -25.7642_dp, -25.7520_dp, -25.7409_dp, &
+    -25.7310_dp, -25.7222_dp, -25.7142_dp, -25.7071_dp, -25.7005_dp, &
+    -25.6942_dp, -25.6878_dp, -25.6811_dp, -25.6733_dp, -25.6641_dp, &
+    -25.6525_dp, -25.6325_dp, -25.6080_dp, -25.5367_dp, -25.4806_dp &
+  /
+
+end module mod_data_dalgarno2
diff --git a/tests/unit_tests/mod_test_cooling_tables.pf b/tests/unit_tests/mod_test_cooling_tables.pf
@@ -6,8 +6,10 @@ module mod_test_cooling_tables
 
   use mod_data_jccorona
   use mod_data_dalgarno
+  use mod_data_dalgarno2
   use mod_data_mlsolar
   use mod_data_spex
+  use mod_data_colgan
 
   implicit none
 
@@ -143,6 +145,22 @@ contains
   end function get_spex_dm_L_table
 
 
+  function get_colgan_dm_T_table()  result(Ttable)
+    real(dp), allocatable :: Ttable(:)
+    real(dp), allocatable :: dummy(:)
+    call get_cooling_table("colgan_dm", Ttable, dummy)
+    deallocate(dummy)
+  end function get_colgan_dm_T_table
+
+
+  function get_colgan_dm_L_table()  result(Ltable)
+    real(dp), allocatable :: Ltable(:)
+    real(dp), allocatable :: dummy(:)
+    call get_cooling_table("colgan_dm", dummy, Ltable)
+    deallocate(dummy)
+  end function get_colgan_dm_L_table
+
+
   @test
   subroutine test_cooling_curve_unknown()
     call set_name("cooling curve: unknown")
@@ -234,6 +252,47 @@ contains
   end subroutine test_cooling_curve_dalgarno
 
 
+  @test
+  subroutine test_cooling_curve_dalgarno2_below()
+    call set_name("cooling curve: dalgarno2 (below Tmin)")
+    call set_cooling_curve("dalgarno2")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(1.0d0, 9.5d0, 100) / unit_temperature
+    expected = 0.0_dp
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_dalgarno2_below
+
+
+  @test
+  subroutine test_cooling_curve_dalgarno2_above()
+    call set_name("cooling curve: dalgarno2 (above Tmax)")
+    call set_cooling_curve("dalgarno2")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(1.1d9, 1.0d10, 100) / unit_temperature
+    expected = get_brehmstrahlung_lambda(T0vals, logT_dalgarno2, logL_dalgarno2)
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+    expected = get_brehmstrahlung_dlambda(T0vals, logT_dalgarno2, logL_dalgarno2)
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_dalgarno2_above
+
+
+  @test
+  subroutine test_cooling_curve_dalgarno2()
+    call set_name("cooling curve: dalgarno2")
+    call set_cooling_curve("dalgarno2")
+    call allocate_arrays(pts=41)
+    T0vals = 10.0d0 ** logT_dalgarno2(20:60) / unit_temperature
+    expected = 10.0d0 ** logL_dalgarno2(20:60) / unit_lambdaT
+    actual = get_actual(physics%heatloss%cooling%lambdaT)
+    @assertEqual(expected, actual, tolerance=2.0d-1)
+  end subroutine test_cooling_curve_dalgarno2
+
+
   @test
   subroutine test_cooling_curve_ml_solar_below()
     call set_name("cooling curve: ml_solar (below Tmin)")
@@ -322,6 +381,7 @@ contains
     call set_cooling_curve("spex_dalgarno")
     call allocate_arrays(pts=100)
     T0vals = linspace(1.0d0, 9.5d0, 100) / unit_temperature
+    expected = 0.0_dp
     actual = get_dLdT()
     @assertEqual(expected, actual, tolerance=TOL)
     actual = get_dLdrho()
@@ -363,4 +423,94 @@ contains
     @assertEqual(expected, actual, tolerance=2.0d-1)
   end subroutine test_cooling_curve_spex_dalgarno
 
+
+  @test
+  subroutine test_cooling_curve_colgan_below()
+    call set_name("cooling curve: colgan (below Tmin)")
+    call set_cooling_curve("colgan")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(1.0d0, 95.0d0, 100) / unit_temperature
+    expected = 0.0_dp
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_colgan_below
+
+
+  @test
+  subroutine test_cooling_curve_colgan_above()
+    call set_name("cooling curve: colgan (above Tmax)")
+    call set_cooling_curve("colgan")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(2.0d9, 1.0d10, 100) / unit_temperature
+    expected = get_brehmstrahlung_lambda(T0vals, logT_colgan, logL_colgan)
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+    expected = get_brehmstrahlung_dlambda(T0vals, logT_colgan, logL_colgan)
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_colgan_above
+
+
+  @test
+  subroutine test_cooling_curve_colgan()
+    call set_name("cooling curve: colgan")
+    call set_cooling_curve("colgan")
+    call allocate_arrays(pts=41)
+    T0vals = 10.0d0 ** logT_colgan(10:50) / unit_temperature
+    expected = 10.0d0 ** logL_colgan(10:50) / unit_lambdaT
+    actual = get_actual(physics%heatloss%cooling%lambdaT)
+    @assertEqual(expected, actual, tolerance=2.0d-1)
+  end subroutine test_cooling_curve_colgan
+
+
+  @test
+  subroutine test_cooling_curve_colgan_dm_below()
+    call set_name("cooling curve: colgan_dm (below Tmin)")
+    call set_cooling_curve("colgan_dm")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(1.0d0, 9.5d0, 100) / unit_temperature
+    expected = 0.0_dp
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_colgan_dm_below
+
+
+  @test
+  subroutine test_cooling_curve_colgan_dm_above()
+    call set_name("cooling curve: colgan_dm (above Tmax)")
+    call set_cooling_curve("colgan_dm")
+    call allocate_arrays(pts=100)
+    T0vals = linspace(2.0d9, 1.0d10, 100) / unit_temperature
+    expected = get_brehmstrahlung_lambda( &
+      T0vals, get_colgan_dm_T_table(), get_colgan_dm_L_table() &
+    )
+    actual = get_dLdrho()
+    @assertEqual(expected, actual, tolerance=TOL)
+    expected = get_brehmstrahlung_dlambda( &
+      T0vals, get_colgan_dm_T_table(), get_colgan_dm_L_table() &
+    )
+    actual = get_dLdT()
+    @assertEqual(expected, actual, tolerance=TOL)
+  end subroutine test_cooling_curve_colgan_dm_above
+
+
+  @test
+  subroutine test_cooling_curve_colgan_dm()
+    real(dp), allocatable :: colgan_dm_T_table(:)
+    real(dp), allocatable :: colgan_dm_L_table(:)
+
+    call set_name("cooling curve: colgan_dm")
+    call set_cooling_curve("colgan_dm")
+    call allocate_arrays(pts=91)
+    call get_cooling_table("colgan_dm", colgan_dm_T_table, colgan_dm_L_table)
+    T0vals = 10.0d0 ** colgan_dm_T_table(30:120) / unit_temperature
+    expected = 10.0_dp ** colgan_dm_L_table(30:120) / unit_lambdaT
+    actual = get_actual(physics%heatloss%cooling%lambdaT)
+    @assertEqual(expected, actual, tolerance=3.5d-1)
+  end subroutine test_cooling_curve_colgan_dm
+
 end module mod_test_cooling_tables