fix: Mark nodes with too large scale radii as unsolvable

During differential evolution is the dark matter profile scale radius temporarily becomes larger than the virial radius (which is unphysical) mark the node as unsolvable to avoid attempts to evaluate the properties of the node in this unphysical regime.
galacticusorg · Aug 10, 2021 · 4dfca95 · 4dfca95
1 parent 4e849b1
commit 4dfca95
Showing 1 changed file with 23 additions and 6 deletions.
diff --git a/source/objects.nodes.components.dark_matter_profile.scale.F90 b/source/objects.nodes.components.dark_matter_profile.scale.F90
@@ -25,6 +25,7 @@ module Node_Component_Dark_Matter_Profile_Scale
   !!{
   Implements a dark matter profile method that provides a scale radius.
   !!}
+  use :: Dark_Matter_Halo_Scales, only : darkMatterHaloScaleClass
   implicit none
   private
   public :: Node_Component_Dark_Matter_Profile_Scale_Scale_Set        , Node_Component_Dark_Matter_Profile_Scale_Plausibility       , &
@@ -54,6 +55,10 @@ module Node_Component_Dark_Matter_Profile_Scale
   </component>
   !!]
 
+  ! Objects used by this component.
+  class(darkMatterHaloScaleClass), pointer :: darkMatterHaloScale_
+  !$omp threadprivate(darkMatterHaloScale_)
+
 contains
 
   !![
@@ -72,9 +77,13 @@ subroutine Node_Component_Dark_Matter_Profile_Scale_Thread_Initialize(parameters
     type(inputParameters), intent(inout) :: parameters_
     !$GLC attributes unused :: parameters_
 
-    if (defaultDarkMatterProfileComponent%scaleIsActive()) &
-         & call nodePromotionEvent%attach(defaultDarkMatterProfileComponent,nodePromotion,openMPThreadBindingAtLevel,label="nodeComponentDarkMatterProfileScale")
-    return
+    if (defaultDarkMatterProfileComponent%scaleIsActive()) then
+       !![
+       <objectBuilder class="darkMatterHaloScale" name="darkMatterHaloScale_" source="parameters_"/>
+       !!]
+       call nodePromotionEvent%attach(defaultDarkMatterProfileComponent,nodePromotion,openMPThreadBindingAtLevel,label="nodeComponentDarkMatterProfileScale")
+     end if
+     return
   end subroutine Node_Component_Dark_Matter_Profile_Scale_Thread_Initialize
 
   !![
@@ -90,8 +99,12 @@ subroutine Node_Component_Dark_Matter_Profile_Scale_Thread_Uninitialize()
     use :: Galacticus_Nodes, only : defaultDarkMatterProfileComponent
     implicit none
 
-    if (defaultDarkMatterProfileComponent%scaleIsActive()) &
-         & call nodePromotionEvent%detach(defaultDarkMatterProfileComponent,nodePromotion)
+    if (defaultDarkMatterProfileComponent%scaleIsActive()) then
+       !![
+       <objectDestructor name="darkMatterHaloScale_"/>
+       !!]
+       call nodePromotionEvent%detach(defaultDarkMatterProfileComponent,nodePromotion)
+    end if
     return
   end subroutine Node_Component_Dark_Matter_Profile_Scale_Thread_Uninitialize
 
@@ -117,7 +130,11 @@ subroutine Node_Component_Dark_Matter_Profile_Scale_Plausibility(node)
     ! Ensure that it is of the scale class.
     select type (darkMatterProfile)
     class is (nodeComponentDarkMatterProfileScale)
-       if (darkMatterProfile%scale() <= 0.0d0) then
+       if     (                                                                      &
+            &   darkMatterProfile%scale() <= 0.0d0                                   &
+            &  .or.                                                                  &
+            &   darkMatterProfile%scale() >  darkMatterHaloScale_%virialRadius(node) &
+            & ) then
           node%isPhysicallyPlausible=.false.
           node%isSolvable           =.false.
        end if