Other pressure implicit (#551)

 Changed other_pressure hook to be implicit.

---------

Co-authored-by: Adam Jermyn <[email protected]>

Author: rjfarmer

docs/source/changelog.rst

@@ -10,6 +10,19 @@ Changes in main
 
 .. _Backwards-incompatible changes main:
 
+Backwards-incompatible changes
+------------------------------
+
+Hooks
+-----
+
+The ``other_pressure`` hook has been converted to use ``auto_diff``
+thus the variable ``s% extra_pressure`` is now an ``auto_diff``
+and allows for the setting of the partial derivatives the
+pressure with respect to other variables.
+
+
+
 Changes in r23.05.1
 ===================
 

star/other/other_pressure.f90

@@ -38,14 +38,18 @@ module other_pressure
 
       subroutine default_other_pressure(id, ierr)
          use star_def
+         use auto_diff
          integer, intent(in) :: id
          integer, intent(out) :: ierr
          type (star_info), pointer :: s
          integer :: k
          ierr = 0
          call star_ptr(id, s, ierr)
          if (ierr /= 0) return
-         s% extra_pressure(:) = 0d0
+         do k=1,s%nz
+            s% extra_pressure(k) = 0d0
+         end do
+         ! note that extra_pressure is type(auto_diff_real_star_order1) so includes partials.
          return
       end subroutine default_other_pressure
 

star/private/alloc.f90

@@ -934,7 +934,7 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             call do1(s% extra_opacity_factor, c% extra_opacity_factor)
             if (failed('extra_opacity_factor')) exit
 
-            call do1(s% extra_pressure, c% extra_pressure)
+            call do1_ad(s% extra_pressure, c% extra_pressure)
             if (failed('extra_pressure')) exit
             call do1(s% eps_heat, c% eps_heat)
             if (failed('eps_heat')) exit

star/private/hydro_energy.f90

@@ -573,12 +573,12 @@ subroutine eval1_work(s, k, skip_Peos, &
          real(dp), dimension(s% species), intent(out) :: &
             d_work_dxa00, d_work_dxam1
          integer, intent(out) :: ierr
-         real(dp) :: alfa, beta, P_theta, extra_P, Peos_face, Av_face
+         real(dp) :: alfa, beta, P_theta, Peos_face, Av_face
          real(dp), dimension(s% species) :: d_Pface_dxa00, d_Pface_dxam1
          type(auto_diff_real_star_order1) :: &
             P_face_ad, A_times_v_face_ad, mlt_Pturb_ad, &
             PtrbR_ad, PtrbL_ad, PvscL_ad, PvscR_ad, Ptrb_div_etrb, PL_ad, PR_ad, &
-            Peos_ad, Ptrb_ad, Pvsc_ad, inv_R2
+            Peos_ad, Ptrb_ad, Pvsc_ad, inv_R2, extra_P
          logical :: test_partials
          integer :: j
          include 'formats'
@@ -694,7 +694,8 @@ subroutine eval1_work(s, k, skip_Peos, &
             if (.not. s% use_other_pressure) then
                extra_P = 0d0
             else if (k > 1) then 
-               extra_P = alfa*s% extra_pressure(k) + beta*s% extra_pressure(k-1) 
+               ! my_val_m1 = shift_m1(get_my_val(s,k-1)) for use in terms going into equation at k
+               extra_P = alfa*s% extra_pressure(k) + beta * shift_m1(s%extra_pressure(k-1))
             else
                extra_P = s% extra_pressure(k)
             end if

star/private/star_utils.f90

@@ -3284,7 +3284,7 @@ subroutine calc_Ptot_ad_tw( &
 
          Ptot_ad = Peos_ad + Pvsc_ad + Ptrb_ad + mlt_Pturb_ad
 
-         if (s% use_other_pressure) Ptot_ad%val = Ptot_ad%val + s% extra_pressure(k)
+         if (s% use_other_pressure) Ptot_ad = Ptot_ad + s% extra_pressure(k)
 
       end subroutine calc_Ptot_ad_tw
 

star/test_suite/other_physics_hooks/src/run_star_extras.f90

@@ -80,6 +80,9 @@ subroutine extras_controls(id, ierr)
          s% use_other_close_gaps = .true.
          s% other_close_gaps => my_close_gaps
 
+         s% use_other_pressure = .true.
+         s% other_pressure => my_other_pressure
+
          s% extras_startup => extras_startup
          s% extras_check_model => extras_check_model
          s% extras_finish_step => extras_finish_step
@@ -418,6 +421,24 @@ subroutine my_close_gaps(id, mix_type, min_gap, ierr)
       end subroutine my_close_gaps
 
 
+      subroutine my_other_pressure(id, ierr)
+         use star_def
+         use auto_diff
+         integer, intent(in) :: id
+         integer, intent(out) :: ierr
+         type (star_info), pointer :: s
+         integer :: k
+         ierr = 0
+         call star_ptr(id, s, ierr)
+         if (ierr /= 0) return
+         do k=1,s%nz
+            s% extra_pressure(k) = 0d0
+         end do
+         ! note that extra_pressure is type(auto_diff_real_star_order1) so includes partials.
+         return
+      end subroutine my_other_pressure
+
+
 
       end module run_star_extras
 

star_data/public/star_data_step_work.inc

@@ -276,7 +276,7 @@
       real(dp), pointer, dimension(:) :: extra_opacity_factor
 
       ! extra pressure profile (can be set by user)
-      real(dp), pointer :: extra_pressure(:) ! (ergs/gm/sec)
+      type(auto_diff_real_star_order1), pointer, dimension(:) :: extra_pressure(:) ! (ergs/gm/sec)
 
       ! extra heat profile (can be set by user)
       type(auto_diff_real_star_order1), pointer, dimension(:) :: extra_heat