Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2026 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module pgstar_equation_residuals
21 : ! Composition equations are not included.
22 : use star_private_def
23 : use const_def, only: dp
24 : use pgstar_support
25 : use star_pgstar
26 : use pgstar_colors
27 :
28 : implicit none
29 :
30 : private
31 : public :: Max_eq_resid_plot, do_Max_eq_resid_plot, Max_eq_resid_record
32 :
33 : ! same line colors as pgstar_abundance
34 : integer, parameter :: line_colors(14) = [ &
35 : clr_Gold, clr_LightSkyBlue, clr_Crimson, clr_LightSkyGreen, clr_MediumSlateBlue, &
36 : clr_Coral, clr_Goldenrod, clr_DarkGray, clr_Lilac, &
37 : clr_Tan, clr_IndianRed, clr_Teal, clr_Silver, clr_BrightBlue ]
38 :
39 : real(dp), parameter :: resid_floor = 1d-20, resid_ceil = 1d30
40 : integer, parameter :: default_capacity = 128
41 :
42 : ! per-star history; indexed by s% id so multiple stars (binaries) don't mix
43 : type resid_hist_t
44 : integer :: nvar = 0
45 : integer :: count = 0
46 : integer, allocatable :: model(:)
47 : real, allocatable :: vals(:, :) ! (nvar, capacity)
48 : character(len=name_len), allocatable :: names(:)
49 : end type resid_hist_t
50 :
51 : type(resid_hist_t), target, save :: hist(max_star_handles)
52 :
53 : contains
54 :
55 0 : subroutine Max_eq_resid_plot(id, device_id, ierr)
56 : integer, intent(in) :: id, device_id
57 : integer, intent(out) :: ierr
58 : type (star_info), pointer :: s
59 :
60 : ierr = 0
61 0 : call get_star_ptr(id, s, ierr)
62 0 : if (ierr /= 0) return
63 :
64 0 : call pgslct(device_id)
65 0 : call pgbbuf()
66 0 : call pgeras()
67 :
68 : call do_Max_eq_resid_plot(s, id, &
69 : s% pg% Max_eq_resid_xleft, s% pg% Max_eq_resid_xright, &
70 : s% pg% Max_eq_resid_ybot, s% pg% Max_eq_resid_ytop, .false., &
71 : s% pg% Max_eq_resid_title, s% pg% Max_eq_resid_txt_scale, &
72 0 : s% pg% Max_eq_resid_max_width, ierr)
73 0 : if (ierr /= 0) return
74 :
75 0 : call pgebuf()
76 :
77 : end subroutine Max_eq_resid_plot
78 :
79 :
80 0 : subroutine Max_eq_resid_record(s, max_width, ierr)
81 : type (star_info), pointer :: s
82 : integer, intent(in) :: max_width ! <= 0: unlimited
83 : integer, intent(out) :: ierr
84 :
85 : type(resid_hist_t), pointer :: h
86 : integer :: nvar, nz, i_eq, k, c
87 :
88 0 : ierr = 0
89 0 : if (s% id < 1 .or. s% id > max_star_handles) then
90 0 : ierr = -1
91 0 : return
92 : end if
93 0 : h => hist(s% id)
94 :
95 0 : nvar = s% nvar_hydro
96 0 : nz = s% nz
97 0 : if (nvar <= 0 .or. nz <= 0 .or. .not. associated(s% equ)) return
98 0 : if (size(s% equ, 1) < nvar .or. size(s% equ, 2) < nz) return
99 :
100 : ! (re)initialize on first call or if the set of equations changed
101 0 : if (.not. allocated(h% model) .or. h% nvar /= nvar) then
102 0 : if (max_width > 0) then
103 0 : c = max_width
104 : else
105 0 : c = default_capacity
106 : end if
107 0 : call reset_hist(h, nvar, c, ierr)
108 0 : if (ierr /= 0) return
109 0 : do i_eq = 1, nvar
110 0 : h% names(i_eq) = s% nameofequ(i_eq)
111 : end do
112 : end if
113 :
114 : ! drop entries that are not older than the current model
115 0 : do while (h% count > 0)
116 0 : if (h% model(h% count) < s% model_number) exit
117 0 : h% count = h% count - 1
118 : end do
119 :
120 0 : if (max_width > 0 .and. h% count >= max_width) then
121 : ! roll: keep the newest max_width-1 entries
122 0 : k = h% count - max_width + 1
123 0 : h% model(1:h% count - k) = h% model(k + 1:h% count)
124 0 : h% vals(:, 1:h% count - k) = h% vals(:, k + 1:h% count)
125 0 : h% count = h% count - k
126 : end if
127 :
128 0 : if (h% count >= size(h% model)) then
129 0 : call grow_hist(h, 2 * size(h% model), ierr)
130 0 : if (ierr /= 0) return
131 : end if
132 :
133 0 : h% count = h% count + 1
134 0 : h% model(h% count) = s% model_number
135 0 : do i_eq = 1, nvar
136 : h% vals(i_eq, h% count) = real(min(resid_ceil, &
137 0 : max(resid_floor, maxval(abs(s% equ(i_eq, 1:nz))))))
138 : end do
139 :
140 : end subroutine Max_eq_resid_record
141 :
142 :
143 0 : subroutine do_Max_eq_resid_plot(s, id, &
144 : win_xleft, win_xright, win_ybot, win_ytop, subplot, title, txt_scale, &
145 : max_width, ierr)
146 :
147 : type (star_info), pointer :: s
148 : integer, intent(in) :: id, max_width
149 : logical, intent(in) :: subplot
150 : real, intent(in) :: win_xleft, win_xright, win_ybot, win_ytop, txt_scale
151 : character (len=*), intent(in) :: title
152 : integer, intent(out) :: ierr
153 :
154 : type(resid_hist_t), pointer :: h
155 0 : real, allocatable :: xvec(:), yvec(:)
156 : real :: xmin, xmax, ymin, ymax, dy_leg
157 : integer :: i_eq, n
158 :
159 : ierr = 0
160 0 : call Max_eq_resid_record(s, max_width, ierr)
161 0 : if (ierr /= 0) return
162 :
163 0 : h => hist(s% id)
164 0 : n = h% count
165 0 : if (n < 1) return
166 :
167 0 : allocate(xvec(n), yvec(n))
168 0 : xvec = real(h% model(1:n))
169 :
170 0 : xmin = xvec(1)
171 0 : xmax = xvec(n)
172 0 : if (xmin >= xmax) xmax = xmin + 1.0
173 :
174 0 : ymin = log10(minval(h% vals(:, 1:n)))
175 0 : ymax = log10(maxval(h% vals(:, 1:n)))
176 0 : if (ymax - ymin < 1.0) then
177 0 : ymin = ymin - 0.5
178 0 : ymax = ymax + 0.5
179 : end if
180 0 : ymin = ymin - 0.05 * (ymax - ymin)
181 0 : ymax = ymax + 0.05 * (ymax - ymin)
182 :
183 0 : call pgsave
184 0 : call pgsci(clr_Background) ! background
185 0 : call pgsvp(win_xleft, win_xright, win_ybot, win_ytop)
186 0 : call pgrect(0.0, 1.0, 0.0, 1.0) ! pseudo erase
187 0 : call pgswin(xmin, xmax, ymin, ymax)
188 0 : call pgsch(txt_scale)
189 0 : call pgscf(1)
190 0 : call pgsci(clr_Foreground)
191 0 : call pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0)
192 0 : call pgmtxt('B', 2.5, 0.5, 0.5, 'Model Number')
193 0 : call pgmtxt('L', 3.0, 0.5, 0.5, 'log\d10\u max\dk\u |residual|')
194 0 : call pgmtxt('T', 1.0, 0.5, 0.5, title)
195 :
196 0 : do i_eq = 1, h% nvar
197 0 : yvec = log10(h% vals(i_eq, 1:n))
198 0 : call pgsci(line_colors(mod(i_eq - 1, size(line_colors)) + 1))
199 0 : if (n == 1) then
200 0 : call pgpt(1, xvec, yvec, 17)
201 : else
202 0 : call pgline(n, xvec, yvec)
203 : end if
204 : end do
205 :
206 : ! legend
207 0 : call pgsch(0.75 * txt_scale)
208 0 : dy_leg = min(0.08, 1.0 / real(h% nvar))
209 0 : do i_eq = 1, h% nvar
210 0 : call pgsci(line_colors(mod(i_eq - 1, size(line_colors)) + 1))
211 : call pgmtxt('RV', 0.5, 1.0 - (real(i_eq) - 0.5) * dy_leg, 0.0, &
212 0 : trim(h% names(i_eq)))
213 : end do
214 :
215 0 : call pgunsa
216 0 : deallocate(xvec, yvec)
217 :
218 0 : end subroutine do_Max_eq_resid_plot
219 :
220 :
221 0 : subroutine reset_hist(h, nvar, capacity, ierr)
222 : type(resid_hist_t), intent(inout) :: h
223 : integer, intent(in) :: nvar, capacity
224 : integer, intent(out) :: ierr
225 :
226 0 : if (allocated(h% model)) deallocate(h% model)
227 0 : if (allocated(h% vals)) deallocate(h% vals)
228 0 : if (allocated(h% names)) deallocate(h% names)
229 0 : allocate(h% model(capacity), h% vals(nvar, capacity), h% names(nvar), stat=ierr)
230 0 : if (ierr /= 0) then
231 0 : write(*, *) 'Max_eq_resid: allocate failed'
232 0 : return
233 : end if
234 0 : h% nvar = nvar
235 0 : h% count = 0
236 : end subroutine reset_hist
237 :
238 :
239 0 : subroutine grow_hist(h, new_capacity, ierr)
240 : type(resid_hist_t), intent(inout) :: h
241 : integer, intent(in) :: new_capacity
242 : integer, intent(out) :: ierr
243 0 : integer, allocatable :: tmp_model(:)
244 0 : real, allocatable :: tmp_vals(:, :)
245 :
246 0 : call move_alloc(h% model, tmp_model)
247 0 : call move_alloc(h% vals, tmp_vals)
248 0 : allocate(h% model(new_capacity), h% vals(h% nvar, new_capacity), stat=ierr)
249 0 : if (ierr /= 0) then
250 0 : write(*, *) 'Max_eq_resid: allocate failed'
251 0 : return
252 : end if
253 0 : h% model(1:h% count) = tmp_model(1:h% count)
254 0 : h% vals(:, 1:h% count) = tmp_vals(:, 1:h% count)
255 :
256 : ! explicitly free temporaries
257 0 : if (allocated(tmp_model)) deallocate(tmp_model)
258 0 : if (allocated(tmp_vals)) deallocate(tmp_vals)
259 :
260 0 : end subroutine grow_hist
261 :
262 0 : end module pgstar_equation_residuals
263 :
|
