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
606 10 : real(dp) :: maxval_delta_xa, next_dq_max, beta_limit, &
607 10 : remaining_dq_old, min_dr
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 :
620 : do ! pick next point location
621 :
622 10495 : dbg = plan_dbg .or. (k_new == kdbg) !.or. (s% mesh_call_number == 2005)
623 :
624 : ! when reach this point,
625 : ! have set xq_new(k) for k = 1 to k_new, and dq_new(k) for k = 1 to k_new-1.
626 : ! and have finished using old points from 1 to k_old
627 : ! i.e., xq_old(k_old+1) > xq_new(k_new) >= xq_old(k_old)
628 :
629 10495 : k_old_init = k_old
630 :
631 10495 : if (k_new == 1) then
632 10 : next_dq_max = max_surface_cell_dq
633 : else
634 10485 : next_dq_max = dq_new(k_new-1)*mesh_max_allowed_ratio
635 : end if
636 :
637 : ! make initial choice based on gradients. may reduce later.
638 10495 : if (dbg) then
639 0 : write(*,'(A)')
640 0 : write(*,3) 'call pick_next_dq', k_old, k_new, next_dq_max
641 : end if
642 :
643 10495 : if (s% gradr(k_old) > s% grada(k_old) .and. &
644 : s% min_dq_for_xa_convective > 0d0) then
645 1325 : min_dq_for_xa = s% min_dq_for_xa_convective
646 : else
647 9170 : min_dq_for_xa = s% min_dq_for_xa
648 : end if
649 :
650 10495 : min_dq_for_logT = s% min_dq_for_logT
651 :
652 : next_dq = pick_next_dq(s, &
653 : dbg, next_dq_max, k_old, k_new, nz_old, num_gvals, &
654 : xq_new, dq_new, xq_old, dq_old, min_dq, min_dq_for_split, &
655 : min_dq_for_xa, min_dq_for_logT, &
656 : max_surface_cell_dq, max_dq_cntr, max_num_subcells, &
657 : gval_is_xa_function, gval_is_logT_function, gvals, &
658 10495 : delta_gval_max, gval_names, which_gval, ierr)
659 10495 : if (ierr /= 0) return
660 :
661 10495 : if (dbg) &
662 0 : write(*,3) 'pick_next_dq next_dq prev_dq', k_old, k_new, next_dq, dq_old(k_old)
663 :
664 10495 : if (k_new == 1) then
665 10 : max_merge = max_num_merge_surface_cells
666 : else
667 10485 : max_merge = max_num_merge_cells
668 : end if
669 :
670 10495 : if (k_old < nz_old) then
671 10485 : remaining_dq_old = xq_old(k_old+1)-xq_new(k_new)
672 : else
673 10 : remaining_dq_old = 1d0-xq_new(k_new)
674 : end if
675 :
676 10495 : if (next_dq < remaining_dq_old) then
677 0 : if (.not. okay_to_split1(k_old, next_dq, remaining_dq_old)) then
678 0 : next_dq = remaining_dq_old ! dq_old(k_old)
679 : end if
680 : end if
681 :
682 10495 : if (next_dq > max_dq) then
683 21 : if (xq_new(k_new) == xq_old(k_old)) then
684 2623 : do i = nz_old, k_old+1, -1
685 2623 : if (xq_old(i) - xq_new(k_new) <= max_dq) then
686 21 : next_dq = xq_old(i) - xq_new(k_new)
687 21 : exit
688 : end if
689 : end do
690 : end if
691 : end if
692 :
693 10495 : if (next_dq > max_dq) then
694 0 : next_dq = xq_old(k_old) + dq_old(k_old) - xq_new(k_new)
695 : end if
696 :
697 10495 : if (k_new == 1 .and. next_dq > max_surface_cell_dq) then
698 0 : if (k_old_init == -1 .and. .true.) write(*,*) 'next_dq > max_surface_cell_dq'
699 0 : next_dq = max_surface_cell_dq
700 : end if
701 :
702 : ! Enforce minimum surface‐cell dq
703 10495 : if (k_old == 1 .and. next_dq < min_surface_cell_dq) then
704 : !write(*,*) 'next_dq < min_surface_cell_dq'
705 0 : next_dq = min_surface_cell_dq
706 : end if
707 :
708 10495 : next_xq = xq_new(k_new) + next_dq
709 10495 : if (next_xq > 1d0 - min_dq) then
710 10 : next_xq = (1d0 + xq_new(k_new))/2d0
711 10 : if (k_old < nz_old) then ! make sure don't split current k_old for this case
712 0 : if (xq_old(k_old+1) > next_xq) next_xq = xq_old(k_old+1)
713 : end if
714 10 : next_dq = next_xq - xq_new(k_new)
715 : end if
716 :
717 10495 : if (k_old < nz_old) then
718 :
719 10485 : if (xq_new(k_new) == xq_old(k_old) .and. &
720 : next_dq > dq_old(k_old) - min_dq/2d0) then
721 :
722 10485 : if (.not. okay_to_merge) then
723 :
724 0 : k_old_next = k_old + 1
725 :
726 10485 : else if (k_old < min_k_old_for_split .or. &
727 : k_old > max_k_old_for_split) then
728 :
729 0 : k_old_next = k_old + 1
730 :
731 : else ! consider doing merge
732 :
733 10485 : if (next_dq > 1.5d0*dq_old(k_old)) then
734 7240 : next_dq = 0.9d0*next_dq ! to avoid split-merge flip-flops
735 7240 : next_xq = xq_new(k_new) + next_dq
736 : end if
737 :
738 10485 : k_old_next_max = min(nz_old, k_old + max_merge)
739 10485 : k_old_next = k_old_next_max ! will cut this back as necessary
740 22526 : do kk=k_old+1,k_old_next_max
741 209260 : maxval_delta_xa = maxval(abs(s% xa(:,kk)-s% xa(:,kk-1)))
742 188334 : j00 = maxloc(s% xa(:,kk),dim=1)
743 188334 : jm1 = maxloc(s% xa(:,kk-1),dim=1)
744 : if (maxval_delta_xa > s% max_delta_x_for_merge .or. &
745 20926 : j00 /= jm1 .or. is_convective_boundary(kk) .or. &
746 1600 : is_crystal_boundary(kk)) then
747 : ! don't merge across convective or crystal boundary
748 58 : k_old_next = kk-1
749 58 : exit
750 20868 : else if (next_xq <= xq_old(kk) + min_dq/2d0) then
751 8827 : k_old_next = max(k_old+1,kk-1)
752 8827 : exit
753 : end if
754 : end do
755 :
756 : end if
757 :
758 10485 : k_old_next = max(k_old_next, k_old+1)
759 10485 : next_xq = xq_old(k_old_next)
760 10485 : next_dq = next_xq - xq_new(k_new)
761 :
762 0 : else if (next_xq >= xq_old(k_old+1) - min_dq/2d0) then
763 : ! this is final subcell of a split, so adjust to finish the parent cell
764 :
765 0 : k_old_next = k_old+1
766 0 : next_xq = xq_old(k_old_next)
767 0 : next_dq = next_xq - xq_new(k_new)
768 :
769 : else ! non-final subcell of split
770 :
771 0 : k_old_next = k_old
772 :
773 : end if
774 :
775 10485 : if (next_xq == xq_old(k_old_next)) then
776 : ! finishing 1 or more old cells and not already at max merge.
777 : ! consider forcing a merge with the next cell to make this one larger.
778 10485 : force_merge_with_one_more = .false.
779 :
780 10485 : if (dq_old(k_old) < min_dq) then
781 : force_merge_with_one_more = .true.
782 10485 : else if (s% merge_if_dlnR_too_small) then
783 0 : if (xq_new(k_new) <= xq_old(k_old) .and. &
784 : s% lnR(k_old) - s% lnR(k_old_next) < s% mesh_min_dlnR) then
785 : force_merge_with_one_more = .true.
786 0 : else if (k_old_next == nz_old .and. s% R_center > 0) then
787 0 : force_merge_with_one_more = s% lnR(k_old_next) - log(s% R_center) < s% mesh_min_dlnR
788 : end if
789 : end if
790 :
791 10485 : if ((.not. force_merge_with_one_more) .and. s% merge_if_dr_div_cs_too_small) then
792 10485 : if (xq_new(k_new) <= xq_old(k_old) .and. &
793 : s% r(k_old) - s% r(k_old_next) < s% mesh_min_dr_div_cs*s% csound(k_old)) then
794 : force_merge_with_one_more = .true.
795 10485 : else if (k_old_next == nz_old) then
796 : force_merge_with_one_more = (s% r(k_old_next) - s% R_center) < &
797 10 : s% mesh_min_dr_div_cs*s% csound(k_old_next)
798 10 : if (force_merge_with_one_more .and. dbg) &
799 0 : write(*,3) 'do merge for k_old_next == nz_old', k_old, k_old_next, &
800 0 : s% r(k_old) - s% R_center, &
801 0 : s% mesh_min_dr_div_cs*s% csound(k_old_next), &
802 0 : s% mesh_min_dr_div_cs, s% csound(k_old_next)
803 : end if
804 : end if
805 :
806 10485 : if ((.not. force_merge_with_one_more) .and. &
807 : s% merge_if_dr_div_dRstar_too_small) then
808 10485 : if (xq_new(k_new) <= xq_old(k_old) .and. &
809 : s% r(k_old) - s% r(k_old_next) < min_dr) then
810 : force_merge_with_one_more = .true.
811 10485 : else if (k_old_next == nz_old) then
812 10 : force_merge_with_one_more = (s% r(k_old_next) - s% R_center) < min_dr
813 10 : if (force_merge_with_one_more .and. dbg) &
814 0 : write(*,3) 'do merge for k_old_next == nz_old', k_old, k_old_next, &
815 0 : s% r(k_old) - s% R_center, min_dr
816 : end if
817 : end if
818 :
819 10485 : if (force_merge_with_one_more) then
820 0 : k_old_next = k_old_next + 1
821 0 : if (k_old_next < nz_old) then
822 0 : next_xq = xq_old(k_old_next)
823 : else
824 0 : next_xq = 1d0
825 : end if
826 0 : next_dq = next_xq - xq_new(k_new)
827 0 : if (dbg) then
828 0 : write(*,3) 'force merge', k_old, k_new, next_xq, next_dq
829 0 : write(*,'(A)')
830 : end if
831 : end if
832 :
833 : end if
834 :
835 10485 : k_old = k_old_next
836 :
837 : end if
838 :
839 10495 : comes_from(k_new) = k_old_init
840 :
841 : ! check if we're done
842 10495 : if (1 - xq_new(k_new) < max_dq_cntr .or. 1 - next_xq < min_dq) then
843 10 : dq_new(k_new) = 1 - xq_new(k_new)
844 : exit
845 : end if
846 :
847 10485 : dq_new(k_new) = next_dq
848 :
849 10485 : if (k_new == new_capacity) then ! increase allocated size
850 0 : new_capacity = (new_capacity*5)/4 + 10
851 0 : call realloc(s, k_new, new_capacity, xq_new, dq_new, which_gval, comes_from, ierr)
852 0 : if (ierr /= 0) return
853 : end if
854 :
855 10485 : if (next_xq < xq_new(k_new)) then
856 0 : write(*,*) 'nz_old', nz_old
857 0 : write(*,*) 'k_new', k_new
858 0 : write(*,1) 'next_xq', next_xq
859 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
860 0 : write(*,*) 'pick_new_points: next_xq < xq_new(k_new)'
861 0 : ierr = -1
862 0 : return
863 : end if
864 :
865 5656552 : dq_sum = sum(dq_new(1:k_new))
866 :
867 10485 : k_new = k_new + 1
868 10485 : if (max_allowed_nz > 0 .and. k_new > max_allowed_nz) then
869 0 : write(*,*) 'tried to increase number of mesh points beyond max allowed nz', max_allowed_nz
870 0 : ierr = -1
871 0 : return
872 : end if
873 :
874 10485 : xq_new(k_new) = next_xq
875 10485 : if (abs(xq_new(k_new) - dq_sum) > 1d-6) then
876 0 : write(*,'(A)')
877 0 : write(*,*) 'k_new', k_new
878 0 : write(*,1) 'xq_new(k_new) - dq_sum', xq_new(k_new) - dq_sum
879 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
880 0 : write(*,1) 'dq_sum', dq_sum
881 0 : write(*,*) 'pick_new_points: abs(xq_new(k_new) - dq_sum) > 1d-6'
882 0 : ierr = -1
883 0 : return
884 : end if
885 : !write(*,2) 'xq_new(k_new) - dq_sum', k_new, xq_new(k_new) - dq_sum
886 :
887 : ! increment k_old if necessary
888 10485 : do while (k_old < nz_old)
889 10475 : if (xq_old(k_old+1) > next_xq) exit
890 10485 : k_old = k_old + 1
891 : end do
892 :
893 : end do
894 :
895 10 : nz_new = k_new
896 :
897 : if (plan_dbg) write(*,2) 'after pick_new_points: nz_new', nz_new
898 :
899 : end subroutine pick_new_points
900 :
901 :
902 20926 : logical function is_convective_boundary(kk)
903 : integer, intent(in) :: kk
904 20926 : is_convective_boundary = .false.
905 20926 : if (kk == nz) return
906 : is_convective_boundary = &
907 : (s% mixing_type(kk) == convective_mixing .and. &
908 : s% mixing_type(kk+1) /= convective_mixing) .or. &
909 : (s% mixing_type(kk+1) == convective_mixing .and. &
910 20916 : s% mixing_type(kk) /= convective_mixing)
911 : end function is_convective_boundary
912 :
913 20868 : logical function is_crystal_boundary(kk)
914 : integer, intent(in) :: kk
915 : if(s% do_phase_separation .and. & ! only need this protection when phase separation is on
916 20868 : s% m(kk) <= s% crystal_core_boundary_mass .and. &
917 : s% m(kk-1) >= s% crystal_core_boundary_mass) then
918 : is_crystal_boundary = .true.
919 : else
920 20868 : is_crystal_boundary = .false.
921 : end if
922 20868 : end function is_crystal_boundary
923 :
924 : subroutine open_debug_file
925 : include 'formats'
926 : open(newunit=iounit, file=trim('plan_debug.data'), action='write', iostat=ierr)
927 : if (ierr /= 0) then
928 : write(*, *) 'open plan_debug.data failed'
929 : call mesa_error(__FILE__,__LINE__,'debug do_mesh_plan')
930 : end if
931 : write(*,*) 'write plan_debug.data'
932 : end subroutine open_debug_file
933 :
934 :
935 : end subroutine do_mesh_plan
936 :
937 :
938 0 : subroutine realloc(s, old_size, new_capacity, xq_new, dq_new, which_gval, comes_from, ierr)
939 : use alloc
940 : type (star_info), pointer :: s
941 : integer, intent(in) :: old_size, new_capacity
942 : real(dp), pointer :: xq_new(:), dq_new(:)
943 : integer, pointer :: which_gval(:), comes_from(:)
944 : integer, intent(out) :: ierr
945 : integer, parameter :: extra = 100
946 0 : call realloc_integer_work_array(s, which_gval, old_size, new_capacity, extra, ierr)
947 0 : if (ierr /= 0) return
948 0 : call realloc_integer_work_array(s, comes_from, old_size, new_capacity, extra, ierr)
949 0 : if (ierr /= 0) return
950 0 : call realloc_work_array(s, .false., xq_new, old_size, new_capacity, extra, 'mesh_plan', ierr)
951 0 : if (ierr /= 0) return
952 0 : call realloc_work_array(s, .false., dq_new, old_size, new_capacity, extra, 'mesh_plan', ierr)
953 0 : if (ierr /= 0) return
954 0 : end subroutine realloc
955 :
956 :
957 10495 : real(dp) function pick_next_dq(s, &
958 : dbg, next_dq_max, k_old, k_new, nz_old, num_gvals, xq_new, dq_new, &
959 : xq_old, dq_old, min_dq, min_dq_for_split, min_dq_for_xa, min_dq_for_logT, &
960 : max_surface_cell_dq, max_dq_cntr, max_num_subcells, &
961 : gval_is_xa_function, gval_is_logT_function, gvals, delta_gval_max, &
962 10495 : gval_names, which_gval, ierr)
963 0 : use mesh_functions, only: max_allowed_gvals
964 : type (star_info), pointer :: s
965 : logical, intent(in) :: dbg
966 : integer, intent(in) :: k_old, k_new, nz_old, num_gvals, max_num_subcells
967 : real(dp), pointer :: xq_new(:), dq_new(:) ! (nz)
968 : real(dp), pointer :: xq_old(:) ! (nz_old)
969 : real(dp), pointer :: dq_old(:) ! (nz_old)
970 : logical, dimension(max_allowed_gvals) :: gval_is_xa_function, gval_is_logT_function
971 : real(dp), pointer :: gvals(:,:) ! (nz_old, num_gvals)
972 : real(dp), intent(in) :: &
973 : next_dq_max, min_dq, min_dq_for_split, &
974 : min_dq_for_xa, min_dq_for_logT, max_surface_cell_dq, max_dq_cntr
975 : real(dp), pointer :: delta_gval_max(:) ! (nz_old, num_gvals)
976 : character (len=32) :: gval_names(:) ! (num_gvals) for debugging.
977 : integer, pointer :: which_gval(:) ! (nz_new) for debugging.
978 : integer, intent(out) :: ierr
979 :
980 41980 : real(dp) :: nxt_dqs(num_gvals), default
981 : integer :: i, j, jmin, op_err
982 : logical :: pkdbg
983 :
984 : include 'formats'
985 :
986 10495 : pkdbg = dbg !.or. k_old == 4000
987 10495 : ierr = 0
988 :
989 10495 : if (pkdbg) write(*,*)
990 0 : if (pkdbg) write(*,2) 'dq_old(k_old)', k_old, dq_old(k_old)
991 :
992 10495 : if (k_new == 1) then
993 10 : pick_next_dq = min(1d0, sqrt(min_dq))
994 10485 : else if (k_new <= 20) then
995 190 : pick_next_dq = min(1-xq_new(k_new), 10*dq_new(k_new-1))
996 : else
997 10295 : pick_next_dq = 1-xq_new(k_new)
998 : end if
999 :
1000 41980 : nxt_dqs(:) = pick_next_dq
1001 10495 : if (k_old == nz_old) then
1002 10 : which_gval(k_new) = 0
1003 10 : pick_next_dq = max(min_dq, (1-xq_new(k_new)) - max_dq_cntr)
1004 10 : do i=1,10
1005 10 : if (1-xq_new(k_new) <= max_dq_cntr*i) then
1006 10 : pick_next_dq = (1-xq_new(k_new))/i
1007 10 : exit
1008 : end if
1009 : end do
1010 10 : return
1011 : end if
1012 :
1013 10485 : default = pick_next_dq ! default size. can be reduced according to gradients of gvals
1014 41940 : do j=1,num_gvals
1015 : nxt_dqs(j) = pick1_dq(s, &
1016 : j, next_dq_max, default, .false., k_old, k_new, nz_old, &
1017 : xq_new, xq_old, dq_old, min_dq, min_dq_for_xa, min_dq_for_logT, max_num_subcells, &
1018 : gval_is_xa_function(j), gval_is_logT_function(j), &
1019 31455 : gvals, delta_gval_max, gval_names, op_err)
1020 41940 : if (op_err /= 0) ierr = op_err
1021 : end do
1022 :
1023 52425 : jmin = minloc(nxt_dqs(:),dim=1)
1024 10485 : if (pkdbg) write(*,3) 'jmin, k_new, init pick_next_dq', jmin, k_new, pick_next_dq
1025 10485 : which_gval(k_new) = jmin
1026 10485 : pick_next_dq = max(min_dq, min(pick_next_dq, nxt_dqs(jmin)))
1027 :
1028 10495 : end function pick_next_dq
1029 :
1030 :
1031 31455 : real(dp) function pick1_dq(s, &
1032 : j, next_dq_max, default, dbg, k_old, k_new, nz_old, &
1033 : xq_new, xq_old, dq_old, min_dq, min_dq_for_xa, min_dq_for_logT, max_num_subcells, &
1034 31455 : is_xa_function, is_logT_function, gvals, delta_gval_max, gval_names, ierr)
1035 10495 : use num_lib, only: binary_search
1036 : type (star_info), pointer :: s
1037 : integer, intent(in) :: j
1038 : real(dp), intent(in) :: next_dq_max, default, min_dq, min_dq_for_xa, min_dq_for_logT
1039 : logical, intent(in) :: dbg
1040 : integer, intent(in) :: k_old, k_new, nz_old, max_num_subcells
1041 : real(dp), pointer :: xq_new(:) ! (nz)
1042 : real(dp), pointer :: xq_old(:) ! (nz_old)
1043 : real(dp), pointer :: dq_old(:) ! (nz_old)
1044 : logical :: is_xa_function, is_logT_function
1045 : real(dp), pointer :: gvals(:,:) ! (nz_old, num_gvals)
1046 : real(dp), pointer :: delta_gval_max(:) ! (nz_old, num_gvals)
1047 : character (len=32) :: gval_names(:) ! (num_gvals) for debugging.
1048 : integer, intent(out) :: ierr
1049 :
1050 31455 : real(dp) :: gnew, dgnew, gmax, gmin, xq, dq_next, dq_sum, sz, dval
1051 : integer :: k
1052 : logical :: dbg1
1053 :
1054 : include 'formats'
1055 :
1056 31455 : ierr = 0
1057 31455 : dbg1 = dbg !.or. (k_old == -1)
1058 :
1059 31455 : pick1_dq = default
1060 31455 : dq_next = -1
1061 :
1062 31455 : if (dbg1) write(*,*)
1063 0 : if (dbg1) write(*,*) 'find max allowed dq for gvals(:,j)', j
1064 : ! find max allowed dq for gvals(:,j)
1065 :
1066 : ! linear interpolate to estimate gvals(:,j) at q_new(k_new)
1067 31455 : dval = gvals(k_old,j) - gvals(k_old+1,j)
1068 :
1069 31455 : gnew = gvals(k_old+1,j) + dval*(xq_old(k_old+1) - xq_new(k_new))/dq_old(k_old)
1070 31455 : if (dbg1) write(*,*) trim(gval_names(j))
1071 :
1072 31455 : dgnew = min(delta_gval_max(k_old), delta_gval_max(k_old+1))
1073 31455 : gmax = gnew + dgnew
1074 31455 : gmin = gnew - dgnew
1075 :
1076 31455 : dq_sum = 0
1077 91915 : do k=k_old+1, nz_old
1078 :
1079 91915 : if (dbg1) write(*,2) 'gvals(k,j)', k, gvals(k,j), dq_sum, dq_old(k-1)
1080 :
1081 91915 : if (gvals(k,j) <= gmax .and. gvals(k,j) >= gmin) then
1082 75436 : if (xq_old(k-1) >= xq_new(k_new)) then
1083 75436 : dq_sum = dq_sum + dq_old(k-1)
1084 : else
1085 0 : if (dq_sum /= 0) then
1086 0 : write(*,*) '(dq_sum /= 0)'
1087 0 : write(*,*) 'pick1_dq'
1088 0 : ierr = -1
1089 0 : return
1090 : end if
1091 0 : dq_sum = xq_old(k) - xq_new(k_new)
1092 : end if
1093 75436 : if (is_bad(dq_sum)) then
1094 0 : write(*,2) 'dq_sum', k, dq_sum
1095 0 : write(*,*) 'pick1_dq'
1096 0 : ierr = -1
1097 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'mesh plan')
1098 0 : return
1099 : end if
1100 75436 : if (dq_sum >= next_dq_max) exit
1101 60690 : if (dq_sum >= default) then
1102 0 : dq_sum = default; exit
1103 : end if
1104 60690 : if (k < nz_old) cycle
1105 : ! pick location inside center zone
1106 230 : if (dbg1) write(*,*) 'pick location inside center zone'
1107 230 : if (gvals(k-1,j) < gvals(k,j)) then ! see where reach gmax in center cell
1108 230 : dq_next = find0(0d0, gvals(k-1,j)-gmax, dq_old(k-1), gvals(k,j)-gmax)
1109 230 : if (dbg1) write(*,1) 'gvals(k-1,j) < gvals(k,j)', dq_next
1110 0 : else if (gvals(k-1,j) > gvals(k,j)) then ! see where reach gmin
1111 0 : dq_next = find0(0d0, gvals(k-1,j)-gmin, dq_old(k-1), gvals(k,j)-gmin)
1112 0 : if (dbg1) then
1113 0 : write(*,1) 'gvals(k-1,j)-gmin', gvals(k-1,j)-gmin
1114 0 : write(*,1) 'gvals(k,j)-gmin', gvals(k,j)-gmin
1115 0 : write(*,1) 'dq_old(k-1)', dq_old(k-1)
1116 0 : write(*,1) 'gvals(k-1,j) > gvals(k,j)', dq_next
1117 0 : call mesa_error(__FILE__,__LINE__,'debug pick1_dq')
1118 : end if
1119 : else ! we're done -- don't need another point for this gval
1120 0 : dq_sum = default; exit ! just return the default
1121 : end if
1122 230 : if (dq_next > 1 - (xq_new(k_new) + min_dq)) then
1123 190 : dq_sum = default
1124 : else
1125 40 : dq_sum = dq_sum + dq_next
1126 : end if
1127 : exit
1128 : end if
1129 :
1130 16479 : if (gvals(k,j) > gmax) then ! estimate where = gmax
1131 16479 : dq_next = find0(0d0, gvals(k-1,j)-gmax, dq_old(k-1), gvals(k,j)-gmax)
1132 0 : else if (gvals(k,j) < gmin) then ! estimate where = gmin
1133 0 : dq_next = find0(0d0, gvals(k-1,j)-gmin, dq_old(k-1), gvals(k,j)-gmin)
1134 : end if
1135 16479 : if (is_bad(dq_next)) then
1136 0 : write(*,2) 'dq_next', k, dq_next
1137 0 : write(*,*) 'gvals(k,j) > gmax', gvals(k,j) > gmax
1138 0 : write(*,*) 'gvals(k,j) < gmin', gvals(k,j) < gmin
1139 0 : write(*,3) 'gvals(k-1,j)', k-1, j, gvals(k-1,j)
1140 0 : write(*,2) 'gmax', k, gmax
1141 0 : write(*,3) 'gvals(k,j)', k, j, gvals(k,j)
1142 0 : write(*,2) 'gmin', k, gmin
1143 0 : write(*,2) 'dq_old(k-1)', k, dq_old(k-1)
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 16479 : if (xq_old(k-1) >= xq_new(k_new)) then
1150 16479 : dq_sum = dq_sum + dq_next
1151 : else
1152 0 : if (dq_sum /= 0) then
1153 0 : write(*,*) '(dq_sum /= 0)'
1154 0 : write(*,*) 'pick1_dq'
1155 0 : ierr = -1
1156 0 : return
1157 : end if
1158 0 : dq_sum = dq_next - (xq_new(k_new) - xq_old(k-1))
1159 : end if
1160 75206 : exit
1161 :
1162 : end do
1163 :
1164 31455 : if (dbg1) then
1165 0 : write(*,1) 'after loop: dq_sum', dq_sum, min_dq, default
1166 : end if
1167 :
1168 31455 : dq_sum = max(min_dq, dq_sum)
1169 31455 : xq = xq_new(k_new) + dq_sum
1170 31455 : if (dbg1) write(*,1) 'before round_off_xq', xq
1171 31455 : xq = round_off_xq(k_old, nz_old, max_num_subcells, xq, xq_old, dq_old, sz, ierr)
1172 31455 : if (ierr /= 0) return
1173 31455 : if (dbg1) then
1174 0 : write(*,1) 'after round_off_xq: xq', xq
1175 0 : write(*,1) 'xq_new(k_new)', xq_new(k_new)
1176 0 : write(*,1) 'xq - xq_new(k_new)', xq - xq_new(k_new)
1177 0 : write(*,1) 'sz', sz
1178 0 : write(*,1) 'min_dq', min_dq
1179 0 : write(*,1) 'default', default
1180 : end if
1181 31455 : pick1_dq = max(min_dq, sz, xq - xq_new(k_new))
1182 :
1183 31455 : if (is_xa_function .and. pick1_dq < min_dq_for_xa) &
1184 310 : pick1_dq = min_dq_for_xa
1185 :
1186 31455 : if (is_logT_function .and. pick1_dq < min_dq_for_logT) &
1187 0 : pick1_dq = min_dq_for_logT
1188 :
1189 31455 : if (dbg1) then
1190 0 : write(*,2) 'dq_sum', k_new, dq_sum
1191 0 : write(*,2) 'xq', k_new, xq
1192 0 : write(*,2) 'pick1_dq', k_new, pick1_dq
1193 0 : write(*,2) 'log pick1_dq', k_new, log10(pick1_dq)
1194 : end if
1195 :
1196 31455 : end function pick1_dq
1197 :
1198 :
1199 31455 : real(dp) function round_off_xq(k_old, nz_old, n, xq, xq_old, dq_old, sz, ierr)
1200 : ! adjust to match one of the candidate subcell locations
1201 : ! this prevents generating too many candidate new points
1202 : integer, intent(in) :: k_old, nz_old, n ! n is number of subcells
1203 : real(dp), intent(in) :: xq, xq_old(:), dq_old(:)
1204 : real(dp), intent(out) :: sz ! subcell size at new location
1205 : integer, intent(out) :: ierr
1206 :
1207 31455 : real(dp) :: dq, tmp
1208 : integer :: i, k, j, knxt
1209 : include 'formats'
1210 31455 : ierr = 0
1211 31455 : knxt = 0
1212 31455 : round_off_xq = -1
1213 31455 : if (xq >= xq_old(nz_old)) then
1214 230 : knxt = nz_old+1
1215 : else
1216 31225 : j = -1
1217 31225 : do k=k_old,1,-1
1218 31225 : if (xq_old(k) <= xq) then
1219 31225 : j = k+1; exit
1220 : end if
1221 : end do
1222 31225 : if (j <= 1) then
1223 0 : write(*,*) 'logic error in mesh plan: round_off_xq'
1224 0 : ierr = -1
1225 0 : return
1226 : end if
1227 104641 : do k=j,nz_old
1228 104641 : if (xq_old(k) > xq) then
1229 : knxt = k; exit
1230 : end if
1231 : end do
1232 : end if
1233 : ! xq is in old cell knxt-1
1234 : ! move location to next subcell boundary
1235 31455 : dq = xq - xq_old(knxt-1)
1236 31455 : sz = dq_old(knxt-1)/dble(n) ! size of subcells
1237 31455 : tmp = dq/sz
1238 31455 : if(tmp>huge(n)) tmp=huge(n)
1239 31455 : i = max(1,floor(tmp))
1240 31455 : if (knxt > nz_old) i = min(i,n/2) ! limit extrapolation into center
1241 31455 : dq = i*sz
1242 31455 : round_off_xq = xq_old(knxt-1) + dq
1243 31455 : if (dq <= 0) then
1244 0 : write(*,2) 'dq', knxt-1, dq, xq, xq_old(knxt-1), dq_old(knxt-1)
1245 0 : write(*,3) 'i n sz', i, n, sz
1246 0 : write(*,*) 'round_off_xq'
1247 0 : ierr = -1
1248 0 : return
1249 : end if
1250 31455 : end function round_off_xq
1251 :
1252 : end module mesh_plan
|
