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 momentum 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 :
379 0 : if (s% split_amr_ignore_core_cells .and. &
380 : s%lnT(k)/ln10 >= s% split_amr_logT_for_ignore_core_cells) then
381 : cycle
382 : end if
383 :
384 : ! first check for cells that are too big and need to be split
385 0 : oversize_ratio = dx_actual/dx_baseline
386 0 : if (TooBig < oversize_ratio .and. s% dq(k) > 5d0*dq_min) then
387 0 : if (k < nz .or. s% split_merge_amr_okay_to_split_nz) then
388 0 : if (k > 1 .or. s% split_merge_amr_okay_to_split_1) then
389 0 : TooBig = oversize_ratio; iTooBig = k
390 : end if
391 : end if
392 : end if
393 :
394 : ! next check for cells that are too small and need to be merged
395 :
396 : ! surface cells
397 0 : if (s% merge_amr_ignore_surface_cells .and. &
398 : k<=s% merge_amr_k_for_ignore_surface_cells) cycle
399 :
400 : ! core cells
401 0 : if (s% merge_amr_ignore_core_cells .and. &
402 : s%lnT(k)/ln10>= s% merge_amr_logT_for_ignore_core_cells) cycle
403 :
404 0 : if (abs(dx_actual)>0d0) then
405 0 : undersize_ratio = max(dx_baseline/dx_actual, dq_min/s% dq(k))
406 : else
407 0 : undersize_ratio = dq_min/s% dq(k)
408 : end if
409 :
410 0 : if (s% merge_amr_max_abs_du_div_cs >= 0d0) then
411 0 : call check_merge_limits
412 0 : else if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then
413 0 : TooSmall = undersize_ratio; iTooSmall = k
414 : end if
415 :
416 : end do
417 :
418 :
419 : contains
420 :
421 0 : subroutine check_merge_limits
422 : ! Pablo's additions to modify when merge
423 : ! merge_amr_max_abs_du_div_cs
424 : ! merge_amr_du_div_cs_limit_only_for_compression
425 : ! merge_amr_inhibit_at_jumps
426 :
427 0 : du_div_cs_limit_flag = .false.
428 :
429 0 : if (.not. s% merge_amr_du_div_cs_limit_only_for_compression) then
430 0 : du_div_cs_limit_flag = .true.
431 0 : else if (associated(v)) then
432 0 : if (k < nz) then
433 0 : if (v(k+1)*pow2(r_for_v(k+1)) > v(k)*pow2(r_for_v(k))) then
434 0 : du_div_cs_limit_flag = .true.
435 : end if
436 : end if
437 0 : if (.not. du_div_cs_limit_flag .and. k > 1) then
438 0 : if (v(k)*pow2(r_for_v(k)) > v(k-1)*pow2(r_for_v(k-1))) then
439 0 : du_div_cs_limit_flag = .true.
440 : end if
441 : end if
442 : end if
443 :
444 0 : if (du_div_cs_limit_flag .and. associated(v)) then
445 0 : if (k == 1) then
446 0 : abs_du_div_cs = abs(v(k) - v(k+1))/s% csound(k)
447 0 : else if (k == nz) then
448 0 : abs_du_div_cs = abs(v(nz-1) - v(nz))/s% csound(nz)
449 : else
450 : abs_du_div_cs = max(abs(v(k) - v(k+1)), &
451 0 : abs(v(k) - v(k-1)))/s% csound(k)
452 : end if
453 : else
454 0 : abs_du_div_cs = 0d0
455 : end if
456 :
457 0 : if (du_div_cs_limit_flag) then
458 0 : if (s% merge_amr_inhibit_at_jumps) then
459 : ! reduce undersize_ratio for large jumps
460 : ! i.e. large jumps inhibit merges but don't prohibit completely
461 0 : if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) &
462 : undersize_ratio = undersize_ratio * &
463 0 : s% merge_amr_max_abs_du_div_cs/abs_du_div_cs
464 0 : if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then ! switch
465 0 : TooSmall = undersize_ratio; iTooSmall = k
466 : end if
467 : else if (TooSmall < undersize_ratio .and. &
468 0 : abs_du_div_cs <= s% merge_amr_max_abs_du_div_cs .and. &
469 : s% dq(k) < dq_max/5d0) then
470 0 : TooSmall = undersize_ratio; iTooSmall = k
471 : end if
472 : else
473 0 : if (TooSmall < undersize_ratio .and. s% dq(k) < dq_max/5d0) then
474 0 : TooSmall = undersize_ratio; iTooSmall = k
475 : end if
476 : end if
477 0 : end subroutine check_merge_limits
478 :
479 :
480 : end subroutine biggest_smallest
481 :
482 :
483 0 : subroutine do_merge(s, i_merge, species, new_xa, ierr)
484 : use mesh_adjust, only: set_lnT_for_energy
485 : use star_utils, only: set_rmid
486 : type (star_info), pointer :: s
487 : integer, intent(in) :: i_merge, species
488 : real(dp), intent(inout) :: new_xa(species)
489 : integer, intent(out) :: ierr
490 : logical :: merge_center
491 : integer :: i, ip, i0, im, q, nz, qi_max, qim_max
492 : real(dp) :: &
493 0 : drR, drL, v, &
494 0 : dm, dm_i, dm_ip, star_PE0, star_PE1, &
495 0 : cell_ie, cell_etrb, &
496 0 : Esum_i, KE_i, PE_i, IE_i, Etrb_i, &
497 0 : Esum_ip, KE_ip, PE_ip, IE_ip, Etrb_ip, &
498 0 : KE, &
499 0 : j_rot_new, j_rot_p1_new, J_old, &
500 0 : dmbar_old, dmbar_p1_old, dmbar_p2_old, &
501 0 : dmbar_new, dmbar_p1_new
502 : include 'formats'
503 :
504 0 : ierr = 0
505 0 : s% need_to_setvars = .true.
506 0 : star_PE0 = get_star_PE(s)
507 0 : nz = s% nz
508 :
509 0 : i = i_merge
510 :
511 0 : s% num_hydro_merges = s% num_hydro_merges+1
512 0 : if (i > 1 .and. i < s% nz) then
513 : ! don't merge across change in most abundance species
514 0 : qi_max = maxloc(s% xa(1:species,i), dim=1)
515 0 : qim_max = maxloc(s% xa(1:species,i-1), dim=1)
516 0 : if (qi_max == qim_max) then ! merge with smaller neighbor
517 0 : if (i+1 == nz) then
518 0 : drL = s% r(nz) - s% R_center
519 : else
520 0 : drL = s% r(i) - s% r(i+1)
521 : end if
522 0 : drR = s% r(i-1) - s% r(i)
523 0 : if (drR < drL) i = i-1
524 : ! else i-1 has different most abundant species,
525 : ! so don't consider merging with it
526 : end if
527 : end if
528 :
529 0 : merge_center = (i == nz)
530 0 : if (merge_center) i = i-1
531 0 : ip = i+1
532 0 : if (s% split_merge_amr_avoid_repeated_remesh .and. &
533 : (s% amr_split_merge_has_undergone_remesh(i) .or. &
534 : s% amr_split_merge_has_undergone_remesh(ip))) then
535 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
536 0 : s% amr_split_merge_has_undergone_remesh(ip) = .true.
537 :
538 0 : return
539 : end if
540 :
541 : ! merge contents of zone ip into zone i; remove zone ip
542 : ! get energies for i and ip before merge
543 0 : call get_cell_energies(s, i, Esum_i, KE_i, PE_i, IE_i, Etrb_i)
544 0 : call get_cell_energies(s, ip, Esum_ip, KE_ip, PE_ip, IE_ip, Etrb_ip)
545 :
546 0 : if (s% rotation_flag) then
547 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
548 0 : j_rot_new = 0d0
549 0 : j_rot_p1_new = 0d0
550 0 : if (i==1) then
551 0 : dmbar_old = 0.5d0*s% dm(i)
552 : else
553 0 : dmbar_old = 0.5d0*(s% dm(i)+s% dm(i-1))
554 : end if
555 0 : if (ip /= nz) then
556 0 : dmbar_p1_old = 0.5d0*(s% dm(ip)+s% dm(ip-1))
557 0 : if (ip+1 /= nz) then
558 0 : dmbar_p2_old = 0.5d0*(s% dm(ip+1)+s% dm(ip))
559 : else
560 0 : dmbar_p2_old = s% dm(nz) + 0.5d0*s% dm(nz-1)
561 : end if
562 : ! before merge we have (for i/=nz)
563 : ! dmbar_old + dmbar_p1_old = 0.5*(m(i-1)+m(i))-0.5*(m(i+1)+m(i+2))
564 : ! after merge we have the merged dm_bar_new
565 : ! dmbar_new = 0.5*(m(i-1)+m(i)) - 0.5*(m(i)+m(i+2)) = 0.5*(m(i-1)-m(i+2))
566 : ! where m is the old mass coordinate. From this we have dmbar_new < dmbar_old,
567 : ! dmbar_old + dmbar_p1_old = dmbar_new + 0.5*(m(i) - m(i+1)) = dmbar_new + 0.5*dm(i)
568 : ! this last expression also holds if i=nz
569 0 : dmbar_new = dmbar_old + dmbar_p1_old - 0.5d0*s% dm(i)
570 : ! this difference in mass goes to the dm_bar below
571 0 : dmbar_p1_new = dmbar_p2_old + 0.5d0*s% dm(i)
572 : ! for the merged cell we take the specific angular momentum of both merged cells
573 0 : J_old = dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip)
574 0 : j_rot_new = J_old/(dmbar_old + dmbar_p1_old)
575 : ! and the j_rot of the cell downwards is adjusted to preserve angular momentum
576 0 : J_old = J_old + dmbar_p2_old*s% j_rot(ip+1)
577 0 : j_rot_p1_new = (J_old - j_rot_new*dmbar_new)/dmbar_p1_new
578 : ! which satisfies:
579 : ! J_old = j_rot_p1_new*dmbar_p1_new + j_rot_new*dmbar_new
580 : else
581 0 : dmbar_p1_old = s% dm(nz) + 0.5d0*s% dm(nz-1)
582 : ! simple case, dmbar_new is equal to dmbar_old plus dmbar_p1_old
583 : ! no need to adjust j_rot of another cell downwards
584 0 : dmbar_new = dmbar_old + dmbar_p1_old
585 0 : j_rot_new = (dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip))/dmbar_new
586 0 : j_rot_p1_new = 0d0
587 : !write(*,*) "check center", i, ip, j_rot_new, j_rot_p1_new
588 : end if
589 : end if
590 :
591 0 : dm_i = s% dm(i)
592 0 : dm_ip = s% dm(ip)
593 0 : dm = dm_i + dm_ip
594 0 : s% dm(i) = dm
595 0 : s% dq(i) = dm/s% xmstar
596 :
597 0 : if (s% RTI_flag) &
598 0 : s% alpha_RTI(i) = (dm_i*s% alpha_RTI(i) + dm_ip*s% alpha_RTI(ip))/dm
599 0 : do q=1,species
600 0 : s% xa(q,i) = (dm_i*s% xa(q,i) + dm_ip*s% xa(q,ip))/dm
601 : end do
602 :
603 0 : KE = KE_i + KE_ip
604 0 : if (s% u_flag) then
605 0 : v = sqrt(KE/(0.5d0*dm))
606 0 : if (s% u(i) < 0d0) v = -v
607 0 : s% u(i) = v
608 : else if (s% v_flag) then
609 : ! there's no good solution for this.
610 : ! so just leave s% v(i) unchanged.
611 : end if
612 :
613 0 : cell_ie = IE_i + IE_ip
614 0 : s% energy(i) = cell_ie/dm
615 :
616 0 : if (s% RSP2_flag) then
617 0 : cell_etrb = Etrb_i + Etrb_ip
618 0 : s% w(i) = sqrt(cell_etrb/dm)
619 : end if
620 :
621 0 : if (s% rotation_flag) then
622 0 : s% j_rot(i) = j_rot_new
623 0 : if (ip/=nz) then
624 0 : s% j_rot(ip+1) = j_rot_p1_new
625 : end if
626 : end if
627 :
628 : ! shift ip+1...nz down by 1 to ip...nz-1
629 0 : do i0 = ip+1, nz
630 0 : im = i0-1
631 : !write(*,3) 'shift to im from i0', im, i0
632 0 : do q=1,species
633 0 : s% xa(q,im) = s% xa(q,i0)
634 : end do
635 0 : do q=1,s% nvar_hydro
636 0 : s% xh(q,im) = s% xh(q,i0)
637 : end do
638 0 : s% r(im) = s% r(i0)
639 0 : s% rmid(im) = s% rmid(i0)
640 0 : s% dm(im) = s% dm(i0)
641 0 : s% m(im) = s% m(i0)
642 0 : s% dq(im) = s% dq(i0)
643 0 : s% q(im) = s% q(i0)
644 0 : if (s% u_flag) then
645 0 : s% u(im) = s% u(i0)
646 0 : else if (s% v_flag) then
647 0 : s% v(im) = s% v(i0)
648 : end if
649 0 : if (s% RSP2_flag) then
650 0 : s% w(im) = s% w(i0)
651 0 : s% Hp_face(im) = s% Hp_face(i0)
652 : end if
653 0 : s% energy(im) = s% energy(i0)
654 0 : s% dPdr_dRhodr_info(im) = s% dPdr_dRhodr_info(i0)
655 0 : s% cgrav(im) = s% cgrav(i0)
656 0 : s% alpha_mlt(im) = s% alpha_mlt(i0)
657 0 : s% lnT(im) = s% lnT(i0)
658 0 : s% D_mix(im) = s% D_mix(i0)
659 0 : s% mlt_vc(im) = s% mlt_vc(i0)
660 0 : s% csound(im) = s% csound(i0)
661 0 : s% tau(im) = s% tau(i0)
662 0 : s% opacity(im) = s% opacity(i0)
663 0 : s% amr_split_merge_has_undergone_remesh(im) = s% amr_split_merge_has_undergone_remesh(i0)
664 0 : if (s% RTI_flag) s% alpha_RTI(im) = s% alpha_RTI(i0)
665 0 : if (s% rotation_flag) s% j_rot(im) = s% j_rot(i0)
666 : end do
667 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
668 :
669 0 : nz = nz - 1
670 0 : s% nz = nz
671 :
672 0 : if (s% u_flag) then
673 0 : s% xh(s% i_u,i) = s% u(i)
674 0 : else if (s% v_flag) then
675 0 : s% xh(s% i_v,i) = s% v(i)
676 : end if
677 :
678 0 : if (s% RTI_flag) s% xh(s% i_alpha_RTI,i) = s% alpha_RTI(i)
679 :
680 0 : if (s% RSP2_flag) then
681 0 : s% xh(s% i_w,i) = s% w(i)
682 0 : s% xh(s% i_Hp,i) = s% Hp_face(i)
683 : end if
684 :
685 : ! do this after move cells since need new r(ip) to calc new rho(i).
686 0 : call update_xh_eos_and_kap(s,i,species,new_xa,ierr)
687 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_merge')
688 :
689 0 : s% rmid_start(i) = -1
690 0 : call set_rmid(s, i, i, ierr)
691 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_merge')
692 :
693 0 : star_PE1 = get_star_PE(s)
694 0 : call revise_star_radius(s, star_PE0, star_PE1)
695 :
696 0 : end subroutine do_merge
697 :
698 :
699 0 : subroutine revise_star_radius(s, star_PE0, star_PE1)
700 0 : use star_utils, only: store_r_in_xh, get_lnR_from_xh
701 : type (star_info), pointer :: s
702 : real(dp), intent(in) :: star_PE0, star_PE1
703 : integer :: k
704 : real(dp) :: frac
705 : include 'formats'
706 0 : if (star_PE1 == 0d0 .or. star_PE0 == star_PE1) return
707 0 : frac = star_PE1/star_PE0
708 0 : if (s% model_number == -6918) write(*,1) 'frac', frac
709 0 : if (s% model_number == -6918) write(*,1) 'star_PE0', star_PE0
710 0 : if (s% model_number == -6918) write(*,1) 'star_PE1', star_PE1
711 0 : do k=1,s% nz
712 0 : s% r(k) = s% r(k)*frac
713 0 : if (s% model_number == -6918) write(*,2) 's% r(k)', k, s% r(k)
714 0 : call store_r_in_xh(s, k, s% r(k))
715 0 : s% lnR(k) = get_lnR_from_xh(s,k)
716 : end do
717 0 : s% r_center = frac*s% r_center
718 0 : end subroutine revise_star_radius
719 :
720 :
721 0 : real(dp) function get_star_PE(s) result(totPE)
722 : type (star_info), pointer :: s
723 : integer :: k
724 0 : real(dp) :: PE, rL, rC, dm, mC
725 0 : totPE = 0d0
726 0 : do k=1,s% nz
727 0 : if (k == s% nz) then
728 0 : rL = s% R_center
729 : else
730 0 : rL = s% r(k+1)
731 : end if
732 0 : rC = 0.5d0*(rL + s% r(k))
733 0 : dm = s% dm(k)
734 0 : mC = s% m(k) - 0.5d0*dm
735 0 : PE = -s% cgrav(k)*mC*dm/rC
736 0 : totPE = totPE + PE
737 : end do
738 0 : end function get_star_PE
739 :
740 :
741 0 : subroutine get_cell_energies(s, k, Etot, KE, PE, IE, Etrb)
742 : type (star_info), pointer :: s
743 : integer, intent(in) :: k
744 : real(dp), intent(out) :: Etot, KE, PE, IE, Etrb
745 0 : real(dp) :: dm, mC, v0, v1, rL, rC
746 : include 'formats'
747 0 : dm = s% dm(k)
748 0 : if (s% u_flag) then
749 0 : KE = 0.5d0*dm*s% u(k)**2
750 0 : else if (s% v_flag) then
751 0 : v0 = s% v(k)
752 0 : if (k < s% nz) then
753 0 : v1 = s% v(k+1)
754 : else
755 0 : v1 = s% v_center
756 : end if
757 0 : KE = 0.25d0*dm*(v0**2 + v1**2)
758 : else
759 0 : KE = 0d0
760 : end if
761 0 : IE = s% energy(k)*dm
762 0 : if (s% RSP2_flag) then
763 0 : Etrb = pow2(s% w(k))*dm
764 : else
765 0 : Etrb = 0d0
766 : end if
767 0 : if (s% cgrav(k) == 0) then
768 0 : PE = 0
769 : else
770 0 : if (k == s% nz) then
771 0 : rL = s% R_center
772 : else
773 0 : rL = s% r(k+1)
774 : end if
775 0 : rC = 0.5d0*(rL + s% r(k))
776 0 : mC = s% m(k) - 0.5d0*dm
777 0 : PE = -s% cgrav(k)*mC*dm/rC
778 : end if
779 0 : Etot = IE + KE + PE
780 : if (is_bad(Etot + IE + KE + PE) .or. &
781 0 : IE <= 0 .or. KE < 0) then
782 0 : write(*,2) 'nz', s% nz
783 0 : write(*,2) 'Etot', k, Etot
784 0 : write(*,2) 'IE', k, IE
785 0 : write(*,2) 'KE', k, KE
786 0 : write(*,2) 'PE', k, PE
787 0 : call mesa_error(__FILE__,__LINE__,'get_cell_energies')
788 : end if
789 0 : end subroutine get_cell_energies
790 :
791 :
792 0 : subroutine split1_non_negative( &
793 : val, grad, dr, dV, dVR, dVL, new_valL, new_valR)
794 : real(dp), intent(in) :: val, grad, dr, dV, dVR, dVL
795 : real(dp), intent(out) :: new_valL, new_valR
796 0 : real(dp) :: xL, xR, xML, xMR, xM, f
797 : include 'formats'
798 0 : new_valR = val
799 0 : new_valL = val
800 0 : if (val < 1d-99 .or. grad == 0d0) return
801 0 : xL = max(0d0, min(1d0, val - grad*dr/4))
802 0 : xR = max(0d0, min(1d0, val + grad*dr/4))
803 0 : xML = xL*dVL
804 0 : xMR = xR*dVR
805 0 : xM = val*dV
806 0 : if (xM < 1d-99 .or. xML + xMR < 1d-99) then
807 : new_valR = val
808 : new_valL = val
809 : else
810 0 : f = xM/(xML + xMR)
811 0 : new_valR = f*xR
812 0 : new_valL = f*xL
813 : end if
814 : end subroutine split1_non_negative
815 :
816 :
817 0 : subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr)
818 : use alloc, only: reallocate_star_info_arrays
819 : use star_utils, only: set_rmid, store_r_in_xh
820 : type (star_info), pointer :: s
821 : integer, intent(in) :: i_split, species
822 : real(dp) :: tau_center, grad_xa(species), new_xa(species)
823 : integer, intent(out) :: ierr
824 : integer :: i, ip, j, jp, q, nz, nz_old, &
825 : iR, iC, iL
826 : real(dp) :: &
827 : cell_Esum_old, cell_KE_old, cell_PE_old, cell_IE_old, cell_Etrb_old, &
828 0 : rho_RR, rho_iR, rR, rL, dr, dr_old, rC, dV, dVR, dVL, dM, dML, dMR, rho, &
829 0 : v, v2, energy, v2_R, energy_R, rho_R, v2_C, energy_C, rho_C, v2_L, energy_L, rho_L, &
830 0 : dLeft, dRght, dCntr, grad_rho, grad_energy, grad_v2, &
831 0 : sumx, sumxp, new_xaL, new_xaR, star_PE0, star_PE1, &
832 0 : grad_alpha, f, new_alphaL, new_alphaR, v_R, v_C, v_L, min_dm, &
833 0 : mlt_vcL, mlt_vcR, tauL, tauR, etrb, etrb_L, etrb_C, etrb_R, grad_etrb, &
834 0 : j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old
835 : logical :: done, use_new_grad_rho
836 : include 'formats'
837 :
838 0 : ierr = 0
839 0 : star_PE0 = get_star_PE(s)
840 0 : s% need_to_setvars = .true.
841 0 : nz = s% nz
842 0 : s% num_hydro_splits = s% num_hydro_splits + 1
843 0 : done = .false.
844 0 : nz_old = nz
845 :
846 0 : i = i_split
847 0 : ip = i+1
848 :
849 : call get_cell_energies( &
850 0 : s, i, cell_Esum_old, cell_KE_old, cell_PE_old, cell_IE_old, cell_Etrb_old)
851 :
852 0 : if (s% rotation_flag .and. i<nz) then
853 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
854 0 : if (i==1) then
855 0 : dmbar_old = 0.5d0*s% dm(i)
856 : else
857 0 : dmbar_old = 0.5d0*(s% dm(i)+s% dm(i-1))
858 : end if
859 0 : if (ip /= nz) then
860 0 : dmbar_p1_old = 0.5d0*(s% dm(ip)+s% dm(ip-1))
861 : else
862 0 : dmbar_p1_old = s% dm(ip)+0.5d0*s% dm(ip-1)
863 : end if
864 :
865 : ! We need to conserve the angular momentum in these two cells
866 0 : J_old = dmbar_old*s% j_rot(i) + dmbar_p1_old*s% j_rot(ip)
867 : end if
868 :
869 0 : rR = s% r(i)
870 0 : mlt_vcR = s% mlt_vc(i)
871 0 : if (i == nz) then
872 0 : rL = s% R_center
873 0 : mlt_vcL = 0d0
874 : tauL = tau_center
875 : else
876 0 : rL = s% r(ip)
877 0 : mlt_vcL = s% mlt_vc(ip)
878 : tauL = s% tau(ip)
879 : end if
880 :
881 0 : tauR = s% tau(i)
882 0 : if (i == nz) then
883 0 : tauL = tau_center
884 : else
885 0 : tauL = s% tau(ip)
886 : end if
887 0 : if (is_bad(tauL+tauR)) then
888 0 : !$omp critical (adjust_mesh_split_merge_crit1)
889 0 : write(*,2) 'tauL', ip, tauL
890 0 : write(*,2) 'tauR', i, tauR
891 0 : write(*,2) 'nz', nz
892 0 : call mesa_error(__FILE__,__LINE__,'do_split')
893 : !$omp end critical (adjust_mesh_split_merge_crit1)
894 : end if
895 :
896 0 : dr = rR - rL
897 0 : dr_old = dr
898 0 : rC = 0.5d0*(rR + rL) ! split at center by radius
899 :
900 0 : dV = four_thirds_pi*(rR*rR*rR - rL*rL*rL)
901 0 : dVR = four_thirds_pi*(rR*rR*rR - rC*rC*rC)
902 0 : dVL = dV - dVR
903 :
904 0 : dm = s% dm(i)
905 0 : rho = dm/dV
906 :
907 0 : if (s% u_flag) then
908 0 : v = s% u(i)
909 0 : v2 = v*v
910 : else ! not used
911 : v = 0
912 : v2 = 0
913 : end if
914 :
915 0 : energy = s% energy(i)
916 0 : if (s% RSP2_flag) etrb = pow2(s% w(i))
917 :
918 : ! use iR, iC, and iL for getting values to determine slopes
919 0 : if (i > 1 .and. i < nz_old) then
920 0 : iR = i-1
921 0 : iC = i
922 0 : iL = i+1
923 0 : else if (i == 1) then
924 0 : iR = 1
925 0 : iC = 2
926 0 : iL = 3
927 : else ! i == nz_old
928 0 : iR = nz_old-2
929 0 : iC = nz_old-1
930 0 : iL = nz_old
931 : end if
932 :
933 0 : energy_R = s% energy(iR)
934 0 : rho_R = s% dm(iR)/get_dV(s,iR)
935 :
936 0 : energy_C = s% energy(iC)
937 0 : rho_C = s% dm(iC)/get_dV(s,iC)
938 :
939 0 : energy_L = s% energy(iL)
940 0 : rho_L = s% dm(iL)/get_dV(s,iL)
941 :
942 : ! get gradients before move cell contents
943 :
944 0 : if (iL == nz_old) then
945 0 : dLeft = 0.5d0*(s% r(iC) - s% R_center)
946 : else
947 0 : dLeft = 0.5d0*(s% r(iC) - s% r(iL+1))
948 : end if
949 0 : dRght = 0.5d0*(s% r(iR) - s% r(iL))
950 0 : dCntr = dLeft + dRght
951 :
952 0 : if (s% equal_split_density_amr) then ! same density in both parts
953 0 : rho_RR = 0
954 : grad_rho = 0d0
955 : use_new_grad_rho = .false.
956 : else if (.false.) then
957 : rho_RR = s% dm(iR-1)/get_dV(s,iR-1)
958 : grad_rho = 2*(rho_RR - rho_R)/(s% r(iR-1) - s% r(iR+1))
959 : use_new_grad_rho = .true.
960 : else
961 0 : rho_RR = 0
962 0 : grad_rho = get1_grad(rho_L, rho_C, rho_R, dLeft, dCntr, dRght)
963 0 : use_new_grad_rho = .false.
964 : end if
965 :
966 0 : grad_energy = get1_grad(energy_L, energy_C, energy_R, dLeft, dCntr, dRght)
967 :
968 0 : if (s% RTI_flag) grad_alpha = get1_grad( &
969 0 : s% alpha_RTI(iL), s% alpha_RTI(iC), s% alpha_RTI(iR), dLeft, dCntr, dRght)
970 :
971 0 : if (s% RSP2_flag) then
972 0 : etrb_R = pow2(s% w(iR))
973 0 : etrb_C = pow2(s% w(iC))
974 0 : etrb_L = pow2(s% w(iL))
975 0 : grad_etrb = get1_grad(etrb_L, etrb_C, etrb_R, dLeft, dCntr, dRght)
976 : end if
977 :
978 0 : if (s% u_flag) then
979 0 : v_R = s% u(iR)
980 0 : v2_R = v_R*v_R
981 0 : v_C = s% u(iC)
982 0 : v2_C = v_C*v_C
983 0 : v_L = s% u(iL)
984 0 : v2_L = v_L*v_L
985 0 : if ((v_L - v_C)*(v_C - v_R) <= 0) then ! not strictly monotonic velocities
986 : grad_v2 = 0d0
987 : else
988 0 : grad_v2 = get1_grad(v2_L, v2_C, v2_R, dLeft, dCntr, dRght)
989 : end if
990 0 : else if (s% v_flag) then
991 0 : if (iL == s% nz) then
992 0 : v_L = s% v_center
993 : else
994 0 : v_L = s% v(ip)
995 : end if
996 0 : v2_L = v_L*v_L
997 0 : v_R = s% v(i)
998 0 : v2_R = v_R*v_R
999 : end if
1000 :
1001 0 : do q = 1, species
1002 : grad_xa(q) = get1_grad( &
1003 0 : s% xa(q,iL), s% xa(q,iC), s% xa(q,iR), dLeft, dCntr, dRght)
1004 : end do
1005 :
1006 0 : nz = nz + 1
1007 0 : s% nz = nz
1008 0 : call reallocate_star_info_arrays(s, ierr)
1009 0 : if (ierr /= 0) then
1010 0 : write(*,2) 'reallocate_star_info_arrays ierr', ierr
1011 0 : call mesa_error(__FILE__,__LINE__,'split failed')
1012 : end if
1013 :
1014 0 : if (i_split < nz_old) then ! move i_split..nz-1 to i_split+1..nz
1015 0 : do j=nz-1,i_split,-1
1016 0 : jp = j+1
1017 0 : do q=1,species
1018 0 : s% xa(q,jp) = s% xa(q,j)
1019 : end do
1020 0 : do q=1,s% nvar_hydro
1021 0 : s% xh(q,jp) = s% xh(q,j)
1022 : end do
1023 0 : s% r(jp) = s% r(j)
1024 0 : s% rmid(jp) = s% rmid(j)
1025 0 : s% dm(jp) = s% dm(j)
1026 0 : s% m(jp) = s% m(j)
1027 0 : s% dq(jp) = s% dq(j)
1028 0 : s% q(jp) = s% q(j)
1029 0 : if (s% u_flag) then
1030 0 : s% u(jp) = s% u(j)
1031 0 : else if (s% v_flag) then
1032 0 : s% v(jp) = s% v(j)
1033 : end if
1034 0 : s% energy(jp) = s% energy(j)
1035 0 : s% dPdr_dRhodr_info(jp) = s% dPdr_dRhodr_info(j)
1036 0 : s% cgrav(jp) = s% cgrav(j)
1037 0 : s% alpha_mlt(jp) = s% alpha_mlt(j)
1038 0 : s% lnT(jp) = s% lnT(j)
1039 0 : s% D_mix(jp) = s% D_mix(j)
1040 0 : s% mlt_vc(jp) = s% mlt_vc(j)
1041 0 : s% csound(jp) = s% csound(j)
1042 0 : s% tau(jp) = s% tau(j)
1043 0 : s% opacity(jp) = s% opacity(j)
1044 0 : s% amr_split_merge_has_undergone_remesh(jp) = s% amr_split_merge_has_undergone_remesh(j)
1045 0 : if (s% RTI_flag) s% alpha_RTI(jp) = s% alpha_RTI(j)
1046 0 : if (s% rotation_flag) s% j_rot(jp) = s% j_rot(j)
1047 : end do
1048 : end if
1049 0 : s% amr_split_merge_has_undergone_remesh(i) = .true.
1050 0 : s% amr_split_merge_has_undergone_remesh(ip) = .true.
1051 :
1052 0 : dM = rho*dV
1053 : if (.not. use_new_grad_rho) then ! do it the old way
1054 :
1055 0 : rho_L = rho - grad_rho*dr/4
1056 0 : rho_R = rho + grad_rho*dr/4
1057 : if (grad_rho == 0d0 .or. &
1058 0 : rho_L <= 0d0 .or. &
1059 : rho_R <= 0d0) then
1060 0 : rho_R = rho
1061 0 : rho_L = rho
1062 0 : dML = rho_L*dVL
1063 0 : dMR = dM - dML ! should = rhoR*dVR
1064 : else
1065 0 : dML = rho_L*dVL
1066 0 : dMR = rho_R*dVR
1067 0 : f = dM/(dML + dMR)
1068 0 : rho_L = f*rho_L
1069 0 : rho_R = f*rho_R
1070 0 : dML = rho_L*dVL
1071 0 : dMR = rho_R*dVR
1072 0 : if (abs(dML + dMR - dM) > 1d-14*dM) then
1073 0 : write(*,2) '(dML + dMR - dM)/dM', i, (dML + dMR - dM)/dM
1074 0 : call mesa_error(__FILE__,__LINE__,'split')
1075 : end if
1076 0 : dMR = dM - dML
1077 : end if
1078 :
1079 : else
1080 :
1081 : ! at this point, rho_R is the density of the cell iR
1082 : ! we are about to redefine it as the density of the right 1/2 of the split
1083 : ! similarly for rho_L
1084 : rho_iR = rho_R
1085 : dR = -(dRght/4 + (s% r(iR) - s% r(iC))/2)
1086 : rho_R = rho_iR + grad_rho*dR
1087 : dMR = rho_R*dVR
1088 : dML = dM - dMR
1089 : rho_L = dML/dVL
1090 : if (rho_R <= 1d-16 .or. rho_L <= 1d-16) then
1091 : !$omp critical (adjust_mesh_split_merge_crit2)
1092 : write(*,'(A)')
1093 : write(*,2) 'nz', nz
1094 : write(*,2) 'rho_RR', iR-1, rho_RR
1095 : write(*,2) 'rho_iR', iR, rho_iR
1096 : write(*,2) 'rho', iC, rho
1097 : write(*,2) 's% r(iR-1)', iR-1, s% r(iR-1)
1098 : write(*,2) 's% r(iR)', iR, s% r(iR)
1099 : write(*,2) 's% r(iR+1)', iR+1, s% r(iR+1)
1100 : write(*,2) 'rho_RR', iR-1, rho_RR
1101 : write(*,2) 'rho_iR', iR, rho_iR
1102 : write(*,2) 'rho_RR - rho_iR', iR, rho_RR - rho_iR
1103 : write(*,2) 'dr for right', iR, (s% r(iR-1) - s% r(iR+1))/2
1104 : write(*,2) 'grad_rho', iR, grad_rho
1105 : write(*,'(A)')
1106 : write(*,2) 's% r(iL)', iL, s% r(iL)
1107 : write(*,2) 'dR', iC, dR
1108 : write(*,2) 'dRho', iC, grad_rho*dR
1109 : write(*,2) 'rho_R', iC, rho_R
1110 : write(*,2) 'rho_L', iC, rho_L
1111 : write(*,'(A)')
1112 : call mesa_error(__FILE__,__LINE__,'failed in do_split extrapolation of density from above')
1113 : !$omp end critical (adjust_mesh_split_merge_crit2)
1114 : end if
1115 :
1116 : end if
1117 :
1118 0 : min_dm = s% split_merge_amr_dq_min*s% xmstar
1119 0 : if (dML < min_dm .or. dMR < min_dm) then
1120 0 : rho_R = rho
1121 0 : rho_L = rho
1122 0 : dM = rho*dV
1123 0 : dML = rho_L*dVL
1124 0 : dMR = dM - dML ! should = rhoR*dVR
1125 : end if
1126 :
1127 0 : energy_R = energy + grad_energy*dr/4
1128 0 : energy_L = (dm*energy - dmR*energy_R)/dmL
1129 0 : if (energy_R < 0d0 .or. energy_L < 0d0) then
1130 0 : energy_R = energy
1131 0 : energy_L = energy
1132 : end if
1133 :
1134 0 : s% energy(i) = energy_R
1135 0 : s% energy(ip) = energy_L
1136 :
1137 0 : if (s% RSP2_flag) then
1138 0 : etrb_R = etrb + grad_etrb*dr/4
1139 0 : etrb_L = (dm*etrb - dmR*etrb_R)/dmL
1140 0 : if (etrb_R < 0d0 .or. etrb_L < 0d0) then
1141 0 : etrb_R = etrb
1142 0 : etrb_L = etrb
1143 : end if
1144 0 : s% w(i) = sqrt(max(0d0,etrb_R))
1145 0 : s% w(ip) = sqrt(max(0d0,etrb_L))
1146 : end if
1147 :
1148 0 : if (s% u_flag) then
1149 0 : v2_R = v2 + grad_v2*dr/4
1150 0 : v2_L = (dm*v2 - dmR*v2_R)/dmL
1151 0 : if (v2_R < 0d0 .or. v2_L < 0d0) then
1152 0 : v2_R = v2
1153 0 : v2_L = v2
1154 : end if
1155 0 : s% u(i) = sqrt(v2_R)
1156 0 : s% u(ip) = sqrt(v2_L)
1157 0 : if (v < 0d0) then
1158 0 : s% u(i) = -s% u(i)
1159 0 : s% u(ip) = -s% u(ip)
1160 : end if
1161 0 : else if (s% v_flag) then ! just make a rough approximation.
1162 0 : s% v(ip) = sqrt(0.5d0*(v2_L + v2_R))
1163 : end if
1164 :
1165 0 : if (s% RTI_flag) then ! set new alpha
1166 0 : if (i == 1) then
1167 0 : s% alpha_RTI(ip) = s% alpha_RTI(i)
1168 : else
1169 : call split1_non_negative( &
1170 : s% alpha_RTI(i), grad_alpha, &
1171 0 : dr, dV, dVR, dVL, new_alphaL, new_alphaR)
1172 0 : s% alpha_RTI(i) = new_alphaR
1173 0 : s% alpha_RTI(ip) = new_alphaL
1174 : end if
1175 0 : s% dPdr_dRhodr_info(ip) = s% dPdr_dRhodr_info(i)
1176 : end if
1177 :
1178 0 : if (i == 1) then
1179 0 : s% mlt_vc(ip) = s% mlt_vc(i)
1180 : else
1181 0 : s% mlt_vc(ip) = (mlt_vcL*dML + mlt_vcR*dMR)/dM
1182 : end if
1183 :
1184 0 : s% tau(ip) = tauR + (tauL - tauR)*dMR/dM
1185 0 : if (is_bad(s% tau(ip))) then
1186 0 : write(*,2) 'tau', ip, s% tau(ip), tauL, tauR, dMR/dM
1187 0 : call mesa_error(__FILE__,__LINE__,'split')
1188 : end if
1189 :
1190 0 : if (i == 1) then
1191 0 : do q=1,species
1192 0 : s% xa(q,ip) = s% xa(q,i)
1193 : end do
1194 : else ! split mass fractions
1195 : ! check that sum of grads for mass fractions is 0
1196 0 : sumx = sum(grad_xa)
1197 0 : if (sumx > 0) then ! reduce largest grad
1198 0 : j = maxloc(grad_xa, dim=1)
1199 : else ! increase smallest grad
1200 0 : j = minloc(grad_xa, dim=1)
1201 : end if
1202 0 : grad_xa(j) = 0
1203 0 : grad_xa(j) = -sum(grad_xa)
1204 : ! set new mass fractions
1205 0 : do q = 1, species
1206 : call split1_non_negative( &
1207 : s% xa(q,i), grad_xa(q), &
1208 0 : dr, dV, dVR, dVL, new_xaL, new_xaR)
1209 0 : s% xa(q,i) = new_xaR
1210 0 : s% xa(q,ip) = new_xaL
1211 : end do
1212 : !check mass fractions >= 0 and <= 1 and sum to 1.0
1213 0 : do q = 1, species
1214 0 : s% xa(q,i) = min(1d0, max(0d0, s% xa(q,i)))
1215 0 : s% xa(q,ip) = min(1d0, max(0d0, s% xa(q,ip)))
1216 : end do
1217 0 : sumx = sum(s% xa(1:species,i))
1218 0 : sumxp = sum(s% xa(1:species,ip))
1219 0 : do q = 1, species
1220 0 : s% xa(q,i) = s% xa(q,i)/sumx
1221 0 : s% xa(q,ip) = s% xa(q,ip)/sumxp
1222 : end do
1223 : end if
1224 :
1225 0 : if (s% u_flag) s% u_face_ad(ip)%val = 0.5d0*(s% u(i) + s% u(ip))
1226 : ! just for setting u_face_start so don't need partials
1227 :
1228 : ! r, q, m, u_face unchanged for face i
1229 0 : dVR = dV - dVL
1230 0 : dmR = s% dm(i) - dmL
1231 0 : s% dm(i) = dmR
1232 0 : s% dq(i) = s% dm(i)/s% xmstar
1233 :
1234 0 : s% r(ip) = rC
1235 0 : s% m(ip) = s% m(i) - s% dm(i) ! <<< using new value dm(i)
1236 0 : s% q(ip) = s% q(i) - s% dq(i) ! <<< using new value dq(i)
1237 0 : s% dm(ip) = dmL
1238 0 : s% dq(ip) = s% dm(ip)/s% xmstar
1239 :
1240 0 : s% cgrav(ip) = s% cgrav(i)
1241 0 : s% alpha_mlt(ip) = s% alpha_mlt(i)
1242 0 : s% lnT(ip) = s% lnT(i)
1243 0 : s% T(ip) = s% T(i)
1244 0 : s% D_mix(ip) = s% D_mix(i)
1245 0 : s% mlt_vc(ip) = s% mlt_vc(i)
1246 0 : s% csound(ip) = s% csound(i)
1247 0 : s% opacity(ip) = s% opacity(i)
1248 0 : if (s% RTI_flag) s% alpha_RTI(ip) = s% alpha_RTI(i)
1249 :
1250 0 : if (s% rotation_flag) then
1251 : ! WARNING! this is designed to conserve angular momentum, but not to explicitly conserve energy
1252 0 : j_rot_new = 0d0 ! specific angular momentum for the newly created cell
1253 0 : if (i==nz_old) then
1254 : ! simple case, dm_bar contained in cells i and i+1 after merge is conserved,
1255 : ! so use same j_rot
1256 0 : j_rot_new = s% j_rot(i)
1257 : else
1258 0 : if (i==1) then
1259 0 : dmbar_new = 0.5d0*s% dm(i)
1260 : else
1261 0 : dmbar_new = 0.5d0*(s% dm(i)+s% dm(i-1))
1262 : end if
1263 0 : if (ip /= nz) then
1264 0 : dmbar_p1_new = 0.5d0*(s% dm(ip)+s% dm(ip-1))
1265 : else
1266 0 : dmbar_p1_new = s% dm(ip)+0.5d0*s% dm(ip-1)
1267 : end if
1268 0 : if (ip+1 /= nz) then
1269 0 : dmbar_p2_new = 0.5d0*(s% dm(ip+1)+s% dm(ip))
1270 : else
1271 0 : dmbar_p2_new = s% dm(ip+1)+0.5d0*s% dm(ip)
1272 : end if
1273 :
1274 : ! we keep the same j_rot for cells i and ip (cells i and ip+1 after merge),
1275 : ! we compute the j_rot needed for the new cell to preserve angular momentum.
1276 0 : j_rot_new = (J_old - dmbar_new*s% j_rot(i) - dmbar_p2_new*s% j_rot(ip+1))/(dmbar_p1_new)
1277 : ! this satisfies:
1278 : ! 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
1279 : end if
1280 0 : s% j_rot(ip) = j_rot_new
1281 : end if
1282 :
1283 0 : if (s% i_lum /= 0) then
1284 0 : if (ip < nz_old) then
1285 : s% xh(s% i_lum,ip) = &
1286 0 : 0.5d0*(s% xh(s% i_lum,i) + s% xh(s% i_lum,ip+1))
1287 : else
1288 : s% xh(s% i_lum,ip) = &
1289 0 : 0.5d0*(s% xh(s% i_lum,i) + s% L_center)
1290 : end if
1291 : end if
1292 :
1293 0 : call store_r_in_xh(s, ip, s% r(ip))
1294 0 : if (s% u_flag) then
1295 0 : s% xh(s% i_u,i) = s% u(i)
1296 0 : s% xh(s% i_u,ip) = s% u(ip)
1297 0 : else if (s% v_flag) then
1298 0 : s% xh(s% i_v,i) = s% v(i)
1299 0 : s% xh(s% i_v,ip) = s% v(ip)
1300 : end if
1301 0 : if (s% RTI_flag) then
1302 0 : s% xh(s% i_alpha_RTI,i) = s% alpha_RTI(i)
1303 0 : s% xh(s% i_alpha_RTI,ip) = s% alpha_RTI(ip)
1304 : end if
1305 :
1306 0 : call update_xh_eos_and_kap(s,i,species,new_xa,ierr)
1307 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1308 :
1309 0 : call update_xh_eos_and_kap(s,ip,species,new_xa,ierr)
1310 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1311 :
1312 0 : s% rmid_start(i) = -1
1313 0 : s% rmid_start(ip) = -1
1314 0 : call set_rmid(s, i, ip, ierr)
1315 0 : if (ierr /= 0) return ! call mesa_error(__FILE__,__LINE__,'update_xh_eos_and_kap failed in do_split')
1316 :
1317 0 : star_PE1 = get_star_PE(s)
1318 0 : call revise_star_radius(s, star_PE0, star_PE1)
1319 :
1320 0 : end subroutine do_split
1321 :
1322 :
1323 0 : subroutine update_xh_eos_and_kap(s,i,species,new_xa,ierr)
1324 0 : use mesh_adjust, only: set_lnT_for_energy
1325 : use micro, only: do_kap_for_cell
1326 : use eos_lib, only: eos_gamma_DE_get_PT
1327 : use chem_lib, only: basic_composition_info
1328 : use star_utils, only: store_lnT_in_xh, get_T_and_lnT_from_xh, &
1329 : store_rho_in_xh, get_rho_and_lnd_from_xh
1330 : type (star_info), pointer :: s
1331 : integer, intent(in) :: i, species
1332 : real(dp) :: new_xa(species)
1333 : integer, intent(out) :: ierr
1334 : real(dp) :: rho, logRho, new_lnT, revised_energy
1335 : integer :: q
1336 : include 'formats'
1337 0 : ierr = 0
1338 0 : rho = s% dm(i)/get_dV(s,i)
1339 0 : call store_rho_in_xh(s, i, rho)
1340 0 : call get_rho_and_lnd_from_xh(s, i, s% rho(i), s% lnd(i))
1341 0 : logRho = s% lnd(i)/ln10
1342 0 : do q=1,species
1343 0 : new_xa(q) = s% xa(q,i)
1344 : end do
1345 : call set_lnT_for_energy(s, i, &
1346 : s% net_iso(ih1), s% net_iso(ihe3), s% net_iso(ihe4), &
1347 : species, new_xa, rho, logRho, s% energy(i), s% lnT(i), &
1348 0 : new_lnT, revised_energy, ierr)
1349 0 : if (ierr /= 0) return
1350 0 : call store_lnT_in_xh(s, i, new_lnT)
1351 0 : call get_T_and_lnT_from_xh(s, i, s% T(i), s% lnT(i))
1352 0 : end subroutine update_xh_eos_and_kap
1353 :
1354 :
1355 0 : real(dp) function get_dV(s,k)
1356 : type (star_info), pointer :: s
1357 : integer, intent(in) :: k
1358 0 : if (k == s% nz) then
1359 0 : get_dV = four_thirds_pi*(s% r(k)*s% r(k)*s% r(k) - s% R_center*s% R_center*s% R_center)
1360 : else
1361 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))
1362 : end if
1363 0 : end function get_dV
1364 :
1365 :
1366 0 : real(dp) function minmod(a, b, c) result(m)
1367 : real(dp), intent(in) :: a, b, c
1368 0 : m = a
1369 0 : if (a*b < 0d0) then
1370 0 : m = 0d0
1371 : return
1372 : end if
1373 0 : if (abs(b) < abs(m)) m = b
1374 0 : if (b*c < 0d0) then
1375 0 : m = 0d0
1376 : return
1377 : end if
1378 0 : if (abs(c) < abs(m)) m = c
1379 : end function minmod
1380 :
1381 :
1382 0 : real(dp) function get1_grad(vL, vC, vR, dLeft, dCntr, dRght) &
1383 0 : result(grad)
1384 : real(dp), intent(in) :: vL, vC, vR, dLeft, dCntr, dRght
1385 0 : real(dp) :: sL, sR, sC
1386 : include 'formats'
1387 0 : if (dLeft <= 0d0 .or. dCntr <= 0d0 .or. dRght <= 0d0) then
1388 0 : grad = 0d0
1389 : return
1390 : end if
1391 0 : sL = (vC - vL)/dLeft
1392 0 : sR = (vR - vC)/dRght
1393 0 : sC = (vR - vL)/dCntr
1394 0 : grad = minmod(sL, sC, sR)
1395 0 : end function get1_grad
1396 :
1397 :
1398 : real(dp) function total_KE(s)
1399 : type (star_info), pointer :: s
1400 : integer :: k
1401 : include 'formats'
1402 : total_KE = 0
1403 : if (s% u_flag) then
1404 : do k=1,s% nz
1405 : total_KE = total_KE + 0.5d0*s% dm(k)*s% u(k)**2
1406 : end do
1407 : else if (s% v_flag) then
1408 : do k=1,s% nz-1
1409 : total_KE = total_KE + &
1410 : 0.25d0*s% dm(k)*(s% v(k)**2 + s% v(k+1)**2)
1411 : end do
1412 : k = s% nz
1413 : total_KE = total_KE + &
1414 : 0.25d0*s% dm(k)*(s% v(k)**2 + s% v_center**2)
1415 : end if
1416 : end function total_KE
1417 :
1418 :
1419 : real(dp) function total_PE(s)
1420 : type (star_info), pointer :: s
1421 : integer :: k
1422 : real(dp) :: rL, rR, rC, mC
1423 : total_PE = 0d0
1424 : rR = s% R_center
1425 : do k = s% nz,1,-1
1426 : rL = rR
1427 : rR = s% r(k)
1428 : rC = 0.5d0*(rL + rR)
1429 : mC = s% m(k) - 0.5d0*s% dm(k)
1430 : total_PE = total_PE - s% cgrav(k)*mC*s% dm(k)/rC
1431 : end do
1432 : end function total_PE
1433 :
1434 :
1435 : real(dp) function total_IE(s)
1436 : type (star_info), pointer :: s
1437 : integer :: k
1438 : real(dp) :: specific_ie
1439 : total_IE = 0
1440 : do k=1,s% nz
1441 : specific_ie = s% energy(k)
1442 : total_IE = total_IE + specific_ie*s% dm(k)
1443 : end do
1444 : end function total_IE
1445 :
1446 :
1447 : subroutine report_energies(s, str)
1448 : type (star_info), pointer :: s
1449 : character (len=*), intent(in) :: str
1450 0 : real(dp) :: KE, IE, PE, Etot
1451 : include 'formats'
1452 :
1453 : return
1454 :
1455 :
1456 : KE = total_KE(s)
1457 : IE = total_IE(s)
1458 : PE = total_PE(s)
1459 : Etot = KE + IE + PE
1460 : write(*,1) 'Etot, KE, IE, PE ' // trim(str), Etot, KE, IE, PE
1461 : end subroutine report_energies
1462 :
1463 :
1464 : end module adjust_mesh_split_merge
|
