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 adjust_mesh
21 :
22 : use star_private_def
23 : use const_def, only: dp, lsun
24 : use adjust_mesh_support
25 :
26 : implicit none
27 :
28 : private
29 : public :: remesh
30 :
31 : logical, parameter :: dbg_remesh = .false.
32 :
33 : contains
34 :
35 : ! makes new mesh and sets new values for xh and xa.
36 10 : integer function remesh(s)
37 : use alloc
38 : use mesh_functions, only: max_allowed_gvals
39 : use mesh_plan, only: do_mesh_plan
40 : use mesh_adjust, only: do_mesh_adjust
41 : use rates_def, only: i_rate
42 : use net_lib, only: clean_up_fractions
43 : use utils_lib
44 : use star_utils
45 : use chem_def
46 : use chem_lib, only: chem_get_iso_id
47 :
48 : type (star_info), pointer :: s
49 :
50 : integer :: alloc_level, i, k, ierr, species, nvar, nz, nz_new, nz_old, &
51 : num_gvals, j, cid, cid_max, unchanged, split, merged, revised
52 980 : type (star_info), target :: copy_info
53 : type (star_info), pointer :: c, prv
54 10 : real(dp) :: delta_coeff, sum_L_other, sum_L_other_limit, A_max, &
55 10 : mesh_max_allowed_ratio, tmp, J_tot1, J_tot2, center_logT, alfa, beta, &
56 10 : d_dlnR00, d_dlnRp1, d_dv00, d_dvp1
57 : real(dp), pointer, dimension(:) :: &
58 10 : delta_gval_max, specific_PE, specific_KE, &
59 10 : xq_old, xq_new, dq_new
60 10 : real(dp), pointer :: gvals(:,:), gvals1(:)
61 10 : integer, pointer :: which_gval(:) , comes_from(:), cell_type(:)
62 10 : logical, pointer :: do_not_split(:)
63 : character (len=32) :: gval_names(max_allowed_gvals)
64 : logical, dimension(max_allowed_gvals) :: &
65 : gval_is_xa_function, gval_is_logT_function
66 : logical, parameter :: dbg = .false.
67 :
68 : real(dp), parameter :: max_sum_abs = 10d0
69 : real(dp), parameter :: xsum_tol = 1d-2
70 : real(dp), parameter :: h_cntr_limit = 0.5d0 ! for pre-MS decision
71 :
72 : include 'formats'
73 :
74 10 : ierr = 0
75 10 : remesh = keep_going
76 10 : alloc_level = 0
77 :
78 10 : s% mesh_adjust_KE_conservation = 0
79 10 : s% mesh_adjust_IE_conservation = 0
80 10 : s% mesh_adjust_PE_conservation = 0
81 :
82 : if (dbg_remesh) write(*,*) 'enter adjust mesh'
83 :
84 10 : if (.not. s% okay_to_remesh) then
85 : if (dbg_remesh) write(*,*) 's% okay_to_remesh', s% okay_to_remesh
86 : return
87 : end if
88 :
89 10 : if (s% remesh_max_allowed_logT < 1d2) then ! check it
90 0 : if (maxval(s% lnT(1:s% nz))/ln10 > s% remesh_max_allowed_logT) then
91 : if (dbg_remesh) write(*,2) &
92 : 's% remesh_max_allowed_logT', s% model_number, s% remesh_max_allowed_logT
93 : return
94 : end if
95 : end if
96 :
97 : call nullify_work_ptrs
98 :
99 10 : if (s% dt < s% remesh_dt_limit) then
100 : if (dbg_remesh) write(*,*) 's% dt < s% remesh_dt_limit'
101 : return
102 : end if
103 :
104 10 : species = s% species
105 10 : nz_old = s% nz
106 10 : nz = nz_old
107 :
108 10 : J_tot1 = total_angular_momentum(s)
109 10 : call clean_up_fractions(1, nz, species, nz, s% xa, max_sum_abs, xsum_tol, ierr)
110 10 : if (ierr /= 0) then
111 : if (dbg_remesh) write(*,*) 'clean_up_fractions failed during adjust mesh'
112 0 : s% result_reason = adjust_mesh_failed
113 0 : s% retry_message = 'clean_up_fractions failed during adjust mesh'
114 0 : remesh = retry
115 0 : return
116 : end if
117 :
118 : call do_alloc1(ierr)
119 10 : if (ierr /= 0) then
120 : if (dbg_remesh) write(*,*) 'do_alloc1 failed for adjust mesh'
121 0 : s% result_reason = adjust_mesh_failed
122 0 : s% termination_code = t_adjust_mesh_failed
123 0 : remesh = terminate
124 0 : call dealloc
125 0 : return
126 : end if
127 :
128 12105 : do k=1,nz
129 12095 : specific_PE(k) = cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
130 12105 : specific_KE(k) = cell_specific_KE(s,k,d_dv00,d_dvp1)
131 : end do
132 :
133 10 : s% mesh_call_number = s% mesh_call_number + 1
134 :
135 10 : sum_L_other = 0
136 10 : sum_L_other_limit = Lsun
137 10 : do j=1,num_categories
138 : if (j == ipp .or. j == icno .or. j == i3alf .or. j == iphoto) cycle
139 : sum_L_other = sum_L_other + dot_product(s% dq(1:nz), s% eps_nuc_categories(j,1:nz))
140 : if (sum_L_other > sum_L_other_limit) exit
141 : end do
142 :
143 10 : center_logT = s% lnT(s% nz)/ln10
144 10 : if (center_logT <= s% logT_max_for_standard_mesh_delta_coeff) then
145 10 : delta_coeff = s% mesh_delta_coeff
146 0 : else if (center_logT >= s% logT_min_for_highT_mesh_delta_coeff) then
147 0 : delta_coeff = s% mesh_delta_coeff_for_highT
148 : else
149 : alfa = (center_logT - s% logT_max_for_standard_mesh_delta_coeff)/ &
150 0 : (s% logT_min_for_highT_mesh_delta_coeff - s% logT_max_for_standard_mesh_delta_coeff)
151 0 : beta = 1d0 - alfa
152 0 : delta_coeff = alfa*s% mesh_delta_coeff_for_highT + beta*s% mesh_delta_coeff
153 : end if
154 :
155 10 : mesh_max_allowed_ratio = s% mesh_max_allowed_ratio
156 10 : if (mesh_max_allowed_ratio < 2.5d0) then
157 0 : write(*,*) 'WARNING: increasing mesh_max_allowed_ratio to 2.5'
158 0 : s% mesh_max_allowed_ratio = 2.5d0
159 : end if
160 :
161 : ! find the heaviest species in the current net
162 10 : cid_max = 0; A_max = 0
163 :
164 90 : do i = 1, s% species
165 80 : cid = s% chem_id(i)
166 90 : if (chem_isos% W(cid) > A_max) then
167 80 : A_max = chem_isos% W(cid); cid_max = cid
168 : end if
169 : end do
170 :
171 100 : do j = 1, num_xa_function
172 90 : if (s% xa_mesh_delta_coeff(j) == 1) cycle
173 0 : if (len_trim(s% xa_function_species(j)) == 0) cycle
174 0 : cid = chem_get_iso_id(s% xa_function_species(j))
175 0 : if (cid <= 0) cycle
176 0 : if (s% net_iso(cid) == 0) cycle
177 0 : if (chem_isos% W(cid) /= A_max) cycle
178 100 : delta_coeff = delta_coeff*s% xa_mesh_delta_coeff(j)
179 : end do
180 :
181 : call check_validity(s, ierr)
182 10 : if (ierr /= 0) then
183 0 : write(*,*) 'check_validity failed at entry to adjust mesh'
184 0 : s% termination_code = t_adjust_mesh_failed
185 0 : remesh = terminate
186 0 : s% result_reason = adjust_mesh_failed
187 0 : call dealloc
188 0 : return
189 : end if
190 :
191 10 : call get_num_gvals
192 10 : if (num_gvals == 0) then
193 0 : s% termination_code = t_adjust_mesh_failed
194 0 : remesh = terminate
195 0 : write(*,*) 'must have at least 1 mesh gradient function with nonzero weight'
196 0 : s% result_reason = adjust_mesh_failed
197 0 : call dealloc
198 0 : return
199 : end if
200 :
201 : call do_alloc2(ierr)
202 10 : if (ierr /= 0) then
203 : if (dbg_remesh) write(*,*) 'do_alloc2 failed for adjust mesh'
204 0 : s% termination_code = t_adjust_mesh_failed
205 0 : remesh = terminate
206 0 : s% result_reason = adjust_mesh_failed
207 0 : call dealloc
208 0 : return
209 : end if
210 :
211 : if (dbg_remesh) write(*,*) 'call mesh_plan'
212 :
213 10 : if (s% show_mesh_changes) then
214 0 : write(*,*) 'doing mesh_call_number', s% mesh_call_number
215 0 : write(*,*) ' mesh_plan nz_old', nz_old
216 : end if
217 :
218 : call get_gval_info( &
219 : s, delta_gval_max, gvals1, nz, &
220 : num_gvals, gval_names, &
221 : gval_is_xa_function, gval_is_logT_function, ierr)
222 10 : if (ierr /= 0) then
223 0 : s% termination_code = t_adjust_mesh_failed
224 0 : remesh = terminate
225 : if (dbg_remesh) write(*,*) 'get_gval_info failed for adjust mesh'
226 0 : s% result_reason = adjust_mesh_failed
227 0 : call dealloc
228 0 : return
229 : end if
230 :
231 12095 : do k=1,nz-1
232 : do_not_split(k) = &
233 : (s% lnR(k) - s% lnR(k+1) < 2*s% mesh_min_dlnR) .or. &
234 12095 : (s% r(k) - s% r(k+1) < 2*s% mesh_min_dr_div_cs*s% csound(k))
235 : end do
236 : do_not_split(nz) = &
237 10 : (s% r(nz) - s% R_center < 2*s% mesh_min_dr_div_cs*s% csound(nz))
238 :
239 : if (dbg_remesh) write(*,*) 'call do_mesh_plan'
240 :
241 : if (dbg) write(*,1) 'sum(s% dq(1:nz))', sum(s% dq(1:nz))
242 : if (dbg) write(*,1) 's% dq(nz)', s% dq(nz)
243 : if (dbg) write(*,1) 'sum(s% dq(1:nz-1))', sum(s% dq(1:nz-1))
244 10 : tmp = s% dq(nz)
245 12095 : s% dq(nz) = 1d0 - sum(s% dq(1:nz-1))
246 10 : if (s% dq(nz) <= 0) then
247 : if (dbg) write(*,1) 'fix starting s% dq(nz) <= 0', s% dq(nz), sum(s% dq(1:nz))
248 0 : s% dq(nz) = tmp
249 0 : s% dq(1:nz) = s% dq(1:nz)/sum(s% dq(1:nz))
250 0 : s% dq(nz) = 1d0 - sum(s% dq(1:nz-1))
251 0 : if (s% dq(nz) <= 0) then
252 0 : ierr = -1
253 : if (dbg) write(*,*) 's% dq(nz) <= 0'
254 : if (dbg) call mesa_error(__FILE__,__LINE__,'debug adjust mesh')
255 0 : s% result_reason = adjust_mesh_failed
256 0 : s% termination_code = t_adjust_mesh_failed
257 0 : remesh = terminate
258 0 : call dealloc
259 0 : return
260 : end if
261 : end if
262 :
263 10 : call set_xqs(nz, xq_old, s% dq, ierr)
264 10 : if (ierr /= 0) then
265 : if (dbg) write(*,*) 'set_xqs ierr for xq_old in adjust mesh'
266 0 : s% result_reason = adjust_mesh_failed
267 0 : s% termination_code = t_adjust_mesh_failed
268 0 : remesh = terminate
269 0 : call dealloc
270 0 : return
271 : end if
272 :
273 12105 : do k=1,nz
274 12105 : delta_gval_max(k) = delta_gval_max(k)*delta_coeff
275 : end do
276 :
277 : call do_mesh_plan( &
278 : s, nz, s% max_allowed_nz, s% mesh_ok_to_merge, s% D_mix, &
279 : s% mesh_max_k_old_for_split, s% mesh_min_k_old_for_split, xq_old, s% dq, &
280 : s% min_dq, s% max_dq*s% mesh_delta_coeff, s% min_dq_for_split, mesh_max_allowed_ratio, &
281 : do_not_split, num_gvals, gval_names, &
282 : gval_is_xa_function, gval_is_logT_function, gvals, delta_gval_max, &
283 : s% max_center_cell_dq*s% mesh_delta_coeff, s% max_surface_cell_dq*s% mesh_delta_coeff, &
284 : s% max_num_subcells, s% max_num_merge_cells, &
285 10 : nz_new, xq_new, dq_new, which_gval, comes_from, ierr)
286 :
287 : if (dbg_remesh .or. dbg) write(*,*) 'back from mesh_plan'
288 :
289 10 : if (ierr /= 0) then
290 0 : write(*,'(A)')
291 0 : write(*,*) 'mesh_plan problem'
292 0 : write(*,*) 'doing mesh_call_number', s% mesh_call_number
293 0 : write(*,*) 's% model_number', s% model_number
294 0 : write(*,'(A)')
295 0 : s% termination_code = t_adjust_mesh_failed
296 0 : remesh = terminate
297 0 : s% result_reason = adjust_mesh_failed
298 0 : call dealloc
299 0 : return
300 : end if
301 :
302 10 : nz = nz_new
303 10 : s% nz = nz
304 10 : nvar = s% nvar_total
305 :
306 10 : if (.not. associated(s% other_star_info)) then
307 98 : allocate(s% other_star_info)
308 1 : prv => s% other_star_info
309 1 : c => null()
310 : else
311 9 : prv => s% other_star_info
312 9 : c => copy_info
313 9 : c = prv
314 : end if
315 :
316 10 : prv = s ! this makes copies of pointers and scalars
317 :
318 : if (dbg_remesh .or. dbg) write(*,*) 'call resize_star_info_arrays'
319 : call resize_star_info_arrays(s, c, ierr)
320 10 : if (ierr /= 0) then
321 : if (dbg) write(*,*) 'resize_star_info_arrays ierr'
322 0 : s% result_reason = adjust_mesh_failed
323 0 : s% termination_code = t_adjust_mesh_failed
324 0 : remesh = terminate
325 0 : call dealloc
326 0 : return
327 : end if
328 :
329 : if (dbg_remesh .or. dbg) write(*,*) 'call do_alloc3'
330 : call do_alloc3(ierr)
331 10 : if (ierr /= 0) then
332 : if (dbg) write(*,*) 'do_alloc3 failed for adjust mesh'
333 0 : s% result_reason = adjust_mesh_failed
334 0 : s% termination_code = t_adjust_mesh_failed
335 0 : remesh = terminate
336 0 : call dealloc
337 0 : return
338 : end if
339 :
340 : if (dbg_remesh .or. dbg) write(*,*) 'call set_types_of_new_cells'
341 : call set_types_of_new_cells(cell_type, ierr)
342 10 : if (ierr /= 0) then
343 : if (dbg) write(*,*) 'set_types_of_new_cells ierr'
344 0 : s% result_reason = adjust_mesh_failed
345 0 : s% termination_code = t_adjust_mesh_failed
346 0 : remesh = terminate
347 0 : call dealloc
348 0 : return
349 : end if
350 :
351 10633 : do k=1,nz
352 10633 : s% dq(k) = dq_new(k)
353 : end do
354 10623 : s% dq(nz) = 1d0 - sum(s% dq(1:nz-1))
355 10 : if (s% dq(nz) <= 0) then
356 0 : write(*,1) 'fix s% dq(nz) <= 0', s% dq(nz), sum(s% dq(1:nz))
357 0 : s% dq(nz) = dq_new(nz)
358 0 : s% dq(1:nz) = s% dq(1:nz)/sum(s% dq(1:nz))
359 0 : s% dq(nz) = 1d0 - sum(s% dq(1:nz-1))
360 0 : if (s% dq(nz) <= 0) then
361 0 : ierr = -1
362 : if (dbg) write(*,*) 's% dq(nz) <= 0 in adjust mesh'
363 0 : s% result_reason = adjust_mesh_failed
364 0 : s% termination_code = t_adjust_mesh_failed
365 0 : remesh = terminate
366 0 : call dealloc
367 0 : return
368 : end if
369 : end if
370 :
371 10 : call set_qs(s, nz, s% q, s% dq, ierr)
372 10 : if (ierr /= 0) then
373 : if (dbg) write(*,*) 'set_qs ierr in adjust mesh'
374 0 : s% result_reason = adjust_mesh_failed
375 0 : s% termination_code = t_adjust_mesh_failed
376 0 : remesh = terminate
377 0 : call dealloc
378 0 : return
379 : end if
380 10 : call set_xqs(nz, xq_new, s% dq, ierr)
381 10 : if (ierr /= 0) then
382 : if (dbg) write(*,*) 'set_qs ierr in adjust mesh'
383 0 : s% result_reason = adjust_mesh_failed
384 0 : s% termination_code = t_adjust_mesh_failed
385 0 : remesh = terminate
386 0 : call dealloc
387 0 : return
388 : end if
389 :
390 : ! testing -- check for q strictly decreasing
391 10623 : do k = 2, nz
392 10623 : if (xq_new(k) <= xq_new(k-1)) then
393 0 : write(*,3) 'bad xq_new before call do_mesh_adjust', &
394 0 : k, nz, xq_new(k), xq_new(k-1), dq_new(k-1), xq_new(k-1) + dq_new(k-1)
395 0 : ierr = -1
396 0 : call dealloc
397 0 : return
398 : call mesa_error(__FILE__,__LINE__,'debug adjust mesh')
399 : end if
400 : end do
401 :
402 10 : call set_m_and_dm(s)
403 10 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
404 :
405 : if (dbg_remesh) write(*,*) 'call do_mesh_adjust'
406 :
407 : call do_mesh_adjust( &
408 : s, nz, nz_old, prv% xh, prv% xa, &
409 : prv% energy, prv% eta, prv% lnd, prv% lnPgas, &
410 : prv% j_rot, prv% i_rot, prv% omega, prv% D_omega, &
411 : prv% mlt_vc, prv% lnT, prv% w, specific_PE, specific_KE, &
412 : prv% m, prv% r, prv% rho, prv% dPdr_dRhodr_info, prv% D_mix, &
413 10 : cell_type, comes_from, prv% dq, xq_old, s% xh, s% xa, s% dq, xq_new, ierr)
414 10 : if (ierr /= 0) then
415 0 : s% retry_message = 'do_mesh_adjust failed in mesh_adjust'
416 0 : if (s% report_ierr) write(*, *) s% retry_message
417 0 : s% termination_code = t_adjust_mesh_failed
418 0 : remesh = terminate
419 0 : call dealloc
420 0 : return
421 : end if
422 : if (dbg_remesh) write(*,*) 'back from do_mesh_adjust'
423 :
424 : ! restore prev_mesh info
425 12105 : do k=1, s% prev_mesh_nz
426 205615 : s% prev_mesh_xa(:,k) = prv% prev_mesh_xa(:,k)
427 108855 : s% prev_mesh_xh(:,k) = prv% prev_mesh_xh(:,k)
428 12095 : s% prev_mesh_j_rot(k) = prv% prev_mesh_j_rot(k)
429 12095 : s% prev_mesh_omega(k) = prv% prev_mesh_omega(k)
430 12095 : s% prev_mesh_dq(k) = prv% prev_mesh_dq(k)
431 12105 : s% prev_mesh_mlt_vc(k) = prv% prev_mesh_mlt_vc(k)
432 : end do
433 :
434 : ! restore ST info (for time smoothing)
435 12105 : do k=1, s% prev_mesh_nz
436 12095 : s% prev_mesh_D_ST_start(k) = prv% prev_mesh_D_ST_start(k)
437 12105 : s% prev_mesh_nu_ST_start(k) = prv% prev_mesh_nu_ST_start(k)
438 : end do
439 :
440 10 : if (s% show_mesh_changes) then
441 : ! note: do_mesh_adjust can change cell_type from unchanged to revised
442 : ! so need to recount
443 0 : unchanged=0; split=0; merged=0; revised=0
444 0 : do k=1,nz
445 0 : select case(cell_type(k))
446 : case (unchanged_type)
447 0 : unchanged = unchanged + 1
448 : case (split_type)
449 0 : split = split + 1
450 : case (merged_type)
451 0 : merged = merged + 1
452 : case (revised_type)
453 0 : revised = revised + 1
454 : case default
455 0 : write(*,3) 'bad value for cell_type(k)', k, cell_type(k)
456 0 : call mesa_error(__FILE__,__LINE__,'adjust_mesh')
457 : end select
458 : end do
459 0 : write(*,*) ' mesh_plan nz_new', nz_new
460 0 : write(*,*) ' unchanged', unchanged
461 0 : write(*,*) ' split', split
462 0 : write(*,*) ' merged', merged
463 0 : write(*,*) ' revised', revised
464 0 : write(*,'(A)')
465 :
466 : end if
467 :
468 : ! testing
469 10623 : do k = 2, nz
470 10623 : if (xq_new(k) <= xq_new(k-1)) then
471 0 : write(*,3) 'bad xq_new after call do_mesh_adjust', k, nz, xq_new(k), xq_new(k-1)
472 0 : ierr = -1
473 0 : call dealloc
474 0 : return
475 : call mesa_error(__FILE__,__LINE__,'debug: adjust mesh')
476 : end if
477 : end do
478 :
479 10 : if (ierr /= 0 .and. s% report_ierr) then
480 0 : write(*,*) 'mesh_adjust problem'
481 0 : write(*,*) 'doing mesh_call_number', s% mesh_call_number
482 0 : write(*,*) 's% model_number', s% model_number
483 0 : write(*,*) 's% nz', s% nz
484 0 : write(*,*) 's% num_retries', s% num_retries
485 0 : write(*,'(A)')
486 : end if
487 :
488 : if (remesh /= keep_going) then
489 : s% result_reason = adjust_mesh_failed
490 : s% retry_message = 'mesh_adjust failed'
491 : if (s% report_ierr) write(*, *) s% retry_message
492 : call dealloc
493 : return
494 : end if
495 :
496 10 : if (s% rotation_flag) then
497 0 : J_tot2 = total_angular_momentum(s)
498 0 : if (abs(J_tot2 - J_tot1) > 1d-2*max(J_tot1,J_tot2)) then
499 0 : s% retry_message = 'adjust_mesh failed'
500 0 : if (s% report_ierr) then
501 0 : write(*,2) 'adjust_mesh J_tot', &
502 0 : s% model_number, (J_tot2 - J_tot1)/J_tot2, J_tot2, J_tot1
503 0 : write(*,*) 'failure to conserve angular momentum in adjust_mesh'
504 0 : write(*,'(A)')
505 : end if
506 0 : ierr = -1
507 0 : call dealloc
508 0 : return
509 : !call mesa_error(__FILE__,__LINE__,'adjust_mesh J_tot conservation error')
510 : end if
511 : end if
512 :
513 10 : call dealloc
514 :
515 20 : if (dbg_remesh .or. dbg) write(*,*) 'finished adjust mesh'
516 :
517 : contains
518 :
519 10 : subroutine set_types_of_new_cells(cell_type, ierr)
520 : integer, pointer :: cell_type(:)
521 : integer, intent(out) :: ierr
522 : integer :: k, k_old, new_type
523 :
524 : include 'formats'
525 10 : ierr = 0
526 10 : unchanged=0; split=0; merged=0; revised=0
527 :
528 10633 : do k=1,nz_new
529 10623 : k_old = comes_from(k)
530 10623 : new_type = -111
531 10623 : if (xq_new(k) < xq_old(k_old)) then
532 0 : write(*,*) 'xq_new(k) < xq_old(k_old)'
533 0 : write(*,2) 'xq_new(k)', k, xq_new(k)
534 0 : write(*,2) 'xq_old(k_old)', k_old, xq_old(k_old)
535 0 : write(*,*) 'adjust mesh set_types_of_new_cells'
536 0 : ierr = -1
537 0 : return
538 10623 : else if (xq_new(k) > xq_old(k_old)) then
539 : new_type = split_type
540 10623 : else if (k_old == nz_old) then
541 10 : if (k == nz_new) then
542 : new_type = unchanged_type
543 : else
544 0 : new_type = split_type
545 : end if
546 10613 : else if (k == nz_new) then
547 : new_type = split_type
548 : else ! k_old < nz_old .and. k < nz .and. xq_new(k) == xq_old(k_old)
549 10613 : if (xq_new(k+1) == xq_old(k_old+1)) then
550 : new_type = unchanged_type
551 1472 : else if (xq_new(k+1) > xq_old(k_old+1)) then
552 : new_type = merged_type
553 : else
554 0 : new_type = split_type
555 : end if
556 : end if
557 10623 : cell_type(k) = new_type
558 10 : select case (new_type)
559 : case (split_type)
560 0 : split = split + 1
561 : case (unchanged_type)
562 9151 : unchanged = unchanged + 1
563 : case (merged_type)
564 1472 : merged = merged + 1
565 : case default
566 : write(*,*) 'failed to set new_type in adjust mesh set_types_of_new_cells'
567 : ierr = -1
568 10623 : return
569 : end select
570 : end do
571 :
572 10 : if (unchanged + split + merged /= nz_new) then
573 0 : write(*,2) 'unchanged + split + merged', unchanged + split + merged
574 0 : write(*,2) 'nz_new', nz_new
575 0 : ierr = -1
576 0 : return
577 : end if
578 :
579 10 : end subroutine set_types_of_new_cells
580 :
581 :
582 10 : subroutine nullify_work_ptrs
583 : nullify( &
584 10 : which_gval, xq_old, xq_new, dq_new, comes_from, &
585 10 : delta_gval_max, do_not_split, gvals, cell_type)
586 : end subroutine nullify_work_ptrs
587 :
588 :
589 10 : subroutine do_alloc1(ierr)
590 : integer, intent(out) :: ierr
591 : ierr = 0
592 10 : call get_integer_work_array(s, comes_from, nz, nz_alloc_extra, ierr)
593 10 : if (ierr /= 0) return
594 10 : call get_integer_work_array(s, which_gval, nz, nz_alloc_extra, ierr)
595 10 : if (ierr /= 0) return
596 10 : call do_work_arrays1(.true., ierr)
597 10 : alloc_level = 1
598 : end subroutine do_alloc1
599 :
600 :
601 10 : subroutine do_dealloc1
602 10 : call return_integer_work_array(s, comes_from)
603 10 : call return_integer_work_array(s, which_gval)
604 10 : call do_work_arrays1(.false., ierr)
605 10 : end subroutine do_dealloc1
606 :
607 :
608 20 : subroutine do_work_arrays1(alloc_flag, ierr)
609 : logical, intent(in) :: alloc_flag
610 : integer, intent(out) :: ierr
611 : logical, parameter :: crit = .false.
612 : ierr = 0
613 : call work_array(s, alloc_flag, crit, &
614 20 : specific_PE, nz, nz_alloc_extra, 'adjust_mesh', ierr)
615 20 : if (ierr /= 0) return
616 : call work_array(s, alloc_flag, crit, &
617 20 : specific_KE, nz, nz_alloc_extra, 'adjust_mesh', ierr)
618 20 : if (ierr /= 0) return
619 : call work_array(s, alloc_flag, crit, &
620 20 : xq_old, nz, nz_alloc_extra, 'adjust_mesh', ierr)
621 20 : if (ierr /= 0) return
622 : call work_array(s, alloc_flag, crit, &
623 20 : xq_new, nz, nz_alloc_extra, 'adjust_mesh', ierr)
624 20 : if (ierr /= 0) return
625 : call work_array(s, alloc_flag, crit, &
626 20 : dq_new, nz, nz_alloc_extra, 'adjust_mesh', ierr)
627 : if (ierr /= 0) return
628 : end subroutine do_work_arrays1
629 :
630 :
631 10 : subroutine do_alloc2(ierr)
632 : integer, intent(out) :: ierr
633 : ierr = 0
634 10 : call do_work_arrays2(.true., ierr)
635 10 : if (ierr /= 0) return
636 10 : call get_logical_work_array(s, do_not_split, nz, nz_alloc_extra, ierr)
637 10 : if (ierr /= 0) return
638 10 : gvals(1:nz,1:num_gvals) => gvals1(1:nz*num_gvals)
639 10 : alloc_level = 2
640 : end subroutine do_alloc2
641 :
642 :
643 10 : subroutine do_dealloc2
644 10 : call do_work_arrays2(.false., ierr)
645 10 : if (ierr /= 0) return
646 10 : call return_logical_work_array(s, do_not_split)
647 : end subroutine do_dealloc2
648 :
649 :
650 20 : subroutine do_work_arrays2(alloc_flag, ierr)
651 : logical, intent(in) :: alloc_flag
652 : integer, intent(out) :: ierr
653 : logical, parameter :: crit = .false.
654 : ierr = 0
655 : call work_array(s, alloc_flag, crit, &
656 20 : delta_gval_max, nz, nz_alloc_extra, 'adjust_mesh', ierr)
657 20 : if (ierr /= 0) return
658 : call work_array(s, alloc_flag, crit, &
659 20 : gvals1, nz*num_gvals, nz_alloc_extra, 'adjust_mesh', ierr)
660 : if (ierr /= 0) return
661 : end subroutine do_work_arrays2
662 :
663 :
664 10 : subroutine do_alloc3(ierr)
665 : integer, intent(out) :: ierr
666 : ierr = 0
667 10 : call get_integer_work_array(s, cell_type, nz, nz_alloc_extra, ierr)
668 10 : alloc_level = 3
669 10 : end subroutine do_alloc3
670 :
671 :
672 10 : subroutine do_dealloc3
673 10 : call return_integer_work_array(s, cell_type)
674 10 : end subroutine do_dealloc3
675 :
676 :
677 10 : subroutine dealloc
678 10 : if (alloc_level >= 3) call do_dealloc3
679 10 : if (alloc_level >= 2) call do_dealloc2
680 10 : if (alloc_level >= 1) call do_dealloc1
681 10 : end subroutine dealloc
682 :
683 :
684 10 : subroutine get_num_gvals
685 : use mesh_functions, only: num_mesh_functions
686 10 : num_gvals = num_mesh_functions(s)
687 10 : end subroutine get_num_gvals
688 :
689 :
690 : subroutine end_dump
691 : write(*,*) ' model_num', s% model_number
692 : write(*,*) ' nz_old', nz_old
693 : write(*,*) ' nz_new', nz_new
694 : write(*,*) 'finished dump_mesh'
695 : call mesa_error(__FILE__,__LINE__,'debugging: end_dump remesh')
696 10 : end subroutine end_dump
697 :
698 : end function remesh
699 :
700 : end module adjust_mesh
|
