fix rti bugs and add documentation for drag

Author: Eb

star/defaults/controls.defaults

@@ -2983,7 +2983,7 @@
     conv_premix_dump_snapshots = .false.
 
 
-   ! Rayleigh Taylor Instability
+   ! Rayleigh-Taylor Instability
    ! ____________________________
 
       ! derived from Paul Duffell's code RT1D.
@@ -8196,8 +8196,11 @@
 
     use_dPrad_dm_form_of_T_gradient_eqn = .false.
     use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.    
-    
-    
+  
+  
+      ! Hydrodynamic drag
+      ! =================
+      
       ! drag_coefficient
       ! ~~~~~~~~~~~~~~~~
       ! min_q_for_drag
@@ -8219,12 +8222,32 @@
     drag_coefficient = 0d0
     min_q_for_drag = 0d0
 
+ 
+      ! v_drag_factor
+      ! ~~~~~~~~~~~~~
+      ! v_drag
+      ! ~~~~~~
+      ! q_for_v_drag_full_off
+      ! ~~~~~~~~~~~~~~~~~~~~~
+      ! q_for_v_drag_full_on
+      ! ~~~~~~~~~~~~~~~~~~~~
 
-      ! for hydro comparison tests (e.g., Sedov)
+      ! Implemented only for u_flag right now. Adds a pseudo drag term of the form
+      ! -v_drag_factor*(v-v_drag)^2/r, can be used damp velocities in outer layers
+      ! of a star. Effect is full on for q>q_for_v_drag_full_on and full off for
+      ! q < q_for_v_drag_full_off.
+
+      ! ::
+
+    v_drag_factor = 0d0
+    v_drag = 0d0
+    q_for_v_drag_full_off = 0.95d0
+    q_for_v_drag_full_on = 0.96d0
 
 
       ! Rayleigh-Taylor Instability
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! for hydro comparison tests (e.g., Sedov)
       ! RTI_A
       ! ~~~~~
       ! RTI_B

star/private/mix_info.f90

@@ -2016,7 +2016,7 @@ subroutine set_RTI_mixing_info(s, ierr)
                r = s% r(k)
                s% eta_RTI(k) = C*alpha_face*cs*r
 
-               if (is_bad(s% eta_RTI(k))) then
+               if (is_bad(s% eta_RTI(k)) .and. s% q(k) <= s% alpha_RTI_src_max_q) then
                   ierr = -1
                   return
                   if (s% stop_for_bad_nums) then

star/private/remove_shells.f90

@@ -1138,7 +1138,7 @@ subroutine do_relax_to_star_cut( &
          use interp_1d_def, only: pm_work_size
          use interp_1d_lib, only: interp_pm, interp_values, interp_value
          use adjust_xyz, only: change_net
-         use set_flags, only: set_v_flag, set_u_flag, set_rotation_flag
+         use set_flags, only: set_v_flag, set_u_flag, set_RTI_flag, set_rotation_flag
          use rotation_mix_info, only: set_rotation_mixing_info
          use hydro_rotation, only: set_i_rot, set_rotation_info
          use relax, only: do_relax_composition, do_relax_angular_momentum, do_relax_entropy
@@ -1150,7 +1150,7 @@ subroutine do_relax_to_star_cut( &
             ! determines if we turn off non_nuc_neu and eps_nuc for entropy relax
          integer, intent(out) :: ierr
 
-         logical :: v_flag, u_flag, rotation_flag
+         logical :: v_flag, u_flag, RTI_flag, rotation_flag
          type (star_info), pointer :: s
          character (len=net_name_len) :: net_name
          integer :: model_number, num_trace_history_values, photo_interval
@@ -1240,13 +1240,21 @@ subroutine do_relax_to_star_cut( &
             if (dbg) write(*,*) "set_v_flag ierr", ierr
             v_flag = .true.
          end if
+
          u_flag = .false.
          if (s% u_flag) then
             call set_u_flag(id, .false., ierr)
             if (dbg) write(*,*) "set_u_flag ierr", ierr
             u_flag = .true.
          end if
 
+         RTI_flag = .false.
+         if (s% RTI_flag) then
+            call set_RTI_flag(id, .false., ierr)
+            if (dbg) write(*,*) "set_RTI_flag ierr", ierr
+            RTI_flag = .true.
+         end if
+
          if (s% rotation_flag) then
             call set_rotation_flag(id, .false., ierr)
             if (dbg) write(*,*) "set_rotation_flag ierr", ierr