Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 mesh_plan
21 :
22 : use const_def, only: dp, convective_mixing
23 : use num_lib
24 : use utils_lib
25 : use star_private_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: do_mesh_plan
31 :
32 : logical, parameter :: plan_dbg = .false.
33 : integer, parameter :: kdbg = -1
34 :
35 : contains
36 :
37 10 : subroutine do_mesh_plan( &
38 : s, nz_old, max_allowed_nz, okay_to_merge, D_mix, &
39 : max_k_old_for_split_in, min_k_old_for_split_in, xq_old, dq_old, &
40 : min_dq_in, max_dq, min_dq_for_split, mesh_max_allowed_ratio, &
41 : do_not_split, num_gvals, gval_names, &
42 : gval_is_xa_function, gval_is_logT_function, gvals, &
43 : delta_gval_max, max_center_cell_dq, max_surface_cell_dq, min_surface_cell_dq, &
44 : max_num_subcells, max_num_merge_cells, max_num_merge_surface_cells, &
45 : nz_new, xq_new, dq_new, which_gval, comes_from, ierr)
46 : ! return keep_going, or terminate
47 : use mesh_functions, only: max_allowed_gvals
48 : ! inputs
49 : type (star_info), pointer :: s
50 : integer, intent(in) :: nz_old, max_allowed_nz, max_num_subcells, &
51 : max_num_merge_cells, max_num_merge_surface_cells, max_k_old_for_split_in, &
52 : min_k_old_for_split_in
53 : logical, intent(in) :: okay_to_merge
54 : real(dp), pointer :: D_mix(:) ! (nz_old)
55 : real(dp), pointer :: xq_old(:) ! (nz_old)
56 : real(dp), pointer :: dq_old(:) ! (nz_old)
57 : real(dp), intent(in) :: min_dq_in, max_dq, min_dq_for_split, mesh_max_allowed_ratio
58 : logical, pointer :: do_not_split(:)
59 : integer, intent(in) :: num_gvals
60 : character (len=32) :: gval_names(max_allowed_gvals)
61 : logical, dimension(max_allowed_gvals) :: gval_is_xa_function, gval_is_logT_function
62 : real(dp), pointer :: gvals(:,:) ! (nz_old, num_gvals)
63 : real(dp), pointer :: delta_gval_max(:) ! (nz_old)
64 : real(dp), intent(in) :: max_center_cell_dq, max_surface_cell_dq, &
65 : min_surface_cell_dq
66 : ! outputs
67 : integer, intent(out) :: nz_new
68 : real(dp), pointer :: xq_new(:), dq_new(:) ! (nz_new)
69 : ! must be allocated on entry; suggested size >= nz_old.
70 : ! reallocated as necessary if need to enlarge.
71 : ! size on return is >= nz_new.
72 : integer, pointer :: which_gval(:) ! (nz_new) for debugging.
73 : ! which_gval(k) = gval number that set the size for new cell k.
74 : ! size may have been reduced below gradient setting by other restrictions.
75 : integer, pointer :: comes_from(:) ! (nz_new)
76 : ! xq_old(comes_from(k)+1) > xq_new(k) >= xq_old(comes_from(k)), if comes_from(k) < nz_old.
77 : integer, intent(out) :: ierr
78 :
79 : integer :: j, k, k_old, k_new, nz, new_capacity, iounit, species, &
80 : max_k_old_for_split, min_k_old_for_split
81 10 : real(dp) :: D_mix_cutoff, next_xq, next_dq, max_dq_cntr, &
82 10 : dq_sum, min_dq, min_dq_for_xa, min_dq_for_logT
83 :
84 : logical, parameter :: write_plan_debug = .false.
85 :
86 : include 'formats'
87 :
88 10 : ierr = 0
89 :
90 10 : min_dq = min_dq_in
91 :
92 10 : if (max_k_old_for_split_in < 0) then
93 0 : max_k_old_for_split = nz_old + max_k_old_for_split_in
94 : else
95 10 : max_k_old_for_split = max_k_old_for_split_in
96 : end if
97 :
98 10 : if (min_k_old_for_split_in < 0) then
99 0 : min_k_old_for_split = nz_old + min_k_old_for_split_in
100 : else
101 10 : min_k_old_for_split = min_k_old_for_split_in
102 : end if
103 :
104 10 : if (max_dq < min_dq) then
105 0 : write(*,1) 'ERROR in controls: max_dq < min_dq', max_dq, min_dq
106 0 : ierr = -1
107 0 : return
108 : end if
109 :
110 10 : if (max_center_cell_dq < min_dq) then
111 0 : write(*,1) 'ERROR in controls: max_center_cell_dq < min_dq', max_center_cell_dq, min_dq
112 0 : ierr = -1
113 0 : return
114 : end if
115 :
116 10 : if (max_surface_cell_dq < min_dq) then
117 0 : write(*,1) 'ERROR in controls: max_surface_cell_dq < min_dq', max_surface_cell_dq, min_dq
118 0 : ierr = -1
119 0 : return
120 : end if
121 :
122 12095 : do k = 2, nz_old
123 12095 : if (xq_old(k) <= xq_old(k-1)) then
124 0 : write(*,3) 'bad xq_old', k, nz_old, xq_old(k), xq_old(k-1)
125 0 : ierr = -1
126 0 : return
127 : end if
128 : end do
129 :
130 10 : if (s% set_min_D_mix) then
131 : D_mix_cutoff = s% min_D_mix
132 : else
133 10 : D_mix_cutoff = 0
134 : end if
135 :
136 : if (write_plan_debug) call open_debug_file
137 :
138 : if (plan_dbg) then
139 : write(*,*) 'num_gvals', num_gvals
140 : do j=1,num_gvals
141 : write(*,*) j, trim(gval_names(j))
142 : end do
143 : write(*,'(A)')
144 : end if
145 :
146 10 : nz = nz_old
147 10 : species = s% species
148 :
149 : new_capacity = min(size(xq_new,dim=1), size(dq_new,dim=1), &
150 10 : size(which_gval,dim=1), size(comes_from,dim=1))
151 10 : max_dq_cntr = max(1d-20, max_center_cell_dq, 0.75d0*dq_old(nz_old))
152 :
153 12105 : comes_from(1:nz) = 0
154 10 : comes_from(1) = 1
155 :
156 10 : call pick_new_points(s, ierr)
157 10 : if (ierr /= 0) return
158 :
159 10 : if (min_k_old_for_split <= 1) then
160 10 : do while (dq_new(1) > max(max_surface_cell_dq,2*min_dq,min_dq_for_split))
161 0 : call split1(1, ierr)
162 10 : if (ierr /= 0) then
163 0 : write(*,*) 'failed trying to split surface cell'
164 0 : stop
165 : return
166 : end if
167 : end do
168 : end if
169 :
170 10 : call smooth_new_points(ierr) ! split as necessary
171 10 : if (ierr /= 0) return
172 :
173 : if (write_plan_debug) then
174 : close(iounit)
175 : end if
176 :
177 10 : if (ierr /= 0) return
178 :
179 :
180 : contains
181 :
182 :
183 : subroutine check_before_smooth(ierr)
184 : integer, intent(out) :: ierr
185 : integer :: k, k_old
186 : real(dp) :: xq0, xqold_start, xqold_end
187 :
188 : include 'formats'
189 : ierr = 0
190 : if (comes_from(1) /= 1 .or. xq_old(1) /= xq_new(1)) then
191 : write(*,*) 'mesh plan check_before_smooth: bad value from k=1'
192 : ierr = -1
193 : return
194 : end if
195 :
196 : do k = 2, nz_new
197 : xq0 = xq_new(k)
198 : k_old = comes_from(k)
199 : xqold_start = xq_old(k_old)
200 : if (k_old == nz_old) then
201 : xqold_end = 1d0
202 : else
203 : xqold_end = xq_old(k_old+1)
204 : end if
205 : if (k_old > comes_from(k-1)) then
206 : if (xq0 /= xqold_start) then
207 : write(*,'(A)')
208 : write(*,2) 'nz_new', nz_new
209 : write(*,2) 'k', k
210 : write(*,2) 'comes_from(k-1)', comes_from(k-1)
211 : write(*,2) 'k_old', k_old
212 : write(*,2) 'comes_from(k)', comes_from(k)
213 : write(*,2) 'nz_old', nz_old
214 : write(*,'(A)')
215 : write(*,2) 'xq0', k, xq0
216 : write(*,2) 'xqold_start', k_old, xqold_start
217 : write(*,2) 'xqold_end', k_old, xqold_end
218 : write(*,2) 'xqold_end - xqold_start', k_old, xqold_end - xqold_start
219 : write(*,2) 'xq0 - xqold_start', k_old, xq0 - xqold_start
220 : write(*,'(A)')
221 : write(*,2) 'xq_old(comes_from(k))', k, xq_old(comes_from(k))
222 : write(*,2) 'xq_old(comes_from(k-1))', k, xq_old(comes_from(k-1))
223 : write(*,'(A)')
224 : write(*,*) 'k_old > comes_from(k-1)', k_old > comes_from(k-1)
225 : write(*,*) 'xq0 /= xqold_start', xq0 /= xqold_start
226 : write(*,*) 'mesh plan check_before_smooth'
227 : ierr = -1
228 : return
229 : end if
230 : else if (k_old == comes_from(k-1)) then
231 : if (.not. (xq0 > xqold_start .and. xq0 < xqold_end)) then
232 : write(*,'(A)')
233 : write(*,*) '(.not. (xq0 > xqold_start .and. xq0 < xqold_end))', k_old
234 : write(*,2) 'xq_new(k)', k, xq_new(k)
235 : write(*,2) 'xq_old(k_old)', k, xq_old(k_old)
236 : write(*,2) 'xq_old(k_old+1)', k, xq_old(k_old+1)
237 : write(*,2) 'xqold_end', k, xqold_end
238 : write(*,'(A)')
239 : write(*,*) 'k_old == comes_from(k-1)', k_old == comes_from(k-1)
240 : write(*,*) '.not. (xq0 > xqold_start .and. xq0 < xqold_end)', &
241 : .not. (xq0 > xqold_start .and. xq0 < xqold_end)
242 : write(*,*) 'mesh plan check_before_smooth'
243 : ierr = -1
244 : return
245 : end if
246 : else
247 : write(*,*) 'comes_from(k) > comes_from(k-1)', k, comes_from(k), comes_from(k-1)
248 : write(*,*) 'mesh plan check_before_smooth'
249 : ierr = -1
250 : return
251 : end if
252 :
253 : cycle ! use the following for debugging
254 : if (dq_new(k-1) < 1d-6*dq_new(k) .or. dq_new(k) < 1d-6*dq_new(k-1)) then
255 : write(*,3) 'bad dq_new ratio', k, nz_new, dq_new(k)/dq_new(k-1), dq_new(k), dq_new(k-1)
256 : write(*,*) 'check_before_smooth'
257 : ierr = -1
258 : return
259 : end if
260 :
261 : end do
262 :
263 10 : end subroutine check_before_smooth
264 :
265 :
266 : subroutine test_new(ierr)
267 : integer, intent(out) :: ierr
268 : integer :: k, k_old
269 : include 'formats'
270 : ierr = 0
271 : do k = 1, nz_new
272 : if (dq_new(k) <= 0) then
273 : write(*,3) 'bad dq_new'
274 : write(*,2) 'dq_new', k, dq_new(k)
275 : write(*,2) 'dq_new', k-1, dq_new(k-1)
276 : write(*,2) 'xq_new', k, xq_new(k)
277 : write(*,2) 'xq_new', k-1, xq_new(k-1)
278 : write(*,3) 'comes_from', k, comes_from(k)
279 : write(*,3) 'comes_from', k-1, comes_from(k-1)
280 : write(*,2) 'nz_new', nz_new
281 : write(*,2) 'nz_old', nz_old
282 : write(*,2) 'dq_old(comes_from(k))', comes_from(k), dq_old(comes_from(k))
283 : write(*,2) 'dq_old(comes_from(k)-1)', comes_from(k)-1, dq_old(comes_from(k)-1)
284 : write(*,2) 'dq_old(comes_from(k)-2)', comes_from(k)-2, dq_old(comes_from(k)-2)
285 : write(*,2) 'xq_old(comes_from(k))', comes_from(k), xq_old(comes_from(k))
286 : write(*,2) 'xq_old(comes_from(k)-1)', comes_from(k)-1, xq_old(comes_from(k)-1)
287 : write(*,2) 'xq_old(comes_from(k)-2)', comes_from(k)-2, xq_old(comes_from(k)-2)
288 : write(*,*) 'test_new: mesh plan'
289 : ierr = -1
290 : return
291 : end if
292 : k_old = comes_from(k)
293 : if (k_old < 1 .or. k_old > nz_old) then
294 : write(*,*) 'bad value for comes_from', k, k_old
295 : write(*,*) 'test_new: mesh plan'
296 : ierr = -1
297 : return
298 : end if
299 : if (xq_new(k) < xq_old(k_old)) then
300 : write(*,3) '(xq_new(k) < xq_old(k_old))', k, k_old, &
301 : xq_new(k) - xq_old(k_old), xq_new(k), xq_old(k_old)
302 : write(*,*) 'test_new: mesh plan'
303 : ierr = -1
304 : return
305 : end if
306 : if (k_old < nz_old) then
307 : if (xq_new(k) > xq_old(k_old+1)) then
308 : write(*,3) '(xq_new(k) > xq_old(k_old+1))', k, k_old+1, &
309 : xq_new(k) - xq_old(k_old+1), xq_new(k), xq_old(k_old+1)
310 : write(*,*) 'test_new: mesh plan'
311 : ierr = -1
312 : return
313 : end if
314 : end if
315 : end do
316 : end subroutine test_new
317 :
318 :
319 128 : subroutine split1(k, ierr)
320 : integer, intent(in) :: k
321 : integer, intent(out) :: ierr
322 : integer :: kk, k_old, k_old_last, split_at_k_old
323 128 : real(dp) :: xq_mid, xq_end
324 : logical :: from_merger, dbg
325 :
326 : include 'formats'
327 128 : ierr = 0
328 :
329 128 : k_old = comes_from(k)
330 128 : k_old_last = -1
331 128 : xq_end = -1
332 128 : from_merger = (xq_new(k) == xq_old(k_old) .and. dq_new(k) > dq_old(k_old) + min_dq*1d-3)
333 :
334 128 : dbg = (k_old == -1)
335 :
336 128 : if (dbg) then
337 0 : write(*,2) 'start split1 k', k
338 0 : write(*,2) 'nz_old', nz_old
339 0 : write(*,2) 'nz_new', nz_new
340 0 : do kk=1400,nz_new
341 0 : if (dq_new(kk) < 1d-12) then
342 0 : write(*,2) 'dq_new(kk)', kk, dq_new(kk)
343 0 : call mesa_error(__FILE__,__LINE__,'debug: split1')
344 : end if
345 : end do
346 0 : do kk = 2, nz_new
347 0 : if (dq_new(kk-1) < 1d-6*dq_new(kk) .or. dq_new(kk) < 1d-6*dq_new(kk-1)) then
348 0 : write(*,3) 'bad dq_new ratio', kk, nz_new, dq_new(kk), dq_new(kk-1)
349 0 : call mesa_error(__FILE__,__LINE__,'debug: split1')
350 : end if
351 : end do
352 : end if
353 :
354 128 : if (from_merger) then ! find range of old cells that were merged to form k
355 128 : if (k == nz_new) then
356 0 : xq_end = 1
357 0 : k_old_last = nz_old
358 : else
359 128 : xq_end = xq_new(k+1)
360 128 : k_old_last = 0
361 256 : do kk = k_old+1, nz_old ! find last old cell included in k
362 256 : if (xq_old(kk) == xq_end) then
363 128 : k_old_last = kk-1; exit
364 : end if
365 128 : if (xq_old(kk) > xq_end) then
366 0 : write(*,*) 'oops'
367 0 : write(*,2) 'xq_old(kk)', kk, xq_old(kk)
368 0 : write(*,2) 'xq_old(kk-1)', kk-1, xq_old(kk-1)
369 0 : write(*,2) 'xq_end', k, xq_end
370 0 : write(*,*) 'split1'
371 0 : ierr = -1
372 0 : return
373 : end if
374 : end do
375 128 : if (k_old_last == k_old) then
376 : from_merger = .false.
377 128 : else if (k_old_last < k_old) then
378 0 : write(*,*) 'confusion in split1 for k_old_last'
379 0 : write(*,2) 'k_old', k_old
380 0 : write(*,2) 'k_old_last', k_old_last
381 0 : write(*,2) 'nz_old', nz_old
382 0 : write(*,2) 'k', k
383 0 : write(*,2) 'nz_new', nz_new
384 0 : write(*,'(A)')
385 0 : write(*,2) 'dq_new(k)', k, dq_new(k)
386 0 : write(*,2) 'dq_old(k_old)', k_old, dq_old(k_old)
387 0 : write(*,2) 'dq_new(k)-dq_old(k_old)', k_old, dq_new(k)-dq_old(k_old)
388 0 : write(*,'(A)')
389 0 : write(*,2) 'xq_new(k)', k, xq_new(k)
390 0 : write(*,2) 'xq_new(k)+dq_new(k)', k, xq_new(k)+dq_new(k)
391 0 : write(*,2) 'xq_end', k, xq_end
392 0 : write(*,*) 'split1'
393 0 : ierr = -1
394 0 : return
395 : end if
396 : end if
397 : end if
398 :
399 128 : if (nz_new == new_capacity) then ! increase allocated size
400 0 : new_capacity = (new_capacity*5)/4 + 10
401 0 : call realloc(s, nz_new, new_capacity, xq_new, dq_new, which_gval, comes_from, ierr)
402 0 : if (ierr /= 0) return
403 : end if
404 128 : nz_new = nz_new + 1
405 128 : if (max_allowed_nz > 0 .and. nz_new > max_allowed_nz) then
406 0 : write(*,*) 'tried to increase number of mesh points beyond max allowed nz', max_allowed_nz
407 0 : ierr = -1
408 0 : return
409 : end if
410 17004 : do kk = nz_new, k+1, -1
411 16876 : xq_new(kk) = xq_new(kk-1)
412 16876 : dq_new(kk) = dq_new(kk-1)
413 16876 : which_gval(kk) = which_gval(kk-1)
414 17004 : comes_from(kk) = comes_from(kk-1)
415 : end do
416 :
417 128 : if (from_merger) then ! split by breaking up the merger
418 128 : xq_mid = xq_new(k) + dq_new(k)*0.5d0
419 128 : split_at_k_old = k_old_last
420 146 : do kk = k_old+1, k_old_last ! check interior boundaries for closest to xq_mid
421 146 : if (xq_old(kk) >= xq_mid) then
422 110 : if (xq_old(kk) - xq_mid < xq_mid - xq_old(kk-1)) then ! xq_mid closer to xq_old(kk)
423 : split_at_k_old = kk
424 : else ! xq_mid closer to xq_old(kk-1)
425 0 : split_at_k_old = kk-1
426 : end if
427 : exit
428 : end if
429 : end do
430 128 : comes_from(k+1) = split_at_k_old
431 128 : xq_new(k+1) = xq_old(split_at_k_old)
432 256 : dq_new(k) = sum(dq_old(k_old:split_at_k_old-1))
433 256 : dq_new(k+1) = sum(dq_old(split_at_k_old:k_old_last))
434 : else
435 0 : dq_new(k:k+1) = dq_new(k)*0.5d0
436 0 : xq_new(k+1) = xq_new(k) + dq_new(k)
437 : ! fix up comes_from(k+1)
438 0 : comes_from(k+1) = nz_old ! just in case
439 0 : do k_old = comes_from(k), nz_old-1
440 0 : if (xq_new(k+1) <= xq_old(k_old+1)) then
441 0 : comes_from(k+1) = k_old
442 0 : exit
443 : end if
444 : end do
445 : end if
446 :
447 128 : if (dbg) then
448 0 : write(*,2) 'split1 dq_new k', k, dq_new(k)
449 0 : write(*,2) 'split1 dq_new k+1', k+1, dq_new(k+1)
450 0 : write(*,2) 'comes_from', k_old
451 0 : write(*,2) 'nz_old', nz_old
452 0 : write(*,2) 'nz_new', nz_new
453 0 : do kk=1400,nz_new
454 0 : if (dq_new(kk) < 1d-12) then
455 0 : write(*,2) 'dq_new(kk)', kk, dq_new(kk)
456 0 : call mesa_error(__FILE__,__LINE__,'debug: split1')
457 : end if
458 : end do
459 : end if
460 :
461 : end subroutine split1
462 :
463 :
464 31907 : logical function okay_to_split1(k_old, dq_new, remaining_dq_old)
465 : integer, intent(in) :: k_old
466 : real(dp), intent(in) :: dq_new, remaining_dq_old
467 31907 : real(dp) :: dr_old, min_dr, rR, rL, dlnR_new
468 : logical :: dbg
469 :
470 : include 'formats'
471 :
472 31907 : dbg = .false.
473 31907 : okay_to_split1 = .false.
474 31907 : if (do_not_split(k_old)) return
475 :
476 31907 : if (max(dq_new,remaining_dq_old) < max(min_dq_for_split,2*min_dq)) return
477 :
478 : if (dbg) then
479 : write(*,2) 'remaining_dq_old', k_old, remaining_dq_old
480 : write(*,1) 'dq_new', dq_new
481 : write(*,1) 'min_dq_for_split', min_dq_for_split
482 : write(*,1) 'min_dq', min_dq
483 : write(*,'(A)')
484 : end if
485 :
486 31907 : if (0d0 < remaining_dq_old .and. dq_new > 0.99d0*remaining_dq_old) then
487 : if (dbg) then
488 : write(*,2) 'dq_old', k_old, remaining_dq_old
489 : write(*,1) 'dq_new', dq_new
490 : write(*,1) 'dq_new/remaining_dq_old', dq_new/remaining_dq_old
491 : write(*,'(A)')
492 : end if
493 : return
494 : end if
495 :
496 31907 : if (k_old < nz_old) then
497 :
498 31907 : rR = s% r(k_old)
499 31907 : rL = s% r(k_old+1)
500 31907 : dr_old = rR - rL
501 31907 : min_dr = s% csound(k_old)*s% mesh_min_dr_div_cs
502 31907 : if (dr_old*dq_new/dq_old(k_old) < 2*min_dr) then
503 : return ! sound crossing time would be too small
504 : end if
505 :
506 31907 : min_dr = s% mesh_min_dr_div_dRstar*(s% r(1) - s% R_center)
507 31907 : if (dr_old*dq_new/dq_old(k_old) < 2*min_dr) then
508 : return ! new dr would be too small
509 : end if
510 :
511 31907 : if (s% mesh_min_dlnR > 0d0) then
512 31907 : dlnR_new = log((rL + dr_old*dq_new/dq_old(k_old))/rL)
513 31907 : if (dlnR_new < 2*s% mesh_min_dlnR) then
514 : ! the factor of 2 is a safety margin.
515 : return
516 : end if
517 : end if
518 :
519 : end if
520 31907 : okay_to_split1 = .true.
521 :
522 : end function okay_to_split1
523 :
524 :
525 10 : subroutine smooth_new_points(ierr)
526 : integer, intent(out) :: ierr
527 :
528 : logical :: dbg, done
529 : integer :: k, k_old
530 10 : real(dp) :: alfa
531 :
532 : include 'formats'
533 :
534 10 : ierr = 0
535 10 : dbg = .false.
536 10 : alfa = mesh_max_allowed_ratio
537 :
538 : do ! repeat until nothing left to do
539 : if (min_k_old_for_split <= 1 .and. &
540 20 : dq_new(1) > alfa*dq_new(2) .and. &
541 : okay_to_split1(1,dq_new(1),0d0)) then
542 0 : call split1(1,ierr)
543 0 : if (ierr /= 0) return
544 : cycle
545 : end if
546 20 : done = .true.
547 20 : k = 2
548 : do ! check for cell that is too large
549 21294 : if (k == nz_new) exit
550 21274 : k_old = comes_from(k)
551 : if (k_old <= max_k_old_for_split .and. &
552 : k_old >= min_k_old_for_split .and. &
553 21274 : okay_to_split1(k_old,dq_new(k),0d0) .and. &
554 : (dq_new(k) > max_dq .or. &
555 : dq_new(k) > alfa*dq_new(k+1) .or. &
556 20 : dq_new(k) > alfa*dq_new(k-1))) then
557 98 : call split1(k,ierr)
558 98 : if (ierr /= 0) return
559 : done = .false.
560 : ! don't increment k; want to recheck the cell in case need to split again
561 : else
562 21176 : k = k + 1
563 : end if
564 : end do
565 :
566 : if (dbg) then
567 : call test_new(ierr)
568 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'debug: mesh_plan, smooth_new_points')
569 : end if
570 :
571 20 : if (done) exit
572 :
573 : ! now go in opposite direction
574 10 : done = .true.
575 10 : k = nz_new-1
576 : do
577 10643 : if (k == 1) exit
578 10633 : if (okay_to_split1(comes_from(k),dq_new(k),0d0) .and. &
579 10 : (dq_new(k) > alfa*dq_new(k+1) .or. dq_new(k) > alfa*dq_new(k-1))) then
580 30 : call split1(k,ierr)
581 30 : if (ierr /= 0) return
582 30 : done = .false.
583 30 : k = k + 1 ! recheck the same cell
584 : else
585 10603 : k = k - 1
586 : end if
587 : end do
588 :
589 : if (dbg) then
590 : call test_new(ierr)
591 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'debug: mesh_plan, smooth_new_points')
592 : end if
593 :
594 20 : if (done) exit
595 :
596 : end do
597 :
598 : end subroutine smooth_new_points
599 :
600 :
601 10 : subroutine pick_new_points(s, ierr)
602 : type (star_info), pointer :: s
603 : integer, intent(out) :: ierr
604 :
605 : logical :: dbg, force_merge_with_one_more, du_div_cs_limit_flag
606 10 : real(dp) :: maxval_delta_xa, next_dq_max, beta_limit, &
607 10 : remaining_dq_old, min_dr, abs_du_div_cs
608 : integer :: kk, k_old_init, k_old_next, k_old_next_max, j00, jm1, i, max_merge
609 :
610 : include 'formats'
611 :
612 10 : beta_limit = 0.1d0
613 :
614 10 : ierr = 0
615 10 : k_old = 1
616 10 : k_new = 1
617 10 : xq_new(1) = 0
618 10 : min_dr = s% mesh_min_dr_div_dRstar*(s% r(1) - s% R_center)
619 10 : abs_du_div_cs = 0d0
620 :
621 : do ! pick next point location
622 :
623 10495 : dbg = plan_dbg .or. (k_new == kdbg) !.or. (s% mesh_call_number == 2005)
624 :
625 : ! when reach this point,
626 : ! have set xq_new(k) for k = 1 to k_new, and dq_new(k) for k = 1 to k_new-1.
627 : ! and have finished using old points from 1 to k_old
628 : ! i.e., xq_old(k_old+1) > xq_new(k_new) >= xq_old(k_old)
629 :
630 10495 : k_old_init = k_old
631 :
632 10495 : if (k_new == 1) then
633 10 : next_dq_max = max_surface_cell_dq
634 : else
635 10485 : next_dq_max = dq_new(k_new-1)*mesh_max_allowed_ratio
636 : end if
637 :
638 : ! make initial choice based on gradients. may reduce later.
639 10495 : if (dbg) then
640 0 : write(*,'(A)')
641 0 : write(*,3) 'call pick_next_dq', k_old, k_new, next_dq_max
642 : end if
643 :
644 10495 : if (s% gradr(k_old) > s% grada(k_old) .and. &
645 : s% min_dq_for_xa_convective > 0d0) then
646 1325 : min_dq_for_xa = s% min_dq_for_xa_convective
647 : else
648 9170 : min_dq_for_xa = s% min_dq_for_xa
649 : end if
650 :
651 10495 : min_dq_for_logT = s% min_dq_for_logT
652 :
653 : next_dq = pick_next_dq(s, &
654 : dbg, next_dq_max, k_old, k_new, nz_old, num_gvals, &
655 : xq_new, dq_new, xq_old, dq_old, min_dq, min_dq_for_split, &
656 : min_dq_for_xa, min_dq_for_logT, &
657 : max_surface_cell_dq, max_dq_cntr, max_num_subcells, &
658 : gval_is_xa_function, gval_is_logT_function, gvals, &
659 10495 : delta_gval_max, gval_names, which_gval, ierr)
660 10495 : if (ierr /= 0) return
661 :
662 10495 : if (dbg) &
663 0 : write(*,3) 'pick_next_dq next_dq prev_dq', k_old, k_new, next_dq, dq_old(k_old)
664 :
665 10495 : if (k_new == 1) then
666 10 : max_merge = max_num_merge_surface_cells
667 : else
668 10485 : max_merge = max_num_merge_cells
669 : end if
670 :
671 10495 : if (k_old < nz_old) then
672 10485 : remaining_dq_old = xq_old(k_old+1)-xq_new(k_new)
673 : else
674 10 : remaining_dq_old = 1d0-xq_new(k_new)
675 : end if
676 :
677 10495 : if (next_dq < remaining_dq_old) then
678 0 : if (.not. okay_to_split1(k_old, next_dq, remaining_dq_old)) then
679 0 : next_dq = remaining_dq_old ! dq_old(k_old)
680 : end if
681 : end if
682 :
683 10495 : if (next_dq > max_dq) then
684 21 : if (xq_new(k_new) == xq_old(k_old)) then
685 2623 : do i = nz_old, k_old+1, -1
686 2623 : if (xq_old(i) - xq_new(k_new) <= max_dq) then
687 21 : next_dq = xq_old(i) - xq_new(k_new)
688 21 : exit
689 : end if
690 : end do
691 : end if
692 : end if
693 :
694 10495 : if (next_dq > max_dq) then
695 0 : next_dq = xq_old(k_old) + dq_old(k_old) - xq_new(k_new)
696 : end if
697 :
698 10495 : if (k_new == 1 .and. next_dq > max_surface_cell_dq) then
699 0 : if (k_old_init == -1 .and. .true.) write(*,*) 'next_dq > max_surface_cell_dq'
700 0 : next_dq = max_surface_cell_dq
701 : end if
702 :
703 : ! Enforce minimum surface‐cell dq
704 10495 : if (k_old == 1 .and. next_dq < min_surface_cell_dq) then
705 : !write(*,*) 'next_dq < min_surface_cell_dq'
706 0 : next_dq = min_surface_cell_dq
707 : end if
708 :
709 10495 : next_xq = xq_new(k_new) + next_dq
710 10495 : if (next_xq > 1d0 - min_dq) then
711 10 : next_xq = (1d0 + xq_new(k_new))/2d0
712 10 : if (k_old < nz_old) then ! make sure don't split current k_old for this case
713 0 : if (xq_old(k_old+1) > next_xq) next_xq = xq_old(k_old+1)
714 : end if
715 10 : next_dq = next_xq - xq_new(k_new)
716 : end if
717 :
718 10495 : if (k_old < nz_old) then
719 :
720 10485 : if (xq_new(k_new) == xq_old(k_old) .and. &
721 : next_dq > dq_old(k_old) - min_dq/2d0) then
722 :
723 10485 : if (.not. okay_to_merge) then
724 :
725 0 : k_old_next = k_old + 1
726 :
727 10485 : else if (k_old < min_k_old_for_split .or. &
728 : k_old > max_k_old_for_split) then
729 :
730 0 : k_old_next = k_old + 1
731 :
732 : else ! consider doing merge
733 :
734 10485 : if (next_dq > 1.5d0*dq_old(k_old)) then
735 7240 : next_dq = 0.9d0*next_dq ! to avoid split-merge flip-flops
736 7240 : next_xq = xq_new(k_new) + next_dq
737 : end if
738 :
739 10485 : k_old_next_max = min(nz_old, k_old + max_merge)
740 10485 : k_old_next = k_old_next_max ! will cut this back as necessary
741 22526 : do kk=k_old+1,k_old_next_max
742 :
743 20926 : if (s% use_hydro_merge_limits_in_mesh_plan) then ! limit merges over steep velocity gradients
744 : ! begin hydro check_merge_limits section
745 0 : du_div_cs_limit_flag = .false.
746 0 : if (.not. s% merge_amr_du_div_cs_limit_only_for_compression) then
747 : du_div_cs_limit_flag = .true.
748 0 : else if (associated(s% v)) then
749 : ! Only set flag for compressive flow across interface (kk-1, kk)
750 0 : if (kk <= nz_old .and. s% v(kk)*pow2(s% r(kk)) > s% v(kk-1)*pow2(s% r(kk-1))) &
751 0 : du_div_cs_limit_flag = .true.
752 0 : if (.not. du_div_cs_limit_flag .and. kk-1 > 1) then
753 0 : if (s% v(kk-1)*pow2(s% r(kk-1)) > s% v(kk-2)*pow2(s% r(kk-2))) du_div_cs_limit_flag = .true.
754 : end if
755 : end if
756 :
757 0 : if (du_div_cs_limit_flag .and. associated(s% v)) then
758 0 : if (kk-1 == 1) then
759 0 : abs_du_div_cs = abs(s% v(1) - s% v(2)) / s% csound(1)
760 0 : else if (kk == nz_old) then
761 0 : abs_du_div_cs = abs(s% v(nz_old-1) - s% v(nz_old)) / s% csound(nz_old)
762 : else
763 0 : abs_du_div_cs = abs(s% v(kk-1) - s% v(kk)) / s% csound(kk-1)
764 : end if
765 : else
766 : abs_du_div_cs = 0.0_dp
767 : end if
768 :
769 : if (du_div_cs_limit_flag) then
770 0 : if (.not. s% merge_amr_inhibit_at_jumps) then
771 0 : if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) then
772 : ! Shock jump too large, so block merge at this interface
773 0 : k_old_next = kk-1
774 0 : exit
775 : end if
776 : else ! inhibit_at_jumps = .true.
777 0 : if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) then
778 0 : if (dq_old(k_old) >= min_dq) then
779 : ! Shock is large and zone is not extremely small, so inhibit merge
780 0 : k_old_next = kk-1
781 0 : exit
782 : end if
783 : ! else: if zone is very small, allow merge despite the shock
784 : end if
785 : end if
786 : end if
787 : end if
788 : ! end hydro check_merge_limits section
789 :
790 209260 : maxval_delta_xa = maxval(abs(s% xa(:,kk)-s% xa(:,kk-1)))
791 188334 : j00 = maxloc(s% xa(:,kk),dim=1)
792 188334 : jm1 = maxloc(s% xa(:,kk-1),dim=1)
793 : if (maxval_delta_xa > s% max_delta_x_for_merge .or. &
794 20926 : j00 /= jm1 .or. is_convective_boundary(kk) .or. &
795 1600 : is_crystal_boundary(kk)) then
796 : ! don't merge across convective or crystal boundary
797 58 : k_old_next = kk-1
798 58 : exit
799 20868 : else if (next_xq <= xq_old(kk) + min_dq/2d0) then
800 8827 : k_old_next = max(k_old+1,kk-1)
801 8827 : exit
802 : end if
803 : end do
804 :
805 : end if
806 :
807 10485 : k_old_next = max(k_old_next, k_old+1)
808 10485 : next_xq = xq_old(k_old_next)
809 10485 : next_dq = next_xq - xq_new(k_new)
810 :
811 0 : else if (next_xq >= xq_old(k_old+1) - min_dq/2d0) then
812 : ! this is final subcell of a split, so adjust to finish the parent cell
813 :
814 0 : k_old_next = k_old+1
815 0 : next_xq = xq_old(k_old_next)
816 0 : next_dq = next_xq - xq_new(k_new)
817 :
818 : else ! non-final subcell of split
819 :
820 0 : k_old_next = k_old
821 :
822 : end if
823 :
824 10485 : if (next_xq == xq_old(k_old_next)) then
825 : ! finishing 1 or more old cells and not already at max merge.
826 : ! consider forcing a merge with the next cell to make this one larger.
827 10485 : force_merge_with_one_more = .false.
828 :
829 10485 : if (dq_old(k_old) < min_dq) then
830 : force_merge_with_one_more = .true.
831 10485 : else if (s% merge_if_dlnR_too_small) then
832 0 : if (xq_new(k_new) <= xq_old(k_old) .and. &
833 : s% lnR(k_old) - s% lnR(k_old_next) < s% mesh_min_dlnR) then
834 : force_merge_with_one_more = .true.
835 0 : else if (k_old_next == nz_old .and. s% R_center > 0) then
836 0 : force_merge_with_one_more = s% lnR(k_old_next) - log(s% R_center) < s% mesh_min_dlnR
837 : end if
838 : end if
839 :
840 10485 : if ((.not. force_merge_with_one_more) .and. s% merge_if_dr_div_cs_too_small) then
841 10485 : if (xq_new(k_new) <= xq_old(k_old) .and. &
842 : s% r(k_old) - s% r(k_old_next) < s% mesh_min_dr_div_cs*s% csound(k_old)) then
843 : force_merge_with_one_more = .true.
844 10485 : else if (k_old_next == nz_old) then
845 : force_merge_with_one_more = (s% r(k_old_next) - s% R_center) < &
846 10 : s% mesh_min_dr_div_cs*s% csound(k_old_next)
847 10 : if (force_merge_with_one_more .and. dbg) &
848 0 : write(*,3) 'do merge for k_old_next == nz_old', k_old, k_old_next, &
849 0 : s% r(k_old) - s% R_center, &
850 0 : s% mesh_min_dr_div_cs*s% csound(k_old_next), &
851 0 : s% mesh_min_dr_div_cs, s% csound(k_old_next)
852 : end if
853 : end if
854 :
855 10485 : if ((.not. force_merge_with_one_more) .and. &
856 : s% merge_if_dr_div_dRstar_too_small) then
857 10485 : if (xq_new(k_new) <= xq_old(k_old) .and. &
858 : s% r(k_old) - s% r(k_old_next) < min_dr) then
859 : force_merge_with_one_more = .true.
860 10485 : else if (k_old_next == nz_old) then
861 10 : force_merge_with_one_more = (s% r(k_old_next) - s% R_center) < min_dr
862 10 : if (force_merge_with_one_more .and. dbg) &
863 0 : write(*,3) 'do merge for k_old_next == nz_old', k_old, k_old_next, &
864 0 : s% r(k_old) - s% R_center, min_dr
865 : end if
866 : end if
867 :
868 10485 : if (force_merge_with_one_more) then
869 0 : k_old_next = k_old_next + 1
870 0 : if (k_old_next < nz_old) then
871 0 : next_xq = xq_old(k_old_next)
872 : else
873 0 : next_xq = 1d0
874 : end if
875 0 : next_dq = next_xq - xq_new(k_new)
876 0 : if (dbg) then
877 0 : write(*,3) 'force merge', k_old, k_new, next_xq, next_dq
878 0 : write(*,'(A)')
879 : end if
880 : end if
881 :
882 : end if
883 :
884 10485 : k_old = k_old_next
885 :
886 : end if
887 :
888 10495 : comes_from(k_new) = k_old_init
889 :
890 : ! check if we're done
891 10495 : if (1 - xq_new(k_new) < max_dq_cntr .or. 1 - next_xq < min_dq) then
892 10 : dq_new(k_new) = 1 - xq_new(k_new)
893 : exit
894 : end if
895 :
896 10485 : dq_new(k_new) = next_dq
897 :
898 10485 : if (k_new == new_capacity) then ! increase allocated size
899 0 : new_capacity = (new_capacity*5)/4 + 10
900 0 : call realloc(s, k_new, new_capacity, xq_new, dq_new, which_gval, comes_from, ierr)
901 0 : if (ierr /= 0) return
902 : end if
903 :
904 10485 : if (next_xq < xq_new(k_new)) then
905 0 : write(*,*) 'nz_old', nz_old
906 0 : write(*,*) 'k_new', k_new
907 0 : write(*,1) 'next_xq', next_xq
908 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
909 0 : write(*,*) 'pick_new_points: next_xq < xq_new(k_new)'
910 0 : ierr = -1
911 0 : return
912 : end if
913 :
914 5656552 : dq_sum = sum(dq_new(1:k_new))
915 :
916 10485 : k_new = k_new + 1
917 10485 : if (max_allowed_nz > 0 .and. k_new > max_allowed_nz) then
918 0 : write(*,*) 'tried to increase number of mesh points beyond max allowed nz', max_allowed_nz
919 0 : ierr = -1
920 0 : return
921 : end if
922 :
923 10485 : xq_new(k_new) = next_xq
924 10485 : if (abs(xq_new(k_new) - dq_sum) > 1d-6) then
925 0 : write(*,'(A)')
926 0 : write(*,*) 'k_new', k_new
927 0 : write(*,1) 'xq_new(k_new) - dq_sum', xq_new(k_new) - dq_sum
928 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
929 0 : write(*,1) 'dq_sum', dq_sum
930 0 : write(*,*) 'pick_new_points: abs(xq_new(k_new) - dq_sum) > 1d-6'
931 0 : ierr = -1
932 0 : return
933 : end if
934 : !write(*,2) 'xq_new(k_new) - dq_sum', k_new, xq_new(k_new) - dq_sum
935 :
936 : ! increment k_old if necessary
937 10485 : do while (k_old < nz_old)
938 10475 : if (xq_old(k_old+1) > next_xq) exit
939 10485 : k_old = k_old + 1
940 : end do
941 :
942 : end do
943 :
944 10 : nz_new = k_new
945 :
946 : if (plan_dbg) write(*,2) 'after pick_new_points: nz_new', nz_new
947 :
948 : end subroutine pick_new_points
949 :
950 :
951 20926 : logical function is_convective_boundary(kk)
952 : integer, intent(in) :: kk
953 20926 : is_convective_boundary = .false.
954 20926 : if (kk == nz) return
955 : is_convective_boundary = &
956 : (s% mixing_type(kk) == convective_mixing .and. &
957 : s% mixing_type(kk+1) /= convective_mixing) .or. &
958 : (s% mixing_type(kk+1) == convective_mixing .and. &
959 20916 : s% mixing_type(kk) /= convective_mixing)
960 : end function is_convective_boundary
961 :
962 20868 : logical function is_crystal_boundary(kk)
963 : integer, intent(in) :: kk
964 : if(s% do_phase_separation .and. & ! only need this protection when phase separation is on
965 20868 : s% m(kk) <= s% crystal_core_boundary_mass .and. &
966 : s% m(kk-1) >= s% crystal_core_boundary_mass) then
967 : is_crystal_boundary = .true.
968 : else
969 20868 : is_crystal_boundary = .false.
970 : end if
971 20868 : end function is_crystal_boundary
972 :
973 : subroutine open_debug_file
974 : include 'formats'
975 : open(newunit=iounit, file=trim('plan_debug.data'), action='write', iostat=ierr)
976 : if (ierr /= 0) then
977 : write(*, *) 'open plan_debug.data failed'
978 : call mesa_error(__FILE__,__LINE__,'debug do_mesh_plan')
979 : end if
980 : write(*,*) 'write plan_debug.data'
981 : end subroutine open_debug_file
982 :
983 :
984 : end subroutine do_mesh_plan
985 :
986 :
987 0 : subroutine realloc(s, old_size, new_capacity, xq_new, dq_new, which_gval, comes_from, ierr)
988 : use alloc
989 : type (star_info), pointer :: s
990 : integer, intent(in) :: old_size, new_capacity
991 : real(dp), pointer :: xq_new(:), dq_new(:)
992 : integer, pointer :: which_gval(:), comes_from(:)
993 : integer, intent(out) :: ierr
994 : integer, parameter :: extra = 100
995 0 : call realloc_integer_work_array(s, which_gval, old_size, new_capacity, extra, ierr)
996 0 : if (ierr /= 0) return
997 0 : call realloc_integer_work_array(s, comes_from, old_size, new_capacity, extra, ierr)
998 0 : if (ierr /= 0) return
999 0 : call realloc_work_array(s, .false., xq_new, old_size, new_capacity, extra, 'mesh_plan', ierr)
1000 0 : if (ierr /= 0) return
1001 0 : call realloc_work_array(s, .false., dq_new, old_size, new_capacity, extra, 'mesh_plan', ierr)
1002 0 : if (ierr /= 0) return
1003 0 : end subroutine realloc
1004 :
1005 :
1006 10495 : real(dp) function pick_next_dq(s, &
1007 : dbg, next_dq_max, k_old, k_new, nz_old, num_gvals, xq_new, dq_new, &
1008 : xq_old, dq_old, min_dq, min_dq_for_split, min_dq_for_xa, min_dq_for_logT, &
1009 : max_surface_cell_dq, max_dq_cntr, max_num_subcells, &
1010 : gval_is_xa_function, gval_is_logT_function, gvals, delta_gval_max, &
1011 10495 : gval_names, which_gval, ierr)
1012 0 : use mesh_functions, only: max_allowed_gvals
1013 : type (star_info), pointer :: s
1014 : logical, intent(in) :: dbg
1015 : integer, intent(in) :: k_old, k_new, nz_old, num_gvals, max_num_subcells
1016 : real(dp), pointer :: xq_new(:), dq_new(:) ! (nz)
1017 : real(dp), pointer :: xq_old(:) ! (nz_old)
1018 : real(dp), pointer :: dq_old(:) ! (nz_old)
1019 : logical, dimension(max_allowed_gvals) :: gval_is_xa_function, gval_is_logT_function
1020 : real(dp), pointer :: gvals(:,:) ! (nz_old, num_gvals)
1021 : real(dp), intent(in) :: &
1022 : next_dq_max, min_dq, min_dq_for_split, &
1023 : min_dq_for_xa, min_dq_for_logT, max_surface_cell_dq, max_dq_cntr
1024 : real(dp), pointer :: delta_gval_max(:) ! (nz_old, num_gvals)
1025 : character (len=32) :: gval_names(:) ! (num_gvals) for debugging.
1026 : integer, pointer :: which_gval(:) ! (nz_new) for debugging.
1027 : integer, intent(out) :: ierr
1028 :
1029 41980 : real(dp) :: nxt_dqs(num_gvals), default
1030 : integer :: i, j, jmin, op_err
1031 : logical :: pkdbg
1032 :
1033 : include 'formats'
1034 :
1035 10495 : pkdbg = dbg !.or. k_old == 4000
1036 10495 : ierr = 0
1037 :
1038 10495 : if (pkdbg) write(*,*)
1039 0 : if (pkdbg) write(*,2) 'dq_old(k_old)', k_old, dq_old(k_old)
1040 :
1041 10495 : if (k_new == 1) then
1042 10 : pick_next_dq = min(1d0, sqrt(min_dq))
1043 10485 : else if (k_new <= 20) then
1044 190 : pick_next_dq = min(1-xq_new(k_new), 10*dq_new(k_new-1))
1045 : else
1046 10295 : pick_next_dq = 1-xq_new(k_new)
1047 : end if
1048 :
1049 41980 : nxt_dqs(:) = pick_next_dq
1050 10495 : if (k_old == nz_old) then
1051 10 : which_gval(k_new) = 0
1052 10 : pick_next_dq = max(min_dq, (1-xq_new(k_new)) - max_dq_cntr)
1053 10 : do i=1,10
1054 10 : if (1-xq_new(k_new) <= max_dq_cntr*i) then
1055 10 : pick_next_dq = (1-xq_new(k_new))/i
1056 10 : exit
1057 : end if
1058 : end do
1059 10 : return
1060 : end if
1061 :
1062 10485 : default = pick_next_dq ! default size. can be reduced according to gradients of gvals
1063 41940 : do j=1,num_gvals
1064 : nxt_dqs(j) = pick1_dq(s, &
1065 : j, next_dq_max, default, .false., k_old, k_new, nz_old, &
1066 : xq_new, xq_old, dq_old, min_dq, min_dq_for_xa, min_dq_for_logT, max_num_subcells, &
1067 : gval_is_xa_function(j), gval_is_logT_function(j), &
1068 31455 : gvals, delta_gval_max, gval_names, op_err)
1069 41940 : if (op_err /= 0) ierr = op_err
1070 : end do
1071 :
1072 52425 : jmin = minloc(nxt_dqs(:),dim=1)
1073 10485 : if (pkdbg) write(*,3) 'jmin, k_new, init pick_next_dq', jmin, k_new, pick_next_dq
1074 10485 : which_gval(k_new) = jmin
1075 10485 : pick_next_dq = max(min_dq, min(pick_next_dq, nxt_dqs(jmin)))
1076 :
1077 10495 : end function pick_next_dq
1078 :
1079 :
1080 31455 : real(dp) function pick1_dq(s, &
1081 : j, next_dq_max, default, dbg, k_old, k_new, nz_old, &
1082 : xq_new, xq_old, dq_old, min_dq, min_dq_for_xa, min_dq_for_logT, max_num_subcells, &
1083 31455 : is_xa_function, is_logT_function, gvals, delta_gval_max, gval_names, ierr)
1084 10495 : use num_lib, only: binary_search
1085 : type (star_info), pointer :: s
1086 : integer, intent(in) :: j
1087 : real(dp), intent(in) :: next_dq_max, default, min_dq, min_dq_for_xa, min_dq_for_logT
1088 : logical, intent(in) :: dbg
1089 : integer, intent(in) :: k_old, k_new, nz_old, max_num_subcells
1090 : real(dp), pointer :: xq_new(:) ! (nz)
1091 : real(dp), pointer :: xq_old(:) ! (nz_old)
1092 : real(dp), pointer :: dq_old(:) ! (nz_old)
1093 : logical :: is_xa_function, is_logT_function
1094 : real(dp), pointer :: gvals(:,:) ! (nz_old, num_gvals)
1095 : real(dp), pointer :: delta_gval_max(:) ! (nz_old, num_gvals)
1096 : character (len=32) :: gval_names(:) ! (num_gvals) for debugging.
1097 : integer, intent(out) :: ierr
1098 :
1099 31455 : real(dp) :: gnew, dgnew, gmax, gmin, xq, dq_next, dq_sum, sz, dval
1100 : integer :: k
1101 : logical :: dbg1
1102 :
1103 : include 'formats'
1104 :
1105 31455 : ierr = 0
1106 31455 : dbg1 = dbg !.or. (k_old == -1)
1107 :
1108 31455 : pick1_dq = default
1109 31455 : dq_next = -1
1110 :
1111 31455 : if (dbg1) write(*,*)
1112 0 : if (dbg1) write(*,*) 'find max allowed dq for gvals(:,j)', j
1113 : ! find max allowed dq for gvals(:,j)
1114 :
1115 : ! linear interpolate to estimate gvals(:,j) at q_new(k_new)
1116 31455 : dval = gvals(k_old,j) - gvals(k_old+1,j)
1117 :
1118 31455 : gnew = gvals(k_old+1,j) + dval*(xq_old(k_old+1) - xq_new(k_new))/dq_old(k_old)
1119 31455 : if (dbg1) write(*,*) trim(gval_names(j))
1120 :
1121 31455 : dgnew = min(delta_gval_max(k_old), delta_gval_max(k_old+1))
1122 31455 : gmax = gnew + dgnew
1123 31455 : gmin = gnew - dgnew
1124 :
1125 31455 : dq_sum = 0
1126 91915 : do k=k_old+1, nz_old
1127 :
1128 91915 : if (dbg1) write(*,2) 'gvals(k,j)', k, gvals(k,j), dq_sum, dq_old(k-1)
1129 :
1130 91915 : if (gvals(k,j) <= gmax .and. gvals(k,j) >= gmin) then
1131 75436 : if (xq_old(k-1) >= xq_new(k_new)) then
1132 75436 : dq_sum = dq_sum + dq_old(k-1)
1133 : else
1134 0 : if (dq_sum /= 0) then
1135 0 : write(*,*) '(dq_sum /= 0)'
1136 0 : write(*,*) 'pick1_dq'
1137 0 : ierr = -1
1138 0 : return
1139 : end if
1140 0 : dq_sum = xq_old(k) - xq_new(k_new)
1141 : end if
1142 75436 : if (is_bad(dq_sum)) then
1143 0 : write(*,2) 'dq_sum', k, dq_sum
1144 0 : write(*,*) 'pick1_dq'
1145 0 : ierr = -1
1146 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'mesh plan')
1147 0 : return
1148 : end if
1149 75436 : if (dq_sum >= next_dq_max) exit
1150 60690 : if (dq_sum >= default) then
1151 0 : dq_sum = default; exit
1152 : end if
1153 60690 : if (k < nz_old) cycle
1154 : ! pick location inside center zone
1155 230 : if (dbg1) write(*,*) 'pick location inside center zone'
1156 230 : if (gvals(k-1,j) < gvals(k,j)) then ! see where reach gmax in center cell
1157 230 : dq_next = find0(0d0, gvals(k-1,j)-gmax, dq_old(k-1), gvals(k,j)-gmax)
1158 230 : if (dbg1) write(*,1) 'gvals(k-1,j) < gvals(k,j)', dq_next
1159 0 : else if (gvals(k-1,j) > gvals(k,j)) then ! see where reach gmin
1160 0 : dq_next = find0(0d0, gvals(k-1,j)-gmin, dq_old(k-1), gvals(k,j)-gmin)
1161 0 : if (dbg1) then
1162 0 : write(*,1) 'gvals(k-1,j)-gmin', gvals(k-1,j)-gmin
1163 0 : write(*,1) 'gvals(k,j)-gmin', gvals(k,j)-gmin
1164 0 : write(*,1) 'dq_old(k-1)', dq_old(k-1)
1165 0 : write(*,1) 'gvals(k-1,j) > gvals(k,j)', dq_next
1166 0 : call mesa_error(__FILE__,__LINE__,'debug pick1_dq')
1167 : end if
1168 : else ! we're done -- don't need another point for this gval
1169 0 : dq_sum = default; exit ! just return the default
1170 : end if
1171 230 : if (dq_next > 1 - (xq_new(k_new) + min_dq)) then
1172 190 : dq_sum = default
1173 : else
1174 40 : dq_sum = dq_sum + dq_next
1175 : end if
1176 : exit
1177 : end if
1178 :
1179 16479 : if (gvals(k,j) > gmax) then ! estimate where = gmax
1180 16479 : dq_next = find0(0d0, gvals(k-1,j)-gmax, dq_old(k-1), gvals(k,j)-gmax)
1181 0 : else if (gvals(k,j) < gmin) then ! estimate where = gmin
1182 0 : dq_next = find0(0d0, gvals(k-1,j)-gmin, dq_old(k-1), gvals(k,j)-gmin)
1183 : end if
1184 16479 : if (is_bad(dq_next)) then
1185 0 : write(*,2) 'dq_next', k, dq_next
1186 0 : write(*,*) 'gvals(k,j) > gmax', gvals(k,j) > gmax
1187 0 : write(*,*) 'gvals(k,j) < gmin', gvals(k,j) < gmin
1188 0 : write(*,3) 'gvals(k-1,j)', k-1, j, gvals(k-1,j)
1189 0 : write(*,2) 'gmax', k, gmax
1190 0 : write(*,3) 'gvals(k,j)', k, j, gvals(k,j)
1191 0 : write(*,2) 'gmin', k, gmin
1192 0 : write(*,2) 'dq_old(k-1)', k, dq_old(k-1)
1193 0 : write(*,*) 'pick1_dq'
1194 0 : ierr = -1
1195 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'mesh plan')
1196 0 : return
1197 : end if
1198 16479 : if (xq_old(k-1) >= xq_new(k_new)) then
1199 16479 : dq_sum = dq_sum + dq_next
1200 : else
1201 0 : if (dq_sum /= 0) then
1202 0 : write(*,*) '(dq_sum /= 0)'
1203 0 : write(*,*) 'pick1_dq'
1204 0 : ierr = -1
1205 0 : return
1206 : end if
1207 0 : dq_sum = dq_next - (xq_new(k_new) - xq_old(k-1))
1208 : end if
1209 75206 : exit
1210 :
1211 : end do
1212 :
1213 31455 : if (dbg1) then
1214 0 : write(*,1) 'after loop: dq_sum', dq_sum, min_dq, default
1215 : end if
1216 :
1217 31455 : dq_sum = max(min_dq, dq_sum)
1218 31455 : xq = xq_new(k_new) + dq_sum
1219 31455 : if (dbg1) write(*,1) 'before round_off_xq', xq
1220 31455 : xq = round_off_xq(k_old, nz_old, max_num_subcells, xq, xq_old, dq_old, sz, ierr)
1221 31455 : if (ierr /= 0) return
1222 31455 : if (dbg1) then
1223 0 : write(*,1) 'after round_off_xq: xq', xq
1224 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
1225 0 : write(*,1) 'xq - xq_new(k_new)', xq - xq_new(k_new)
1226 0 : write(*,1) 'sz', sz
1227 0 : write(*,1) 'min_dq', min_dq
1228 0 : write(*,1) 'default', default
1229 : end if
1230 31455 : pick1_dq = max(min_dq, sz, xq - xq_new(k_new))
1231 :
1232 31455 : if (is_xa_function .and. pick1_dq < min_dq_for_xa) &
1233 310 : pick1_dq = min_dq_for_xa
1234 :
1235 31455 : if (is_logT_function .and. pick1_dq < min_dq_for_logT) &
1236 0 : pick1_dq = min_dq_for_logT
1237 :
1238 31455 : if (dbg1) then
1239 0 : write(*,2) 'dq_sum', k_new, dq_sum
1240 0 : write(*,2) 'xq', k_new, xq
1241 0 : write(*,2) 'pick1_dq', k_new, pick1_dq
1242 0 : write(*,2) 'log pick1_dq', k_new, log10(pick1_dq)
1243 : end if
1244 :
1245 31455 : end function pick1_dq
1246 :
1247 :
1248 31455 : real(dp) function round_off_xq(k_old, nz_old, n, xq, xq_old, dq_old, sz, ierr)
1249 : ! adjust to match one of the candidate subcell locations
1250 : ! this prevents generating too many candidate new points
1251 : integer, intent(in) :: k_old, nz_old, n ! n is number of subcells
1252 : real(dp), intent(in) :: xq, xq_old(:), dq_old(:)
1253 : real(dp), intent(out) :: sz ! subcell size at new location
1254 : integer, intent(out) :: ierr
1255 :
1256 31455 : real(dp) :: dq, tmp
1257 : integer :: i, k, j, knxt
1258 : include 'formats'
1259 31455 : ierr = 0
1260 31455 : knxt = 0
1261 31455 : round_off_xq = -1
1262 31455 : if (xq >= xq_old(nz_old)) then
1263 230 : knxt = nz_old+1
1264 : else
1265 31225 : j = -1
1266 31225 : do k=k_old,1,-1
1267 31225 : if (xq_old(k) <= xq) then
1268 31225 : j = k+1; exit
1269 : end if
1270 : end do
1271 31225 : if (j <= 1) then
1272 0 : write(*,*) 'logic error in mesh plan: round_off_xq'
1273 0 : ierr = -1
1274 0 : return
1275 : end if
1276 104641 : do k=j,nz_old
1277 104641 : if (xq_old(k) > xq) then
1278 : knxt = k; exit
1279 : end if
1280 : end do
1281 : end if
1282 : ! xq is in old cell knxt-1
1283 : ! move location to next subcell boundary
1284 31455 : dq = xq - xq_old(knxt-1)
1285 31455 : sz = dq_old(knxt-1)/dble(n) ! size of subcells
1286 31455 : tmp = dq/sz
1287 31455 : if(tmp>huge(n)) tmp=huge(n)
1288 31455 : i = max(1,floor(tmp))
1289 31455 : if (knxt > nz_old) i = min(i,n/2) ! limit extrapolation into center
1290 31455 : dq = i*sz
1291 31455 : round_off_xq = xq_old(knxt-1) + dq
1292 31455 : if (dq <= 0) then
1293 0 : write(*,2) 'dq', knxt-1, dq, xq, xq_old(knxt-1), dq_old(knxt-1)
1294 0 : write(*,3) 'i n sz', i, n, sz
1295 0 : write(*,*) 'round_off_xq'
1296 0 : ierr = -1
1297 0 : return
1298 : end if
1299 31455 : end function round_off_xq
1300 :
1301 : end module mesh_plan
|
