Reintroduce velocity drag (#541)

Reincorporated velocity drag into energy and momentum equations, and added controls to toggle

Added a new test suite make_pre_ccsn_13BVN for testing velocity drag, for v_flag (#535)
---------

Co-authored-by: Rob Farmer <[email protected]>
Co-authored-by: EbF <[email protected]>
Co-authored-by: Evan Bauer <[email protected]>
Co-authored-by: Ebraheem Farag <[email protected]>

Author: 3 people

docs/source/test_suite/make_pre_ccsn_13bvn

@@ -0,0 +1 @@
+../../../star/test_suite/make_pre_ccsn_13bvn/README.rst

star/defaults/controls.defaults

@@ -8136,10 +8136,29 @@
       ! ::
 
     use_dPrad_dm_form_of_T_gradient_eqn = .false.
-    use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.
-    
+    use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.    
 
 
+      ! drag_coefficient
+      ! ~~~~~~~~~~~~~~~~
+      ! min_q_for_drag
+      ! ~~~~~~~~~~~~~~
+
+      ! only when v_flag.  adjusts both v and energy transfer from kinetic to thermal.
+      ! only for v(k) when q(k) > min_q_for_drag.
+      ! kill off fraction of v = drag_coefficient (i.e. set to 1 to keep v near 0)
+      ! useful for preventing the development radial pulsations during advanced
+      ! burning in massive stars and AGB stars.
+
+      ! Under certain circumstances we will not have drag in the surface k=1 zone
+      ! To force the drag term to be on in the outer zone you must enable one of the 
+      ! following surface boundary conditions:
+      ! ``use_momentum_outer_BC``, ``use_zero_Pgas_outer_BC``, or ``use_fixed_Psurf_outer_BC``
+
+      ! ::
+
+    drag_coefficient = 0d0
+    min_q_for_drag = 0d0
 
 
       ! for hydro comparison tests (e.g., Sedov)

star/private/alloc.f90

@@ -985,6 +985,12 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
                   action == do_copy_pointers_and_resize) &
                s% alpha_mlt(1:nz) = s% mixing_length_alpha
 
+            call do1(s% dvdt_drag, c% dvdt_drag)
+            if (failed('dvdt_drag')) exit
+
+            call do1(s% FdotV_drag_energy, c% FdotV_drag_energy)
+            if (failed('FdotV_drag_energy')) exit
+
             call do1(s% vc, c% vc)
             if (failed('vc')) exit
 

star/private/ctrls_io.f90

@@ -350,6 +350,7 @@ module ctrls_io
     P_theta_for_velocity_time_centering, L_theta_for_velocity_time_centering, &
     steps_before_use_TDC, use_P_d_1_div_rho_form_of_work_when_time_centering_velocity, compare_TDC_to_MLT, &
     velocity_logT_lower_bound, max_dt_yrs_for_velocity_logT_lower_bound, velocity_q_upper_bound, &
+    drag_coefficient, min_q_for_drag, &
     retry_for_v_above_clight, &
 
     ! hydro solver
@@ -1856,6 +1857,9 @@ subroutine store_controls(s, ierr)
  s% RTI_m_full_boost = RTI_m_full_boost
  s% RTI_m_no_boost = RTI_m_no_boost
 
+ s% drag_coefficient = drag_coefficient
+ s% min_q_for_drag = min_q_for_drag
+
  s% velocity_logT_lower_bound = velocity_logT_lower_bound
  s% max_dt_yrs_for_velocity_logT_lower_bound = max_dt_yrs_for_velocity_logT_lower_bound
  s% velocity_q_upper_bound = velocity_q_upper_bound
@@ -3520,6 +3524,10 @@ subroutine set_controls_for_writing(s, ierr)
  RTI_m_full_boost = s% RTI_m_full_boost
  RTI_m_no_boost = s% RTI_m_no_boost
 
+
+ drag_coefficient = s% drag_coefficient
+ min_q_for_drag = s% min_q_for_drag
+
  velocity_logT_lower_bound = s% velocity_logT_lower_bound
  max_dt_yrs_for_velocity_logT_lower_bound = s% max_dt_yrs_for_velocity_logT_lower_bound
  velocity_q_upper_bound = s% velocity_q_upper_bound

star/private/hydro_energy.f90

@@ -224,7 +224,8 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
             !use hydro_rsp2, only: compute_Eq_cell
             integer, intent(out) :: ierr
             type(auto_diff_real_star_order1) :: &
-               eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad
+               eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad, &
+               v_00, v_p1, drag_force, drag_energy
             include 'formats'
             ierr = 0
 
@@ -272,7 +273,26 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
 
             call setup_RTI_diffusion(RTI_diffusion_ad)
 
-            sources_ad = eps_nuc_ad - non_nuc_neu_ad + extra_heat_ad + Eq_ad + RTI_diffusion_ad
+            drag_energy = 0d0
+            s% FdotV_drag_energy(k) = 0
+            if (k .ne. s% nz) then
+               if (s% q(k) > s% min_q_for_drag .and. s% drag_coefficient > 0) then
+                  v_00 = wrap_v_00(s,k)
+                  drag_force = s% drag_coefficient*v_00/s% dt
+                  drag_energy = 0.5d0*v_00*drag_force
+                  s% FdotV_drag_energy(k) = drag_energy%val
+               ! drag energy for outer half-cell.   the 0.5d0 is for dm/2
+               end if
+               if (s% q(k+1) > s% min_q_for_drag .and. s% drag_coefficient > 0) then
+                  v_p1 = wrap_v_p1(s,k)
+                  drag_force = s% drag_coefficient*v_p1/s% dt
+                  drag_energy = drag_energy + 0.5d0*v_p1*drag_force
+                  s% FdotV_drag_energy(k) = drag_energy%val
+               ! drag energy for inner half-cell.   the 0.5d0 is for dm/2
+               end if
+            end if
+
+            sources_ad = eps_nuc_ad - non_nuc_neu_ad + extra_heat_ad + Eq_ad + RTI_diffusion_ad + drag_energy
 
             sources_ad%val = sources_ad%val + s% irradiation_heat(k)
 

star/private/hydro_eqns.f90

@@ -745,6 +745,14 @@ subroutine PT_eqns_surf(s, nvar, do_du_dt, do_dv_dt, do_equL, ierr)
             i_P_eqn = s% i_dv_dt
          end if
 
+
+         if(s% drag_coefficient > 0) then
+            ! We dont call expected_non_HSE_term with k==1 unless we call set_momentum_BC
+            ! so lets initilize this to zero, then if we dont call set_momentum_BC we have a 
+            ! sensible value here.
+            s% dvdt_drag(1) = 0 
+         end if
+         
          need_P_surf = .false.
          if (s% use_compression_outer_BC) then
             call set_compression_BC(ierr)

star/private/hydro_momentum.f90

@@ -387,7 +387,8 @@ subroutine expected_non_HSE_term( &
             other_ad, accel_ad,Uq_ad
          real(dp), intent(out) :: other
          integer, intent(out) :: ierr
-         type(auto_diff_real_star_order1) :: extra_ad, v_00
+         type(auto_diff_real_star_order1) :: extra_ad, v_00, &
+            drag
          real(dp) :: accel, d_accel_dv, fraction_on
          logical :: test_partials, local_v_flag
 
@@ -401,6 +402,8 @@ subroutine expected_non_HSE_term( &
          end if
 
          accel_ad = 0d0
+         drag = 0d0
+         s% dvdt_drag(k) = 0d0
          if (s% v_flag) then
 
             if (s% i_lnT == 0) then
@@ -423,6 +426,12 @@ subroutine expected_non_HSE_term( &
             end if
             accel_ad%val = accel
             accel_ad%d1Array(i_v_00) = d_accel_dv
+
+            if (s% q(k) > s% min_q_for_drag .and. s% drag_coefficient > 0) then
+               v_00 = wrap_v_00(s,k)
+               drag = -s% drag_coefficient*v_00/s% dt
+               s% dvdt_drag(k) = drag%val
+            end if
 
          end if ! v_flag
 
@@ -432,7 +441,7 @@ subroutine expected_non_HSE_term( &
             if (ierr /= 0) return
          end if
 
-         other_ad = extra_ad - accel_ad + Uq_ad
+         other_ad = extra_ad - accel_ad + drag + Uq_ad
          other = other_ad%val
 
          !test_partials = (k == s% solver_test_partials_k)

star/test_suite/do1_test_source

@@ -12,6 +12,7 @@ do_one c13_pocket "all values are within tolerance" "final.mod" auto
 do_one split_burn_big_net "stop because have dropped below central lower limit for si28" "final.mod" skip
 do_one 12M_pre_ms_to_core_collapse "all values are within tolerance" "final.mod" auto
 do_one 20M_pre_ms_to_core_collapse "all values are within tolerance" "final.mod" auto
+do_one make_pre_ccsn_13bvn "termination code: fe_core_infall_limit" "final.mod" auto
 do_one zams_to_cc_80 "all values are within tolerance" "final.mod" auto
 do_one pisn "termination code: Star is unbound" "final.mod" skip
 do_one 15M_dynamo "all values are within tolerances" "final.mod" auto

star/test_suite/each_test_run

@@ -92,6 +92,11 @@ function check_restart {
         fi
 
         # check that final model matches
+        if [[ ! -f ./ck ]]; then
+            failure_msg "$1 restart failed: No ck script"
+            return 1
+        fi
+
 
         if ! ./ck
         then

star/test_suite/make_pre_ccsn_13bvn/README.rst

@@ -0,0 +1,37 @@
+.. _make_pre_ccsn_13bvn:
+
+*******************
+make_pre_ccsn_13bvn
+*******************
+
+This test suite evolves a solar metalicity 12 |MSun| model from the pre-ms to helium depletion, strips the hydrogen envelope, and then evolves to core collapse.
+This test suite is made to be similar (almost identical) to the 13BVN model from MESA IV. (This is a sensitive test_suite)
+
+Physical checks
+===============
+
+Tests the functionality of velocity drag for v_flag for damping large surfaces velocities from He shell (N-alpha) flashes during advanced burning. 
+Tests energy conservation for regions where v_drag is applied.
+Tests the use of OP_split_burn as split burn <4d9 leads to large temperature swings and single zone burning in the core, and split burn >4.5d9 needs smaller timesteps. 
+This case illustrates strong coupling between burning and hydro: OP_split_burn lower and higher thresholds can only be exceeded with very tight timesteps.
+
+Inlists
+=======
+
+This test case has five parts.
+
+* Part 1 (``inlist_to_zams``) creates a 12 |Msun|, Z=2*10^-2 metallicity, pre-main sequence model and evolves it the zero age main sequence.
+
+* Part 2 (``inlist_before_remove``) takes the model to core helium depletion.
+
+* Part 3 (``inlist_remove``) removes the hydrogen envelope, leaving a stripped, evolved helium core.
+
+* Part 4 (``inlist_after_remove``) replaces the remaining hydrogen in the envelope with helium.
+
+* Part 5 (``inlist_to_pre_si_burn``) evolves the model until log_center_T = 9.45d0, just before Si-core burning.
+
+* Part 6 (``inlist_to_cc``) evolves the model until core collapse (300 km/s infall).
+
+
+Last-Updated: 24May2023 by EbF
+