[Draft]: reviving adding eddy viscous damping to TDC

star/private/hydro_energy.f90

@@ -221,7 +221,7 @@ subroutine setup_dL_dm(ierr)
          end subroutine setup_dL_dm
 
          subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
-            !use hydro_rsp2, only: compute_Eq_cell
+            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, &
@@ -264,12 +264,9 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
                others_ad%val = others_ad%val + s% eps_pre_mix(k)
             if (s% do_phase_separation .and. s% do_phase_separation_heating) &
                others_ad%val = others_ad%val + s% eps_phase_separation(k)
-
-            Eq_ad = 0d0
-            if (s% RSP2_flag) then             
-               Eq_ad = s% Eq_ad(k) ! compute_Eq_cell(s, k, ierr)
-               if (ierr /= 0) return
-            end if   
+
+            Eq_ad = compute_Eq_cell(s, k, ierr) ! s% Eq_ad(k) XXX 
+            if (ierr /= 0) return
 
             call setup_RTI_diffusion(RTI_diffusion_ad)
 

star/private/hydro_momentum.f90

@@ -434,11 +434,8 @@ subroutine expected_non_HSE_term( &
 
          end if ! v_flag
 
-         Uq_ad = 0d0
-         if (s% RSP2_flag) then ! Uq(k) is turbulent viscosity drag at face k
-            Uq_ad = compute_Uq_face(s, k, ierr)
-            if (ierr /= 0) return
-         end if
+         Uq_ad = compute_Uq_face(s, k, ierr)
+         if (ierr /= 0) return
 
          other_ad = extra_ad - accel_ad + drag + Uq_ad
          other = other_ad%val

star/private/hydro_riemann.f90

@@ -410,11 +410,10 @@ subroutine do1_uface_and_Pface(s, k, ierr)
              end if
          end if
 
-         if (s% RSP2_flag) then ! include Uq in u_face
-            Uq_ad = compute_Uq_face(s, k, ierr)
-            if (ierr /= 0) return
-            s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
-         end if
+         ! include Uq in u_face
+         Uq_ad = compute_Uq_face(s, k, ierr)
+         if (ierr /= 0) return
+         s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
 
          s% u_face_val(k) = s% u_face_ad(k)%val
 

star/private/hydro_rsp2.f90

@@ -664,7 +664,7 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell)
          real(dp) :: f, ALFAM_ALFA
          include 'formats'
          ierr = 0
-         ALFAM_ALFA = s% RSP2_alfam*s% mixing_length_alpha
+         ALFAM_ALFA = s% RSP2_alfam * s% mixing_length_alpha
          if (ALFAM_ALFA == 0d0 .or. &
                k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
                k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
@@ -735,10 +735,10 @@ function compute_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration
             Uq_face = 0d0
          else
             r_00 = wrap_opt_time_center_r_00(s,k)
-            Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
+            Chi_00 = compute_Chi_cell(s,k,ierr)  ! s% Chi_ad(k) XXX
             if (k > 1) then
-               !Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr))
-               Chi_m1 = shift_m1(s% Chi_ad(k-1))
+               Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr))
+               !Chi_m1 = shift_m1(s% Chi_ad(k-1)) XXX
                if (ierr /= 0) return
             else
                Chi_m1 = 0d0

star/private/turb_support.f90

@@ -170,6 +170,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
          alpha_semiconvection, thermohaline_coeff, &
          mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr)
       use star_utils
+      use hydro_rsp2, only: compute_Eq_cell
       type (star_info), pointer :: s
       integer, intent(in) :: k
       character (len=*), intent(in) :: MLT_option
@@ -189,13 +190,27 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
       ! these are used by use_superad_reduction
       real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit
       type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, &
-         diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled
+         diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq
 
       logical ::  test_partials, using_TDC
       logical, parameter :: report = .false.
       include 'formats'
 
+      ! check if this particular k can be done with TDC
+      using_TDC = .false.
+      if (s% MLT_option == 'TDC') using_TDC = .true.
+      if (.not. s% have_mlt_vc) using_TDC = .false.
+      if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false.
+      if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k)
+
       ! Pre-calculate some things. 
+      Eq_div_w = 0d0
+      if (using_TDC) then
+         if (s% mlt_vc_old(k) > 0d0) then
+            check_Eq = compute_Eq_cell(s, k, ierr)
+            Eq_div_w = check_Eq/(s% mlt_vc_old(k)/sqrt_2_div_3) ! maybe should be using conv_vel???
+         end if
+      end if
       Pr = crad*pow4(T)/3d0
       Pg = P - Pr
       beta = Pg / P
@@ -234,14 +249,6 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
             k, s% solver_iter, s% model_number, gradr%val, grada%val, scale_height%val
       end if
 
-
-      ! check if this particular k can be done with TDC
-      using_TDC = .false.
-      if (s% MLT_option == 'TDC') using_TDC = .true.
-      if (.not. s% have_mlt_vc) using_TDC = .false.
-      if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false.
-      if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k)
-
       if (k >= 1) then
          s% dvc_dt_TDC(k) = 0d0
       end if
@@ -264,7 +271,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
          call set_TDC(&
             conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
             mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr)
+            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), Eq_div_w, ierr)
          s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
 
             if (ierr /= 0) then
@@ -282,7 +289,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
                call set_TDC(&
                   conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
                   mixing_type, scale, chiT, chiRho, gradr_scaled, r, P, T, rho, dV, Cp, opacity, &
-                  scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr)
+                  scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), Eq_div_w, ierr)
                s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
                if (ierr /= 0) then
                   if (s% report_ierr) write(*,*) 'ierr from set_TDC when using superad_reduction'

turb/private/tdc_support.f90

@@ -65,7 +65,7 @@ module tdc_support
       logical :: report
       real(dp) :: mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt
       type(auto_diff_real_tdc) :: A0, c0, L, L0, gradL, grada
-      type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma
+      type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w
    end type tdc_info
 
 contains
@@ -522,7 +522,7 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2)
       real(dp), parameter :: x_GAMMAR = 2.d0*sqrt(3.d0)
 
       S0 = convert(x_ALFAS*info%mixing_length_alpha*info%Cp*info%T/info%Hp)*info%grada
-      S0 = S0*Y
+      S0 = S0*Y + convert(info%Eq_div_w)
       D0 = convert(info%alpha_TDC_DAMP*x_CEDE/(info%mixing_length_alpha*info%Hp))
       gammar_div_alfa = info%alpha_TDC_DAMPR*x_GAMMAR/(info%mixing_length_alpha*info%Hp)
       DR0 = convert(4d0*boltz_sigma*pow2(gammar_div_alfa)*pow3(info%T)/(pow2(info%rho)*info%Cp*info%kap))

turb/public/turb.f90

@@ -92,12 +92,12 @@ end subroutine set_thermohaline
    subroutine set_TDC( &
             conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
             mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr)
+            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr)
       use tdc
       use tdc_support
       real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale
       type(auto_diff_real_star_order1), intent(in) :: &
-         chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada
+         chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w
       logical, intent(in) :: report
       type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D
       integer, intent(out) :: tdc_num_iters, mixing_type, ierr
@@ -139,6 +139,7 @@ subroutine set_TDC( &
       info%kap = opacity
       info%Hp = scale_height
       info%Gamma = Gamma
+      info%Eq_div_w = Eq_div_w
 
       ! Get solution
       Zub = upper_bound_Z

turb/test/src/test_turb.f90

@@ -79,7 +79,7 @@ subroutine compare_TDC_and_Cox_MLT()
       real(dp) :: mixing_length_alpha, conv_vel_start, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale
       type(auto_diff_real_star_order1) :: &
          r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL, grav, Lambda
-      type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma
+      type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma, Eq_div_w
       real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param
       character(len=3) :: MLT_option
       integer :: mixing_type, ierr, tdc_num_iters
@@ -125,6 +125,7 @@ subroutine compare_TDC_and_Cox_MLT()
       scale = L%val*1d-3
       report = .false.
       dt = 1d40 ! Long time-step so we get into equilibrium
+      Eq_div_w = 0d0
 
       ! MLT
       MLT_option = 'Cox'
@@ -136,7 +137,7 @@ subroutine compare_TDC_and_Cox_MLT()
       call set_TDC(&
          conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
          mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-         scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr)
+         scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr)
 
       write(*,1) 'TDC: Y, conv_vel_start, conv_vel, dt   ', Y_face%val, conv_vel_start, conv_vel% val, dt
 
@@ -154,7 +155,7 @@ subroutine check_TDC()
       real(dp) :: mixing_length_alpha, conv_vel_start, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale
       type(auto_diff_real_star_order1) :: &
          r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL
-      type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D
+      type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w
       integer :: mixing_type, ierr, tdc_num_iters
       logical :: report
       integer :: j
@@ -187,6 +188,7 @@ subroutine check_TDC()
       chiT = 1d0
       chiRho = 1d0
       gradr = 3d0 * P * opacity * L / (64 * pi * boltz_sigma * pow4(T) * cgrav * m)
+      Eq_div_w = 0d0
 
       write(*,*) "####################################"
       write(*,*) "Running dt test"
@@ -196,7 +198,7 @@ subroutine check_TDC()
          call set_TDC(&
             conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
             mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr)
+            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr)
 
          write(*,1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel% val
          if (report) stop