ensure proper behaviour of forced thermal balance

src/dataIO/mod_input.f08

@@ -259,7 +259,8 @@ subroutine read_physicslist(unit, settings)
     if (incompressible) call settings%physics%set_incompressible()
     if (flow) call settings%physics%flow%enable()
     if (radiative_cooling) call settings%physics%enable_cooling(cooling_curve, ncool)
-    if (heating) call settings%physics%enable_heating(force_thermal_balance)
+    settings%physics%heating%force_thermal_balance = force_thermal_balance
+    if (heating) call settings%physics%enable_heating()
     if (external_gravity) call settings%physics%enable_gravity()
     if (parallel_conduction) then
       call settings%physics%enable_parallel_conduction(fixed_tc_para_value)

src/mod_version.f08

@@ -14,6 +14,6 @@ module mod_version
   implicit none
 
   !> legolas version number
-  character(len=10), parameter    :: LEGOLAS_VERSION = "2.0.3"
+  character(len=10), parameter    :: LEGOLAS_VERSION = "2.0.4"
 
 end module mod_version

src/physics/mod_heatloss.f08

@@ -17,6 +17,11 @@ module mod_heatloss
   type(conduction_t), pointer :: conduction => null()
   type(grid_t), pointer :: grid => null()
 
+  character(len=*), parameter :: msg_heatloss_balance = ( &
+    "unable to enforce thermal balance when heating is disabled. " // &
+    "Set 'force_thermal_balance = .false.' in the parfile or enable heating." &
+  )
+
   type, public :: heatloss_t
     type(cooling_t) :: cooling
     type(heating_t) :: heating
@@ -30,6 +35,7 @@ module mod_heatloss
   end type heatloss_t
 
   public :: new_heatloss
+  public :: msg_heatloss_balance
 
 contains
 
@@ -86,9 +92,12 @@ subroutine check_if_thermal_balance_needs_enforcing(this, conduction_tgt, grid_t
     type(conduction_t), intent(in) :: conduction_tgt
     type(grid_t), intent(in) :: grid_tgt
 
-    if (.not. settings%physics%heating%is_enabled()) return
+    if (is_adiabatic()) return
     if (.not. settings%physics%heating%force_thermal_balance) return
-
+    if (.not. settings%physics%heating%is_enabled()) then
+      call logger%error(msg_heatloss_balance)
+      return
+    end if
     if (.not. associated(this%heating%H, zero_func)) call log_usr_H_func_warning()
 
     call set_module_pointers(conduction_tgt, grid_tgt, this%cooling)
@@ -145,6 +154,15 @@ subroutine log_usr_H_func_warning() ! LCOV_EXCL_START
   end subroutine log_usr_H_func_warning ! LCOV_EXCL_STOP
 
 
+  logical function is_adiabatic()
+    logical :: no_cooling, no_heating, no_conduction
+    no_cooling = .not. settings%physics%cooling%is_enabled()
+    no_heating = .not. settings%physics%heating%is_enabled()
+    no_conduction = .not. settings%physics%conduction%is_enabled()
+    is_adiabatic = no_cooling .and. no_heating .and. no_conduction
+  end function is_adiabatic
+
+
   subroutine delete(this)
     class(heatloss_t), intent(inout) :: this
     nullify(settings)

src/settings/mod_physics_settings.f08

@@ -116,11 +116,10 @@ end subroutine enable_cooling
   pure subroutine enable_heating(this, force_thermal_balance)
     class(physics_settings_t), intent(inout) :: this
     logical, intent(in), optional :: force_thermal_balance
-    logical :: force_balance
 
-    force_balance = .true.
-    if (present(force_thermal_balance)) force_balance = force_thermal_balance
-    this%heating%force_thermal_balance = force_balance
+    if (present(force_thermal_balance)) then
+      this%heating%force_thermal_balance = force_thermal_balance
+    end if
     call this%heating%enable()
   end subroutine enable_heating
 

tests/unit_tests/mod_test_heatloss.pf

@@ -1,6 +1,7 @@
 module mod_test_heatloss
   use mod_suite_utils
   use mod_matrix_structure, only: matrix_t
+  use mod_heatloss, only: msg_heatloss_balance
   use funit
   implicit none
 
@@ -11,11 +12,6 @@ module mod_test_heatloss
   type(matrix_t) :: matrix_A, matrix_B
   complex(dp), allocatable :: eigenvals(:)
 
-  character(len=*), parameter :: msg_heatloss_balance = ( &
-    "forced thermal balance needs cooling + heating. " // &
-    "Set 'force_thermal_balance = .false.' in the parfile or enable heating/cooling" &
-  )
-
 contains
 
   @before
@@ -183,7 +179,13 @@ contains
     call set_name("heatloss: only heating, force thermal balance")
     call enable_heating(force_balance=.true.)
     call assemble_and_solve()
-    @assertExceptionRaised(msg_heatloss_balance)
+    @assertTrue(settings%physics%heating%is_enabled())
+    @assertFalse(settings%physics%cooling%is_enabled())
+    @assertEqual(0.0_dp, get_actual_lambdaT(), tolerance=TOL)
+    ! heating overridden through forced thermal balance, rho = 1 so H0 = lambdaT = 0
+    @assertEqual(0.0_dp, get_actual_H0(), tolerance=TOL)
+    @assertTrue(settings%physics%heating%force_thermal_balance)
+    @assertTrue(eigenvals_all_real())
   end subroutine test_H_only_force_balance