Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012 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 solve_omega_mix
21 :
22 : use star_private_def
23 : use const_def, only: qp, dp, i8
24 :
25 : implicit none
26 :
27 : private
28 : public :: do_solve_omega_mix
29 :
30 : contains
31 :
32 0 : integer function do_solve_omega_mix(s, dt_total)
33 : use star_utils, only: start_time, update_time, total_angular_momentum
34 : use mix_info, only: update_rotation_mixing_info
35 : use hydro_rotation, only: get_rotation_sigmas, set_omega, set_i_rot
36 :
37 : type (star_info), pointer :: s
38 : real(dp), intent(in) :: dt_total
39 :
40 : integer :: ierr, nz, k, max_iters_per_substep, &
41 : max_iters_total, total_num_iters, num_iters
42 : integer(i8) :: time0
43 : integer :: steps_used, max_steps, min_steps
44 0 : real(qp) :: remaining_time, total_time, time, dt, &
45 0 : J_tot0, J_tot1, max_del, avg_del, &
46 0 : tol_correction_max, tol_correction_norm
47 0 : real(dp) :: total
48 0 : real(dp), pointer, dimension(:) :: am_sig_omega, am_sig_j
49 0 : real(qp), pointer, dimension(:) :: du, d, dl, x, b, bp, vp, xp, dX, X_0, X_1, rhs, del
50 : logical :: okay, recalc_mixing_info_each_substep
51 : logical, parameter :: dbg = .false.
52 :
53 : include 'formats'
54 :
55 0 : do_solve_omega_mix = keep_going
56 :
57 0 : ierr = 0
58 :
59 0 : nz = s% nz
60 0 : total_time = dt_total
61 0 : time = 0
62 0 : max_steps = 20
63 0 : min_steps = 4
64 0 : max_iters_per_substep = 4
65 0 : max_iters_total = 40
66 0 : total_num_iters = 0
67 0 : tol_correction_max = 1d-4
68 0 : tol_correction_norm = 1d-7
69 :
70 : ! update omega for new i_rot and previous j_rot to conserve angular momentum
71 0 : call set_i_rot(s, .false.)
72 0 : call set_omega(s, 'solve_omega_mix')
73 :
74 0 : if (dt_total <= 0d0) return
75 :
76 0 : if (s% doing_timing) call start_time(s, time0, total)
77 :
78 0 : s% extra_jdot(1:nz) = 0
79 0 : s% extra_omegadot(1:nz) = 0
80 :
81 0 : if (s% use_other_torque) then
82 0 : call s% other_torque(s% id, ierr)
83 0 : if (ierr /= 0) then
84 0 : if (s% report_ierr .or. dbg) &
85 0 : write(*, *) 'solve_omega_mix: other_torque returned ierr', ierr
86 0 : return
87 : end if
88 : end if
89 :
90 0 : if (associated(s% binary_other_torque)) then
91 : call s% binary_other_torque(s% id, ierr)
92 0 : if (ierr /= 0) then
93 0 : if (s% report_ierr .or. dbg) &
94 0 : write(*, *) 'solve_omega_mix: binary_other_torque returned ierr', ierr
95 0 : return
96 : end if
97 : end if
98 :
99 : call do_alloc(ierr)
100 0 : if (ierr /= 0) then
101 0 : do_solve_omega_mix = terminate
102 0 : s% termination_code = t_solve_omega_mix
103 0 : s% result_reason = nonzero_ierr
104 0 : if (s% report_ierr) write(*,*) 'allocate failed in do_solve_omega_mix'
105 0 : return
106 : end if
107 :
108 0 : J_tot0 = total_angular_momentum(s)
109 :
110 0 : okay = .true.
111 0 : do k=1,nz
112 0 : if (is_bad_num(s% omega(k)) .or. abs(s% omega(k)) > 1d50) then
113 0 : if (s% stop_for_bad_nums) then
114 0 : write(*,2) 's% omega(k)', k, s% omega(k)
115 0 : call mesa_error(__FILE__,__LINE__,'solve omega')
116 : end if
117 : okay = .false.
118 : exit
119 : end if
120 : end do
121 : if (.not. okay) then
122 0 : write(*,2) 'model_number', s% model_number
123 0 : call mesa_error(__FILE__,__LINE__,'start solve omega: bad num omega')
124 : end if
125 :
126 0 : steps_used = 0
127 0 : recalc_mixing_info_each_substep = s% recalc_mixing_info_each_substep
128 :
129 : step_loop: do while &
130 0 : (total_time - time > 1d-10*total_time .and. &
131 0 : steps_used < max_steps)
132 :
133 0 : if (steps_used > 0 .and. recalc_mixing_info_each_substep) then
134 : call update_rotation_mixing_info(s,ierr)
135 0 : if (ierr /= 0) then
136 0 : do_solve_omega_mix = terminate
137 0 : s% termination_code = t_solve_omega_mix
138 0 : s% result_reason = nonzero_ierr
139 0 : if (s% report_ierr) write(*,*) 'update_rotation_mixing_info failed in do_solve_omega_mix'
140 0 : return
141 : end if
142 : end if
143 :
144 0 : if (steps_used == 0 .or. recalc_mixing_info_each_substep) then
145 0 : if (steps_used > 0) then
146 0 : do k=1,nz
147 0 : am_sig_omega(k) = s% am_sig_omega(k)
148 0 : am_sig_j(k) = s% am_sig_j(k)
149 : end do
150 : end if
151 0 : call get_rotation_sigmas(s, 1, nz, dt_total, ierr)
152 0 : if (ierr /= 0) then
153 0 : do_solve_omega_mix = terminate
154 0 : s% termination_code = t_solve_omega_mix
155 0 : s% result_reason = nonzero_ierr
156 0 : if (s% report_ierr) write(*,*) 'get_rotation_sigmas failed in do_solve_omega_mix'
157 0 : return
158 : end if
159 0 : if (steps_used > 0) then
160 0 : do k=1,nz
161 0 : s% am_sig_omega(k) = 0.5d0*(s% am_sig_omega(k) + am_sig_omega(k))
162 0 : s% am_sig_j(k) = 0.5d0*(s% am_sig_j(k) + am_sig_j(k))
163 : end do
164 : end if
165 : end if
166 :
167 0 : steps_used = steps_used + 1
168 :
169 0 : dt = 0.5d0*min_mixing_timescale()
170 0 : remaining_time = total_time - time
171 0 : dt = max(dt, 1d-6*remaining_time)
172 0 : dt = min(dt, real(dt_total/min_steps,kind=qp))
173 0 : if (dt >= remaining_time) then
174 0 : dt = remaining_time
175 : else
176 0 : dt = min(dt, 0.5d0*remaining_time)
177 : end if
178 0 : if (steps_used >= max_steps) dt = remaining_time ! just go for it
179 : if (dbg) write(*,3) 'mix dt', &
180 : s% model_number, steps_used, dt, dt/remaining_time
181 :
182 : ! X_0 is omega at start of substep
183 : ! X_1 is current candidate for omega at end of substep
184 : ! dX = X_1 - X_0
185 0 : do k=1,nz
186 0 : X_0(k) = s% omega(k)
187 0 : X_1(k) = X_0(k)
188 0 : dX(k) = 0d0
189 : end do
190 :
191 0 : solve_loop: do num_iters = 1, max_iters_per_substep
192 :
193 0 : if (total_num_iters >= max_iters_total) then
194 0 : s% retry_message = 'solve omega mix failed to converge in allowed number of steps'
195 0 : do_solve_omega_mix = retry
196 0 : exit step_loop
197 : end if
198 :
199 0 : total_num_iters = total_num_iters+1
200 :
201 0 : if (s% use_other_torque_implicit) then
202 0 : call s% other_torque_implicit(s% id, ierr)
203 0 : if (ierr /= 0) then
204 0 : s% retry_message = 'other_torque_implicit returned ierr'
205 0 : do_solve_omega_mix = retry
206 0 : exit step_loop
207 : end if
208 : end if
209 :
210 0 : if (associated(s% binary_other_torque_implicit)) then
211 : call s% binary_other_torque_implicit(s% id, ierr)
212 0 : if (ierr /= 0) then
213 0 : s% retry_message = 'binary_other_torque_implicit returned ierr'
214 0 : do_solve_omega_mix = retry
215 0 : exit step_loop
216 : end if
217 : end if
218 :
219 0 : call create_matrix_and_rhs(dt)
220 :
221 : ! solve for del
222 0 : call solve_tridiag(dl, d, du, rhs(1:nz), del(1:nz), nz, ierr)
223 0 : if (ierr /= 0) then
224 0 : s% retry_message = 'matrix solve failed in solve mix'
225 0 : do_solve_omega_mix = retry
226 0 : exit step_loop
227 : end if
228 :
229 : ! apply the correction dX = dX + del
230 : ! X_1 = X_0 + dX
231 : ! X_0 is omega at start of substep
232 : ! X_1 is candidate for omega at end of substep
233 0 : do k=2,nz
234 0 : dX(k) = dX(k) + del(k)
235 0 : X_1(k) = X_0(k) + dX(k)
236 0 : s% omega(k) = X_1(k)
237 : end do
238 0 : s% omega(1) = s% omega(2)
239 :
240 : ! if correction small enough, exit solve_loop
241 0 : max_del = maxval(abs(del(1:nz)))
242 0 : avg_del = sum(abs(del(1:nz)))/nz
243 0 : if (max_del <= tol_correction_max .and. avg_del <= tol_correction_norm) then
244 : if (dbg) &
245 : write(*,3) 'substep converged: iters max_del avg_del dt/total', &
246 : steps_used, num_iters, max_del, avg_del, dt/total_time
247 : exit solve_loop ! this substep is done
248 : end if
249 :
250 0 : if (num_iters == max_iters_per_substep) then
251 0 : s% retry_message = 'num_iters == max_iters_per_substep in solve mix'
252 0 : do_solve_omega_mix = retry
253 0 : exit step_loop
254 : end if
255 :
256 : end do solve_loop
257 :
258 0 : time = time + dt
259 :
260 : end do step_loop
261 :
262 : !if (recalc_mixing_info_each_substep) &
263 : ! write(*,3) 'omega mix steps_used', steps_used, s% model_number
264 :
265 : if (dbg) write(*,2) 'omega mix steps_used', steps_used
266 :
267 0 : s% num_rotation_solver_steps = max(steps_used, s% num_rotation_solver_steps)
268 :
269 0 : if (do_solve_omega_mix == keep_going .and. total_time - time > 1d-10*total_time) then
270 0 : do_solve_omega_mix = retry
271 0 : s% retry_message = 'failed in mixing angular momentum'
272 : end if
273 :
274 : if (do_solve_omega_mix == keep_going) then
275 :
276 0 : okay = .true.
277 0 : do k=1,nz
278 0 : if (is_bad_num(s% omega(k)) .or. abs(s% omega(k)) > 1d50) then
279 0 : write(*,2) 's% omega(k)', k, s% omega(k)
280 : okay = .false.
281 : exit
282 : end if
283 : end do
284 : if (.not. okay) then
285 0 : write(*,2) 'model_number', s% model_number
286 0 : call mesa_error(__FILE__,__LINE__,'end solve omega')
287 : end if
288 :
289 0 : do k=1,nz
290 0 : s% j_rot(k) = s% i_rot(k)% val*s% omega(k)
291 : end do
292 :
293 0 : if (.not. (s% use_other_torque .or. s% use_other_torque_implicit .or. &
294 : associated(s% binary_other_torque))) then
295 :
296 : ! check conservation for cases with no extra torque
297 0 : J_tot1 = total_angular_momentum(s) ! what we have
298 :
299 0 : if (abs(J_tot0 - J_tot1) > s% angular_momentum_error_retry*abs(J_tot0)) then
300 0 : s% retry_message = 'retry: failed to conserve angular momentum in mixing'
301 0 : write(*,*) "angular momentum error larger than angular_momentum_error_retry", abs(J_tot0 - J_tot1)/abs(J_tot0)
302 0 : do_solve_omega_mix = retry
303 0 : else if (abs(J_tot0 - J_tot1) > s% angular_momentum_error_warn*abs(J_tot0)) then
304 0 : write(*,*) "angular momentum error larger than angular_momentum_error_warn", abs(J_tot0 - J_tot1)/abs(J_tot0)
305 : end if
306 : if (dbg) then
307 : write(*,2) 'final J_tot1', s% model_number, J_tot1
308 : write(*,2) '(J_tot1 - J_tot0)/J_tot0', &
309 : steps_used, (J_tot1 - J_tot0)/J_tot0, J_tot0, J_tot1
310 : end if
311 : end if
312 :
313 : end if
314 :
315 : if (dbg) write(*,*)
316 :
317 0 : call dealloc
318 :
319 0 : if (s% doing_timing) &
320 0 : call update_time(s, time0, total, s% time_solve_omega_mix)
321 :
322 :
323 : contains
324 :
325 :
326 0 : subroutine do_alloc(ierr)
327 0 : use alloc, only: non_crit_get_quad_array
328 : integer, intent(out) :: ierr
329 0 : call do_work_arrays(.true.,ierr)
330 :
331 0 : call non_crit_get_quad_array(s, du, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
332 0 : if (ierr /= 0) return
333 0 : call non_crit_get_quad_array(s, d, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
334 0 : if (ierr /= 0) return
335 0 : call non_crit_get_quad_array(s, dl, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
336 0 : if (ierr /= 0) return
337 0 : call non_crit_get_quad_array(s, x, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
338 0 : if (ierr /= 0) return
339 0 : call non_crit_get_quad_array(s, b, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
340 0 : if (ierr /= 0) return
341 0 : call non_crit_get_quad_array(s, bp, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
342 0 : if (ierr /= 0) return
343 0 : call non_crit_get_quad_array(s, vp, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
344 0 : if (ierr /= 0) return
345 0 : call non_crit_get_quad_array(s, xp, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
346 0 : if (ierr /= 0) return
347 0 : call non_crit_get_quad_array(s, dX, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
348 0 : if (ierr /= 0) return
349 0 : call non_crit_get_quad_array(s, X_0, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
350 0 : if (ierr /= 0) return
351 0 : call non_crit_get_quad_array(s, X_1, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
352 0 : if (ierr /= 0) return
353 0 : call non_crit_get_quad_array(s, rhs, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
354 0 : if (ierr /= 0) return
355 0 : call non_crit_get_quad_array(s, del, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
356 0 : if (ierr /= 0) return
357 :
358 0 : end subroutine do_alloc
359 :
360 :
361 0 : subroutine dealloc
362 0 : use alloc, only: non_crit_return_quad_array
363 0 : call do_work_arrays(.false.,ierr)
364 :
365 0 : call non_crit_return_quad_array(s, du, 'solve_omega_mix')
366 0 : call non_crit_return_quad_array(s, d, 'solve_omega_mix')
367 0 : call non_crit_return_quad_array(s, dl, 'solve_omega_mix')
368 0 : call non_crit_return_quad_array(s, x, 'solve_omega_mix')
369 0 : call non_crit_return_quad_array(s, b, 'solve_omega_mix')
370 0 : call non_crit_return_quad_array(s, bp, 'solve_omega_mix')
371 0 : call non_crit_return_quad_array(s, vp, 'solve_omega_mix')
372 0 : call non_crit_return_quad_array(s, xp, 'solve_omega_mix')
373 :
374 0 : call non_crit_return_quad_array(s, dX, 'solve_omega_mix')
375 0 : call non_crit_return_quad_array(s, X_0, 'solve_omega_mix')
376 0 : call non_crit_return_quad_array(s, X_1, 'solve_omega_mix')
377 0 : call non_crit_return_quad_array(s, rhs, 'solve_omega_mix')
378 0 : call non_crit_return_quad_array(s, del, 'solve_omega_mix')
379 :
380 0 : end subroutine dealloc
381 :
382 :
383 0 : subroutine do_work_arrays(alloc_flag, ierr)
384 0 : use alloc, only: work_array
385 : logical, intent(in) :: alloc_flag
386 : integer, intent(out) :: ierr
387 : logical, parameter :: crit = .false.
388 : ierr = 0
389 : call work_array(s, alloc_flag, crit, &
390 0 : am_sig_omega, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
391 0 : if (ierr /= 0) return
392 : call work_array(s, alloc_flag, crit, &
393 0 : am_sig_j, nz, nz_alloc_extra, 'solve_omega_mix', ierr)
394 0 : if (ierr /= 0) return
395 0 : end subroutine do_work_arrays
396 :
397 :
398 0 : subroutine solve_tridiag(sub, diag, sup, rhs, x, n, ierr)
399 : ! sub - sub-diagonal
400 : ! diag - the main diagonal
401 : ! sup - sup-diagonal
402 : ! rhs - right hand side
403 : ! x - the answer
404 : ! n - number of equations
405 :
406 : integer, intent(in) :: n
407 : real(qp), dimension(:), intent(in) :: sup, diag, sub
408 : real(qp), dimension(:), intent(in) :: rhs
409 : real(qp), dimension(:), intent(out) :: x
410 : integer, intent(out) :: ierr
411 :
412 0 : real(qp) :: m
413 : integer :: i
414 :
415 0 : ierr = 0
416 :
417 0 : bp(1) = diag(1)
418 0 : vp(1) = rhs(1)
419 :
420 0 : do i = 2,n
421 0 : m = sub(i-1)/bp(i-1)
422 0 : bp(i) = diag(i) - m*sup(i-1)
423 0 : vp(i) = rhs(i) - m*vp(i-1)
424 : end do
425 :
426 0 : xp(n) = vp(n)/bp(n)
427 0 : x(n) = xp(n)
428 0 : do i = n-1, 1, -1
429 0 : xp(i) = (vp(i) - sup(i)*xp(i+1))/bp(i)
430 0 : x(i) = xp(i)
431 : end do
432 :
433 0 : end subroutine solve_tridiag
434 :
435 :
436 0 : real(dp) function min_mixing_timescale() result(dt)
437 : integer :: k
438 : real(dp) :: & ! use dp instead of qp to get same answer in ifort and gfortran
439 0 : omega, irot, irot_mid_00, am_sig_omega_00, c_omega_00, del00_omega, &
440 0 : omega_mid_00, am_sig_irot_00, c_irot_00, del00_irot, &
441 0 : dmbar, irot_mid_m1, am_sig_omega_m1, c_omega_m1, delm1_omega, &
442 0 : omega_mid_m1, am_sig_irot_m1, c_irot_m1, delm1_irot, &
443 0 : d2omega, d2irot, dt00
444 :
445 : include 'formats'
446 :
447 0 : dt = 1d99
448 :
449 0 : do k = 1, nz
450 :
451 0 : omega = s% omega(k)
452 0 : irot = s% i_rot(k)% val
453 :
454 0 : if (k < nz) then
455 :
456 0 : irot_mid_00 = 0.5d0*(irot + s% i_rot(k+1)% val)
457 0 : am_sig_omega_00 = s% am_sig_omega(k) + s% am_sig_j(k)
458 0 : c_omega_00 = am_sig_omega_00*irot_mid_00
459 0 : del00_omega = omega - s% omega(k+1)
460 :
461 0 : omega_mid_00 = 0.5d0*(omega + s% omega(k+1))
462 0 : am_sig_irot_00 = s% am_sig_j(k)
463 0 : c_irot_00 = am_sig_irot_00*omega_mid_00
464 0 : del00_irot = irot - s% i_rot(k+1)% val
465 :
466 : else
467 :
468 : c_omega_00 = 0
469 : del00_omega = 0
470 : c_irot_00 = 0
471 : del00_irot = 0
472 :
473 : end if
474 :
475 0 : if (k > 1) then
476 :
477 0 : if (k < nz) then
478 0 : dmbar = 0.5d0*(s% dm(k-1) + s% dm(k))
479 : else
480 0 : dmbar = 0.5d0*s% dm(k-1) + s% dm(k)
481 : end if
482 :
483 0 : irot_mid_m1 = 0.5d0*(s% i_rot(k-1)% val + irot)
484 0 : am_sig_omega_m1 = s% am_sig_omega(k-1) + s% am_sig_j(k-1)
485 0 : c_omega_m1 = am_sig_omega_m1*irot_mid_m1
486 0 : delm1_omega = s% omega(k-1) - omega
487 :
488 0 : omega_mid_m1 = 0.5d0*(s% omega(k-1) + omega)
489 0 : am_sig_irot_m1 = s% am_sig_j(k-1)
490 0 : c_irot_m1 = am_sig_irot_m1*omega_mid_m1
491 0 : delm1_irot = s% i_rot(k-1)% val - irot
492 :
493 : else
494 :
495 0 : dmbar = 0.5d0*s% dm(k)
496 0 : c_omega_m1 = 0
497 0 : delm1_omega = 0
498 0 : c_irot_m1 = 0
499 0 : delm1_irot = 0
500 :
501 : end if
502 :
503 0 : if (k == 1) then
504 0 : d2omega = -c_omega_00*del00_omega
505 0 : d2irot = -c_irot_00*del00_irot
506 0 : else if (k == nz) then
507 0 : d2omega = c_omega_m1*delm1_omega
508 0 : d2irot = c_irot_m1*delm1_irot
509 : else
510 0 : d2omega = c_omega_m1*delm1_omega - c_omega_00*del00_omega
511 0 : d2irot = c_irot_m1*delm1_irot - c_irot_00*del00_irot
512 : end if
513 :
514 : dt00 = max(1d-12,abs(omega))*irot/ &
515 : max(1d-50, abs((d2omega + d2irot)/dmbar + &
516 0 : s% extra_omegadot(k)*irot + s% extra_jdot(k)))
517 0 : if (dt00 < dt) dt = dt00
518 :
519 : end do
520 :
521 0 : end function min_mixing_timescale
522 :
523 :
524 0 : subroutine create_matrix_and_rhs(dt)
525 : ! basic equation from Heger, Langer, & Woosley, 2000, eqn 46.
526 : ! with source terms added.
527 : ! and term for j curvature as well as omega curvature
528 : real(qp), intent(in) :: dt
529 : integer :: k
530 : real(qp) :: &
531 0 : dmbar, f, &
532 0 : omega, omega_mid_00, omega_mid_m1, &
533 0 : irot, irot_mid_00, irot_mid_m1, &
534 0 : am_sig_omega_00, am_sig_omega_m1, c_omega_00, c_omega_m1, &
535 0 : am_sig_irot_00, am_sig_irot_m1, c_irot_00, c_irot_m1, &
536 0 : d_c_irot_00_domega_p1, d_c_irot_00_domega_00, &
537 0 : d_c_irot_m1_domega_00, d_c_irot_m1_domega_m1, &
538 0 : del00_omega, delm1_omega, &
539 0 : del00_irot, delm1_irot, &
540 0 : d2omega, d_d2omega_domega_p1, d_d2omega_domega_m1, d_d2omega_domega_00, &
541 0 : d2irot, d_d2irot_domega_p1, d_d2irot_domega_m1, d_d2irot_domega_00
542 :
543 : include 'formats'
544 :
545 0 : do k = 1, nz
546 :
547 0 : omega = s% omega(k)
548 0 : irot = s% i_rot(k)% val
549 :
550 0 : if (k < nz) then
551 :
552 0 : irot_mid_00 = 0.5d0*(irot + s% i_rot(k+1)% val)
553 0 : am_sig_omega_00 = s% am_sig_omega(k) + s% am_sig_j(k)
554 0 : c_omega_00 = am_sig_omega_00*irot_mid_00
555 0 : del00_omega = omega - s% omega(k+1)
556 :
557 0 : omega_mid_00 = 0.5d0*(omega + s% omega(k+1))
558 0 : am_sig_irot_00 = s% am_sig_j(k)
559 0 : c_irot_00 = am_sig_irot_00*omega_mid_00
560 0 : d_c_irot_00_domega_p1 = 0.5d0*am_sig_irot_00
561 0 : d_c_irot_00_domega_00 = 0.5d0*am_sig_irot_00
562 0 : del00_irot = irot - s% i_rot(k+1)% val
563 :
564 : else
565 :
566 : c_omega_00 = 0
567 : del00_omega = 0
568 : c_irot_00 = 0
569 : d_c_irot_00_domega_p1 = 0
570 : d_c_irot_00_domega_00 = 0
571 : del00_irot = 0
572 :
573 : end if
574 :
575 0 : if (k > 1) then
576 :
577 0 : if (k < nz) then
578 0 : dmbar = 0.5d0*(s% dm(k-1) + s% dm(k))
579 : else
580 0 : dmbar = 0.5d0*s% dm(k-1) + s% dm(k)
581 : end if
582 :
583 0 : irot_mid_m1 = 0.5d0*(s% i_rot(k-1)% val + irot)
584 0 : am_sig_omega_m1 = s% am_sig_omega(k-1) + s% am_sig_j(k-1)
585 0 : c_omega_m1 = am_sig_omega_m1*irot_mid_m1
586 0 : delm1_omega = s% omega(k-1) - omega
587 :
588 0 : omega_mid_m1 = 0.5d0*(s% omega(k-1) + omega)
589 0 : am_sig_irot_m1 = s% am_sig_j(k-1)
590 0 : c_irot_m1 = am_sig_irot_m1*omega_mid_m1
591 0 : d_c_irot_m1_domega_00 = 0.5d0*am_sig_irot_m1
592 0 : d_c_irot_m1_domega_m1 = 0.5d0*am_sig_irot_m1
593 0 : delm1_irot = s% i_rot(k-1)% val - irot
594 :
595 : else
596 :
597 0 : dmbar = 0.5d0*s% dm(k)
598 0 : c_omega_m1 = 0
599 0 : delm1_omega = 0
600 0 : c_irot_m1 = 0
601 0 : d_c_irot_m1_domega_00 = 0
602 0 : d_c_irot_m1_domega_m1 = 0
603 0 : delm1_irot = 0
604 :
605 : end if
606 :
607 0 : if (k == 1) then
608 0 : d2omega = -c_omega_00*del00_omega
609 0 : d2irot = -c_irot_00*del00_irot
610 0 : else if (k == nz) then
611 0 : d2omega = c_omega_m1*delm1_omega
612 0 : d2irot = c_irot_m1*delm1_irot
613 : else
614 0 : d2omega = c_omega_m1*delm1_omega - c_omega_00*del00_omega
615 0 : d2irot = c_irot_m1*delm1_irot - c_irot_00*del00_irot
616 : end if
617 0 : d_d2omega_domega_00 = -(c_omega_m1 + c_omega_00)
618 : d_d2irot_domega_00 = &
619 0 : d_c_irot_m1_domega_00*delm1_irot - d_c_irot_00_domega_00*del00_irot
620 :
621 : ! X_1 = X_0 + dX
622 : ! X_0 is omega at start of substep
623 : ! X_1 is candidate for omega at end of substep
624 :
625 : ! residual = dX - dt*(((d2omega+d2irot)/dmbar + extra_jdot)/irot + extra_omegadot)
626 : ! J = d(residual)/d(omega)
627 : ! del is linear estimate of change to dX to make residual = 0
628 : ! solve J*del = -residual == rhs
629 :
630 : rhs(k) = -dX(k) + &
631 : dt*(((d2omega + d2irot)/dmbar + &
632 0 : s% extra_jdot(k))/irot + s% extra_omegadot(k))
633 :
634 0 : f = dt/(dmbar*irot)
635 0 : d(k) = 1d0 - (d_d2omega_domega_00 + d_d2irot_domega_00)*f
636 :
637 0 : if (k < nz) then
638 0 : d_d2omega_domega_p1 = c_omega_00
639 0 : d_d2irot_domega_p1 = -d_c_irot_00_domega_p1*del00_irot
640 0 : du(k) = -(d_d2omega_domega_p1 + d_d2irot_domega_p1)*f
641 : end if
642 :
643 0 : if (k > 1) then
644 0 : d_d2omega_domega_m1 = c_omega_m1
645 0 : d_d2irot_domega_m1 = d_c_irot_m1_domega_m1*delm1_irot
646 0 : dl(k-1) = -(d_d2omega_domega_m1 + d_d2irot_domega_m1)*f
647 : end if
648 :
649 0 : if (s% use_other_torque_implicit) then
650 : d(k) = d(k) - &
651 : dt*(s% d_extra_jdot_domega_00(k)/irot + &
652 0 : s% d_extra_omegadot_domega_00(k))
653 0 : if (k < nz) &
654 : du(k) = du(k) - &
655 : dt*(s% d_extra_jdot_domega_p1(k)/irot + &
656 0 : s% d_extra_omegadot_domega_p1(k))
657 0 : if (k > 1) &
658 : dl(k-1) = dl(k-1) - &
659 : dt*(s% d_extra_jdot_domega_m1(k)/irot + &
660 0 : s% d_extra_omegadot_domega_m1(k))
661 : end if
662 :
663 : end do
664 :
665 0 : end subroutine create_matrix_and_rhs
666 :
667 : end function do_solve_omega_mix
668 :
669 : end module solve_omega_mix
|
