Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2016-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_split_merge
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi4, four_thirds_pi
24 : use chem_def, only: ih1, ihe3, ihe4
25 : use utils_lib
26 : use auto_diff_support
27 :
28 : implicit none
29 :
30 : private
31 : public :: remesh_split_merge
32 :
33 : contains
34 :
35 : ! makes new mesh and sets new values for xh and xa.
36 0 : integer function remesh_split_merge(s)
37 : use star_utils, only: total_angular_momentum, set_dm_bar
38 :
39 : type (star_info), pointer :: s
40 0 : real(dp) :: old_J, new_J
41 : integer :: ierr
42 :
43 : include 'formats'
44 :
45 0 : if (s% RSP2_flag) then
46 0 : call mesa_error(__FILE__,__LINE__,'need to add mlt_vc and Hp_face to remesh_split_merge')
47 : end if
48 :
49 0 : s% amr_split_merge_has_undergone_remesh(:) = .false.
50 :
51 0 : remesh_split_merge = keep_going
52 0 : if (.not. s% okay_to_remesh) return
53 :
54 0 : if (s% rotation_flag) old_J = total_angular_momentum(s)
55 :
56 0 : call amr(s,ierr)
57 0 : if (ierr /= 0) then
58 0 : s% retry_message = 'remesh_split_merge failed'
59 0 : if (s% report_ierr) write(*, *) s% retry_message
60 0 : s% result_reason = adjust_mesh_failed
61 0 : s% termination_code = t_adjust_mesh_failed
62 0 : remesh_split_merge = terminate
63 0 : return
64 : end if
65 :
66 0 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
67 :
68 0 : if (s% rotation_flag) then
69 0 : new_J = total_angular_momentum(s)
70 0 : if (abs((old_J-new_J)/old_J)>1d-14) then
71 0 : write(*,*) "Error in angular momemtum conservation from amr split merge"
72 0 : s% result_reason = adjust_mesh_failed
73 0 : s% termination_code = t_adjust_mesh_failed
74 0 : remesh_split_merge = terminate
75 0 : return
76 : end if
77 : end if
78 :
79 0 : end function remesh_split_merge
80 :
81 :
82 0 : subroutine amr(s,ierr)
83 0 : use hydro_rotation, only: w_div_w_roche_jrot, update1_i_rot_from_xh
84 : use star_utils, only: get_r_from_xh
85 : type (star_info), pointer :: s
86 : integer, intent(out) :: ierr
87 0 : real(dp) :: TooBig, TooSmall, MaxTooBig, MaxTooSmall
88 0 : real(dp) :: grad_xa(s% species), new_xa(s% species), &
89 0 : tau_center, r00
90 : integer :: iTooBig, iTooSmall, iter, k, species, &
91 : nz, num_split, num_merge, nz_old
92 : include 'formats'
93 0 : species = s% species
94 0 : nz_old = s% nz
95 0 : ierr = 0
96 0 : num_split = 0
97 0 : num_merge = 0
98 0 : MaxTooSmall = s% split_merge_amr_MaxShort
99 0 : MaxTooBig = s% split_merge_amr_MaxLong
100 0 : k = nz_old
101 : tau_center = s% tau(k) + &
102 0 : s% dm(k)*s% opacity(k)/(pi4*s% rmid(k)*s% rmid(k))
103 0 : do iter = 1, s% split_merge_amr_max_iters
104 : call biggest_smallest(s, tau_center, TooBig, TooSmall, iTooBig, iTooSmall)
105 0 : if (mod(iter,2) == 0) then
106 0 : if (iTooSmall > 0 .and. TooSmall > MaxTooSmall) then
107 0 : call split1
108 0 : else if (iTooBig > 0 .and. TooBig > MaxTooBig) then
109 0 : call merge1
110 : else
111 : exit
112 : end if
113 : else
114 0 : if (iTooBig > 0 .and. TooBig > MaxTooBig) then
115 0 : call merge1
116 0 : else if (iTooSmall > 0 .and. TooSmall > MaxTooSmall) then
117 0 : call split1
118 : else
119 : exit
120 : end if
121 : end if
122 0 : if (ierr /= 0) return
123 : end do
124 0 : do iter = 1, 100
125 : call emergency_merge(s, iTooSmall)
126 0 : if (iTooSmall == 0) exit
127 0 : if (s% split_merge_amr_avoid_repeated_remesh .and. &
128 : s% amr_split_merge_has_undergone_remesh(iTooSmall)) exit
129 0 : if (s% trace_split_merge_amr) &
130 0 : write(*,2) 'emergency_merge', iTooSmall, s% dq(iTooSmall)
131 0 : call do_merge(s, iTooSmall, species, new_xa, ierr)
132 0 : if (ierr /= 0) return
133 0 : num_merge = num_merge + 1
134 0 : if (s% trace_split_merge_amr) then
135 0 : call report_energies(s,'after emergency_merge')
136 : !write(*,*)
137 : end if
138 : end do
139 0 : do iter = 1, 100
140 : call emergency_split(s, iTooBig)
141 0 : if (iTooBig == 0) exit
142 0 : if (s% split_merge_amr_avoid_repeated_remesh .and. &
143 : s% amr_split_merge_has_undergone_remesh(iTooBig)) exit
144 0 : if (s% trace_split_merge_amr) &
145 0 : write(*,2) 'emergency_split', iTooBig, s% dq(iTooBig)
146 0 : call do_split(s, iTooBig, species, tau_center, grad_xa, new_xa, ierr)
147 0 : if (ierr /= 0) return
148 0 : num_split = num_split + 1
149 0 : if (s% trace_split_merge_amr) then
150 0 : call report_energies(s,'after emergency_split')
151 : !write(*,*)
152 : end if
153 : end do
154 0 : s% mesh_call_number = s% mesh_call_number + 1
155 0 : nz = s% nz
156 0 : s% q(nz) = s% dq(nz)
157 0 : s% m(nz) = s% dm(nz) + s% m_center
158 0 : if (s% show_mesh_changes) then
159 0 : write(*,*) 'doing mesh_call_number', s% mesh_call_number
160 0 : write(*,*) ' nz_old', nz_old
161 0 : write(*,*) ' nz_new', nz
162 0 : write(*,*) ' split', num_split
163 0 : write(*,*) ' merged', num_merge
164 : end if
165 0 : if (s% rotation_flag) then !update moments of inertia and omega
166 0 : do k=1, s% nz
167 0 : r00 = get_r_from_xh(s,k)
168 : s% w_div_w_crit_roche(k) = &
169 : w_div_w_roche_jrot(r00,s% m(k),s% j_rot(k),s% cgrav(k), &
170 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
171 0 : call update1_i_rot_from_xh(s, k)
172 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
173 : end do
174 : end if
175 0 : if (s% model_number == -6918) call mesa_error(__FILE__,__LINE__,'amr')
176 :
177 : contains
178 :
179 0 : subroutine split1 ! ratio of desired/actual is too large
180 : include 'formats'
181 0 : if (s% trace_split_merge_amr) then
182 0 : write(*,4) 'do_merge TooSmall', &
183 0 : iTooSmall, s% nz, s% model_number, &
184 0 : TooSmall, MaxTooSmall, s% dq(iTooSmall)
185 0 : call report_energies(s,'before merge TooSmall')
186 : end if
187 0 : call do_merge(s, iTooSmall, species, new_xa, ierr)
188 0 : if (ierr /= 0) return
189 0 : num_merge = num_merge + 1
190 0 : if (s% trace_split_merge_amr) then
191 0 : call report_energies(s,'after merge')
192 : !write(*,*)
193 : end if
194 0 : end subroutine split1
195 :
196 0 : subroutine merge1 ! ratio of actual/desired is too large
197 : include 'formats'
198 0 : if (s% trace_split_merge_amr) then
199 0 : write(*,4) 'do_split TooBig', iTooBig, &
200 0 : s% nz, s% model_number, TooBig, MaxTooBig, s% dq(iTooBig)
201 0 : call report_energies(s,'before split')
202 : end if
203 0 : call do_split(s, iTooBig, species, tau_center, grad_xa, new_xa, ierr)
204 0 : if (ierr /= 0) return
205 0 : num_split = num_split + 1
206 0 : if (s% trace_split_merge_amr) then
207 0 : call report_energies(s,'after split')
208 : !write(*,*)
209 : end if
210 : end subroutine merge1
211 :
212 : end subroutine amr
213 :
214 :
215 0 : subroutine emergency_merge(s, iTooSmall)
216 : type (star_info), pointer :: s
217 : integer, intent(out) :: iTooSmall
218 : integer :: k_min_dq
219 : include 'formats'
220 0 : k_min_dq = minloc(s% dq(1:s% nz),dim=1)
221 0 : if (s% dq(k_min_dq) < s% split_merge_amr_dq_min) then
222 0 : iTooSmall = k_min_dq
223 : else
224 0 : iTooSmall = 0
225 : end if
226 0 : end subroutine emergency_merge
227 :
228 :
229 0 : subroutine emergency_split(s, iTooBig)
230 : type (star_info), pointer :: s
231 : integer, intent(out) :: iTooBig
232 : integer :: k_max_dq
233 : include 'formats'
234 0 : k_max_dq = maxloc(s% dq(1:s% nz),dim=1)
235 0 : if (s% dq(k_max_dq) > s% split_merge_amr_dq_max) then
236 0 : iTooBig = k_max_dq
237 : else
238 0 : iTooBig = 0
239 : end if
240 0 : end subroutine emergency_split
241 :
242 :
243 0 : subroutine biggest_smallest( &
244 : s, tau_center, TooBig, TooSmall, iTooBig, iTooSmall)
245 : type (star_info), pointer :: s
246 : real(dp), intent(in) :: tau_center
247 : real(dp), intent(out) :: TooBig, TooSmall
248 : integer, intent(out) :: iTooBig, iTooSmall
249 : real(dp) :: &
250 0 : oversize_ratio, undersize_ratio, abs_du_div_cs, &
251 0 : xmin, xmax, dx_actual, xR, xL, dq_min, dq_max, dx_baseline, &
252 0 : outer_dx_baseline, inner_dx_baseline, inner_outer_q, r_core_cm, &
253 0 : target_dr_core, target_dlnR_envelope, target_dlnR_core, target_dr_envelope
254 : logical :: hydrid_zoning, flipped_hydrid_zoning, log_zoning, logtau_zoning, &
255 : du_div_cs_limit_flag
256 : integer :: nz, nz_baseline, k, nz_r_core
257 0 : real(dp), pointer :: v(:), r_for_v(:)
258 :
259 : include 'formats'
260 :
261 0 : nz = s% nz
262 0 : hydrid_zoning = s% split_merge_amr_hybrid_zoning
263 0 : flipped_hydrid_zoning = s% split_merge_amr_flipped_hybrid_zoning
264 0 : log_zoning = s% split_merge_amr_log_zoning
265 0 : logtau_zoning = s% split_merge_amr_logtau_zoning
266 0 : nz_baseline = s% split_merge_amr_nz_baseline
267 0 : nz_r_core = s% split_merge_amr_nz_r_core
268 0 : if (s% split_merge_amr_mesh_delta_coeff /= 1d0) then
269 0 : nz_baseline = int(dble(nz_baseline)/s% split_merge_amr_mesh_delta_coeff)
270 0 : nz_r_core = int(dble(nz_r_core)/s% split_merge_amr_mesh_delta_coeff)
271 : end if
272 0 : if (s% split_merge_amr_r_core_cm > 0d0) then
273 : r_core_cm = s% split_merge_amr_r_core_cm
274 0 : else if (s% split_merge_amr_nz_r_core_fraction > 0d0) then
275 : r_core_cm = s% R_center + &
276 0 : s% split_merge_amr_nz_r_core_fraction*(s% r(1) - s% R_center)
277 : else
278 0 : r_core_cm = s% R_center
279 : end if
280 0 : dq_min = s% split_merge_amr_dq_min
281 0 : dq_max = s% split_merge_amr_dq_max
282 0 : inner_outer_q = 0d0
283 0 : if (s% u_flag) then
284 0 : v => s% u
285 0 : r_for_v => s% rmid
286 0 : else if (s% v_flag) then
287 0 : v => s% v
288 0 : r_for_v => s% r
289 : else
290 0 : nullify(v,r_for_v)
291 : end if
292 :
293 0 : if (hydrid_zoning) then
294 0 : target_dr_core = (r_core_cm - s% R_center)/nz_r_core
295 : target_dlnR_envelope = &
296 0 : (s% lnR(1) - log(max(1d0,r_core_cm)))/(nz_baseline - nz_r_core)
297 0 : else if (flipped_hydrid_zoning) then
298 0 : target_dlnR_core = (log(max(1d0,r_core_cm)) - s% R_center)/(nz_baseline - nz_r_core)
299 0 : target_dr_envelope = (s% r(1) - r_core_cm)/nz_r_core
300 0 : else if (logtau_zoning) then
301 : k = nz
302 0 : xmin = log(tau_center)
303 0 : xmax = log(s% tau(1))
304 0 : inner_dx_baseline = (xmin - xmax)/nz_baseline ! keep it > 0
305 0 : outer_dx_baseline = inner_dx_baseline
306 0 : else if (log_zoning) then
307 0 : xmin = log(max(1d0,s% R_center))
308 0 : xmax = s% lnR(1)
309 0 : inner_dx_baseline = (xmax - xmin)/nz_baseline
310 0 : outer_dx_baseline = inner_dx_baseline
311 : else
312 0 : xmin = s% R_center
313 0 : xmax = s% r(1)
314 0 : inner_dx_baseline = (xmax - xmin)/nz_baseline
315 0 : outer_dx_baseline = inner_dx_baseline
316 : end if
317 0 : dx_baseline = outer_dx_baseline
318 :
319 0 : TooBig = 0d0
320 0 : TooSmall = 0d0
321 0 : iTooBig = -1
322 0 : iTooSmall = -1
323 0 : xR = xmin ! start at center
324 0 : do k = nz, 1, -1
325 :
326 0 : xL = xR
327 0 : dx_baseline = inner_dx_baseline
328 0 : if (hydrid_zoning) then
329 0 : if (s% r(k) < r_core_cm) then
330 0 : xR = s% r(k)
331 0 : if (k == nz) then
332 0 : xL = s% R_center
333 : else
334 0 : xL = s% r(k+1)
335 : end if
336 : dx_baseline = target_dr_core
337 : else
338 0 : xR = log(s% r(k))
339 0 : if (k == nz) then
340 0 : xL = log(max(1d0,s% R_center))
341 : else
342 0 : xL = log(s% r(k+1))
343 : end if
344 : dx_baseline = target_dlnR_envelope
345 : end if
346 0 : else if (flipped_hydrid_zoning) then
347 0 : if (s% r(k) <= r_core_cm) then
348 0 : xR = log(s% r(k))
349 0 : if (k == nz) then
350 0 : xL = log(max(1d0,s% R_center))
351 : else
352 0 : xL = log(s% r(k+1))
353 : end if
354 : dx_baseline = target_dlnR_core
355 : else
356 0 : xR = s% r(k)
357 0 : if (k == nz) then
358 0 : xL = s% R_center
359 : else
360 0 : xL = s% r(k+1)
361 : end if
362 : dx_baseline = target_dr_core
363 : end if
364 0 : else if (logtau_zoning) then
365 0 : xR = log(s% tau(k))
366 0 : else if (log_zoning) then
367 0 : xR = log(s% r(k)) ! s% lnR(k) may not be set since making many changes
368 : else
369 0 : xR = s% r(k)
370 : end if
371 :
372 0 : if (s% split_merge_amr_avoid_repeated_remesh .and. &
373 : (s% split_merge_amr_avoid_repeated_remesh .and. &
374 : s% amr_split_merge_has_undergone_remesh(k))) cycle
375 0 : dx_actual = xR - xL
376 0 : if (logtau_zoning) dx_actual = -dx_actual ! make dx_actual > 0
377 :
378 : ! first check for cells that are too big and need to be split
379 0 : oversize_ratio = dx_actual/dx_baseline
380 0 : if (TooBig < oversize_ratio .and. s% dq(k) > 5d0*dq_min) then
381 0 : if (k < nz .or. s% split_merge_amr_okay_to_split_nz) then
382 0 : if (k > 1 .or. s% split_merge_amr_okay_to_split_1) then
383 0 : TooBig = oversize_ratio; iTooBig = k
384 : end if
385 : end if
386 : end if
387 :
388 : ! next check for cells that are too small and need to be merged
389 :
390 0 : if (s% merge_amr_ignore_surface_cells .and. &
391 : k<=s% merge_amr_k_for_ignore_surface_cells) cycle
392 :
393 0 : if (abs(dx_actual)>0d0) then
394 0 : undersize_ratio = max(dx_baseline/dx_actual, dq_min/s% dq(k))
395 : else
396 0 : undersize_ratio = dq_min/s% dq(k)
397 : end if
398 :
399 0 : if (s% merge_amr_max_abs_du_div_cs >= 0d0) then
400 0 : call check_merge_limits
401 0 : else if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then
402 0 : TooSmall = undersize_ratio; iTooSmall = k
403 : end if
404 :
405 : end do
406 :
407 :
408 : contains
409 :
410 0 : subroutine check_merge_limits
411 : ! Pablo's additions to modify when merge
412 : ! merge_amr_max_abs_du_div_cs
413 : ! merge_amr_du_div_cs_limit_only_for_compression
414 : ! merge_amr_inhibit_at_jumps
415 :
416 0 : du_div_cs_limit_flag = .false.
417 :
418 0 : if (.not. s% merge_amr_du_div_cs_limit_only_for_compression) then
419 0 : du_div_cs_limit_flag = .true.
420 0 : else if (associated(v)) then
421 0 : if (k < nz) then
422 0 : if (v(k+1)*pow2(r_for_v(k+1)) > v(k)*pow2(r_for_v(k))) then
423 0 : du_div_cs_limit_flag = .true.
424 : end if
425 : end if
426 0 : if (.not. du_div_cs_limit_flag .and. k > 1) then
427 0 : if (v(k)*pow2(r_for_v(k)) > v(k-1)*pow2(r_for_v(k-1))) then
428 0 : du_div_cs_limit_flag = .true.
429 : end if
430 : end if
431 : end if
432 :
433 0 : if (du_div_cs_limit_flag .and. associated(v)) then
434 0 : if (k == 1) then
435 0 : abs_du_div_cs = abs(v(k) - v(k+1))/s% csound(k)
436 0 : else if (k == nz) then
437 0 : abs_du_div_cs = abs(v(nz-1) - v(nz))/s% csound(nz)
438 : else
439 : abs_du_div_cs = max(abs(v(k) - v(k+1)), &
440 0 : abs(v(k) - v(k-1)))/s% csound(k)
441 : end if
442 : else
443 0 : abs_du_div_cs = 0d0
444 : end if
445 :
446 0 : if (du_div_cs_limit_flag) then
447 0 : if (s% merge_amr_inhibit_at_jumps) then
448 : ! reduce undersize_ratio for large jumps
449 : ! i.e. large jumps inhibit merges but don't prohibit completely
450 0 : if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) &
451 : undersize_ratio = undersize_ratio * &
452 0 : s% merge_amr_max_abs_du_div_cs/abs_du_div_cs
453 0 : if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then ! switch
454 0 : TooSmall = undersize_ratio; iTooSmall = k
455 : end if
456 : else if (TooSmall < undersize_ratio .and. &
457 0 : abs_du_div_cs <= s% merge_amr_max_abs_du_div_cs .and. &
458 : s% dq(k) < dq_max/5d0) then
459 0 : TooSmall = undersize_ratio; iTooSmall = k
460 : end if
461 : else
462 0 : if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then
463 0 : TooSmall = undersize_ratio; iTooSmall = k
464 : end if
465 : end if
466 0 : end subroutine check_merge_limits
467 :
468 :
469 : end subroutine biggest_smallest
470 :
471 :
472 0 : subroutine do_merge(s, i_merge, species, new_xa, ierr)
473 : use mesh_adjust, only: set_lnT_for_energy
474 : use star_utils, only: set_rmid
475 : type (star_info), pointer :: s
476 : integer, intent(in) :: i_merge, species
477 : real(dp), intent(inout) :: new_xa(species)
478 : integer, intent(out) :: ierr
479 : logical :: merge_center
480 : integer :: i, ip, i0, im, q, nz, qi_max, qim_max
481 : real(dp) :: &
482 0 : drR, drL, v, &
483 0 : dm, dm_i, dm_ip, star_PE0, star_PE1, &
484 0 : cell_ie, cell_etrb, &
485 0 : Esum_i, KE_i, PE_i, IE_i, Etrb_i, &
486 0 : Esum_ip, KE_ip, PE_ip, IE_ip, Etrb_ip, &
487 0 : KE, &
488 0 : j_rot_new, j_rot_p1_new, J_old, &
489 0 : dmbar_old, dmbar_p1_old, dmbar_p2_old, &
490 0 : dmbar_new, dmbar_p1_new
491 : include 'formats'
492 :
493 0 : ierr = 0
494 0 : s% need_to_setvars = .true.
495 0 : star_PE0 = get_star_PE(s)
496 0 : nz = s% nz
497 :
498 0 : i = i_merge
499 :
500 0 : s% num_hydro_merges = s% num_hydro_merges+1
501 0 : if (i > 1 .and. i < s% nz) then
502 : ! don't merge across change in most abundance species
503 0 : qi_max = maxloc(s% xa(1:species,i), dim=1)
504 0 : qim_max = maxloc(s% xa(1:species,i-1), dim=1)
505 0 : if (qi_max == qim_max) then ! merge with smaller neighbor
506 0 : if (i+1 == nz) then
507 0 : drL = s% r(nz) - s% R_center
508 : else
509 0 : drL = s% r(i) - s% r(i+1)
510 : end if
511 0 : drR = s% r(i-1) - s% r(i)
512 0 : if (drR < drL) i = i-1
513 : ! else i-1 has different most abundant species,
514 : ! so don't consider merging with it
515 : end if
516 : end if
517 :
518 0 : merge_center = (i == nz)
519 0 : if (merge_center) i = i-1
520 0 : ip = i+1
521 0 : if (s% split_merge_amr_avoid_repeated_remesh .and. &
522 : (s% amr_split_merge_has_undergone_remesh(i) .or. &
523 : s% amr_split_merge_has_undergone_remesh(ip))) then
524 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
525 0 : s% amr_split_merge_has_undergone_remesh(ip) = .true.
526 :
527 0 : return
528 : end if
529 :
530 : ! merge contents of zone ip into zone i; remove zone ip
531 : ! get energies for i and ip before merge
532 0 : call get_cell_energies(s, i, Esum_i, KE_i, PE_i, IE_i, Etrb_i)
533 0 : call get_cell_energies(s, ip, Esum_ip, KE_ip, PE_ip, IE_ip, Etrb_ip)
534 :
535 0 : if (s% rotation_flag) then
536 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
537 0 : j_rot_new = 0d0
538 0 : j_rot_p1_new = 0d0
539 0 : if (i==1) then
540 0 : dmbar_old = 0.5d0*s% dm(i)
541 : else
542 0 : dmbar_old = 0.5d0*(s% dm(i)+s% dm(i-1))
543 : end if
544 0 : if (ip /= nz) then
545 0 : dmbar_p1_old = 0.5d0*(s% dm(ip)+s% dm(ip-1))
546 0 : if (ip+1 /= nz) then
547 0 : dmbar_p2_old = 0.5d0*(s% dm(ip+1)+s% dm(ip))
548 : else
549 0 : dmbar_p2_old = s% dm(nz) + 0.5d0*s% dm(nz-1)
550 : end if
551 : ! before merge we have (for i/=nz)
552 : ! dmbar_old + dmbar_p1_old = 0.5*(m(i-1)+m(i))-0.5*(m(i+1)+m(i+2))
553 : ! after merge we have the merged dm_bar_new
554 : ! dmbar_new = 0.5*(m(i-1)+m(i)) - 0.5*(m(i)+m(i+2)) = 0.5*(m(i-1)-m(i+2))
555 : ! where m is the old mass coordinate. From this we have dmbar_new < dmbar_old,
556 : ! dmbar_old + dmbar_p1_old = dmbar_new + 0.5*(m(i) - m(i+1)) = dmbar_new + 0.5*dm(i)
557 : ! this last expression also holds if i=nz
558 0 : dmbar_new = dmbar_old + dmbar_p1_old - 0.5d0*s% dm(i)
559 : ! this difference in mass goes to the dm_bar below
560 0 : dmbar_p1_new = dmbar_p2_old + 0.5d0*s% dm(i)
561 : ! for the merged cell we take the specific angular momentum of both merged cells
562 0 : J_old = dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip)
563 0 : j_rot_new = J_old/(dmbar_old + dmbar_p1_old)
564 : ! and the j_rot of the cell downwards is adjusted to preserve angular momentum
565 0 : J_old = J_old + dmbar_p2_old*s% j_rot(ip+1)
566 0 : j_rot_p1_new = (J_old - j_rot_new*dmbar_new)/dmbar_p1_new
567 : ! which satisfies:
568 : ! J_old = j_rot_p1_new*dmbar_p1_new + j_rot_new*dmbar_new
569 : else
570 0 : dmbar_p1_old = s% dm(nz) + 0.5d0*s% dm(nz-1)
571 : ! simple case, dmbar_new is equal to dmbar_old plus dmbar_p1_old
572 : ! no need to adjust j_rot of another cell downwards
573 0 : dmbar_new = dmbar_old + dmbar_p1_old
574 0 : j_rot_new = (dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip))/dmbar_new
575 0 : j_rot_p1_new = 0d0
576 : !write(*,*) "check center", i, ip, j_rot_new, j_rot_p1_new
577 : end if
578 : end if
579 :
580 0 : dm_i = s% dm(i)
581 0 : dm_ip = s% dm(ip)
582 0 : dm = dm_i + dm_ip
583 0 : s% dm(i) = dm
584 0 : s% dq(i) = dm/s% xmstar
585 :
586 0 : if (s% RTI_flag) &
587 0 : s% alpha_RTI(i) = (dm_i*s% alpha_RTI(i) + dm_ip*s% alpha_RTI(ip))/dm
588 0 : do q=1,species
589 0 : s% xa(q,i) = (dm_i*s% xa(q,i) + dm_ip*s% xa(q,ip))/dm
590 : end do
591 :
592 0 : KE = KE_i + KE_ip
593 0 : if (s% u_flag) then
594 0 : v = sqrt(KE/(0.5d0*dm))
595 0 : if (s% u(i) < 0d0) v = -v
596 0 : s% u(i) = v
597 : else if (s% v_flag) then
598 : ! there's no good solution for this.
599 : ! so just leave s% v(i) unchanged.
600 : end if
601 :
602 0 : cell_ie = IE_i + IE_ip
603 0 : s% energy(i) = cell_ie/dm
604 :
605 0 : if (s% RSP2_flag) then
606 0 : cell_etrb = Etrb_i + Etrb_ip
607 0 : s% w(i) = sqrt(cell_etrb/dm)
608 : end if
609 :
610 0 : if (s% rotation_flag) then
611 0 : s% j_rot(i) = j_rot_new
612 0 : if (ip/=nz) then
613 0 : s% j_rot(ip+1) = j_rot_p1_new
614 : end if
615 : end if
616 :
617 : ! shift ip+1...nz down by 1 to ip...nz-1
618 0 : do i0 = ip+1, nz
619 0 : im = i0-1
620 : !write(*,3) 'shift to im from i0', im, i0
621 0 : do q=1,species
622 0 : s% xa(q,im) = s% xa(q,i0)
623 : end do
624 0 : do q=1,s% nvar_hydro
625 0 : s% xh(q,im) = s% xh(q,i0)
626 : end do
627 0 : s% r(im) = s% r(i0)
628 0 : s% rmid(im) = s% rmid(i0)
629 0 : s% dm(im) = s% dm(i0)
630 0 : s% m(im) = s% m(i0)
631 0 : s% dq(im) = s% dq(i0)
632 0 : s% q(im) = s% q(i0)
633 0 : if (s% u_flag) then
634 0 : s% u(im) = s% u(i0)
635 0 : else if (s% v_flag) then
636 0 : s% v(im) = s% v(i0)
637 : end if
638 0 : if (s% RSP2_flag) then
639 0 : s% w(im) = s% w(i0)
640 0 : s% Hp_face(im) = s% Hp_face(i0)
641 : end if
642 0 : s% energy(im) = s% energy(i0)
643 0 : s% dPdr_dRhodr_info(im) = s% dPdr_dRhodr_info(i0)
644 0 : s% cgrav(im) = s% cgrav(i0)
645 0 : s% alpha_mlt(im) = s% alpha_mlt(i0)
646 0 : s% lnT(im) = s% lnT(i0)
647 0 : s% D_mix(im) = s% D_mix(i0)
648 0 : s% mlt_vc(im) = s% mlt_vc(i0)
649 0 : s% csound(im) = s% csound(i0)
650 0 : s% tau(im) = s% tau(i0)
651 0 : s% opacity(im) = s% opacity(i0)
652 0 : s% amr_split_merge_has_undergone_remesh(im) = s% amr_split_merge_has_undergone_remesh(i0)
653 0 : if (s% RTI_flag) s% alpha_RTI(im) = s% alpha_RTI(i0)
654 0 : if (s% rotation_flag) s% j_rot(im) = s% j_rot(i0)
655 : end do
656 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
657 :
658 0 : nz = nz - 1
659 0 : s% nz = nz
660 :
661 0 : if (s% u_flag) then
662 0 : s% xh(s% i_u,i) = s% u(i)
663 0 : else if (s% v_flag) then
664 0 : s% xh(s% i_v,i) = s% v(i)
665 : end if
666 :
667 0 : if (s% RTI_flag) s% xh(s% i_alpha_RTI,i) = s% alpha_RTI(i)
668 :
669 0 : if (s% RSP2_flag) then
670 0 : s% xh(s% i_w,i) = s% w(i)
671 0 : s% xh(s% i_Hp,i) = s% Hp_face(i)
672 : end if
673 :
674 : ! do this after move cells since need new r(ip) to calc new rho(i).
675 0 : call update_xh_eos_and_kap(s,i,species,new_xa,ierr)
676 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_merge')
677 :
678 0 : s% rmid_start(i) = -1
679 0 : call set_rmid(s, i, i, ierr)
680 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_merge')
681 :
682 0 : star_PE1 = get_star_PE(s)
683 0 : call revise_star_radius(s, star_PE0, star_PE1)
684 :
685 0 : end subroutine do_merge
686 :
687 :
688 0 : subroutine revise_star_radius(s, star_PE0, star_PE1)
689 0 : use star_utils, only: store_r_in_xh, get_lnR_from_xh
690 : type (star_info), pointer :: s
691 : real(dp), intent(in) :: star_PE0, star_PE1
692 : integer :: k
693 : real(dp) :: frac
694 : include 'formats'
695 0 : if (star_PE1 == 0d0 .or. star_PE0 == star_PE1) return
696 0 : frac = star_PE1/star_PE0
697 0 : if (s% model_number == -6918) write(*,1) 'frac', frac
698 0 : if (s% model_number == -6918) write(*,1) 'star_PE0', star_PE0
699 0 : if (s% model_number == -6918) write(*,1) 'star_PE1', star_PE1
700 0 : do k=1,s% nz
701 0 : s% r(k) = s% r(k)*frac
702 0 : if (s% model_number == -6918) write(*,2) 's% r(k)', k, s% r(k)
703 0 : call store_r_in_xh(s, k, s% r(k))
704 0 : s% lnR(k) = get_lnR_from_xh(s,k)
705 : end do
706 0 : s% r_center = frac*s% r_center
707 0 : end subroutine revise_star_radius
708 :
709 :
710 0 : real(dp) function get_star_PE(s) result(totPE)
711 : type (star_info), pointer :: s
712 : integer :: k
713 0 : real(dp) :: PE, rL, rC, dm, mC
714 0 : totPE = 0d0
715 0 : do k=1,s% nz
716 0 : if (k == s% nz) then
717 0 : rL = s% R_center
718 : else
719 0 : rL = s% r(k+1)
720 : end if
721 0 : rC = 0.5d0*(rL + s% r(k))
722 0 : dm = s% dm(k)
723 0 : mC = s% m(k) - 0.5d0*dm
724 0 : PE = -s% cgrav(k)*mC*dm/rC
725 0 : totPE = totPE + PE
726 : end do
727 0 : end function get_star_PE
728 :
729 :
730 0 : subroutine get_cell_energies(s, k, Etot, KE, PE, IE, Etrb)
731 : type (star_info), pointer :: s
732 : integer, intent(in) :: k
733 : real(dp), intent(out) :: Etot, KE, PE, IE, Etrb
734 0 : real(dp) :: dm, mC, v0, v1, rL, rC
735 : include 'formats'
736 0 : dm = s% dm(k)
737 0 : if (s% u_flag) then
738 0 : KE = 0.5d0*dm*s% u(k)**2
739 0 : else if (s% v_flag) then
740 0 : v0 = s% v(k)
741 0 : if (k < s% nz) then
742 0 : v1 = s% v(k+1)
743 : else
744 0 : v1 = s% v_center
745 : end if
746 0 : KE = 0.25d0*dm*(v0**2 + v1**2)
747 : else
748 0 : KE = 0d0
749 : end if
750 0 : IE = s% energy(k)*dm
751 0 : if (s% RSP2_flag) then
752 0 : Etrb = pow2(s% w(k))*dm
753 : else
754 0 : Etrb = 0d0
755 : end if
756 0 : if (s% cgrav(k) == 0) then
757 0 : PE = 0
758 : else
759 0 : if (k == s% nz) then
760 0 : rL = s% R_center
761 : else
762 0 : rL = s% r(k+1)
763 : end if
764 0 : rC = 0.5d0*(rL + s% r(k))
765 0 : mC = s% m(k) - 0.5d0*dm
766 0 : PE = -s% cgrav(k)*mC*dm/rC
767 : end if
768 0 : Etot = IE + KE + PE
769 : if (is_bad(Etot + IE + KE + PE) .or. &
770 0 : IE <= 0 .or. KE < 0) then
771 0 : write(*,2) 'nz', s% nz
772 0 : write(*,2) 'Etot', k, Etot
773 0 : write(*,2) 'IE', k, IE
774 0 : write(*,2) 'KE', k, KE
775 0 : write(*,2) 'PE', k, PE
776 0 : call mesa_error(__FILE__,__LINE__,'get_cell_energies')
777 : end if
778 0 : end subroutine get_cell_energies
779 :
780 :
781 0 : subroutine split1_non_negative( &
782 : val, grad, dr, dV, dVR, dVL, new_valL, new_valR)
783 : real(dp), intent(in) :: val, grad, dr, dV, dVR, dVL
784 : real(dp), intent(out) :: new_valL, new_valR
785 0 : real(dp) :: xL, xR, xML, xMR, xM, f
786 : include 'formats'
787 0 : new_valR = val
788 0 : new_valL = val
789 0 : if (val < 1d-99 .or. grad == 0d0) return
790 0 : xL = max(0d0, min(1d0, val - grad*dr/4))
791 0 : xR = max(0d0, min(1d0, val + grad*dr/4))
792 0 : xML = xL*dVL
793 0 : xMR = xR*dVR
794 0 : xM = val*dV
795 0 : if (xM < 1d-99 .or. xML + xMR < 1d-99) then
796 : new_valR = val
797 : new_valL = val
798 : else
799 0 : f = xM/(xML + xMR)
800 0 : new_valR = f*xR
801 0 : new_valL = f*xL
802 : end if
803 : end subroutine split1_non_negative
804 :
805 :
806 0 : subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
807 : use alloc, only: reallocate_star_info_arrays
808 : use star_utils, only: set_rmid, store_r_in_xh
809 : type (star_info), pointer :: s
810 : integer, intent(in) :: i_split, species
811 : real(dp) :: tau_center, grad_xa(species), new_xa(species)
812 : integer, intent(out) :: ierr
813 : integer :: i, ip, j, jp, q, nz, nz_old, &
814 : iR, iC, iL
815 : real(dp) :: &
816 : cell_Esum_old, cell_KE_old, cell_PE_old, cell_IE_old, cell_Etrb_old, &
817 0 : rho_RR, rho_iR, rR, rL, dr, dr_old, rC, dV, dVR, dVL, dM, dML, dMR, rho, &
818 0 : v, v2, energy, v2_R, energy_R, rho_R, v2_C, energy_C, rho_C, v2_L, energy_L, rho_L, &
819 0 : dLeft, dRght, dCntr, grad_rho, grad_energy, grad_v2, &
820 0 : sumx, sumxp, new_xaL, new_xaR, star_PE0, star_PE1, &
821 0 : grad_alpha, f, new_alphaL, new_alphaR, v_R, v_C, v_L, min_dm, &
822 0 : mlt_vcL, mlt_vcR, tauL, tauR, etrb, etrb_L, etrb_C, etrb_R, grad_etrb, &
823 0 : j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old
824 : logical :: done, use_new_grad_rho
825 : include 'formats'
826 :
827 0 : ierr = 0
828 0 : star_PE0 = get_star_PE(s)
829 0 : s% need_to_setvars = .true.
830 0 : nz = s% nz
831 0 : s% num_hydro_splits = s% num_hydro_splits + 1
832 0 : done = .false.
833 0 : nz_old = nz
834 :
835 0 : i = i_split
836 0 : ip = i+1
837 :
838 : call get_cell_energies( &
839 0 : s, i, cell_Esum_old, cell_KE_old, cell_PE_old, cell_IE_old, cell_Etrb_old)
840 :
841 0 : if (s% rotation_flag .and. i<nz) then
842 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
843 0 : if (i==1) then
844 0 : dmbar_old = 0.5d0*s% dm(i)
845 : else
846 0 : dmbar_old = 0.5d0*(s% dm(i)+s% dm(i-1))
847 : end if
848 0 : if (ip /= nz) then
849 0 : dmbar_p1_old = 0.5d0*(s% dm(ip)+s% dm(ip-1))
850 : else
851 0 : dmbar_p1_old = s% dm(ip)+0.5d0*s% dm(ip-1)
852 : end if
853 :
854 : ! We need to conserve the angular momentum in these two cells
855 0 : J_old = dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip)
856 : end if
857 :
858 0 : rR = s% r(i)
859 0 : mlt_vcR = s% mlt_vc(i)
860 0 : if (i == nz) then
861 0 : rL = s% R_center
862 0 : mlt_vcL = 0d0
863 : tauL = tau_center
864 : else
865 0 : rL = s% r(ip)
866 0 : mlt_vcL = s% mlt_vc(ip)
867 : tauL = s% tau(ip)
868 : end if
869 :
870 0 : tauR = s% tau(i)
871 0 : if (i == nz) then
872 0 : tauL = tau_center
873 : else
874 0 : tauL = s% tau(ip)
875 : end if
876 0 : if (is_bad(tauL+tauR)) then
877 0 : !$omp critical (adjust_mesh_split_merge_crit1)
878 0 : write(*,2) 'tauL', ip, tauL
879 0 : write(*,2) 'tauR', i, tauR
880 0 : write(*,2) 'nz', nz
881 0 : call mesa_error(__FILE__,__LINE__,'do_split')
882 : !$omp end critical (adjust_mesh_split_merge_crit1)
883 : end if
884 :
885 0 : dr = rR - rL
886 0 : dr_old = dr
887 0 : rC = 0.5d0*(rR + rL) ! split at center by radius
888 :
889 0 : dV = four_thirds_pi*(rR*rR*rR - rL*rL*rL)
890 0 : dVR = four_thirds_pi*(rR*rR*rR - rC*rC*rC)
891 0 : dVL = dV - dVR
892 :
893 0 : dm = s% dm(i)
894 0 : rho = dm/dV
895 :
896 0 : if (s% u_flag) then
897 0 : v = s% u(i)
898 0 : v2 = v*v
899 : else ! not used
900 : v = 0
901 : v2 = 0
902 : end if
903 :
904 0 : energy = s% energy(i)
905 0 : if (s% RSP2_flag) etrb = pow2(s% w(i))
906 :
907 : ! use iR, iC, and iL for getting values to determine slopes
908 0 : if (i > 1 .and. i < nz_old) then
909 0 : iR = i-1
910 0 : iC = i
911 0 : iL = i+1
912 0 : else if (i == 1) then
913 0 : iR = 1
914 0 : iC = 2
915 0 : iL = 3
916 : else ! i == nz_old
917 0 : iR = nz_old-2
918 0 : iC = nz_old-1
919 0 : iL = nz_old
920 : end if
921 :
922 0 : energy_R = s% energy(iR)
923 0 : rho_R = s% dm(iR)/get_dV(s,iR)
924 :
925 0 : energy_C = s% energy(iC)
926 0 : rho_C = s% dm(iC)/get_dV(s,iC)
927 :
928 0 : energy_L = s% energy(iL)
929 0 : rho_L = s% dm(iL)/get_dV(s,iL)
930 :
931 : ! get gradients before move cell contents
932 :
933 0 : if (iL == nz_old) then
934 0 : dLeft = 0.5d0*(s% r(iC) - s% R_center)
935 : else
936 0 : dLeft = 0.5d0*(s% r(iC) - s% r(iL+1))
937 : end if
938 0 : dRght = 0.5d0*(s% r(iR) - s% r(iL))
939 0 : dCntr = dLeft + dRght
940 :
941 0 : if (s% equal_split_density_amr) then ! same density in both parts
942 0 : rho_RR = 0
943 : grad_rho = 0d0
944 : use_new_grad_rho = .false.
945 : else if (.false.) then
946 : rho_RR = s% dm(iR-1)/get_dV(s,iR-1)
947 : grad_rho = 2*(rho_RR - rho_R)/(s% r(iR-1) - s% r(iR+1))
948 : use_new_grad_rho = .true.
949 : else
950 0 : rho_RR = 0
951 0 : grad_rho = get1_grad(rho_L, rho_C, rho_R, dLeft, dCntr, dRght)
952 0 : use_new_grad_rho = .false.
953 : end if
954 :
955 0 : grad_energy = get1_grad(energy_L, energy_C, energy_R, dLeft, dCntr, dRght)
956 :
957 0 : if (s% RTI_flag) grad_alpha = get1_grad( &
958 0 : s% alpha_RTI(iL), s% alpha_RTI(iC), s% alpha_RTI(iR), dLeft, dCntr, dRght)
959 :
960 0 : if (s% RSP2_flag) then
961 0 : etrb_R = pow2(s% w(iR))
962 0 : etrb_C = pow2(s% w(iC))
963 0 : etrb_L = pow2(s% w(iL))
964 0 : grad_etrb = get1_grad(etrb_L, etrb_C, etrb_R, dLeft, dCntr, dRght)
965 : end if
966 :
967 0 : if (s% u_flag) then
968 0 : v_R = s% u(iR)
969 0 : v2_R = v_R*v_R
970 0 : v_C = s% u(iC)
971 0 : v2_C = v_C*v_C
972 0 : v_L = s% u(iL)
973 0 : v2_L = v_L*v_L
974 0 : if ((v_L - v_C)*(v_C - v_R) <= 0) then ! not strictly monotonic velocities
975 : grad_v2 = 0d0
976 : else
977 0 : grad_v2 = get1_grad(v2_L, v2_C, v2_R, dLeft, dCntr, dRght)
978 : end if
979 0 : else if (s% v_flag) then
980 0 : if (iL == s% nz) then
981 0 : v_L = s% v_center
982 : else
983 0 : v_L = s% v(ip)
984 : end if
985 0 : v2_L = v_L*v_L
986 0 : v_R = s% v(i)
987 0 : v2_R = v_R*v_R
988 : end if
989 :
990 0 : do q = 1, species
991 : grad_xa(q) = get1_grad( &
992 0 : s% xa(q,iL), s% xa(q,iC), s% xa(q,iR), dLeft, dCntr, dRght)
993 : end do
994 :
995 0 : nz = nz + 1
996 0 : s% nz = nz
997 0 : call reallocate_star_info_arrays(s, ierr)
998 0 : if (ierr /= 0) then
999 0 : write(*,2) 'reallocate_star_info_arrays ierr', ierr
1000 0 : call mesa_error(__FILE__,__LINE__,'split failed')
1001 : end if
1002 :
1003 0 : if (i_split < nz_old) then ! move i_split..nz-1 to i_split+1..nz
1004 0 : do j=nz-1,i_split,-1
1005 0 : jp = j+1
1006 0 : do q=1,species
1007 0 : s% xa(q,jp) = s% xa(q,j)
1008 : end do
1009 0 : do q=1,s% nvar_hydro
1010 0 : s% xh(q,jp) = s% xh(q,j)
1011 : end do
1012 0 : s% r(jp) = s% r(j)
1013 0 : s% rmid(jp) = s% rmid(j)
1014 0 : s% dm(jp) = s% dm(j)
1015 0 : s% m(jp) = s% m(j)
1016 0 : s% dq(jp) = s% dq(j)
1017 0 : s% q(jp) = s% q(j)
1018 0 : if (s% u_flag) then
1019 0 : s% u(jp) = s% u(j)
1020 0 : else if (s% v_flag) then
1021 0 : s% v(jp) = s% v(j)
1022 : end if
1023 0 : s% energy(jp) = s% energy(j)
1024 0 : s% dPdr_dRhodr_info(jp) = s% dPdr_dRhodr_info(j)
1025 0 : s% cgrav(jp) = s% cgrav(j)
1026 0 : s% alpha_mlt(jp) = s% alpha_mlt(j)
1027 0 : s% lnT(jp) = s% lnT(j)
1028 0 : s% D_mix(jp) = s% D_mix(j)
1029 0 : s% mlt_vc(jp) = s% mlt_vc(j)
1030 0 : s% csound(jp) = s% csound(j)
1031 0 : s% tau(jp) = s% tau(j)
1032 0 : s% opacity(jp) = s% opacity(j)
1033 0 : s% amr_split_merge_has_undergone_remesh(jp) = s% amr_split_merge_has_undergone_remesh(j)
1034 0 : if (s% RTI_flag) s% alpha_RTI(jp) = s% alpha_RTI(j)
1035 0 : if (s% rotation_flag) s% j_rot(jp) = s% j_rot(j)
1036 : end do
1037 : end if
1038 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
1039 0 : s% amr_split_merge_has_undergone_remesh(ip) = .true.
1040 :
1041 0 : dM = rho*dV
1042 : if (.not. use_new_grad_rho) then ! do it the old way
1043 :
1044 0 : rho_L = rho - grad_rho*dr/4
1045 0 : rho_R = rho + grad_rho*dr/4
1046 : if (grad_rho == 0d0 .or. &
1047 0 : rho_L <= 0d0 .or. &
1048 : rho_R <= 0d0) then
1049 0 : rho_R = rho
1050 0 : rho_L = rho
1051 0 : dML = rho_L*dVL
1052 0 : dMR = dM - dML ! should = rhoR*dVR
1053 : else
1054 0 : dML = rho_L*dVL
1055 0 : dMR = rho_R*dVR
1056 0 : f = dM/(dML + dMR)
1057 0 : rho_L = f*rho_L
1058 0 : rho_R = f*rho_R
1059 0 : dML = rho_L*dVL
1060 0 : dMR = rho_R*dVR
1061 0 : if (abs(dML + dMR - dM) > 1d-14*dM) then
1062 0 : write(*,2) '(dML + dMR - dM)/dM', i, (dML + dMR - dM)/dM
1063 0 : call mesa_error(__FILE__,__LINE__,'split')
1064 : end if
1065 0 : dMR = dM - dML
1066 : end if
1067 :
1068 : else
1069 :
1070 : ! at this point, rho_R is the density of the cell iR
1071 : ! we are about to redefine it as the density of the right 1/2 of the split
1072 : ! similarly for rho_L
1073 : rho_iR = rho_R
1074 : dR = -(dRght/4 + (s% r(iR) - s% r(iC))/2)
1075 : rho_R = rho_iR + grad_rho*dR
1076 : dMR = rho_R*dVR
1077 : dML = dM - dMR
1078 : rho_L = dML/dVL
1079 : if (rho_R <= 1d-16 .or. rho_L <= 1d-16) then
1080 : !$omp critical (adjust_mesh_split_merge_crit2)
1081 : write(*,'(A)')
1082 : write(*,2) 'nz', nz
1083 : write(*,2) 'rho_RR', iR-1, rho_RR
1084 : write(*,2) 'rho_iR', iR, rho_iR
1085 : write(*,2) 'rho', iC, rho
1086 : write(*,2) 's% r(iR-1)', iR-1, s% r(iR-1)
1087 : write(*,2) 's% r(iR)', iR, s% r(iR)
1088 : write(*,2) 's% r(iR+1)', iR+1, s% r(iR+1)
1089 : write(*,2) 'rho_RR', iR-1, rho_RR
1090 : write(*,2) 'rho_iR', iR, rho_iR
1091 : write(*,2) 'rho_RR - rho_iR', iR, rho_RR - rho_iR
1092 : write(*,2) 'dr for right', iR, (s% r(iR-1) - s% r(iR+1))/2
1093 : write(*,2) 'grad_rho', iR, grad_rho
1094 : write(*,'(A)')
1095 : write(*,2) 's% r(iL)', iL, s% r(iL)
1096 : write(*,2) 'dR', iC, dR
1097 : write(*,2) 'dRho', iC, grad_rho*dR
1098 : write(*,2) 'rho_R', iC, rho_R
1099 : write(*,2) 'rho_L', iC, rho_L
1100 : write(*,'(A)')
1101 : call mesa_error(__FILE__,__LINE__,'failed in do_split extrapolation of density from above')
1102 : !$omp end critical (adjust_mesh_split_merge_crit2)
1103 : end if
1104 :
1105 : end if
1106 :
1107 0 : min_dm = s% split_merge_amr_dq_min*s% xmstar
1108 0 : if (dML < min_dm .or. dMR < min_dm) then
1109 0 : rho_R = rho
1110 0 : rho_L = rho
1111 0 : dM = rho*dV
1112 0 : dML = rho_L*dVL
1113 0 : dMR = dM - dML ! should = rhoR*dVR
1114 : end if
1115 :
1116 0 : energy_R = energy + grad_energy*dr/4
1117 0 : energy_L = (dm*energy - dmR*energy_R)/dmL
1118 0 : if (energy_R < 0d0 .or. energy_L < 0d0) then
1119 0 : energy_R = energy
1120 0 : energy_L = energy
1121 : end if
1122 :
1123 0 : s% energy(i) = energy_R
1124 0 : s% energy(ip) = energy_L
1125 :
1126 0 : if (s% RSP2_flag) then
1127 0 : etrb_R = etrb + grad_etrb*dr/4
1128 0 : etrb_L = (dm*etrb - dmR*etrb_R)/dmL
1129 0 : if (etrb_R < 0d0 .or. etrb_L < 0d0) then
1130 0 : etrb_R = etrb
1131 0 : etrb_L = etrb
1132 : end if
1133 0 : s% w(i) = sqrt(max(0d0,etrb_R))
1134 0 : s% w(ip) = sqrt(max(0d0,etrb_L))
1135 : end if
1136 :
1137 0 : if (s% u_flag) then
1138 0 : v2_R = v2 + grad_v2*dr/4
1139 0 : v2_L = (dm*v2 - dmR*v2_R)/dmL
1140 0 : if (v2_R < 0d0 .or. v2_L < 0d0) then
1141 0 : v2_R = v2
1142 0 : v2_L = v2
1143 : end if
1144 0 : s% u(i) = sqrt(v2_R)
1145 0 : s% u(ip) = sqrt(v2_L)
1146 0 : if (v < 0d0) then
1147 0 : s% u(i) = -s% u(i)
1148 0 : s% u(ip) = -s% u(ip)
1149 : end if
1150 0 : else if (s% v_flag) then ! just make a rough approximation.
1151 0 : s% v(ip) = sqrt(0.5d0*(v2_L + v2_R))
1152 : end if
1153 :
1154 0 : if (s% RTI_flag) then ! set new alpha
1155 0 : if (i == 1) then
1156 0 : s% alpha_RTI(ip) = s% alpha_RTI(i)
1157 : else
1158 : call split1_non_negative( &
1159 : s% alpha_RTI(i), grad_alpha, &
1160 0 : dr, dV, dVR, dVL, new_alphaL, new_alphaR)
1161 0 : s% alpha_RTI(i) = new_alphaR
1162 0 : s% alpha_RTI(ip) = new_alphaL
1163 : end if
1164 0 : s% dPdr_dRhodr_info(ip) = s% dPdr_dRhodr_info(i)
1165 : end if
1166 :
1167 0 : if (i == 1) then
1168 0 : s% mlt_vc(ip) = s% mlt_vc(i)
1169 : else
1170 0 : s% mlt_vc(ip) = (mlt_vcL*dML + mlt_vcR*dMR)/dM
1171 : end if
1172 :
1173 0 : s% tau(ip) = tauR + (tauL - tauR)*dMR/dM
1174 0 : if (is_bad(s% tau(ip))) then
1175 0 : write(*,2) 'tau', ip, s% tau(ip), tauL, tauR, dMR/dM
1176 0 : call mesa_error(__FILE__,__LINE__,'split')
1177 : end if
1178 :
1179 0 : if (i == 1) then
1180 0 : do q=1,species
1181 0 : s% xa(q,ip) = s% xa(q,i)
1182 : end do
1183 : else ! split mass fractions
1184 : ! check that sum of grads for mass fractions is 0
1185 0 : sumx = sum(grad_xa)
1186 0 : if (sumx > 0) then ! reduce largest grad
1187 0 : j = maxloc(grad_xa, dim=1)
1188 : else ! increase smallest grad
1189 0 : j = minloc(grad_xa, dim=1)
1190 : end if
1191 0 : grad_xa(j) = 0
1192 0 : grad_xa(j) = -sum(grad_xa)
1193 : ! set new mass fractions
1194 0 : do q = 1, species
1195 : call split1_non_negative( &
1196 : s% xa(q,i), grad_xa(q), &
1197 0 : dr, dV, dVR, dVL, new_xaL, new_xaR)
1198 0 : s% xa(q,i) = new_xaR
1199 0 : s% xa(q,ip) = new_xaL
1200 : end do
1201 : !check mass fractions >= 0 and <= 1 and sum to 1.0
1202 0 : do q = 1, species
1203 0 : s% xa(q,i) = min(1d0, max(0d0, s% xa(q,i)))
1204 0 : s% xa(q,ip) = min(1d0, max(0d0, s% xa(q,ip)))
1205 : end do
1206 0 : sumx = sum(s% xa(1:species,i))
1207 0 : sumxp = sum(s% xa(1:species,ip))
1208 0 : do q = 1, species
1209 0 : s% xa(q,i) = s% xa(q,i)/sumx
1210 0 : s% xa(q,ip) = s% xa(q,ip)/sumxp
1211 : end do
1212 : end if
1213 :
1214 0 : if (s% u_flag) s% u_face_ad(ip)%val = 0.5d0*(s% u(i) + s% u(ip))
1215 : ! just for setting u_face_start so don't need partials
1216 :
1217 : ! r, q, m, u_face unchanged for face i
1218 0 : dVR = dV - dVL
1219 0 : dmR = s% dm(i) - dmL
1220 0 : s% dm(i) = dmR
1221 0 : s% dq(i) = s% dm(i)/s% xmstar
1222 :
1223 0 : s% r(ip) = rC
1224 0 : s% m(ip) = s% m(i) - s% dm(i) ! <<< using new value dm(i)
1225 0 : s% q(ip) = s% q(i) - s% dq(i) ! <<< using new value dq(i)
1226 0 : s% dm(ip) = dmL
1227 0 : s% dq(ip) = s% dm(ip)/s% xmstar
1228 :
1229 0 : s% cgrav(ip) = s% cgrav(i)
1230 0 : s% alpha_mlt(ip) = s% alpha_mlt(i)
1231 0 : s% lnT(ip) = s% lnT(i)
1232 0 : s% T(ip) = s% T(i)
1233 0 : s% D_mix(ip) = s% D_mix(i)
1234 0 : s% mlt_vc(ip) = s% mlt_vc(i)
1235 0 : s% csound(ip) = s% csound(i)
1236 0 : s% opacity(ip) = s% opacity(i)
1237 0 : if (s% RTI_flag) s% alpha_RTI(ip) = s% alpha_RTI(i)
1238 :
1239 0 : if (s% rotation_flag) then
1240 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
1241 0 : j_rot_new = 0d0 ! specific angular momentum for the newly created cell
1242 0 : if (i==nz_old) then
1243 : ! simple case, dm_bar contained in cells i and i+1 after merge is conserved,
1244 : ! so use same j_rot
1245 0 : j_rot_new = s% j_rot(i)
1246 : else
1247 0 : if (i==1) then
1248 0 : dmbar_new = 0.5d0*s% dm(i)
1249 : else
1250 0 : dmbar_new = 0.5d0*(s% dm(i)+s% dm(i-1))
1251 : end if
1252 0 : if (ip /= nz) then
1253 0 : dmbar_p1_new = 0.5d0*(s% dm(ip)+s% dm(ip-1))
1254 : else
1255 0 : dmbar_p1_new = s% dm(ip)+0.5d0*s% dm(ip-1)
1256 : end if
1257 0 : if (ip+1 /= nz) then
1258 0 : dmbar_p2_new = 0.5d0*(s% dm(ip+1)+s% dm(ip))
1259 : else
1260 0 : dmbar_p2_new = s% dm(ip+1)+0.5d0*s% dm(ip)
1261 : end if
1262 :
1263 : ! we keep the same j_rot for cells i and ip (cells i and ip+1 after merge),
1264 : ! we compute the j_rot needed for the new cell to preserve angular momentum.
1265 0 : j_rot_new = (J_old - dmbar_new*s% j_rot(i) - dmbar_p2_new*s% j_rot(ip+1))/(dmbar_p1_new)
1266 : ! this satisfies:
1267 : ! dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip) = dmbar_new*s% j_rot(i) + dmbar_p2_new*s% j_rot(ip) + dmbar_p1_new*j_rot_new
1268 : end if
1269 0 : s% j_rot(ip) = j_rot_new
1270 : end if
1271 :
1272 0 : if (s% i_lum /= 0) then
1273 0 : if (ip < nz_old) then
1274 : s% xh(s% i_lum,ip) = &
1275 0 : 0.5d0*(s% xh(s% i_lum,i) + s% xh(s% i_lum,ip+1))
1276 : else
1277 : s% xh(s% i_lum,ip) = &
1278 0 : 0.5d0*(s% xh(s% i_lum,i) + s% L_center)
1279 : end if
1280 : end if
1281 :
1282 0 : call store_r_in_xh(s, ip, s% r(ip))
1283 0 : if (s% u_flag) then
1284 0 : s% xh(s% i_u,i) = s% u(i)
1285 0 : s% xh(s% i_u,ip) = s% u(ip)
1286 0 : else if (s% v_flag) then
1287 0 : s% xh(s% i_v,i) = s% v(i)
1288 0 : s% xh(s% i_v,ip) = s% v(ip)
1289 : end if
1290 0 : if (s% RTI_flag) then
1291 0 : s% xh(s% i_alpha_RTI,i) = s% alpha_RTI(i)
1292 0 : s% xh(s% i_alpha_RTI,ip) = s% alpha_RTI(ip)
1293 : end if
1294 :
1295 0 : call update_xh_eos_and_kap(s,i,species,new_xa,ierr)
1296 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1297 :
1298 0 : call update_xh_eos_and_kap(s,ip,species,new_xa,ierr)
1299 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1300 :
1301 0 : s% rmid_start(i) = -1
1302 0 : s% rmid_start(ip) = -1
1303 0 : call set_rmid(s, i, ip, ierr)
1304 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1305 :
1306 0 : star_PE1 = get_star_PE(s)
1307 0 : call revise_star_radius(s, star_PE0, star_PE1)
1308 :
1309 0 : end subroutine do_split
1310 :
1311 :
1312 0 : subroutine update_xh_eos_and_kap(s,i,species,new_xa,ierr)
1313 0 : use mesh_adjust, only: set_lnT_for_energy
1314 : use micro, only: do_kap_for_cell
1315 : use eos_lib, only: eos_gamma_DE_get_PT
1316 : use chem_lib, only: basic_composition_info
1317 : use star_utils, only: store_lnT_in_xh, get_T_and_lnT_from_xh, &
1318 : store_rho_in_xh, get_rho_and_lnd_from_xh
1319 : type (star_info), pointer :: s
1320 : integer, intent(in) :: i, species
1321 : real(dp) :: new_xa(species)
1322 : integer, intent(out) :: ierr
1323 : real(dp) :: rho, logRho, new_lnT, revised_energy
1324 : integer :: q
1325 : include 'formats'
1326 0 : ierr = 0
1327 0 : rho = s% dm(i)/get_dV(s,i)
1328 0 : call store_rho_in_xh(s, i, rho)
1329 0 : call get_rho_and_lnd_from_xh(s, i, s% rho(i), s% lnd(i))
1330 0 : logRho = s% lnd(i)/ln10
1331 0 : do q=1,species
1332 0 : new_xa(q) = s% xa(q,i)
1333 : end do
1334 : call set_lnT_for_energy(s, i, &
1335 : s% net_iso(ih1), s% net_iso(ihe3), s% net_iso(ihe4), &
1336 : species, new_xa, rho, logRho, s% energy(i), s% lnT(i), &
1337 0 : new_lnT, revised_energy, ierr)
1338 0 : if (ierr /= 0) return
1339 0 : call store_lnT_in_xh(s, i, new_lnT)
1340 0 : call get_T_and_lnT_from_xh(s, i, s% T(i), s% lnT(i))
1341 0 : end subroutine update_xh_eos_and_kap
1342 :
1343 :
1344 0 : real(dp) function get_dV(s,k)
1345 : type (star_info), pointer :: s
1346 : integer, intent(in) :: k
1347 0 : if (k == s% nz) then
1348 0 : get_dV = four_thirds_pi*(s% r(k)*s% r(k)*s% r(k) - s% R_center*s% R_center*s% R_center)
1349 : else
1350 0 : get_dV = four_thirds_pi*(s% r(k)*s% r(k)*s% r(k) - s% r(k+1)*s% r(k+1)*s% r(k+1))
1351 : end if
1352 0 : end function get_dV
1353 :
1354 :
1355 0 : real(dp) function minmod(a, b, c) result(m)
1356 : real(dp), intent(in) :: a, b, c
1357 0 : m = a
1358 0 : if (a*b < 0d0) then
1359 0 : m = 0d0
1360 : return
1361 : end if
1362 0 : if (abs(b) < abs(m)) m = b
1363 0 : if (b*c < 0d0) then
1364 0 : m = 0d0
1365 : return
1366 : end if
1367 0 : if (abs(c) < abs(m)) m = c
1368 : end function minmod
1369 :
1370 :
1371 0 : real(dp) function get1_grad(vL, vC, vR, dLeft, dCntr, dRght) &
1372 0 : result(grad)
1373 : real(dp), intent(in) :: vL, vC, vR, dLeft, dCntr, dRght
1374 0 : real(dp) :: sL, sR, sC
1375 : include 'formats'
1376 0 : if (dLeft <= 0d0 .or. dCntr <= 0d0 .or. dRght <= 0d0) then
1377 0 : grad = 0d0
1378 : return
1379 : end if
1380 0 : sL = (vC - vL)/dLeft
1381 0 : sR = (vR - vC)/dRght
1382 0 : sC = (vR - vL)/dCntr
1383 0 : grad = minmod(sL, sC, sR)
1384 0 : end function get1_grad
1385 :
1386 :
1387 : real(dp) function total_KE(s)
1388 : type (star_info), pointer :: s
1389 : integer :: k
1390 : include 'formats'
1391 : total_KE = 0
1392 : if (s% u_flag) then
1393 : do k=1,s% nz
1394 : total_KE = total_KE + 0.5d0*s% dm(k)*s% u(k)**2
1395 : end do
1396 : else if (s% v_flag) then
1397 : do k=1,s% nz-1
1398 : total_KE = total_KE + &
1399 : 0.25d0*s% dm(k)*(s% v(k)**2 + s% v(k+1)**2)
1400 : end do
1401 : k = s% nz
1402 : total_KE = total_KE + &
1403 : 0.25d0*s% dm(k)*(s% v(k)**2 + s% v_center**2)
1404 : end if
1405 : end function total_KE
1406 :
1407 :
1408 : real(dp) function total_PE(s)
1409 : type (star_info), pointer :: s
1410 : integer :: k
1411 : real(dp) :: rL, rR, rC, mC
1412 : total_PE = 0d0
1413 : rR = s% R_center
1414 : do k = s% nz,1,-1
1415 : rL = rR
1416 : rR = s% r(k)
1417 : rC = 0.5d0*(rL + rR)
1418 : mC = s% m(k) - 0.5d0*s% dm(k)
1419 : total_PE = total_PE - s% cgrav(k)*mC*s% dm(k)/rC
1420 : end do
1421 : end function total_PE
1422 :
1423 :
1424 : real(dp) function total_IE(s)
1425 : type (star_info), pointer :: s
1426 : integer :: k
1427 : real(dp) :: specific_ie
1428 : total_IE = 0
1429 : do k=1,s% nz
1430 : specific_ie = s% energy(k)
1431 : total_IE = total_IE + specific_ie*s% dm(k)
1432 : end do
1433 : end function total_IE
1434 :
1435 :
1436 : subroutine report_energies(s, str)
1437 : type (star_info), pointer :: s
1438 : character (len=*), intent(in) :: str
1439 0 : real(dp) :: KE, IE, PE, Etot
1440 : include 'formats'
1441 :
1442 : return
1443 :
1444 :
1445 : KE = total_KE(s)
1446 : IE = total_IE(s)
1447 : PE = total_PE(s)
1448 : Etot = KE + IE + PE
1449 : write(*,1) 'Etot, KE, IE, PE ' // trim(str), Etot, KE, IE, PE
1450 : end subroutine report_energies
1451 :
1452 :
1453 : end module adjust_mesh_split_merge
1454 :
1455 :
|
