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