diff --git a/star/Makefile b/star/Makefile
index c248e3d6b..23cde982d 100644
--- a/star/Makefile
+++ b/star/Makefile
@@ -202,6 +202,8 @@ SRCS += \
private/pgstar_profile_panels.f90 \
private/pgstar_grid.f90 \
private/pgstar_network.f90 \
+ private/pgstar_equation_residuals.f90 \
+ private/pgstar_kipp_residuals.f90 \
other/sample_pgstar_plot.f90 \
private/pgstar_full.f90
EXTERNAL_DEPENDS_ON := pgplot
@@ -237,11 +239,13 @@ MODULES += pgstar_abundance.mod \
pgstar_decorator.mod \
pgstar_dpg_dnu.mod \
pgstar_dynamo.mod \
+ pgstar_equation_residuals.mod \
pgstar_grid.mod \
pgstar_history_panels.mod \
pgstar_hist_track.mod \
pgstar_hr.mod \
pgstar_kipp.mod \
+ pgstar_kipp_residuals.mod \
pgstar_logg_logt.mod \
pgstar_logg_teff.mod \
pgstar_logl_r.mod \
diff --git a/star/defaults/pgstar.defaults b/star/defaults/pgstar.defaults
index 346009501..341c90d7b 100644
--- a/star/defaults/pgstar.defaults
+++ b/star/defaults/pgstar.defaults
@@ -7814,6 +7814,66 @@
Grid9_file_width = -1
Grid9_file_aspect_ratio = -1
+
+ ! Max_eq_resid
+ ! ------------
+
+ ! ::
+
+ Max_eq_resid_win_flag = .false.
+ Max_eq_resid_win_width = 7
+ Max_eq_resid_win_aspect_ratio = 5d-1
+
+ Max_eq_resid_xleft = 0.15
+ Max_eq_resid_xright = 0.85
+ Max_eq_resid_ybot = 0.15
+ Max_eq_resid_ytop = 0.85
+ Max_eq_resid_txt_scale = 1
+ Max_eq_resid_title = 'Max equation residuals'
+ Max_eq_resid_max_width = -1
+
+ ! **File output**
+
+ ! ::
+
+ Max_eq_resid_file_flag = .false.
+ Max_eq_resid_file_dir = 'png'
+ Max_eq_resid_file_prefix = 'max_eq_resid_'
+ Max_eq_resid_file_width = -1
+ Max_eq_resid_file_aspect_ratio = -1
+
+
+ ! Kipp_residuals
+ ! ------------
+
+ ! ::
+
+ Kipp_residuals_win_flag = .false.
+ Kipp_residuals_win_width = 9
+ Kipp_residuals_win_aspect_ratio = 6d-1
+
+ ! ::
+
+ Kipp_residuals_xleft = 0.15
+ Kipp_residuals_xright = 0.85
+ Kipp_residuals_ybot = 0.15
+ Kipp_residuals_ytop = 0.85
+ Kipp_residuals_txt_scale = 1
+ Kipp_residuals_title = 'Kippenhahn residuals'
+ Kipp_residuals_max_width = -1
+ Kipp_residuals_yaxis_name = 'mass'
+
+ ! **File output**
+
+ ! ::
+
+ Kipp_residuals_file_flag = .false.
+ Kipp_residuals_file_dir = 'png'
+ Kipp_residuals_file_prefix = 'kipp_residuals_'
+ Kipp_residuals_file_width = -1
+ Kipp_residuals_file_aspect_ratio = -1
+
+
! Annotation
! ==========
diff --git a/star/private/pgstar_ctrls_io.f90 b/star/private/pgstar_ctrls_io.f90
index dcf53870f..071f4561b 100644
--- a/star/private/pgstar_ctrls_io.f90
+++ b/star/private/pgstar_ctrls_io.f90
@@ -3023,6 +3023,38 @@ module pgstar_ctrls_io
Grid9_file_width, &
Grid9_file_aspect_ratio, &
+ Max_eq_resid_win_flag, &
+ Max_eq_resid_win_width, &
+ Max_eq_resid_win_aspect_ratio, &
+ Max_eq_resid_file_flag, &
+ Max_eq_resid_file_width, &
+ Max_eq_resid_file_aspect_ratio, &
+ Max_eq_resid_file_dir, &
+ Max_eq_resid_file_prefix, &
+ Max_eq_resid_xleft, &
+ Max_eq_resid_xright, &
+ Max_eq_resid_ybot, &
+ Max_eq_resid_ytop, &
+ Max_eq_resid_txt_scale, &
+ Max_eq_resid_title, &
+ Max_eq_resid_max_width, &
+
+ Kipp_residuals_win_flag, &
+ Kipp_residuals_win_width, &
+ Kipp_residuals_win_aspect_ratio, &
+ Kipp_residuals_file_flag, &
+ Kipp_residuals_file_width, &
+ Kipp_residuals_file_aspect_ratio, &
+ Kipp_residuals_file_dir, &
+ Kipp_residuals_file_prefix, &
+ Kipp_residuals_xleft, &
+ Kipp_residuals_xright, &
+ Kipp_residuals_ybot, &
+ Kipp_residuals_ytop, &
+ Kipp_residuals_txt_scale, &
+ Kipp_residuals_title, &
+ Kipp_residuals_max_width, &
+ Kipp_residuals_yaxis_name, &
annotation1_ci, &
annotation1_ch, &
@@ -6203,6 +6235,39 @@ subroutine store_pgstar_controls(s, ierr)
s% pg% Grid9_file_width = Grid9_file_width
s% pg% Grid9_file_aspect_ratio = Grid9_file_aspect_ratio
+ s% pg% Max_eq_resid_win_flag = Max_eq_resid_win_flag
+ s% pg% Max_eq_resid_win_width = Max_eq_resid_win_width
+ s% pg% Max_eq_resid_win_aspect_ratio = Max_eq_resid_win_aspect_ratio
+ s% pg% Max_eq_resid_file_flag = Max_eq_resid_file_flag
+ s% pg% Max_eq_resid_file_width = Max_eq_resid_file_width
+ s% pg% Max_eq_resid_file_aspect_ratio = Max_eq_resid_file_aspect_ratio
+ s% pg% Max_eq_resid_file_dir = Max_eq_resid_file_dir
+ s% pg% Max_eq_resid_file_prefix = Max_eq_resid_file_prefix
+ s% pg% Max_eq_resid_xleft = Max_eq_resid_xleft
+ s% pg% Max_eq_resid_xright = Max_eq_resid_xright
+ s% pg% Max_eq_resid_ybot = Max_eq_resid_ybot
+ s% pg% Max_eq_resid_ytop = Max_eq_resid_ytop
+ s% pg% Max_eq_resid_txt_scale = Max_eq_resid_txt_scale
+ s% pg% Max_eq_resid_title = Max_eq_resid_title
+ s% pg% Max_eq_resid_max_width = Max_eq_resid_max_width
+
+ s% pg% Kipp_residuals_win_flag = Kipp_residuals_win_flag
+ s% pg% Kipp_residuals_win_width = Kipp_residuals_win_width
+ s% pg% Kipp_residuals_win_aspect_ratio = Kipp_residuals_win_aspect_ratio
+ s% pg% Kipp_residuals_file_flag = Kipp_residuals_file_flag
+ s% pg% Kipp_residuals_file_width = Kipp_residuals_file_width
+ s% pg% Kipp_residuals_file_aspect_ratio = Kipp_residuals_file_aspect_ratio
+ s% pg% Kipp_residuals_file_dir = Kipp_residuals_file_dir
+ s% pg% Kipp_residuals_file_prefix = Kipp_residuals_file_prefix
+ s% pg% Kipp_residuals_xleft = Kipp_residuals_xleft
+ s% pg% Kipp_residuals_xright = Kipp_residuals_xright
+ s% pg% Kipp_residuals_ybot = Kipp_residuals_ybot
+ s% pg% Kipp_residuals_ytop = Kipp_residuals_ytop
+ s% pg% Kipp_residuals_txt_scale = Kipp_residuals_txt_scale
+ s% pg% Kipp_residuals_title = Kipp_residuals_title
+ s% pg% Kipp_residuals_max_width = Kipp_residuals_max_width
+ s% pg% Kipp_residuals_yaxis_name = Kipp_residuals_yaxis_name
+
s% pg% annotation1_ci = annotation1_ci
s% pg% annotation1_ch = annotation1_ch
diff --git a/star/private/pgstar_equation_residuals.f90 b/star/private/pgstar_equation_residuals.f90
new file mode 100644
index 000000000..11d440364
--- /dev/null
+++ b/star/private/pgstar_equation_residuals.f90
@@ -0,0 +1,209 @@
+! ***********************************************************************
+!
+! Copyright (C) 2026 The MESA Team
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU Lesser General Public License
+! as published by the Free Software Foundation,
+! either version 3 of the License, or (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+! See the GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with this program. If not, see .
+!
+! ***********************************************************************
+
+
+module pgstar_equation_residuals
+
+ use star_private_def
+ use const_def, only: dp
+ use pgstar_support
+ use star_pgstar
+
+ implicit none
+
+contains
+
+ subroutine Max_eq_resid_plot(id, device_id, ierr)
+
+ integer, intent(in) :: id
+ integer, intent(in) :: device_id
+ integer, intent(out) :: ierr
+
+ type (star_info), pointer :: s
+
+ ierr = 0
+ call get_star_ptr(id, s, ierr)
+ if (ierr /= 0) return
+
+ call pgslct(device_id)
+ call pgbbuf()
+ call pgeras()
+
+ call do_Max_eq_resid_plot(s, id, &
+ s% pg% Max_eq_resid_xleft, s% pg% Max_eq_resid_xright, &
+ s% pg% Max_eq_resid_ybot, s% pg% Max_eq_resid_ytop, .false., &
+ s% pg% Max_eq_resid_title, s% pg% Max_eq_resid_txt_scale, &
+ s% pg% Max_eq_resid_max_width, ierr)
+ if (ierr /= 0) return
+
+ call pgebuf()
+
+ end subroutine Max_eq_resid_plot
+
+
+
+ ! history of residuals for each structure equation
+ subroutine do_Max_eq_resid_plot(s, id, &
+ win_xleft, win_xright, win_ybot, win_ytop, subplot, title, txt_scale, &
+ max_width, ierr)
+
+ integer, intent(in) :: id, max_width
+ logical, intent(in) :: subplot
+ real, intent(in) :: &
+ win_xleft, win_xright, win_ybot, win_ytop, txt_scale
+ character (len=*), intent(in) :: title
+ integer, intent(out) :: ierr
+ type (star_info), pointer :: s
+
+ integer :: i_eq, initial_width
+ integer :: n
+ integer, save :: resid_hist_nvar
+ integer, save :: init_model
+ real :: xmin, xmax
+ real :: ymin, ymax
+
+ ! use save, we set these once, and they are reused in later calls of this function
+ real, allocatable, dimension(:), save :: xvec, yvec, resid_hist_model
+ real, allocatable, dimension(:, :), save :: resid_hist_vals
+ character(len=strlen), allocatable, dimension(:), save :: resid_equ_names
+
+ ierr = 0
+
+ call star_ptr(id, s, ierr)
+ if (ierr /= 0) return
+
+ if (.not. allocated(resid_hist_model)) then
+ init_model = s% model_number - 1
+ end if
+
+ if (max_width < 0) then
+ n = s% model_number - init_model
+ else
+ n = min(s% model_number - init_model, max_width)
+ end if
+
+ if (.not. allocated(resid_hist_model)) then
+ resid_hist_nvar = s% nvar_hydro
+ if (max_width < 0) then
+ initial_width = 100
+ else
+ initial_width = max_width
+ end if
+ allocate(resid_hist_model(initial_width))
+ allocate(resid_hist_vals(resid_hist_nvar, initial_width))
+ allocate(resid_equ_names(resid_hist_nvar))
+ allocate(xvec(initial_width), yvec(initial_width), stat=ierr)
+ if (ierr /= 0) then
+ write(*,*) 'allocate failed for PGSTAR'
+ return
+ end if
+
+ do i_eq = 1, resid_hist_nvar ! set names
+ resid_equ_names(i_eq) = trim(s% nameofequ(i_eq))
+ end do
+ else if (s% model_number - init_model >= size(resid_hist_model) .and. max_width < 0) then
+ call realloc_resid_hist(2*size(resid_hist_model)) ! grow the arrays
+ end if
+
+ ! roll if using a max width and you go over the limit
+ if (max_width > 1 .and. s% model_number - init_model > max_width) then
+ resid_hist_model(1:n-1) = resid_hist_model(2:n)
+ resid_hist_vals(:,1:n-1) = resid_hist_vals(:,2:n)
+ end if
+
+ resid_hist_model(n) = s% model_number
+
+ ! get the max resid values for this step
+ do i_eq = 1, resid_hist_nvar
+ resid_hist_vals(i_eq, n) = &
+ real(max(1d-40, maxval(abs(s% equ(i_eq, 1:s% nz)))))
+ end do
+
+ xvec(1:n) = real(resid_hist_model(1:n))
+
+ xmin = xvec(1)
+ xmax = xvec(n)
+
+ if (xmin >= xmax) xmax = xmin + 1.0
+
+ ymin = minval(resid_hist_vals(:, 1:n))
+ ymax = maxval(resid_hist_vals(:, 1:n))
+ ymin = max(ymin, 1e-20)
+ ymax = max(ymax, 10.0*ymin)
+
+ call pgsave
+ call pgsci(0) ! color index 0 = background
+ call pgsvp(win_xleft, win_xright, win_ybot, win_ytop)
+ call pgrect(0.0, 1.0, 0.0, 1.0) ! pseudo erase
+ call pgswin(xmin, xmax, log10(ymin), log10(ymax))
+ call pgsch(txt_scale)
+ call pgscf(1)
+ call pgsci(1)
+ call pgbox('BCNST',0.0,0,'BCNST',0.0,0)
+ call pgmtxt('B',2.5,0.5,0.5,'Model Number')
+ call pgmtxt('L',3.0,0.5,0.5,'log10(max |residual|)')
+ call pgmtxt('T',1.0,0.5,0.5, title)
+
+ ! draw lines
+ do i_eq = 1, resid_hist_nvar
+ yvec(1:n) = log10(resid_hist_vals(i_eq, 1:n))
+ call pgsci(mod(i_eq-1,13) + 2)
+ call pgline(n, xvec(1:n), yvec(1:n))
+ end do
+
+ ! legend
+ call pgsch(0.75 * txt_scale)
+
+ do i_eq = 1, resid_hist_nvar
+ call pgsci(mod(i_eq-1,13)+2)
+ call pgptxt( &
+ xmax + 0.03, &
+ log10(ymax) - (0.03 + 0.06*real(i_eq-1))*(log10(ymax)-log10(ymin)), &
+ 0.0, 0.0, &
+ trim(resid_equ_names(i_eq)))
+ end do
+
+ call pgsci(1)
+ call pgunsa
+
+ contains
+
+ subroutine realloc_resid_hist(new_size)
+ integer, intent(in) :: new_size
+ real, allocatable :: tmp_model(:)
+ real, allocatable :: tmp_vals(:,:)
+
+ call move_alloc(resid_hist_model, tmp_model)
+ call move_alloc(resid_hist_vals, tmp_vals)
+
+ allocate(resid_hist_model(new_size))
+ allocate(resid_hist_vals(resid_hist_nvar, new_size))
+
+ deallocate(xvec, yvec)
+ allocate(xvec(new_size), yvec(new_size))
+
+ resid_hist_model(1:n) = tmp_model(1:n)
+ resid_hist_vals(:,1:n) = tmp_vals(:,1:n)
+
+ deallocate(tmp_model, tmp_vals)
+ end subroutine
+
+ end subroutine do_Max_eq_resid_plot
+
+end module pgstar_equation_residuals
diff --git a/star/private/pgstar_full.f90 b/star/private/pgstar_full.f90
index 1ae79de2a..d65d2d1bb 100644
--- a/star/private/pgstar_full.f90
+++ b/star/private/pgstar_full.f90
@@ -161,6 +161,10 @@ subroutine set_win_file_data(s, ierr)
Color_Magnitude1_plot, Color_Magnitude2_plot, Color_Magnitude3_plot, &
Color_Magnitude4_plot, Color_Magnitude5_plot, Color_Magnitude6_plot, &
Color_Magnitude7_plot, Color_Magnitude8_plot, Color_Magnitude9_plot
+ use pgstar_equation_residuals, only: &
+ Max_eq_resid_plot
+ use pgstar_kipp_residuals, only: &
+ Kipp_residuals_plot
type (star_info), pointer :: s
integer, intent(out) :: ierr
@@ -1261,6 +1265,34 @@ subroutine set_win_file_data(s, ierr)
p% file_width = s% pg% Grid9_file_width
p% file_aspect_ratio = s% pg% Grid9_file_aspect_ratio
+ p => s% pg% pgstar_win_file_ptr(i_max_eq_resid)
+ p% plot => Max_eq_resid_plot
+ p% id = i_max_eq_resid
+ p% name = 'Max_eq_resid'
+ p% win_flag = s% pg% max_eq_resid_win_flag
+ p% win_width = s% pg% max_eq_resid_win_width
+ p% win_aspect_ratio = s% pg% max_eq_resid_win_aspect_ratio
+ p% file_flag = s% pg% max_eq_resid_file_flag
+ p% file_dir = s% pg% max_eq_resid_file_dir
+ p% file_prefix = s% pg% max_eq_resid_file_prefix
+ p% file_interval = s% pg% max_eq_resid_file_interval
+ p% file_width = s% pg% max_eq_resid_file_width
+ p% file_aspect_ratio = s% pg% max_eq_resid_file_aspect_ratio
+
+ p => s% pg% pgstar_win_file_ptr(i_Kipp_residuals)
+ p% plot => Kipp_residuals_plot
+ p% id = i_Kipp_residuals
+ p% name = 'Kipp_resid'
+ p% win_flag = s% pg% Kipp_residuals_win_flag
+ p% win_width = s% pg% Kipp_residuals_win_width
+ p% win_aspect_ratio = s% pg% Kipp_residuals_win_aspect_ratio
+ p% file_flag = s% pg% Kipp_residuals_file_flag
+ p% file_dir = s% pg% Kipp_residuals_file_dir
+ p% file_prefix = s% pg% Kipp_residuals_file_prefix
+ p% file_interval = s% pg% Kipp_residuals_file_interval
+ p% file_width = s% pg% Kipp_residuals_file_width
+ p% file_aspect_ratio = s% pg% Kipp_residuals_file_aspect_ratio
+
do i = 1, max_num_Other_plots
p => s% pg% pgstar_win_file_ptr(i_Other + i - 1)
p% win_flag = .false.
@@ -1343,8 +1375,8 @@ subroutine onScreen_Plots(s, must_write_files_in, ierr)
use utils_lib
use chem_def
use net_def
- use net_lib, only : get_net_reaction_table
- use const_def, only : Msun, Rsun
+ use net_lib, only: get_net_reaction_table
+ use const_def, only: Msun, Rsun
type (star_info), pointer :: s
logical :: must_write_files_in
@@ -1412,7 +1444,7 @@ subroutine onScreen_Plots(s, must_write_files_in, ierr)
do i = 1, num_pgstar_plots
p => s% pg% pgstar_win_file_ptr(i)
- if(show_plot_now) then
+ if (show_plot_now) then
! call to check_window opens device
call check_window(s, p, ierr)
if (failed('check_window')) return
@@ -1568,7 +1600,7 @@ end subroutine read_pgstar_data
subroutine update_pgstar_data(s, ierr)
- use star_utils, only : eval_csound
+ use star_utils, only: eval_csound
type (star_info), pointer :: s
integer, intent(out) :: ierr
@@ -1595,7 +1627,7 @@ subroutine update_pgstar_data(s, ierr)
contains
subroutine get_hist_values(num, ierr)
- use history, only : do_get_data_for_history_columns
+ use history, only: do_get_data_for_history_columns
integer, intent(in) :: num
integer, intent(out) :: ierr
integer :: i
diff --git a/star/private/pgstar_grid.f90 b/star/private/pgstar_grid.f90
index 98a4a2465..0c3735888 100644
--- a/star/private/pgstar_grid.f90
+++ b/star/private/pgstar_grid.f90
@@ -566,6 +566,10 @@ subroutine Grid_plot(s, id, device_id, &
do_Color_Magnitude1_plot, do_Color_Magnitude2_plot, do_Color_Magnitude3_plot, &
do_Color_Magnitude4_plot, do_Color_Magnitude5_plot, do_Color_Magnitude6_plot, &
do_Color_Magnitude7_plot, do_Color_Magnitude8_plot, do_Color_Magnitude9_plot
+ use pgstar_equation_residuals, only : &
+ do_Max_eq_resid_plot
+ use pgstar_kipp_residuals, only : &
+ do_Kipp_residuals_plot
type (star_info), pointer :: s
logical, intent(in) :: subplot
@@ -947,6 +951,15 @@ subroutine Grid_plot(s, id, device_id, &
call do_Text_Summary9_plot(&
s, id, device_id, xleft, xright, ybot, ytop, grid_subplot, s% pg% Text_Summary9_title, &
Grid_txt_scale_factor(i) * s% pg% Text_Summary9_txt_scale, s% pg% Text_Summary9_dxval, ierr)
+ case ('max_eq_resid')
+ call do_Max_eq_resid_plot(&
+ s, id, xleft, xright, ybot, ytop, grid_subplot, s% pg% Max_eq_resid_title, &
+ Grid_txt_scale_factor(i) * s% pg% Max_eq_resid_txt_scale, s% pg% Max_eq_resid_max_width, ierr)
+ case ('kipp_residuals')
+ call do_Kipp_residuals_plot(&
+ s, id, xleft, xright, ybot, ytop, grid_subplot, s% pg% Kipp_residuals_title, &
+ Grid_txt_scale_factor(i) * s% pg% Kipp_residuals_txt_scale, &
+ s% pg% Kipp_residuals_max_width, s% pg% Kipp_residuals_yaxis_name, ierr)
case default
! check for "other" plot
found_it = .false.
@@ -998,6 +1011,7 @@ subroutine Grid_plot(s, id, device_id, &
'History_Panels1,..,9', &
'History_Tracks1,..,9', &
'Color_Magnitude1,..,9', &
+ 'Max_eq_resid', &
'and if you are using star/astero', &
'Echelle', &
'Ratios'
diff --git a/star/private/pgstar_kipp_residuals.f90 b/star/private/pgstar_kipp_residuals.f90
new file mode 100644
index 000000000..7d99fc9eb
--- /dev/null
+++ b/star/private/pgstar_kipp_residuals.f90
@@ -0,0 +1,377 @@
+! ***********************************************************************
+!
+! Copyright (C) 2010-2019 The MESA Team
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU Lesser General Public License
+! as published by the Free Software Foundation,
+! either version 3 of the License, or (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+! See the GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with this program. If not, see .
+!
+! ***********************************************************************
+
+module pgstar_kipp_residuals
+
+ use const_def
+ use star_private_def
+ use star_pgstar
+ use pgstar_colors
+ use const_def, only: dp, Msun
+
+ implicit none
+
+ real, parameter, dimension(9) :: &
+ l = (/ 0.00, 0.12, 0.25, 0.38, 0.50, 0.62, 0.75, 0.88, 1.00 /), &
+ r = (/ 0.00, 0.05, 0.15, 0.32, 0.55, 0.78, 0.93, 0.99, 1.00 /), &
+ g = (/ 0.00, 0.02, 0.04, 0.06, 0.10, 0.18, 0.32, 0.58, 0.95 /), &
+ b = (/ 0.00, 0.08, 0.18, 0.32, 0.36, 0.28, 0.16, 0.08, 0.80 /)
+ real, parameter :: contra = 1.0, bright = 0.5
+
+contains
+
+ ! Kippenhahn-style log10(max(|res|)) across all equations for a
+ ! given space-time location
+ subroutine Kipp_residuals_plot(id, device_id, ierr)
+ integer, intent(in) :: id, device_id
+ integer, intent(out) :: ierr
+
+ type (star_info), pointer :: s
+
+ ierr = 0
+ call get_star_ptr(id, s, ierr)
+ if (ierr /= 0) return
+
+ call pgslct(device_id)
+ call pgbbuf()
+ call pgeras()
+
+ call do_Kipp_residuals_Plot(s, id, &
+ s% pg% Kipp_residuals_xleft, s% pg% Kipp_residuals_xright, &
+ s% pg% Kipp_residuals_ybot, s% pg% Kipp_residuals_ytop, .false., &
+ s% pg% Kipp_residuals_title, s% pg% Kipp_residuals_txt_scale, &
+ s% pg% Kipp_residuals_max_width, s% pg% Kipp_residuals_yaxis_name, ierr)
+ if (ierr /= 0) return
+
+ call pgebuf()
+
+ end subroutine Kipp_residuals_plot
+
+
+ subroutine do_Kipp_residuals_plot(s, id, &
+ win_xleft, win_xright, win_ybot, win_ytop, subplot, title, txt_scale, &
+ max_width, yaxis_name, &
+ ierr)
+ integer, intent(in) :: id, max_width
+ real, intent(in) :: win_xleft, win_xright, win_ybot, win_ytop, txt_scale
+ logical, intent(in) :: subplot
+ character(len=*), intent(in) :: title, yaxis_name
+ integer, intent(out) :: ierr
+
+ type (star_info), pointer :: s
+ integer :: k, neq, initial_width, n, x_size, y_size, nr_n_cells
+ integer, save :: init_model
+
+ real :: resid_hi, resid_lo
+
+ real, allocatable, dimension(:, :), save :: nr_resid_buf, nr_ycoord_buf
+ integer, allocatable, dimension(:), save :: nr_model_buf, nr_zone_buf
+
+ ! colorbar color mapping info
+ real, allocatable, save :: r_save(:), g_save(:), b_save(:)
+ integer, save :: n_colors, c_low, c_high
+
+ real, allocatable :: tmp_resid_buf(:, :), tmp_ycoord_buf(:, :)
+ integer, allocatable :: tmp_model_buf(:), tmp_zone_buf(:)
+
+ neq = size(s% equ, 1) ! number of equations per cell
+
+ if (.not. allocated(nr_resid_buf)) then
+ if (max_width < 0) then
+ initial_width = 10
+ else
+ initial_width = max_width
+ end if
+ ! make 5 times bigger than current nz
+ allocate(nr_resid_buf(initial_width, int(1.5 * s% nz)))
+ allocate(nr_ycoord_buf(initial_width, int(1.5 * s% nz)))
+ allocate(nr_model_buf(initial_width))
+ allocate(nr_zone_buf(initial_width))
+ nr_resid_buf(:, :) = -99.0
+ nr_ycoord_buf(:, :) = 0.0 ! zero-init so unused cells are safely masked
+
+ init_model = s% model_number - 1
+ else if (s% nz > 0.95 * size(nr_resid_buf, 2)) then
+ ! check if remeshing has increased nz beyond 95% of the
+ ! current array size
+ x_size = size(nr_resid_buf, 1)
+ call move_alloc(nr_resid_buf, tmp_resid_buf)
+ call move_alloc(nr_ycoord_buf, tmp_ycoord_buf)
+
+ allocate(nr_resid_buf(x_size, 5 * s% nz))
+ allocate(nr_ycoord_buf(x_size, 5 * s% nz))
+ nr_resid_buf(:, :) = -99.0
+ nr_ycoord_buf(:, :) = 0.0 ! zero-init so unused cells are safely masked
+
+ nr_resid_buf(:, 1:nr_n_cells) = tmp_resid_buf(:, 1:nr_n_cells)
+ nr_ycoord_buf(:, 1:nr_n_cells) = tmp_ycoord_buf(:, 1:nr_n_cells)
+
+ deallocate(tmp_resid_buf, tmp_ycoord_buf)
+ end if
+
+ nr_n_cells = s% nz
+ if (max_width < 0) then
+ n = s% model_number - init_model
+ else
+ n = min(s% model_number - init_model, max_width)
+ end if
+
+ ! resize n_model dimension
+ if (n > size(nr_resid_buf, 1)) then
+ x_size = size(nr_resid_buf, 1)
+ y_size = size(nr_resid_buf, 2)
+
+ call move_alloc(nr_resid_buf, tmp_resid_buf)
+ call move_alloc(nr_ycoord_buf, tmp_ycoord_buf)
+ call move_alloc(nr_zone_buf, tmp_zone_buf)
+ call move_alloc(nr_model_buf, tmp_model_buf)
+
+ allocate(nr_resid_buf(2 * x_size, y_size))
+ allocate(nr_ycoord_buf(2 * x_size, y_size))
+ allocate(nr_zone_buf(2 * x_size))
+ allocate(nr_model_buf(2 * x_size))
+
+ nr_resid_buf(:, :) = -99.0
+ nr_ycoord_buf(:, :) = 0.0 ! zero-init so unused cells are safely masked
+
+ nr_resid_buf(1:n, :) = tmp_resid_buf(1:n, :)
+ nr_ycoord_buf(1:n, :) = tmp_ycoord_buf(1:n, :)
+ nr_zone_buf(1:n) = tmp_zone_buf(1:n)
+ nr_model_buf(1:n) = tmp_model_buf(1:n)
+
+ deallocate(tmp_resid_buf, tmp_ycoord_buf, tmp_zone_buf, tmp_model_buf)
+ end if
+
+ if (max_width > 1 .and. s% model_number - init_model > max_width) then
+ nr_model_buf(1:max_width-1) = nr_model_buf(2:max_width)
+ nr_resid_buf(1:max_width-1, :) = nr_resid_buf(2:max_width, :)
+ nr_ycoord_buf(1:max_width-1, :) = nr_ycoord_buf(2:max_width, :)
+ nr_zone_buf(1:max_width-1) = nr_zone_buf(2:max_width)
+ end if
+
+ nr_model_buf(n) = s% model_number
+ nr_zone_buf(n) = s% nz
+
+ ! fill arrays with new vals
+ do k = 1, s% nz
+ nr_resid_buf(n, k) = real(safe_log10(maxval(abs(s% equ(1:neq, k)))))
+ select case (trim(yaxis_name))
+ case ('mass')
+ nr_ycoord_buf(n, k) = real(s% m(k)/Msun)
+ case ('logR')
+ if (k < s%nz) then
+ nr_ycoord_buf(n, k) = real(safe_log10(0.5d0 * (s% r(k) + s% r(k+1)) / Rsun))
+ else
+ nr_ycoord_buf(n, k) = real(safe_log10(0.5d0 * s% r(k) / Rsun))
+ end if
+ case ('tau')
+ ! optical depth
+ nr_ycoord_buf(n, k) = real(safe_log10(s% tau(k)))
+ case default
+ ! cell-centered mass coordinate
+ nr_ycoord_buf(n, k) = real((s% m(k) - 0.5d0 * s% dm(k))/Msun)
+ end select
+ end do
+
+ ! so actual plotting
+ call pgsave
+ call pgsci(0) ! color index 0 = background
+ call pgsvp(win_xleft, win_xright, win_ybot, win_ytop)
+ call pgrect(0.0, 1.0, 0.0, 1.0) ! pseudo erase
+ call Kipp_residuals_render(ierr, id)
+
+ call pgunsa
+
+ contains
+
+ subroutine Kipp_residuals_render(ierr, id)
+ integer, intent(in) :: id
+ integer, intent(out) :: ierr
+ real, parameter :: missing_val = TINY(1.0)
+ real, allocatable :: img(:,:), coord_grid(:)
+ real :: tr(6), fg, bg, xlo, xhi, dx
+ real :: clo, chi, dcoord, coord_j, frac
+ integer :: nx, ny, i, j, k, nz_i
+ logical :: y_reverse
+
+ if (.not. allocated(r_save)) then
+ call pgqcir(c_low, c_high)
+ n_colors = c_high - c_low + 1
+ allocate(r_save(n_colors), g_save(n_colors), b_save(n_colors))
+
+ ! save existing color table
+ do i = 1, n_colors
+ call pgqcr(c_low + i - 1, r_save(i), g_save(i), b_save(i))
+ end do
+ end if
+
+ call pgqcir(c_low, c_high)
+ call pgscir(16,255)
+ call pgctab(l, r, g, b, 9, contra, bright)
+
+ nx = n
+ ny = nr_n_cells
+
+ allocate(img(nx, ny))
+ allocate(coord_grid(ny))
+
+ chi = -1.0e30
+ clo = 1.0e30
+ do i = 1, nx
+ nz_i = nr_zone_buf(i)
+ do k = 1, nz_i
+ if (nr_ycoord_buf(i,k) > chi) chi = nr_ycoord_buf(i,k)
+ ! m(k) is positive; skip any zero-initialised slots
+ if (nr_ycoord_buf(i,k) > 0.0 .and. nr_ycoord_buf(i,k) < clo) &
+ clo = nr_ycoord_buf(i,k)
+ end do
+ end do
+
+ resid_hi = -1.0e30
+ resid_lo = 1.0e30
+ do i = 1, nx
+ nz_i = nr_zone_buf(i)
+ do k = 1, nz_i
+ if (nr_resid_buf(i,k) > resid_hi) resid_hi = nr_resid_buf(i,k)
+ if (nr_resid_buf(i,k) < resid_lo) resid_lo = nr_resid_buf(i,k)
+ end do
+ end do
+
+ select case (yaxis_name)
+ case ('mass')
+ y_reverse = .false.
+ case ('logR')
+ y_reverse = .false.
+ case ('tau')
+ y_reverse = .true.
+ case default
+ y_reverse = .false.
+ end select
+
+ dcoord = (chi - clo) / real(ny - 1)
+ do j = 1, ny
+ if (.not. y_reverse) then
+ coord_grid(j) = clo + real(j - 1) * dcoord
+ else
+ coord_grid(j) = chi - real(j - 1) * dcoord
+ end if
+ end do
+
+ ! --- interpolate residuals onto the uniform grid ---
+ img = missing_val
+ do i = 1, nx
+ nz_i = nr_zone_buf(i)
+ do j = 1, ny
+ coord_j = coord_grid(j)
+ do k = 1, nz_i - 1
+ ! skip invalid cells
+ if (nr_resid_buf(i,k) == missing_val) cycle
+ if (nr_resid_buf(i,k+1) == missing_val) cycle
+ ! works for either increasing or decreasing coordinates
+ if ((nr_ycoord_buf(i,k) - coord_j) * &
+ (nr_ycoord_buf(i,k+1) - coord_j) <= 0d0) then
+ if (abs(nr_ycoord_buf(i,k) - nr_ycoord_buf(i,k+1)) > 1d-99) then
+ frac = (coord_j - nr_ycoord_buf(i,k+1)) / &
+ (nr_ycoord_buf(i,k) - nr_ycoord_buf(i,k+1))
+ else
+ frac = 0.5d0
+ end if
+ img(i,j) = nr_resid_buf(i,k+1) + &
+ frac * (nr_resid_buf(i,k) - nr_resid_buf(i,k+1))
+ exit
+ end if
+ end do
+ end do
+ end do
+ ! --- x range centred on model numbers ---
+ dx = max(1.0, real(nr_model_buf(nx) - nr_model_buf(1)) / real(nx - 1))
+ xlo = real(nr_model_buf(1)) - 0.5*dx
+ xhi = real(nr_model_buf(nx)) + 0.5*dx
+
+ ! --- transformation matrix for pgimag ---
+ ! world_y(j) = tr(4) + tr(6)*j
+ ! j=1 -> mhi => tr(4) = mhi - tr(6) = mhi + dm
+ ! j=ny -> mlo => confirms tr(6) = -dm
+ tr(1) = xlo - 0.5
+ tr(2) = (xhi - xlo) / real(nx)
+ tr(3) = 0.0
+ tr(4) = chi + dcoord
+ tr(5) = 0.0
+ tr(6) = -dcoord
+
+ ! color scale
+ fg = resid_hi
+ bg = resid_lo
+ if (fg <= -100) fg = maxval(img)
+ if (bg <= -100) bg = minval(img, mask=(img > -19.9))
+ if (bg >= fg) bg = fg - 6.0
+
+ ! --- main panel ---
+ ! y-axis: mlo at bottom (centre), mhi at top (surface)
+ if (.not. y_reverse) then
+ call pgswin(xlo, xhi, clo - 0.5*dcoord, chi + 0.5*dcoord)
+ else
+ call pgswin(xlo, xhi, chi + 0.5*dcoord, clo - 0.5*dcoord)
+ end if
+ call pgimag(img, nx, ny, 1, nx, 1, ny, fg, bg, tr)
+
+ call pgsci(clr_Foreground)
+ call pgsch(txt_scale)
+ call pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0)
+
+ select case (yaxis_name)
+ case ('mass')
+ y_reverse = .false.
+ call pglab('Model Number', 'Mass (M\d\(2281)\u)', &
+ 'Newton Raphson Solver Residuals of accepted step')
+ case ('logR')
+ y_reverse = .false.
+ call pglab('Model Number', 'log(R/R\d\(2281)\u)', &
+ 'Newton Raphson Solver Residuals of accepted step')
+ case ('tau')
+ y_reverse = .true.
+ call pglab('Model Number', 'log(tau)', &
+ 'Newton Raphson Solver Residuals of accepted step')
+ case default
+ y_reverse = .false.
+ call pglab('Model Number', 'Mass (M\d\(2281)\u)', &
+ 'Newton Raphson Solver Residuals of accepted step')
+ end select
+
+ ! --- color wedge ---
+ call pgsvp(win_xright + 0.01 * (win_xright - win_xleft), win_xright + 0.03 * (win_xright - win_xleft), &
+ win_ybot, win_ytop)
+ call pgswin(0.0, 1.0, bg, fg)
+ call pgwedg('RI', 0.0, 4.0, bg, fg, 'log10(max(|res|))')
+
+ ! reset color range
+ call pgscir(c_low, c_high)
+ do i = 1, n_colors
+ call pgscr(c_low + i - 1, r_save(i), g_save(i), b_save(i))
+ end do
+
+ deallocate(img, coord_grid)
+
+ end subroutine Kipp_residuals_render
+
+ end subroutine do_Kipp_residuals_plot
+
+
+
+end module pgstar_kipp_residuals
diff --git a/star/work/inlist_pgstar b/star/work/inlist_pgstar
index 2be909017..2a18cc736 100644
--- a/star/work/inlist_pgstar
+++ b/star/work/inlist_pgstar
@@ -1,37 +1,6 @@
&pgstar
- ! see star/defaults/pgstar.defaults
+Kipp_residuals_win_flag = .true.
- ! MESA uses PGPLOT for live plotting and gives the user a tremendous
- ! amount of control of the presentation of the information.
-
- ! show HR diagram
- ! this plots the history of L,Teff over many timesteps
- HR_win_flag = .true.
-
- ! set static plot bounds
- HR_logT_min = 3.5
- HR_logT_max = 4.6
- HR_logL_min = 2.0
- HR_logL_max = 6.0
-
- ! set window size (aspect_ratio = height/width)
- HR_win_width = 6
- HR_win_aspect_ratio = 1.0
-
-
- ! show temperature/density profile
- ! this plots the internal structure at single timestep
- TRho_Profile_win_flag = .true.
-
- ! add legend explaining colors
- show_TRho_Profile_legend = .true.
-
- ! display numerical info about the star
- show_TRho_Profile_text_info = .true.
-
- ! set window size (aspect_ratio = height/width)
- TRho_Profile_win_width = 8
- TRho_Profile_win_aspect_ratio = 0.75
/ ! end of pgstar namelist
diff --git a/star/work/inlist_project b/star/work/inlist_project
index 6e9a86824..97e3ded9f 100644
--- a/star/work/inlist_project
+++ b/star/work/inlist_project
@@ -6,78 +6,78 @@
&star_job
- ! see star/defaults/star_job.defaults
+! see star/defaults/star_job.defaults
- ! begin with a pre-main sequence model
- create_pre_main_sequence_model = .true.
+! begin with a pre-main sequence model
+create_pre_main_sequence_model = .false.
- ! save a model at the end of the run
- save_model_when_terminate = .false.
- save_model_filename = '15M_at_TAMS.mod'
+! save a model at the end of the run
+save_model_when_terminate = .false.
+save_model_filename = '15M_at_TAMS.mod'
- ! display on-screen plots
- pgstar_flag = .true.
+! display on-screen plots
+pgstar_flag = .true.
/ ! end of star_job namelist
&eos
- ! eos options
- ! see eos/defaults/eos.defaults
+! eos options
+! see eos/defaults/eos.defaults
/ ! end of eos namelist
&kap
- ! kap options
- ! see kap/defaults/kap.defaults
- use_Type2_opacities = .true.
- Zbase = 0.02
+! kap options
+! see kap/defaults/kap.defaults
+use_Type2_opacities = .true.
+Zbase = 0.02
/ ! end of kap namelist
&controls
- ! see star/defaults/controls.defaults
+! see star/defaults/controls.defaults
- ! starting specifications
- initial_mass = 15 ! in Msun units
- initial_z = 0.02
+! starting specifications
+initial_mass = 15 ! in Msun units
+initial_z = 0.02
- ! when to stop
+! when to stop
- ! stop when the star nears ZAMS (Lnuc/L > 0.99)
- Lnuc_div_L_zams_limit = 0.99d0
- stop_near_zams = .true.
+! ! stop when the star nears ZAMS (Lnuc/L > 0.99)
+! Lnuc_div_L_zams_limit = 0.99d0
+! stop_near_zams = .true.
- ! stop when the center mass fraction of h1 drops below this limit
- xa_central_lower_limit_species(1) = 'h1'
- xa_central_lower_limit(1) = 1d-3
+! stop when the center mass fraction of h1 drops below this limit
+xa_central_lower_limit_species(1) = 'h1'
+xa_central_lower_limit(1) = 1d-4
- ! wind
+! wind
- ! atmosphere
+! atmosphere
- ! rotation
+! rotation
- ! element diffusion
+! element diffusion
- ! mlt
+! mlt
- ! mixing
+! mixing
- ! timesteps
+! timesteps
- ! mesh
+! mesh
- ! solver
- ! options for energy conservation (see MESA V, Section 3)
- energy_eqn_option = 'dedt'
- use_gold_tolerances = .true.
+! solver
+! options for energy conservation (see MESA V, Section 3)
+energy_eqn_option = 'dedt'
+use_gold_tolerances = .true.
- ! output
+! output
/ ! end of controls namelist
diff --git a/star_data/private/pgstar_controls.inc b/star_data/private/pgstar_controls.inc
index 6098721b2..9aa9c4cea 100644
--- a/star_data/private/pgstar_controls.inc
+++ b/star_data/private/pgstar_controls.inc
@@ -2047,6 +2047,38 @@
integer, dimension(max_num_pgstar_grid_plots) :: &
Grid9_plot_row, Grid9_plot_rowspan, &
Grid9_plot_col, Grid9_plot_colspan
+
+ logical :: Max_eq_resid_win_flag, Max_eq_resid_file_flag
+ logical :: do_Max_eq_resid_win, do_Max_eq_resid_file
+ integer :: Max_eq_resid_file_interval
+ integer :: id_Max_eq_resid_win, id_Max_eq_resid_file
+ character (len=strlen) :: Max_eq_resid_file_dir, Max_eq_resid_file_prefix
+ real :: &
+ Max_eq_resid_win_width, Max_eq_resid_win_aspect_ratio, &
+ Max_eq_resid_xleft, Max_eq_resid_xright, &
+ Max_eq_resid_ybot, Max_eq_resid_ytop, &
+ prev_Max_eq_resid_win_width, prev_Max_eq_resid_win_ratio, &
+ Max_eq_resid_file_width, Max_eq_resid_file_aspect_ratio, &
+ prev_Max_eq_resid_file_width, prev_Max_eq_resid_file_ratio, &
+ Max_eq_resid_txt_scale
+ character (len=strlen) :: Max_eq_resid_title
+ integer :: Max_eq_resid_max_width
+
+ logical :: Kipp_residuals_win_flag, Kipp_residuals_file_flag
+ logical :: do_Kipp_residuals_win, do_Kipp_residuals_file
+ integer :: Kipp_residuals_file_interval
+ integer :: id_Kipp_residuals_win, id_Kipp_residuals_file
+ character (len=strlen) :: Kipp_residuals_file_dir, Kipp_residuals_file_prefix
+ real :: &
+ Kipp_residuals_win_width, Kipp_residuals_win_aspect_ratio, &
+ Kipp_residuals_xleft, Kipp_residuals_xright, &
+ Kipp_residuals_ybot, Kipp_residuals_ytop, &
+ prev_Kipp_residuals_win_width, prev_Kipp_residuals_win_ratio, &
+ Kipp_residuals_file_width, Kipp_residuals_file_aspect_ratio, &
+ prev_Kipp_residuals_file_width, prev_Kipp_residuals_file_ratio
+ character (len=strlen) :: Kipp_residuals_title, Kipp_residuals_yaxis_name
+ real :: Kipp_residuals_txt_scale
+ integer :: Kipp_residuals_max_width
integer :: annotation1_ci, annotation1_lw, annotation1_cf
@@ -2133,3 +2165,12 @@
trho_use_decorator, &
tmaxrho_use_decorator, &
trho_profile_use_decorator
+
+
+ ! for Kippenhahn-like visualization of max residual across all equations
+ logical :: Kipp_residuals_win_flag, Kipp_residuals_file_flag
+ real :: Kipp_residuals_win_width, Kipp_residuals_win_aspect_ratio, &
+ Kipp_residuals_file_width, Kipp_residuals_file_aspect_ratio, &
+ Kipp_residuals_min, Kipp_residuals_max
+ character (len=strlen) :: Kipp_residuals_file_dir, Kipp_residuals_file_prefix
+ integer :: Kipp_residuals_file_interval
diff --git a/star_data/public/star_pgstar.f90 b/star_data/public/star_pgstar.f90
index a1b7197a0..cb1b4857f 100644
--- a/star_data/public/star_pgstar.f90
+++ b/star_data/public/star_pgstar.f90
@@ -175,7 +175,11 @@ end subroutine pgstar_decorator_interface
integer, parameter :: i_Grid7 = i_Grid6 + 1
integer, parameter :: i_Grid8 = i_Grid7 + 1
integer, parameter :: i_Grid9 = i_Grid8 + 1
- integer, parameter :: i_Other = i_Grid9 + 1
+
+ integer, parameter :: i_Max_eq_resid = i_Grid9 + 1
+ integer, parameter :: i_Kipp_residuals = i_Max_eq_resid + 1
+
+ integer, parameter :: i_Other = i_Kipp_residuals + 1
integer, parameter :: num_pgstar_plots = i_Other + max_num_Other_plots