Minor changes and edits related to 2.0

post_processing/pylbo/visualisation/continua.py

@@ -2,17 +2,28 @@
 from pylbo.utilities.logger import pylboLogger
 from pylbo.visualisation.legend_handler import LegendHandler
 
+SLOW_MIN = "slow-"
+SLOW_PLUS = "slow+"
+ALFVEN_MIN = "alfven-"
+ALFVEN_PLUS = "alfven+"
+THERMAL = "thermal"
+DOPPLER = "doppler"
+
 CONTINUA_NAMES = {
-    "slow-": r"$\Omega_S^-",
-    "slow+": r"$\Omega_S^+",
-    "alfven-": r"$\Omega_A^-",
-    "alfven+": r"$\Omega_A^+",
-    "thermal": r"$\Omega_T",
-    "doppler": r"$\Omega_0",
+    SLOW_MIN: r"$\Omega_S^-",
+    SLOW_PLUS: r"$\Omega_S^+",
+    ALFVEN_MIN: r"$\Omega_A^-",
+    ALFVEN_PLUS: r"$\Omega_A^+",
+    THERMAL: r"$\Omega_T",
+    DOPPLER: r"$\Omega_0",
 }
 CONTINUA_COLORS = ["red", "red", "cyan", "cyan", "green", "grey"]
 
 
+def is_zero(values):
+    return np.all(np.isclose(values, 0, atol=1e-10))
+
+
 def calculate_continua(ds):
     """
     Calculates the different continua for a given dataset.
@@ -47,7 +58,6 @@ def calculate_continua(ds):
     doppler = k2 * v02 / ds.scale_factor + k3 * v03
     slow_min = -np.sqrt(slow2)
     slow_plus = np.sqrt(slow2)
-    slow_is_zero = np.all(np.isclose(slow2, 0))
 
     # Thermal continuum calculation
     # wave vector parallel to magnetic field, uses vector projection and scale factor
@@ -65,17 +75,17 @@ def calculate_continua(ds):
     kappa_perp = ds.equilibria["kappa_perp"]
     # if there is no conduction/cooling, there is no thermal continuum.
     if (
-        all(dLdT == 0)
-        and all(dLdrho == 0)
-        and all(kappa_para == 0)
-        and all(kappa_perp == 0)
+        is_zero(dLdT)
+        and is_zero(dLdrho)
+        and is_zero(kappa_para)
+        and is_zero(kappa_perp)
     ):
         thermal = np.zeros_like(ds.grid_gauss)
     # if temperature is zero (no pressure), set to zero and return
-    elif (T == 0).all():
+    elif is_zero(T):
         thermal = np.zeros_like(ds.grid_gauss)
     # if slow continuum vanishes, thermal continuum is analytical
-    elif slow_is_zero:
+    elif is_zero(slow2):
         thermal = 1j * (gamma - 1) * (rho * dLdrho - dLdT * (ci2 + vA2)) / (cs2 + vA2)
     else:
         sigma_A2 = kpara**2 * vA2  # Alfvén frequency
@@ -139,14 +149,16 @@ def calculate_continua(ds):
             slow_plus[idx] = s_pos
 
     # get doppler-shifted continua and return
-    continua = {
-        list(CONTINUA_NAMES.keys())[0]: doppler + slow_min,  # minus already in slow_min
-        list(CONTINUA_NAMES.keys())[1]: doppler + slow_plus,
-        list(CONTINUA_NAMES.keys())[2]: doppler - np.sqrt(alfven2),
-        list(CONTINUA_NAMES.keys())[3]: doppler + np.sqrt(alfven2),
-        list(CONTINUA_NAMES.keys())[4]: thermal,
-        list(CONTINUA_NAMES.keys())[5]: doppler,
-    }
+    continua = {THERMAL: thermal, DOPPLER: doppler}
+    # minus already in slow_min
+    continua[SLOW_MIN] = doppler + slow_min if not is_zero(slow_min) else slow_min
+    continua[SLOW_PLUS] = doppler + slow_plus if not is_zero(slow_plus) else slow_plus
+    continua[ALFVEN_MIN] = (
+        doppler - np.sqrt(alfven2) if not is_zero(alfven2) else -np.sqrt(alfven2)
+    )
+    continua[ALFVEN_PLUS] = (
+        doppler + np.sqrt(alfven2) if not is_zero(alfven2) else np.sqrt(alfven2)
+    )
     return continua
 
 

src/boundaries/smod_natural_bounds_conduction.f08

@@ -92,7 +92,7 @@ subroutine add_natural_conduction_terms_bfield( &
 
     Gop_min = k3 * B02 - k2 * B03 / eps
     Fop = k2 * B02 / eps + k3 * B03
-    Kp = physics%conduction%tcprefactor(x)
+    Kp = physics%conduction%get_tcprefactor(x)
     Kp_plus = Kp + dkappa_perp_dB2
     Kp_plusplus = dkappa_perp_dB2 - (B01**2 * Kp_plus / B0**2)
 

src/dataIO/mod_console.f08

@@ -81,7 +81,7 @@ subroutine print_whitespace(lines)
     integer, intent(in) :: lines
     integer :: i
 
-    if (logger%get_logging_level() >= 1) then
+    if (logger%get_logging_level() > 1) then
       do i = 1, lines
         write(*, *) ""
       end do
@@ -91,14 +91,15 @@ end subroutine print_whitespace
 
   subroutine log_grid_info(settings)
     type(settings_t), intent(in) :: settings
+    integer :: dims
 
     call logger%info("              << Grid settings >>")
-    call logger%info("geometry           : " // settings%grid%get_geometry())
-    call logger%info("grid start         : " // str(settings%grid%get_grid_start()))
-    call logger%info("grid end           : " // str(settings%grid%get_grid_end()))
-    call logger%info("gridpoints (base)  : " // str(settings%grid%get_gridpts()))
-    call logger%info("gridpoints (Gauss) : " // str(settings%grid%get_gauss_gridpts()))
-    call logger%info("gridpoints (matrix): " // str(settings%dims%get_dim_matrix()))
+    call logger%info("geometry         : " // settings%grid%get_geometry())
+    call logger%info("grid start       : " // str(settings%grid%get_grid_start()))
+    call logger%info("grid end         : " // str(settings%grid%get_grid_end()))
+    call logger%info("points base grid : " // str(settings%grid%get_gridpts()))
+    dims = settings%dims%get_dim_matrix()
+    call logger%info("matrix dimensions: " // str(dims) // " x " // str(dims))
   end subroutine log_grid_info
 
 
@@ -135,20 +136,22 @@ subroutine log_physics_info(settings)
     subroutine log_flow_info()
       logical :: flow
       flow = settings%physics%flow%is_enabled()
+      if (.not. flow) return
       call logger%info("flow                     : " // str(flow))
     end subroutine log_flow_info
 
     subroutine log_gravity_info()
       logical :: gravity
       gravity = settings%physics%gravity%is_enabled()
+      if (.not. gravity) return
       call logger%info("external gravity         : " // str(gravity))
     end subroutine log_gravity_info
 
     subroutine log_cooling_info()
       logical :: cooling
       cooling = settings%physics%cooling%is_enabled()
-      call logger%info("radiative cooling        : " // str(cooling))
       if (.not. cooling) return
+      call logger%info("radiative cooling        : " // str(cooling))
       call logger%info("  cooling curve          : " // &
         trim(adjustl(settings%physics%cooling%get_cooling_curve())) &
       )
@@ -157,8 +160,8 @@ end subroutine log_cooling_info
     subroutine log_heating_info()
       logical :: heating
       heating = settings%physics%heating%is_enabled()
-      call logger%info("heating                  : " // str(heating))
       if (.not. heating) return
+      call logger%info("heating                  : " // str(heating))
       call logger%info("  forcing thermal balance: " // &
         str(settings%physics%heating%force_thermal_balance) &
       )
@@ -167,8 +170,8 @@ end subroutine log_heating_info
     subroutine log_parallel_conduction_info()
       logical :: tc_para
       tc_para = settings%physics%conduction%has_parallel_conduction()
-      call logger%info("parallel conduction      : " // str(tc_para))
       if (.not. tc_para) return
+      call logger%info("parallel conduction      : " // str(tc_para))
       if (settings%physics%conduction%has_fixed_tc_para()) then
         call logger%info("  fixed value            : " // &
           str(settings%physics%conduction%get_fixed_tc_para(), fmt=exp_fmt) &
@@ -179,8 +182,8 @@ end subroutine log_parallel_conduction_info
     subroutine log_perpendicular_conduction_info()
       logical :: tc_perp
       tc_perp = settings%physics%conduction%has_perpendicular_conduction()
-      call logger%info("perpendicular conduction : " // str(tc_perp))
       if (.not. tc_perp) return
+      call logger%info("perpendicular conduction : " // str(tc_perp))
       if (settings%physics%conduction%has_fixed_tc_perp()) then
         call logger%info("  fixed value            : " // &
           str(settings%physics%conduction%get_fixed_tc_perp(), fmt=exp_fmt) &
@@ -191,8 +194,8 @@ end subroutine log_perpendicular_conduction_info
     subroutine log_resistivity_info()
       logical :: resistivity
       resistivity = settings%physics%resistivity%is_enabled()
-      call logger%info("resistivity              : " // str(resistivity))
       if (.not. resistivity) return
+      call logger%info("resistivity              : " // str(resistivity))
       if (settings%physics%resistivity%has_fixed_resistivity()) then
         call logger%info("  fixed eta value        : " // &
           str(settings%physics%resistivity%get_fixed_resistivity(), fmt=exp_fmt) &
@@ -203,8 +206,8 @@ end subroutine log_resistivity_info
     subroutine log_viscosity_info()
       logical :: viscosity
       viscosity = settings%physics%viscosity%is_enabled()
-      call logger%info("viscosity                : " // str(viscosity))
       if (.not. viscosity) return
+      call logger%info("viscosity                : " // str(viscosity))
       call logger%info( &
         "  viscosity value        : " // &
         str(settings%physics%viscosity%get_viscosity_value()) &
@@ -218,8 +221,8 @@ end subroutine log_viscosity_info
     subroutine log_hall_info()
       logical :: hall
       hall = settings%physics%hall%is_enabled()
-      call logger%info("hall mhd                 : " // str(hall))
       if (.not. hall) return
+      call logger%info("hall mhd                 : " // str(hall))
       call logger%info( &
         "  by substitution        : " // &
         str(settings%physics%hall%is_using_substitution()) &
@@ -276,13 +279,12 @@ subroutine log_io_info(settings)
       "write derived eigenfunctions : " &
       // str(settings%io%write_derived_eigenfunctions) &
     )
+    if (.not. settings%io%write_ef_subset) return
     call logger%info( &
       "write eigenfunction subset : " // str(settings%io%write_ef_subset) &
     )
-    if (settings%io%write_ef_subset) then
-      call logger%info("    subset center : " // str(settings%io%ef_subset_center))
-      call logger%info("    subset radius : " // str(settings%io%ef_subset_radius))
-    end if
+    call logger%info("    subset center : " // str(settings%io%ef_subset_center))
+    call logger%info("    subset radius : " // str(settings%io%ef_subset_radius))
   end subroutine log_io_info
 
 end module mod_console

src/equilibria/smod_user_defined.f08

@@ -60,10 +60,8 @@ real(dp) function B03(x)
     B03 = x
   end function B03
 
-  real(dp) function g0(x, settings, background)
+  real(dp) function g0(x)
     real(dp), intent(in) :: x
-    type(settings_t), intent(in) :: settings
-    type(background_t), intent(in) :: background
     g0 = 1.0_dp / x**2
   end function g0
 

src/main.f08

@@ -10,6 +10,7 @@ program legolas
   use mod_global_variables, only: dp, str_len, initialise_globals
   use mod_matrix_structure, only: matrix_t
   use mod_equilibrium, only: set_equilibrium
+  use mod_inspections, only: do_equilibrium_inspections
   use mod_matrix_manager, only: build_matrices
   use mod_solvers, only: solve_evp
   use mod_output, only: datfile_path, create_datfile
@@ -56,6 +57,7 @@ program legolas
   timer%init_time = timer%end_timer()
 
   call print_console_info(settings)
+  call do_equilibrium_inspections(settings, grid, background, physics)
 
   call timer%start_timer()
   call build_matrices(matrix_B, matrix_A, settings, grid, background, physics)

src/matrices/smod_conduction_matrix.f08

@@ -107,8 +107,8 @@ subroutine add_conduction_matrix_terms_bfield( &
     dkappa_perp_dT = physics%conduction%dtcperpdT(x_gauss)
     dkappa_perp_dB2 = physics%conduction%dtcperpdB2(x_gauss)
     ! prefactors
-    Kp = physics%conduction%tcprefactor(x_gauss)
-    diffKp = physics%conduction%dtcprefactordr(x_gauss)
+    Kp = physics%conduction%get_tcprefactor(x_gauss)
+    diffKp = physics%conduction%get_dtcprefactordr(x_gauss)
     Kp_plus = Kp + dkappa_perp_dB2
     Kp_plusplus = dkappa_perp_dB2 - (B01**2 * Kp_plus / B0**2)
     ! operators

src/mod_equilibrium.f08

@@ -222,8 +222,6 @@ end subroutine user_defined_eq
   !!          not balanced, contains NaN or if density/temperature contains
   !!          negative values.
   subroutine set_equilibrium(settings, grid, background, physics)
-    use mod_inspections, only: perform_NaN_and_negative_checks, perform_sanity_checks
-
     type(settings_t), intent(inout) :: settings
     type(grid_t), intent(inout) :: grid
     type(background_t), intent(inout) :: background
@@ -235,19 +233,13 @@ subroutine set_equilibrium(settings, grid, background, physics)
     call set_equilibrium_values(settings, grid, background, physics)
     call grid%initialise()
 
-    ! Do initial checks for NaN and negative density/temperature
-    call perform_NaN_and_negative_checks(settings, grid, background, physics)
-
     if (settings%physics%cooling%is_enabled()) then
       call physics%heatloss%cooling%initialise()
     end if
     call physics%heatloss%check_if_thermal_balance_needs_enforcing( &
       physics%conduction, grid &
     )
     call physics%hall%validate_scale_ratio(grid%gaussian_grid)
-
-    ! Do final sanity checks on values
-    call perform_sanity_checks(settings, grid, background, physics)
   end subroutine set_equilibrium