adding eddy viscous damping to TDC

star/defaults/controls_dev.defaults

@@ -162,7 +162,7 @@
 
 
 ! split burn
-! ================
+! ==========
 
 
       ! op_split_burn_min_T_for_variable_T_solver
@@ -174,3 +174,24 @@
       ! ::
 
     op_split_burn_min_T_for_variable_T_solver = 1d99
+
+
+! TDC
+! ===
+
+
+      ! alpha_TDC_DampM
+      ! ~~~~~~~~~~~~~~~
+
+      ! If alpha_TDC_DampM >0, then include eddy viscous damping in TDC alpha_TDC_DampM
+      ! This control is analogous to `RSP_alfam`, where the default is `RSP_alfam` = 0.25d0. 
+
+      ! This control is only supported when v_flag = .true.
+      ! If hydrostatic (v_flag = .false. , v = 0 ) there are no velocity gradients,
+      ! and thus no shear to drive turbulence. Without shear, the eddy viscosity term becomes zero.
+
+      ! ::  
+
+    alpha_TDC_DampM = 0.25d0 
+
+

star/private/ctrls_io.f90

@@ -111,7 +111,7 @@ module ctrls_io
     mixing_D_limit_for_log, trace_mass_location, min_tau_for_max_abs_v_location, &
     min_q_for_inner_mach1_location, max_q_for_outer_mach1_location, &
     conv_core_gap_dq_limit, &
-    alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, &
+    alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, alpha_TDC_DAMPM, &
 
     ! burn zone eps definitions for use in logs and profiles
     burn_min1, burn_min2, &
@@ -2059,6 +2059,7 @@ subroutine store_controls(s, ierr)
  s% alpha_TDC_DAMP = alpha_TDC_DAMP
  s% alpha_TDC_DAMPR = alpha_TDC_DAMPR
  s% alpha_TDC_PtdVdt = alpha_TDC_PtdVdt
+ s% alpha_TDC_DAMPM = alpha_TDC_DAMPM
  s% compare_TDC_to_MLT = compare_TDC_to_MLT
 
  s% RSP2_alfap = RSP2_alfap
@@ -3738,6 +3739,7 @@ subroutine set_controls_for_writing(s, ierr)
  alpha_TDC_DAMP = s% alpha_TDC_DAMP
  alpha_TDC_DAMPR = s% alpha_TDC_DAMPR
  alpha_TDC_PtdVdt = s% alpha_TDC_PtdVdt
+ alpha_TDC_DAMPM = s% alpha_TDC_DAMPM
  compare_TDC_to_MLT = s% compare_TDC_to_MLT
 
  RSP2_alfap= s% RSP2_alfap

star/private/hydro_energy.f90

@@ -220,8 +220,10 @@ subroutine setup_dL_dm(ierr)
             dL_dm_ad = (L00_ad - Lp1_ad)/dm
          end subroutine setup_dL_dm
 
-         subroutine setup_sources_and_others(ierr)  ! sources_ad, others_ad
-            !use hydro_rsp2, only: compute_Eq_cell
+
+         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, &
@@ -271,6 +273,17 @@ subroutine setup_sources_and_others(ierr)  ! sources_ad, others_ad
                if (ierr /= 0) return
             end if
 
+! we do not need this term for tdc, because Eq is included
+! in the turbulent energy equation solved inside TDC. It's effects
+! are seen in hydro_momentum
+!            if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC') then
+!                Eq_ad = s% Eq_ad(k) !compute_Eq_cell(s, k, ierr)
+!!                if (k==91) then
+!!                write(*,*) 'test Eq, k', Eq_ad %val , k
+!!                end if
+!                if (ierr /= 0) return
+!            end if
+
             call setup_RTI_diffusion(RTI_diffusion_ad)
 
             drag_energy = 0d0

star/private/hydro_momentum.f90

@@ -400,6 +400,7 @@ subroutine expected_non_HSE_term( &
             extra_ad = s% extra_grav(k)
          end if
 
+         Uq_ad = 0d0
          accel_ad = 0d0
          drag = 0d0
          s% dvdt_drag(k) = 0d0
@@ -432,14 +433,20 @@ subroutine expected_non_HSE_term( &
                s% dvdt_drag(k) = drag%val
             end if
 
-         end if  ! v_flag
+            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 = 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
+            ! need to add check for v_flag.
+            if (s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') then ! Uq(k) is turbulent viscosity drag at face k
+               Uq_ad = compute_Uq_face(s, k, ierr)
+               !write (*,*) Uq_ad %val, 'Uq_ad %val'
+               if (ierr /= 0) return
+            end if
 
+         end if  ! v_flag
+
          other_ad = extra_ad - accel_ad + drag + Uq_ad
          other = other_ad%val
 

star/private/hydro_riemann.f90

@@ -433,12 +433,23 @@ 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
 
+        if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC') then ! include Uq in u_face
+            call mesa_error(__FILE__,__LINE__,'Uq is not currently supported by u_flag')
+!           Uq_ad = compute_Uq_face(s, k, ierr)
+!           if (k==63) then
+!               write(*,*) 'test Uq, k', Uq_ad %val, k
+!           end if
+!           if (ierr /= 0) return
+!           s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
+        end if
+
          s% u_face_val(k) = s% u_face_ad(k)%val
 
          if (s% P_face_start(k) < 0d0) then

star/private/hydro_rsp2.f90

@@ -40,7 +40,8 @@ module hydro_rsp2
          do1_rsp2_L_eqn, do1_turbulent_energy_eqn, do1_rsp2_Hp_eqn, &
          compute_Eq_cell, compute_Uq_face, set_RSP2_vars, &
          Hp_face_for_rsp2_val, Hp_face_for_rsp2_eqn, set_etrb_start_vars, &
-         RSP2_adjust_vars_before_call_solver, get_RSP2_alfa_beta_face_weights
+         RSP2_adjust_vars_before_call_solver, get_RSP2_alfa_beta_face_weights, &
+         set_viscosity_vars_TDC
 
       real(dp), parameter :: &
          x_ALFAP = 2.d0/3.d0, &  ! Ptrb
@@ -113,6 +114,53 @@ subroutine set_RSP2_vars(s,ierr)
       end subroutine set_RSP2_vars
 
 
+        ! This routine is called to initialize eq and uq for TDC.
+        subroutine set_viscosity_vars_TDC(s,ierr)
+           type (star_info), pointer :: s
+           integer, intent(out) :: ierr
+           type(auto_diff_real_star_order1) :: x
+           integer :: k, op_err
+           include 'formats'
+           ierr = 0
+           op_err = 0
+
+            !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
+            do k=1,s%nz
+               ! Hp_face(k) <= 0 means it needs to be set.  e.g., after read file
+               if (s% Hp_face(k) <= 0) then
+                   ! this scale height for face is already calculated in TDC
+                   s% Hp_face(k) = s% scale_height(k) !get_scale_height_face_val(s,k)
+               end if
+            end do
+            !$OMP END PARALLEL DO
+           if (ierr /= 0) then
+              if (s% report_ierr) write(*,2) 'failed in set_viscosity_vars_TDC loop 1', s% model_number
+              return
+           end if
+           !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
+           do k=1,s% nz
+              x = compute_Chi_cell(s, k, op_err)
+              if (op_err /= 0) ierr = op_err
+              x = compute_Eq_cell(s, k, op_err)
+              if (op_err /= 0) ierr = op_err
+              x = compute_Uq_face(s, k, op_err)
+              if (op_err /= 0) ierr = op_err
+           end do
+           !$OMP END PARALLEL DO
+           if (ierr /= 0) then
+              if (s% report_ierr) write(*,2) 'failed in set_viscosity_vars_TDC loop 2', s% model_number
+              return
+           end if
+           if (.not. s% v_flag) then ! set values 0 if not using v_flag.
+              do k = 1, s% nz
+                 s% Eq(k) = 0d0; s% Eq_ad(k) = 0d0
+                 s% Chi(k) = 0d0; s% Chi_ad(k) = 0d0
+                 s% Uq(k) = 0d0
+              end do
+           end if
+        end subroutine set_viscosity_vars_TDC
+
+
       subroutine do1_rsp2_L_eqn(s, k, nvar, ierr)
          use star_utils, only: save_eqn_residual_info
          type (star_info), pointer :: s
@@ -585,15 +633,25 @@ function compute_d_v_div_r(s, k, ierr) result(d_v_div_r)  ! s^-1
          integer, intent(in) :: k
          type(auto_diff_real_star_order1) :: d_v_div_r
          integer, intent(out) :: ierr
-         type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
+         type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1, term1, term2
+         logical :: dbg
          include 'formats'
          ierr = 0
+         dbg = .false.
          v_00 = wrap_v_00(s,k)
          v_p1 = wrap_v_p1(s,k)
          r_00 = wrap_r_00(s,k)
          r_p1 = wrap_r_p1(s,k)
          if (r_p1%val == 0d0) r_p1 = 1d0
-         d_v_div_r = v_00/r_00 - v_p1/r_p1  ! units s^-1
+         d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
+
+        ! Debugging output to trace values
+        if (dbg .and. k == -63) then
+            write(*,*) 'test d_v_div_r, k:', k
+            write(*,*) 'v_00:', v_00%val, 'v_p1:', v_p1%val
+            write(*,*) 'r_00:', r_00%val, 'r_p1:', r_p1%val
+            write(*,*) 'd_v_div_r:', d_v_div_r %val
+        end if
       end function compute_d_v_div_r
 
 
@@ -662,9 +720,21 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell)
          type(auto_diff_real_star_order1) :: &
             rho2, r6_cell, d_v_div_r, Hp_cell, w_00, d_00, r_00, r_p1
          real(dp) :: f, ALFAM_ALFA
+         logical :: dbg
          include 'formats'
          ierr = 0
-         ALFAM_ALFA = s% RSP2_alfam*s% mixing_length_alpha
+         dbg = .false.
+
+         ! check where we are getting alfam from.
+         if (s% MLT_option == 'TDC' .and. .not. s% RSP2_flag) then
+            ALFAM_ALFA = s% alpha_TDC_DampM * s% mixing_length_alpha
+         else if (s% RSP2_flag) then
+            ALFAM_ALFA = s% RSP2_alfam * s% mixing_length_alpha
+         else ! this is for safety, but probably is never called.
+            ALFAM_ALFA = 0d0
+         end if
+
+
          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
@@ -678,7 +748,17 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell)
             if (ierr /= 0) return
             d_v_div_r = compute_d_v_div_r(s, k, ierr)
             if (ierr /= 0) return
-            w_00 = wrap_w_00(s,k)
+
+            ! don't need to check if mlt_vc > 0 here.
+            if (s% MLT_option == 'TDC' .and. .not. s% RSP2_flag) then
+                if (s% have_mlt_vc .and. s% okay_to_set_mlt_vc) then
+                   w_00 = s% mlt_vc_old(k)/sqrt_2_div_3! same as info%A0 from TDC
+                else
+                   w_00 = s% mlt_vc(k)/sqrt_2_div_3! same as info%A0 from TDC
+                end if
+            else ! normal RSP2
+                w_00 = wrap_w_00(s,k)
+            end if
             d_00 = wrap_d_00(s,k)
             f = (16d0/3d0)*pi*ALFAM_ALFA/s% dm(k)
             rho2 = pow2(d_00)
@@ -693,6 +773,18 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell)
          s% Chi(k) = Chi_cell%val
          s% Chi_ad(k) = Chi_cell
 
+        if (dbg .and. k==-100) then
+                write(*,*) ' s% ALFAM_ALFA', ALFAM_ALFA
+                write(*,*) 'Hp_cell', Hp_cell %val
+                write(*,*) 'd_v_div_r', d_v_div_r %val
+                write(*,*) ' f',  f
+                write(*,*) 'w_00',w_00 %val
+                write(*,*) 'd_00 ', d_00 %val
+                write(*,*) 'rho2 ', rho2 %val
+                write(*,*) 'r_00',  r_00 %val
+                write(*,*) 'r_p1 ',  r_p1 %val
+                write(*,*) 'r6_cell',  r6_cell %val
+        end if
       end function compute_Chi_cell
 
 
@@ -735,10 +827,14 @@ 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)
+
+            ! which do we adopt?
+            Chi_00 = compute_Chi_cell(s,k,ierr)  ! s% Chi_ad(k) XXX
+            !Chi_00 = s% Chi_ad(k)  ! compute_Chi_cell(s,k,ierr)
+
             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/hydro_vars.f90

@@ -483,7 +483,7 @@ subroutine set_hydro_vars( &
          use hydro_rotation, only: set_rotation_info, compute_j_fluxes_and_extra_jdot
          use brunt, only: do_brunt_B, do_brunt_N2
          use mix_info, only: set_mixing_info
-         use hydro_rsp2, only: set_RSP2_vars
+         use hydro_rsp2, only: set_RSP2_vars, set_viscosity_vars_TDC
 
          type (star_info), pointer :: s
          integer, intent(in) :: nzlo, nzhi
@@ -596,6 +596,15 @@ subroutine set_hydro_vars( &
             else
                s% gradr_factor(nzlo:nzhi) = 1d0
             end if
+
+            if (s% alpha_TDC_DampM > 0) then
+               call set_viscosity_vars_TDC(s,ierr)
+               if (ierr /= 0) then
+                  if (len_trim(s% retry_message) == 0) s% retry_message = 'set_viscosity_vars_TDC failed'
+                  if (s% report_ierr) write(*,*) 'ierr from set_viscosity_vars_TDC'
+                  return
+               end if
+            end if
 
             call set_mlt_vars(s, nzlo, nzhi, ierr)
             if (failed('set_mlt_vars')) return
@@ -634,6 +643,7 @@ subroutine set_hydro_vars( &
             end if
          end if
 
+
          if (s% doing_timing) &
             call update_time(s, time0, total, s% time_set_hydro_vars)
 

star/private/profile_getval.f90

@@ -1886,9 +1886,9 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val)
             case(p_Chi)
                if (rsp_or_w) val = s% Chi(k)
             case(p_Eq)
-               if (rsp_or_w) val = s% Eq(k)
+               val = s% Eq(k)
             case(p_Uq)
-               if (rsp_or_w) val = s% Uq(k)
+               val = s% Uq(k)
             case(p_Lr)
                val = get_Lrad(s,k)
             case(p_Lr_div_L)

star/private/read_model.f90

@@ -80,7 +80,7 @@ subroutine finish_load_model(s, restart, ierr)
             total_angular_momentum, reset_epsnuc_vectors, set_qs
          use hydro_rotation, only: use_xh_to_update_i_rot_and_j_rot, &
             set_i_rot_from_omega_and_j_rot, use_xh_to_update_i_rot, set_rotation_info
-         use hydro_RSP2, only: set_RSP2_vars
+         use hydro_RSP2, only: set_RSP2_vars, set_viscosity_vars_TDC
          use RSP, only: RSP_setup_part1, RSP_setup_part2
          use report, only: do_report
          use alloc, only: fill_ad_with_zeros
@@ -165,6 +165,7 @@ subroutine finish_load_model(s, restart, ierr)
          s% doing_finish_load_model = .true.
          call set_vars(s, s% dt, ierr)
          if (ierr == 0 .and. s% RSP2_flag) call set_RSP2_vars(s,ierr)
+         if (ierr == 0 .and. s% alpha_TDC_DampM > 0) call set_viscosity_vars_TDC(s,ierr)
          s% doing_finish_load_model = .false.
          if (ierr /= 0) then
             write(*,*) 'finish_load_model: failed in set_vars'

star/private/set_flags.f90

@@ -192,6 +192,7 @@ subroutine set_u_flag(id, u_flag, ierr)
             call set_v_flag(id, .false., ierr)
          end if
 
+
          contains
 
          subroutine del(xs)