add tdc uq into HLLC

star/defaults/controls.defaults

@@ -2284,8 +2284,8 @@
       ! mlt_Pturb_factor
       ! ~~~~~~~~~~~~~~~~
 
-      ! include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3
-      ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 3d0
+      ! For MESA's momentume equation, include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3
+      ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 2d0
       ! corresponds to alpha_TDC_Ptdvdt = 1d0, where Pt = (2/3)*rho(k)*(mlt_vc(k)/sqrt(2/3))**2
 
       ! this replaces conv_dP_term_factor.   also see extra_pressure vector and other_pressure routine

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar

@@ -181,9 +181,9 @@ Profile_Panels1_title = ''
 Profile_Panels1_txt_scale = 0.7
 Profile_Panels1_num_panels = 4
 
-Profile_Panels1_yaxis_name(1) = 'Lc_div_L'
+Profile_Panels1_yaxis_name(1) = 'Eq'
 !Profile_Panels1_dymin(1) = 0.14
-Profile_Panels1_other_yaxis_name(1) = 'Frad_div_cUrad'!'lum_conv'
+Profile_Panels1_other_yaxis_name(1) = 'Uq'!'lum_conv'
 !Profile_Panels1_other_dymin(1) = 0.14
 
 Profile_Panels1_yaxis_name(2) = 'logT'!'logdq'!'radius'!'entropy'

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90

@@ -490,6 +490,7 @@ integer function extras_start_step(id)
          ! if v= 0, turn on v so we can kick
          if (.not. s%v_flag .and. .not. s%u_flag) then
             call star_set_v_flag(id, .true., ierr)
+            write (*,*) 'turning on v_flag hydro for kick'
          end if
 
          call gyre_in_mesa_extras_set_velocities(s, .false., ierr)

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc

@@ -158,8 +158,8 @@
                v_surf = s% v(1)
                v_surf_start = s% v_start(1)
             else if (s% u_flag) then
-               v_surf = s% u_face_val(1)
-               v_surf_start = s% u_face_start(1)
+               v_surf = s% u(1) ! s% u_face_val(1)
+               v_surf_start = s% u_start(1) ! s% u_face_start(1)
             else !
                v_surf = 0d0
                v_surf_start = 0d0

star/private/auto_diff_support.f90

@@ -1235,6 +1235,39 @@ function wrap_u_p1(s, k) result(v_p1)
          end if
       end function wrap_u_p1
 
+      function wrap_opt_time_center_u_m1(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = 0d0
+         if (k == 1) return
+         v_tc = wrap_u_m1(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k-1))
+      end function wrap_opt_time_center_u_m1
+
+      function wrap_opt_time_center_u_00(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = 0d0
+         v_tc = wrap_u_00(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k))
+      end function wrap_opt_time_center_u_00
+
+      function wrap_opt_time_center_u_p1(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = 0d0
+         if (k == s%nz) return
+         v_tc = wrap_u_p1(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k+1))
+      end function wrap_opt_time_center_u_p1
+
+
       function wrap_Hp_m1(s, k) result(Hp_m1)
          type (star_info), pointer :: s
          type(auto_diff_real_star_order1) :: Hp_m1

star/private/hydro_energy.f90

@@ -216,7 +216,6 @@ end subroutine setup_dL_dm
 
          subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
             use hydro_rsp2, only: compute_Eq_cell
-            use star_utils, only: get_face_weights
             use tdc_hydro, only: compute_tdc_Eq_div_w_face ! compute_Eq_cell
             real(dp) :: alfa, beta
             integer, intent(out) :: ierr
@@ -267,14 +266,13 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
                Eq_ad = s% Eq_ad(k)  ! compute_Eq_cell(s, k, ierr)
                if (ierr /= 0) return
             else if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC' .and. &
-               s% TDC_include_eturb_in_energy_equation) then ! not checking for v or u flag.
+               s% TDC_include_eturb_in_energy_equation .and. (s% v_flag .or. s% u_flag)) then ! not checking for v or u flag.
                 !Eq_ad = compute_tdc_Eq_cell(s, k, ierr) ! safe to just recompute
-                if (k == 1) then
-                   Eq_ad = compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3)
-                else
-                   call get_face_weights(s, k, alfa, beta)
-                   Eq_ad = alfa*compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
-                      beta*shift_m1(compute_tdc_Eq_div_w_face(s, k-1, ierr))*(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3)
+                if (k < s% nz) then
+                  Eq_ad = 0.5d0*(compute_tdc_Eq_div_w_face(s, k, ierr)*s% mlt_vc_ad(k) + &
+                     shift_p1(compute_tdc_Eq_div_w_face(s, k+1, ierr))*shift_p1(s% mlt_vc_ad(k+1)))/sqrt_2_div_3
+                else ! center cell is 0 at inner face
+                     Eq_ad = 0.5d0*compute_tdc_Eq_div_w_face(s, k, ierr)*s% mlt_vc_ad(k)/sqrt_2_div_3
                 end if
                 if (ierr /= 0) return
             end if
@@ -355,27 +353,25 @@ subroutine setup_RTI_diffusion(diffusion_eps_ad)
          end subroutine setup_RTI_diffusion
 
          subroutine setup_d_turbulent_energy_dt(ierr)
-            use star_utils, only: get_face_weights
             use const_def, only: sqrt_2_div_3
             integer, intent(out) :: ierr
             type(auto_diff_real_star_order1) :: TDC_eturb_cell
-            real (dp) :: TDC_eturb_cell_start, alfa, beta
+            real (dp) :: TDC_eturb_cell_start
             include 'formats'
             ierr = 0
             if (s% RSP2_flag) then
                d_turbulent_energy_dt_ad = (wrap_etrb_00(s,k) - get_etrb_start(s,k))/dt
             else if (s% mlt_vc_old(k) > 0d0 .and. s% MLT_option == 'TDC' .and. &
                s% TDC_include_eturb_in_energy_equation) then
                ! write a wrapper for this.
-               if (k == 1) then
-                  TDC_eturb_cell_start = pow2(s% mlt_vc_old(k)/sqrt_2_div_3)
-                  TDC_eturb_cell = pow2(s% mlt_vc(k)/sqrt_2_div_3)
-               else
-                  call get_face_weights(s, k, alfa, beta)
-                  TDC_eturb_cell_start = alfa*pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + &
-                     beta*pow2(s% mlt_vc_old(k-1)/sqrt_2_div_3)
-                  TDC_eturb_cell = alfa*pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
-                     beta*pow2(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3)
+               if (k < s% nz) then
+                  TDC_eturb_cell_start = 0.5d0*(pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + &
+                     pow2(s% mlt_vc_old(k+1)/sqrt_2_div_3))
+                  TDC_eturb_cell = 0.5d0*(pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
+                     pow2(shift_p1(s% mlt_vc_ad(k+1))/sqrt_2_div_3))
+               else ! center cell averaged with 0 for inner face
+                  TDC_eturb_cell_start = 0.5d0*pow2(s% mlt_vc_old(k)/sqrt_2_div_3)
+                  TDC_eturb_cell = 0.5d0*pow2(s% mlt_vc(k)/sqrt_2_div_3)
                end if
                d_turbulent_energy_dt_ad = (TDC_eturb_cell - TDC_eturb_cell_start) / dt
             else

star/private/hydro_riemann.f90

@@ -73,18 +73,19 @@ subroutine do1_dudt_eqn( &
             s, k, P_surf_ad, nvar, ierr)
          use accurate_sum_auto_diff_star_order1
          use star_utils, only: get_area_info_opt_time_center, save_eqn_residual_info
+         use tdc_hydro, only: compute_tdc_Uq_dm_cell
          type (star_info), pointer :: s
          integer, intent(in) :: k
          type(auto_diff_real_star_order1), intent(in) :: P_surf_ad  ! only for k=1
          integer, intent(in) :: nvar
          integer, intent(out) :: ierr
-
          integer :: nz, i_du_dt
          type(auto_diff_real_star_order1) :: &
             flux_in_ad, flux_out_ad, diffusion_source_ad, &
             geometry_source_ad, gravity_source_ad, &
             area_00, area_p1, inv_R2_00, inv_R2_p1, &
-            dudt_expected_ad, dudt_actual_ad, resid_ad
+            dudt_expected_ad, dudt_actual_ad, resid_ad, &
+            Uq_cell
          type(accurate_auto_diff_real_star_order1) :: sum_ad
          real(dp) :: dt, dm, ie_plus_ke, scal, residual
          logical :: dbg, do_diffusion, test_partials
@@ -123,8 +124,15 @@ subroutine do1_dudt_eqn( &
          call setup_gravity_source
          call setup_diffusion_source
 
+         ! Add turbulent eddy viscous acceleration Uq for TDC as source
+         Uq_cell = 0d0
+         if (s% MLT_option == 'TDC' .and. s%alpha_TDC_DampM > 0d0) then
+            Uq_cell = compute_tdc_Uq_dm_cell(s, k, ierr) ! Uq * dm
+            if (ierr /= 0) return
+         end if
+
          sum_ad = flux_in_ad - flux_out_ad + &
-            geometry_source_ad + gravity_source_ad + diffusion_source_ad
+            geometry_source_ad + gravity_source_ad + diffusion_source_ad + Uq_cell
          dudt_expected_ad = sum_ad
          dudt_expected_ad = dudt_expected_ad/dm
 
@@ -310,7 +318,6 @@ subroutine do1_uface_and_Pface(s, k, ierr)
          use eos_def, only: i_gamma1, i_lnfree_e, i_lnPgas
          use star_utils, only: calc_Ptot_ad_tw, get_face_weights
          use hydro_rsp2, only: compute_Uq_face
-         use tdc_hydro, only: compute_tdc_Uq_face
          type (star_info), pointer :: s
          integer, intent(in) :: k
          integer, intent(out) :: ierr
@@ -427,16 +434,10 @@ subroutine do1_uface_and_Pface(s, k, ierr)
          end if
 
 
-         if (s% RSP2_flag) then  ! include Uq in u_face
+         if (s% RSP2_flag) then  ! include Uq in u_face, To do: implement in sources instead ~ EbF
             Uq_ad = compute_Uq_face(s, k, ierr)
             if (ierr /= 0) return
             s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
-         else if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC') then ! include Uq in u_face
-            Uq_ad = compute_tdc_Uq_face(s, k, ierr)
-            if (ierr /= 0) return
-            s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
-         else
-            Uq_ad = 0d0
          end if
 
          s% u_face_val(k) = s% u_face_ad(k)%val

star/private/hydro_vars.f90

@@ -599,7 +599,7 @@ subroutine set_hydro_vars( &
                s% gradr_factor(nzlo:nzhi) = 1d0
             end if
 
-            if (s% alpha_TDC_DampM > 0) then
+            if (s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') 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'

star/private/read_model.f90

@@ -161,7 +161,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 .and. s% mlt_option == 'TDC') call set_viscosity_vars_TDC(s,ierr)
+         if (ierr == 0 .and. s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') 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'