MAMMA  1.0
Conduitsolver
steady_solver.f90
Go to the documentation of this file.
1 !*********************************************************************
3 !
9 !*********************************************************************
10 
11 MODULE steady_solver
12  ! external variables
13  USE geometry, ONLY : z0 , zn, zeta_exit , radius
14 
15  USE parameters, ONLY : n_eqns , n_vars
16  USE parameters, ONLY : max_nl_iter
17  USE parameters, ONLY : verbose_level
18 
19  USE parameters, ONLY : idx_p1 , idx_p2 , idx_u1 , idx_u2 , idx_t , &
20  idx_xd_first , idx_xd_last , idx_alfa_first , idx_alfa_last , &
21  idx_beta_first , idx_beta_last
22 
23  USE parameters, ONLY : idx_mix_mass_eqn , idx_vol1_eqn , idx_mix_mom_eqn , &
24  idx_rel_vel_eqn , idx_mix_engy_eqn , idx_dis_gas_eqn_first , &
25  idx_dis_gas_eqn_last , idx_ex_gas_eqn_first , idx_ex_gas_eqn_last , &
26  idx_cry_eqn_first , idx_cry_eqn_last
27 
28  USE equations, ONLY : lateral_degassing_flag
29  USE equations, ONLY : lateral_degassing
30  USE equations, ONLY : alfa2_lat_thr
31 
32  USE constitutive, ONLY : explosive
33 
34  USE init, ONLY : p_out
35 
36  USE equations, ONLY : isothermal
37 
38  IMPLICIT NONE
39 
40  REAL*8 :: zeta
41  REAL*8, ALLOCATABLE :: fluxes_old(:)
42  REAL*8, ALLOCATABLE :: nh_terms_old(:)
43 
44  REAL*8 :: zeta_old , dz_max
45 
46  LOGICAL :: increase_flow_rate
47 
48  INTEGER :: nl_iter
49 
50 CONTAINS
51 
52  !******************************************************************************
54  !
62  !******************************************************************************
63 
64  SUBROUTINE steady_shooting
65 
66  ! external procedures
67  USE constitutive, ONLY : eos
68  USE constitutive, ONLY : eval_densities
69  USE equations, ONLY : phys_var_qp
70  USE init, ONLY : init_steady
71 
72  ! external variables
73  USE constitutive, ONLY : rho_m, rho_c, rho_mix, u_mix, u_1, u_2
74  USE constitutive, ONLY : p_1 , p_2 , t , s_1, s_2, e_1, e_2
75  USE constitutive, ONLY : rho_2 , rho_g, rho_1, rho_c, rho_m
76  USE constitutive, ONLY : mu_1, mu_2, s_c, mu_c, s_g, mu_g
77  USE constitutive, ONLY : e_c, e_g, e_m, mu_m, s_m
78  USE init, ONLY : u1_in
79  USE parameters, ONLY : shooting
80  USE parameters, ONLY : eps_conv
81  USE geometry, ONLY : update_radius
82 
83  IMPLICIT NONE
84 
85  REAL*8, ALLOCATABLE :: qp(:)
86 
87  REAL*8 :: V_0 , V_1, V_2
88  REAL*8 :: u_inlet
89 
90  LOGICAL :: check_flow_rate
91 
92  INTEGER :: iter
93 
94  REAL*8 :: extrap_z , extrap_z_p , extrap_z_mach
95  INTEGER :: extrap_flag
96  REAL*8 :: r_p_1 , r_p_2 , r_p
97  REAL*8 :: mach
98 
99  REAL*8 :: extrap_z_0 , extrap_z_2
100 
101  REAL*8 :: V_temp , V_coeff
102  LOGICAL :: check_interval
103 
104  LOGICAL :: flag_output , increase_temp
105 
106  REAL*8 :: initial_mass_flow_rate, final_mass_flow_rate
107 
108 
109  ALLOCATE( qp(n_vars) )
110  ALLOCATE( fluxes_old(n_eqns) )
111  ALLOCATE( nh_terms_old(n_eqns) )
112 
113 
114  IF ( shooting ) THEN
115 
116  IF ( explosive ) THEN
117 
118  v_temp = u1_in
119 
120  ELSE
121 
122  v_temp = u1_in
123 
124  END IF
125 
126  flag_output = .false.
127 
128  ELSE
129 
130  v_temp = u1_in
131 
132  flag_output = .false.
133 
134  END IF
135 
136  ! ------------------- SOLVE with u_inlet = V_0 -----------------------------
137 
138  u_inlet = v_temp
139 
140  CALL init_steady( u_inlet , qp )
141 
142  CALL phys_var_qp(qp)
143  CALL eos
144  CALL eval_densities
145 
146  IF ( verbose_level .GE. 1 ) THEN
147  p_1 = dcmplx( qp(idx_p1) , 0.d0 )
148  p_2 = dcmplx( qp(idx_p2) , 0.d0 )
149  u_1 = dcmplx( qp(idx_u1) , 0.d0 )
150  u_2 = dcmplx( qp(idx_u2) , 0.d0 )
151  t = dcmplx( qp(idx_t) , 0.d0 )
152 
153  WRITE(*,*) ''
154  WRITE(*,*) '<<<<<<<<<<<< ThermoPhysical quantities >>>>>>>>>>>'
155  WRITE(*,*) 'P_1 =', REAL(p_1)
156  WRITE(*,*) 'P_2 =', REAL(p_2)
157  WRITE(*,*) 'T =', REAL(t)
158  WRITE(*,*) ''
159 
160  WRITE(*,*) 'rho_c =', REAL(rho_c)
161  WRITE(*,*) 's_c =', REAL(s_c)
162  WRITE(*,*) 'e_c =', REAL(e_c)
163  WRITE(*,*) 'enth_c =', REAL(e_c + p_1 / rho_c)
164  WRITE(*,*) 'gibbs_c =', REAL(mu_c)
165  WRITE(*,*) ''
166 
167  WRITE(*,*) 'rho_g =', REAL(rho_g)
168  WRITE(*,*) 's_g =', REAL(s_g)
169  WRITE(*,*) 'e_g =', REAL(e_g)
170  WRITE(*,*) 'enth_g =', REAL(e_g + p_2 / rho_g)
171  WRITE(*,*) 'gibbs_g =', REAL(mu_g)
172  WRITE(*,*) ''
173 
174  WRITE(*,*) 'rho_m =', REAL(rho_m)
175  WRITE(*,*) 's_m =', REAL(s_m)
176  WRITE(*,*) 'e_m =', REAL(e_m)
177  WRITE(*,*) 'enth_m =', REAL(e_m + p_1 / rho_m)
178  WRITE(*,*) 'gibbs_m =', REAL(mu_m)
179  WRITE(*,*) ''
180 
181 
182  WRITE(*,*) 'rho_1 =', REAL(rho_1)
183  WRITE(*,*) 's_1 =', REAL(s_1)
184  WRITE(*,*) 'e_1 =', REAL(e_1)
185  WRITE(*,*) 'enth_1 =', REAL(e_1 + p_1 / rho_1)
186  WRITE(*,*) 'gibbs_1 =', REAL(mu_1)
187  WRITE(*,*) ''
188 
189  WRITE(*,*) 'rho_2 =', REAL(rho_2)
190  WRITE(*,*) 's_2 =', REAL(s_2)
191  WRITE(*,*) 'e_2 =', REAL(e_2)
192  WRITE(*,*) 'enth_2 =', REAL(e_2 + p_2 / rho_2)
193  WRITE(*,*) 'gibbs_2 =', REAL(mu_2)
194  WRITE(*,*) ''
195  WRITE(*,*) '<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>'
196 
197  END IF
198 
199  WRITE(*,*) ''
200  WRITE(*,*) '********* V_inlet = ',v_temp,'**********'
201 
202  CALL update_radius(zeta)
203 
204  initial_mass_flow_rate = REAL(rho_mix * u_mix)*3.1415*radius*radius
205 
206  WRITE(*,*) 'Initial mass flow rate = ',initial_mass_flow_rate,' kg/s'
207 
208  CALL integrate_equations(qp , flag_output , extrap_z , extrap_z_p , &
209  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
210 
211  IF ( .NOT. shooting ) THEN
212 
213  flag_output = .true.
214  zeta_exit = zeta
215  u_inlet = v_temp
216 
217  CALL init_steady( u_inlet , qp )
218 
219  CALL phys_var_qp(qp)
220 
221  WRITE(*,*) ''
222  WRITE(*,*) 'Saving solution for V at the inlet',v_temp
223 
224  CALL integrate_equations(qp , flag_output , extrap_z , extrap_z_p , &
225  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
226 
227  END IF
228 
229  CALL update_radius(zeta)
230 
231  final_mass_flow_rate = REAL(rho_mix * u_mix)*3.1415*radius*radius
232 
233  WRITE(*,*) 'Initial mass flow rate = ',initial_mass_flow_rate,' kg/s'
234  WRITE(*,*) 'Final mass flow rate = ',final_mass_flow_rate,' kg/s'
235 
236  extrap_z_0 = extrap_z
237 
238  IF ( increase_flow_rate ) THEN
239 
240  WRITE(*,*) 'Pressure 1 at exit:', r_p_1
241  WRITE(*,*) 'Pressure 2 at exit:', r_p_2
242  WRITE(*,*) 'Mach at exit:', mach
243 
244  IF ( extrap_flag .EQ. -3 ) THEN
245 
246  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
247 
248  ELSE
249 
250  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
251 
252  END IF
253 
254  v_0 = v_temp
255  v_coeff = 1.d1
256  v_2 = v_temp * v_coeff
257 
258  ELSE
259 
260  v_2 = v_temp
261  v_coeff = 5.d-1
262  v_0 = v_temp * v_coeff
263 
264  WRITE(*,*) 'At zeta ',zeta
265  WRITE(*,*) 'Pressure 1 is: ',r_p_1
266  WRITE(*,*) 'Pressure 2 is: ',r_p_2
267  WRITE(*,*) 'Mach is:', mach
268 
269  IF ( extrap_flag .EQ. 3 ) THEN
270 
271  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
272 
273  ELSE
274 
275  WRITE(*,*) 'Extrap_flag:',extrap_flag
276  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
277 
278  END IF
279 
280  END IF
281 
282  IF ( .NOT.shooting ) THEN
283 
284  RETURN
285 
286  END IF
287 
288  IF ( verbose_level .GE. 1 ) READ(*,*)
289 
290 
291  ! -------------------- SOLVE with u_inlet = V_2 ----------------------------
292 
293  check_interval = .false.
294  increase_temp = increase_flow_rate
295 
296  DO WHILE ( .NOT. check_interval )
297 
298  IF (v_coeff .LT. 1.0) THEN
299 
300  v_2 = v_temp
301 
302  ELSE
303 
304  v_0 = v_temp
305 
306  END IF
307 
308  v_temp = v_temp * v_coeff
309 
310  WRITE(*,*) '********* V_inlet = ',v_temp,'**********'
311 
312  u_inlet = v_temp
313 
314  CALL init_steady( u_inlet , qp )
315 
316  CALL phys_var_qp(qp)
317 
318  flag_output = .false.
319 
320  zeta = z0
321  CALL update_radius(zeta)
322  initial_mass_flow_rate = REAL(rho_mix * u_mix)*3.1415*radius*radius
323 
324  CALL integrate_equations(qp , flag_output , extrap_z , extrap_z_p , &
325  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
326 
327  CALL update_radius(zeta)
328  final_mass_flow_rate = REAL(rho_mix * u_mix)*3.1415*radius*radius
329 
330  WRITE(*,*) 'Initial mass flow rate = ',initial_mass_flow_rate,' kg/s'
331  WRITE(*,*) 'Final mass flow rate = ',final_mass_flow_rate,' kg/s'
332 
333  extrap_z_2 = extrap_z
334 
335 
336  IF ( increase_flow_rate ) THEN
337 
338  WRITE(*,*) 'Pressure 1 at exit:', r_p_1
339  WRITE(*,*) 'Pressure 2 at exit:', r_p_2
340  WRITE(*,*) 'Mach at exit:', mach
341 
342  IF ( extrap_flag .EQ. -3 ) THEN
343 
344  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
345 
346  ELSE
347 
348  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
349 
350  END IF
351 
352 
353  IF ( .NOT. increase_temp ) THEN
354 
355  check_interval = .true.
356  v_0 = v_temp
357 
358  END IF
359 
360  ELSE
361 
362  WRITE(*,*) 'At zeta ',zeta
363  WRITE(*,*) 'Pressure 1 is: ',r_p_1
364  WRITE(*,*) 'Pressure 2 is: ',r_p_2
365  WRITE(*,*) 'Mach is:', mach
366 
367  IF ( extrap_flag .EQ. 3 ) THEN
368 
369  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
370 
371  ELSE
372 
373  WRITE(*,*) 'Extrap_flag:',extrap_flag
374  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
375 
376  END IF
377 
378  IF ( increase_temp ) THEN
379 
380  check_interval = .true.
381  v_2 = v_temp
382 
383  END IF
384 
385  END IF
386 
387  END DO
388 
389  IF ( verbose_level .GE. 1 ) READ(*,*)
390 
391  ! ----------------------START LOOP FOR SOLUTION -----------------------------
392 
393  check_flow_rate = .false.
394 
395  iter = 0
396 
397  DO WHILE ( .NOT. check_flow_rate )
398 
399  iter = iter + 1
400 
401  v_1 = 0.5d0 * ( v_0 + v_2 )
402 
403  WRITE(*,*) '********* V_inlet = ',v_1,'********** iter = ', iter
404  WRITE(*,*) 'V_0 = ',v_0,' V_1 = ',v_1,' V_2 = ',v_2
405  WRITE(*,*) 'Delta V = ', 0.5d0 * ( v_2 - v_0 )
406 
407  u_inlet = v_1
408 
409  CALL init_steady( u_inlet , qp )
410 
411  CALL phys_var_qp(qp)
412 
413  flag_output = .false.
414 
415  zeta = z0
416  CALL update_radius(zeta)
417  initial_mass_flow_rate = REAL(rho_mix * u_mix)*3.14*radius*radius
418 
419  CALL integrate_equations( qp , flag_output , extrap_z , extrap_z_p , &
420  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
421 
422  CALL update_radius(zeta)
423  final_mass_flow_rate = REAL(rho_mix * u_mix)*3.14*radius*radius
424 
425  WRITE(*,*) 'Initial mass flow rate = ',initial_mass_flow_rate,' kg/s'
426  WRITE(*,*) 'Final mass flow rate = ',final_mass_flow_rate,' kg/s'
427 
428  WRITE(*,*) 'Extrap_flag:',extrap_flag
429  WRITE(*,*) 'Extrap_z_p',extrap_z_p
430  WRITE(*,*) 'Extrap_z_mach',extrap_z_mach
431 
432  r_p = min( r_p_1 , r_p_2 )
433 
434  IF ( increase_flow_rate ) THEN
435 
436  WRITE(*,*) 'Pressure 1 at exit:', r_p_1
437  WRITE(*,*) 'Pressure 2 at exit:', r_p_2
438  WRITE(*,*) 'Mach at exit:', mach
439 
440 
441  IF ( extrap_z_mach .LE. extrap_z_p ) THEN
442 
443  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
444 
445  ELSE
446 
447  WRITE(*,*) 'Extrap_flag:',extrap_flag
448  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
449 
450  END IF
451 
452  IF ( abs(r_p - p_out) .LT. eps_conv * p_out ) THEN
453 
454  check_flow_rate = .true.
455 
456  END IF
457 
458  ELSE
459 
460  IF ( ( abs( mach - 1.d0 ) .LT. eps_conv ) .AND. ( zeta .EQ. zn ) ) THEN
461 
462  WRITE(*,*) 'Mach at exit:', mach
463  WRITE(*,*) 'Pressure 1 at exit:', r_p_1
464  WRITE(*,*) 'Pressure 2 at exit:', r_p_2
465 
466  check_flow_rate = .true.
467 
468  ELSE
469 
470  IF ( zeta .GE. zn ) THEN
471 
472  WRITE(*,*) 'Mach at exit:', mach
473  WRITE(*,*) 'Pressure 1 at exit:', r_p_1
474  WRITE(*,*) 'Pressure 2 at exit:', r_p_2
475 
476  ELSE
477 
478  WRITE(*,*) 'At zeta ',zeta
479  WRITE(*,*) 'Pressure 1 is: ',r_p_1
480  WRITE(*,*) 'Pressure 2 is: ',r_p_2
481  WRITE(*,*) 'Mach is:', mach
482 
483  IF ( extrap_flag .EQ. 3 ) THEN
484 
485  WRITE(*,*) 'Mach = 1 is reached at:',extrap_z
486 
487  ELSE
488 
489  WRITE(*,*) 'Extrap_flag:',extrap_flag
490  WRITE(*,*) 'Exit pressure is reached at:',extrap_z
491 
492  END IF
493 
494  END IF
495 
496  END IF
497 
498  END IF
499 
500  IF ( increase_flow_rate ) THEN
501 
502  v_0 = v_1
503 
504 
505  extrap_z_0 = extrap_z
506 
507  ELSE
508 
509  v_2 = v_1
510 
511  extrap_z_2 = extrap_z
512 
513  END IF
514 
515 
516  IF ( abs(r_p - p_out) .GT. p_out ) THEN
517 
518  check_flow_rate = .false.
519 
520  END IF
521 
522 
523  IF ( iter .GT. 20 ) THEN
524 
525  IF ( ( ( v_2 - v_0 ) / v_0 .LT. eps_conv ) .AND. &
526  ( zeta .GE. zn - (iter-20)/100.d0 ) .AND. ( extrap_flag .NE. 4 ) &
527  .AND. ( extrap_flag .NE. -3 ) ) THEN
528 
529  WRITE(*,*) 'Relative change in flow rate', ( v_2 - v_0 ) / v_0
530 
531  check_flow_rate = .true.
532 
533  END IF
534 
535  IF ( ( ( v_2 - v_0 ) / v_0 .LT. eps_conv ) .AND. &
536  ( zeta .GE. zn - (iter-20)/100.d0 ) .AND. &
537  ( abs(mach - 1.0) .LT. 0.01) ) THEN
538 
539  WRITE(*,*) 'Relative change in flow rate', ( v_2 - v_0 ) / v_0
540 
541  check_flow_rate = .true.
542 
543  END IF
544 
545  END IF
546 
547  IF ( iter .GT. 100 ) THEN
548 
549  IF ( ( zeta .GE. zn - (iter-20)/100.d0 ) .AND. ( extrap_flag .NE. 4 ) &
550  .AND. ( extrap_flag .NE. -3 ) ) THEN
551 
552  WRITE(*,*) 'Relative change in flow rate', ( v_2 - v_0 ) / v_0
553 
554  check_flow_rate = .true.
555 
556  END IF
557 
558  IF ( ( zeta .GE. zn - (iter-20)/100.d0 ) .AND. ( abs(r_p - p_out) &
559  .LT. (1.0d0 + (iter-50)/10.d0) * p_out) ) THEN
560 
561  WRITE(*,*) 'Relative change in flow rate', ( v_2 - v_0 ) / v_0
562 
563  check_flow_rate = .true.
564 
565  END IF
566 
567  END IF
568 
569  IF ( ( ( v_2 - v_0 ) / v_0 .LT. eps_conv ) .AND. &
570  ( zeta .GE. zn ) .AND. ( extrap_flag .NE. 4 ) &
571  ! .AND. ( extrap_flag .NE. -3 ) &
572  ) THEN
573 
574  WRITE(*,*) 'Relative change in flow rate', ( v_2 - v_0 ) / v_0
575 
576  check_flow_rate = .true.
577 
578  END IF
579 
580  IF ( verbose_level .GE. 1 ) READ(*,*)
581 
582  END DO
583 
584  ! --- Integrate again with the found value of u_inlet saving the output -----
585 
586  flag_output = .true.
587  zeta_exit = zeta
588  u_inlet = v_1
589 
590  CALL init_steady( u_inlet , qp )
591 
592  CALL phys_var_qp(qp)
593 
594  WRITE(*,*) 'Number of iterations',iter
595  WRITE(*,*) 'Saving solution for V at the inlet',v_1
596 
597  CALL integrate_equations(qp , flag_output , extrap_z , extrap_z_p , &
598  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
599 
600  END SUBROUTINE steady_shooting
601 
602  !******************************************************************************
604  !
621  !******************************************************************************
622 
623  SUBROUTINE integrate_equations( qp , flag_output , extrap_z , extrap_z_p , &
624  extrap_z_mach , extrap_flag , r_p_1 , r_p_2 , mach )
625 
626  ! external procedures
627  USE constitutive, ONLY : eos
628  USE constitutive, ONLY : sound_speeds
629 
630  USE constitutive, ONLY : eval_densities
631  USE constitutive, ONLY : f_alfa3, f_alfa
632 
633  USE equations, ONLY : eval_nonhyperbolic_terms_qp
634  USE equations, ONLY : eval_fluxes_qp
635  USE equations, ONLY : phys_var_qp
636 
637 
638  USE geometry, ONLY : update_radius
639 
640  ! external variables
641  USE constitutive, ONLY : frag_thr, frag_eff
642  USE constitutive, ONLY : rho_mix, u_mix
643 
644  USE inpout, ONLY : output_steady
645 
646  ! external variables
647  USE constitutive, ONLY : zeta_lith
648 
649  USE geometry, ONLY : z_comp , comp_cells, radius
650 
651  USE parameters, ONLY : shooting
652 
653  USE parameters, ONLY : tol_abs , tol_rel, n_gas
654 
655  IMPLICIT NONE
656 
657  REAL*8, INTENT(INOUT) :: qp(n_eqns)
658  LOGICAL, INTENT(IN) :: flag_output
659  REAL*8,INTENT(OUT) :: extrap_z , extrap_z_p , extrap_z_mach
660  INTEGER,INTENT(OUT) :: extrap_flag
661  REAL*8, INTENT(OUT) :: r_p_1 , r_p_2
662  REAL*8, INTENT(OUT) :: mach
663 
665  REAL*8 :: dz
666 
668  REAL*8 :: C_mix
669 
671  REAL*8 :: check_error , check_error_old, coeff_z
672 
673  REAL*8 :: r_alfa_2
674 
675  LOGICAL :: fragmentation
676 
677  INTEGER :: idx_zeta
678 
679  REAL*8 :: coeff_zeta
680 
682  REAL*8 :: qp_comp(n_eqns)
683 
685  REAL*8 :: qp_old(n_eqns)
686 
688  REAL*8 :: qp_full(n_eqns)
689 
691  REAL*8 :: qp_half(n_eqns)
692 
694  REAL*8 :: qp_half2(n_eqns)
695 
696  REAL*8 :: fluxes_temp(n_vars) , nh_terms_temp(n_vars)
697 
698  REAL*8 :: dqp_dz(n_vars)
699  ! REAL*8 :: dqp_dz_half(n_vars)
700 
701  LOGICAL :: check_convergence
702 
703  REAL*8 :: qp_guess(n_vars)
704 
705  REAL*8 :: delta_full , delta_half2
706 
707  LOGICAL :: fragmentation_half2, fragmentation_full
708 
709  REAL*8 :: alfa_2_half2, alfa_2_full, alfa_2_qp
710 
711  REAL*8 :: strain_rate_qp
712 
713  REAL*8 :: u_1_old
714 
715  INTEGER :: counter
716 
717  delta_full = 0.d0
718  delta_half2 = 0.d0
719 
720 
721  zeta = z0
722 
723  idx_zeta = 1
724 
725  IF ( flag_output) shooting = .false.
726 
727  CALL update_radius(zeta)
728 
729  IF ( ( .NOT. shooting ) .OR. ( flag_output ) ) CALL output_steady( zeta , &
730  qp , radius )
731 
732  dz_max = ( zn - z0 ) / comp_cells
733 
734  dz = dz_max
735 
736  IF ( verbose_level .GE. 1 ) THEN
737 
738  WRITE(*,*) 'qp = '
739  WRITE(*,*) qp
740 
741  READ(*,*)
742 
743  ENDIF
744 
745  lateral_degassing = .false.
746  fragmentation = .false.
747 
748  check_error = 1.d0
749  check_error_old = 1.d0
750 
751  dqp_dz(1:n_vars) = 0.d0
752 
753  counter = 1
754 
755  frag_eff = 0.d0
756 
757  u_1_old = qp(idx_u1)
758  strain_rate_qp = 0.0
759 
760  zeta_integration:DO
761 
762  counter = counter + 1
763 
764  IF ( counter .GT. 10 * comp_cells ) THEN
765 
766  WRITE(*,*)'Convergence error: Too many iterations'
767  WRITE(*,*)'counter = ', counter
768  stop
769 
770  END IF
771 
772  zeta_old = zeta
773 
774  qp_old = qp
775 
776  u_1_old = qp_old(idx_u1)
777 
778  zeta_lith = zeta
779 
780  CALL update_radius(zeta)
781 
782  IF ( verbose_level .GE. 1 ) THEN
783 
784  WRITE(*,*) 'zeta, radius, counter',zeta, radius, counter
785 
786  END IF
787 
788  CALL eval_fluxes_qp( r_qp = qp_old , r_flux = fluxes_old )
789 
790  CALL eval_nonhyperbolic_terms_qp( r_qp = qp_old , r_nh_term_impl = &
791  nh_terms_old )
792 
793  IF ( verbose_level .GE. 3 ) THEN
794 
795  WRITE(*,*) 'Non hyperbolic term = '
796  WRITE(*,*) nh_terms_old/(radius**2)
797  READ(*,*)
798 
799  END IF
800 
801  IF ( isothermal ) THEN
802 
803  fluxes_old(idx_mix_engy_eqn) = qp_old(idx_t)
804  nh_terms_old(idx_mix_engy_eqn) = 0.d0
805 
806  END IF
807 
808 
809  dz = min( dz , zn - zeta_old )
810 
811  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'z = ',zeta,'dz = ',dz
812 
813  find_step_loop:DO
814 
815  zeta = zeta_old + dz
816  zeta_lith = zeta
817  CALL update_radius(zeta)
818 
819  ! ------------integrate the whole step --------------------------------
820 
821  IF ( verbose_level .GE. 1 ) THEN
822  WRITE(*,*)''
823  WRITE(*,*) '/--------full step---------/'
824  END IF
825 
826  qp_full = qp_old
827 
828  qp_guess = qp_full
829 
830  CALL perturbe_qp(qp_full)
831 
832  CALL advance_dz( qp_full , dz , check_convergence )
833 
834  IF ( .NOT. check_convergence ) THEN
835 
836  qp_full = qp_old + dqp_dz * dz
837 
838  CALL perturbe_qp(qp_full)
839 
840  CALL advance_dz( qp_full , dz , check_convergence )
841 
842  END IF
843 
844  fluxes_temp = fluxes_old
845  nh_terms_temp = nh_terms_old
846 
847  IF ( check_convergence ) THEN
848 
849  ! ------------integrate the first half step ------------------------
850 
851  zeta = zeta_old + 0.5d0 * dz
852  zeta_lith = zeta
853  CALL update_radius(zeta)
854 
855 
856  IF ( verbose_level .GE. 1 ) THEN
857  WRITE(*,*)''
858  WRITE(*,*) '/--------fist half step---------/'
859  END IF
860 
861  qp_half(1:n_vars) = qp_old(1:n_vars)
862 
863  CALL perturbe_qp(qp_half)
864 
865  CALL advance_dz( qp_half , 0.5 * dz , check_convergence )
866 
867  ! if the integration failed try with a different initial guess:
868  ! use an average between qp_old and qp_full
869  IF ( .NOT. check_convergence ) THEN
870 
871  qp_half(1:n_vars) = 0.5d0 * (qp_old(1:n_vars) + &
872  qp_full(1:n_vars) )
873 
874  CALL perturbe_qp(qp_half)
875 
876  CALL advance_dz( qp_half , 0.5 * dz , check_convergence )
877 
878  END IF
879 
880  ! if the integration failed try with a different initial guess:
881  ! use qp_full
882  IF ( .NOT. check_convergence ) THEN
883 
884  qp_half = qp_full
885 
886  CALL perturbe_qp(qp_half)
887 
888  CALL advance_dz( qp_half , 0.5 * dz , check_convergence )
889 
890  END IF
891 
892 
893  END IF
894 
895  IF ( check_convergence ) THEN
896 
897  ! ------------integrate the second half step -----------------------
898 
899  zeta = zeta_old + dz
900 
901  IF ( verbose_level .GE. 1 ) THEN
902  WRITE(*,*)''
903  WRITE(*,*) '/--------second half step---------/'
904  END IF
905 
906  ! use the values from the first step as old values
907  CALL update_radius(zeta_old + 0.5d0 * dz)
908  CALL eval_fluxes_qp( r_qp = qp_half , r_flux = fluxes_old )
909 
910  CALL eval_nonhyperbolic_terms_qp( r_qp = qp_half , r_nh_term_impl =&
911  nh_terms_old )
912 
913  IF ( isothermal ) THEN
914 
915  fluxes_old(idx_mix_engy_eqn) = qp_old(idx_t)
916  nh_terms_old(idx_mix_engy_eqn) = 0.d0
917 
918  END IF
919 
920  zeta_lith = zeta
921  CALL update_radius(zeta)
922 
923  ! use the solution of the whole step as initial guess
924  qp_half2 = qp_half
925 
926  CALL perturbe_qp(qp_half2)
927 
928  CALL advance_dz( qp_half2 , 0.5 * dz , check_convergence )
929 
930  ! if the integration failed try with a different initial guess:
931  ! use the solution from the whole step
932  IF ( .NOT. check_convergence ) THEN
933 
934  qp_half2 = qp_full
935 
936  CALL perturbe_qp(qp_half2)
937 
938  CALL advance_dz( qp_half2 , 0.5 * dz , check_convergence )
939 
940  END IF
941 
942  ! if the integration failed again try with another initial guess:
943  ! use a solution extrapolated from the first half step
944  IF ( .NOT. check_convergence ) THEN
945 
946  qp_half2 = qp_old + 0.50 * dz * ( qp_half - qp_old )
947 
948  CALL perturbe_qp(qp_half2)
949 
950  CALL advance_dz( qp_half2 , 0.5 * dz , check_convergence )
951 
952  END IF
953 
954  END IF
955 
956  IF ( verbose_level .GE. 1 ) THEN
957  WRITE(*,*)''
958  WRITE(*,*) '<<<<<<<<< check_convergence >>>>>>>> ',check_convergence
959  READ(*,*)
960  END IF
961 
962  IF ( check_convergence ) THEN
963 
964  ! -------- compare the solutions obtained with dz and twice dz/2 ---
965 
966  check_error = maxval(abs( qp_full(1:n_vars) - &
967  qp_half2(1:n_vars) ) / ( tol_abs + tol_rel * &
968  max( abs( qp_full(1:n_vars) ) , &
969  abs( qp_half2(1:n_vars) ) ) ) )
970 
971  delta_full = maxval(abs( qp_full(1:n_vars) - &
972  qp_guess(1:n_vars) ) / ( tol_abs + tol_rel * &
973  max( abs( qp_full(1:n_vars) ) , &
974  abs( qp_guess(1:n_vars) ) ) ) )
975 
976  delta_half2 = maxval(abs( qp_half2(1:n_vars) - &
977  qp_guess(1:n_vars) ) / ( tol_abs + tol_rel * &
978  max( abs( qp_half2(1:n_vars) ) , &
979  abs( qp_guess(1:n_vars) ) ) ) )
980 
981  IF ( verbose_level .GE. 1 ) THEN
982 
983  WRITE(*,*) 'check dz and dz/2 error',check_error
984 
985  END IF
986 
987  IF ( max( qp_half2(idx_p1) , qp_half2(idx_p2) ) &
988  .LT. p_out) THEN
989 
990  check_convergence = .false.
991  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'pressure', &
992  max( qp_half2(idx_p1) , qp_half2(idx_p2) )
993 
994  ELSE
995 
996  IF ( verbose_level .GE. 2 ) THEN
997  WRITE(*,*)''
998  WRITE(*,*) 'qp_half2'
999  WRITE(*,*) qp_half2
1000  WRITE(*,*)''
1001  WRITE(*,*) 'qp_full'
1002  WRITE(*,*) qp_full
1003  WRITE(*,*)''
1004  WRITE(*,*) 'alfa_2',sum(qp_half2(1:n_gas))
1005  READ(*,*)
1006 
1007  END IF
1008 
1009  END IF
1010 
1011  END IF
1012 
1013 
1014  ! ---- Check if the fragmentation threshold is reached ---------------
1015 
1016  IF ( explosive ) THEN
1017 
1018  alfa_2_half2 = sum(qp_half2(idx_alfa_first:idx_alfa_last))
1019  alfa_2_full = sum(qp_full(idx_alfa_first:idx_alfa_last))
1020  alfa_2_qp = sum(qp(idx_alfa_first:idx_alfa_last))
1021 
1022  IF ( alfa_2_half2 .GT. frag_thr ) THEN
1023 
1024  fragmentation_half2 = .true.
1025 
1026  ELSE
1027 
1028  fragmentation_half2 = .false.
1029 
1030  END IF
1031 
1032  IF ( alfa_2_full .GT. frag_thr ) THEN
1033 
1034  fragmentation_full = .true.
1035 
1036  ELSE
1037 
1038  fragmentation_full = .false.
1039 
1040  END IF
1041 
1042 
1043  IF ( ( .NOT.fragmentation ) &
1044  .AND. ( fragmentation_half2 .OR. fragmentation_full ) &
1045  .AND. ( frag_thr - alfa_2_qp .GT. 1.d-4 ) &
1046 ! .AND. ( MAX(alfa_2_half2,alfa_2_full) - frag_thr .GT. 1.D-4 ) &
1047  ) THEN
1048 
1049  check_convergence = .false.
1050 
1051  END IF
1052 
1053  END IF
1054 
1055  IF ( ( check_error .LT. 1.d0 ) .AND. ( check_convergence ) ) THEN
1056 
1057  ! --- if the error is small accept the solution obtained with dz/2
1058  ! --- and exit from the search loop
1059 
1060  coeff_z = max( 1.05d0 , min( 1.50d0 , check_error ** ( -0.35d0 ) * &
1061  check_error_old ** ( 0.2d0 ) ) )
1062 
1063  check_error_old = check_error
1064 
1065  ! evaluate the slope of the solution for the guess at the next step
1066 
1067  IF ( delta_half2 .LT. delta_full ) THEN
1068 
1069  qp = qp_half2
1070  dqp_dz = ( qp_half2 - qp_old ) / dz
1071 
1072  ELSE
1073 
1074  qp = qp_full
1075  dqp_dz = ( qp_full - qp_old ) / dz
1076 
1077  END IF
1078 
1079  CALL phys_var_qp(qp)
1080 
1081  IF ( verbose_level .GE. 1 ) THEN
1082 
1083  WRITE(*,*) '::::::::: checks ok :::::::::'
1084  WRITE(*,*) 'zeta = ',zeta,' check_error = ',check_error, &
1085  ' coeff_z = ',coeff_z
1086  WRITE(*,*) ''
1087  WRITE(*,*) 'qp ='
1088  WRITE(*,*) qp
1089  WRITE(*,*) ''
1090  WRITE(*,*) 'alfa_2 =',sum(qp(idx_alfa_first:idx_alfa_last)), &
1091  ' fragmentation = ', frag_eff, 'Mass flow rate = ', &
1092  REAL(rho_mix * u_mix)*3.14*radius*radius
1093  WRITE(*,*)''
1094  READ(*,*)
1095  END IF
1096 
1097  EXIT find_step_loop
1098 
1099  ELSE
1100 
1101 
1102  ! --- if the error is big repeat the step with smaller dz
1103 
1104  dz = 0.950d0 * dz
1105 
1106  fluxes_old = fluxes_temp
1107  nh_terms_old = nh_terms_temp
1108 
1109  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'zeta_old',zeta_old, &
1110  'dz = ',dz
1111 
1112  IF ( dz .LT. 1e-20 ) THEN
1113 
1114  WRITE(*,*)'Convergence Error: dz too small'
1115  WRITE(*,*)'dz =', dz
1116  stop
1117 
1118  END IF
1119 
1120  END IF
1121 
1122  END DO find_step_loop
1123 
1124  ! ---------- Write the output on file when the solution is found ---------
1125 
1126  IF ( flag_output ) THEN
1127 
1128  IF ( shooting ) THEN
1129 
1130  DO WHILE ( ( z_comp(idx_zeta) .LT. zeta ) .AND. &
1131  ( idx_zeta .LE. comp_cells ) )
1132 
1133  coeff_zeta = ( z_comp(idx_zeta) - ( zeta - dz ) ) / dz
1134 
1135  qp_comp = coeff_zeta * qp + ( 1.d0 - coeff_zeta ) * qp_old
1136 
1137  CALL output_steady(z_comp(idx_zeta),qp_comp,radius)
1138 
1139  idx_zeta = idx_zeta + 1
1140 
1141  END DO
1142 
1143  IF ( zeta .EQ. zn ) CALL output_steady(zeta,qp,radius)
1144 
1145  ELSE
1146 
1147  CALL output_steady(zeta,qp,radius)
1148 
1149  END IF
1150 
1151  END IF
1152 
1153  ! ----- Update the integration step with the coefficient coeff_z ---------
1154 
1155  dz = min( coeff_z * dz , dz_max )
1156 
1157  ! ----- Evaluate some phys. variables to pass out of the subroutine ------
1158 
1159  r_alfa_2 = sum(qp_half2(idx_alfa_first:idx_alfa_last))
1160 
1161  r_p_1 = qp(idx_p1)
1162  r_p_2 = qp(idx_p2)
1163 
1164  CALL phys_var_qp( r_qp = qp )
1165  CALL eos
1166  CALL sound_speeds( c_mix , mach )
1167 
1168  ! ------ Check if the extrapolated solution or the solution at zeta ------
1169  ! ------ reach the boundary conditions -----------------------------------
1170 
1171  CALL linear_extrapolation(zeta_old,qp_old,zeta,qp,extrap_z,extrap_z_p, &
1172  extrap_z_mach,extrap_flag)
1173 
1174  IF ( extrap_flag .GT. 0 ) THEN
1175 
1176  increase_flow_rate = .false.
1177 
1178  RETURN
1179 
1180  END IF
1181 
1182  ! ----- Check if the top is reached --------------------------------------
1183 
1184  IF ( zeta .GE. zn ) THEN
1185 
1186  increase_flow_rate = .true.
1187 
1188  RETURN
1189 
1190  END IF
1191 
1192  ! ---- Check if the fragmentation threshold is reached -------------------
1193 
1194  IF ( explosive ) THEN
1195 
1196  alfa_2_qp = sum(qp(idx_alfa_first:idx_alfa_last))
1197 
1198  IF ( ( alfa_2_qp .GT. frag_thr) .AND. &
1199  ( .NOT. fragmentation ) )THEN
1200 
1201  frag_eff = 1.0d0
1202  fragmentation = .true.
1203 
1204  WRITE(*,*) 'Fragmentation at z = ',zeta
1205  ! verbose_level = 3
1206 
1207  END IF
1208 
1209  END IF
1210 
1211  ! ---- Check if the lateral degassing threshold is reached ---------------
1212 
1213  IF ( ( r_alfa_2 .GE. alfa2_lat_thr ) .AND. &
1214  ( lateral_degassing_flag ) .AND. ( .NOT.lateral_degassing ) ) THEN
1215 
1216  lateral_degassing = .true.
1217 
1218  WRITE(*,*) 'Lateral degassing from z = ',zeta
1219 
1220  END IF
1221 
1222  END DO zeta_integration
1223 
1224  !WRITE(*,*) counter
1225 
1226  END SUBROUTINE integrate_equations
1227 
1228  !******************************************************************************
1230  !
1237  !******************************************************************************
1238 
1239  SUBROUTINE perturbe_qp(qp)
1240 
1241  IMPLICIT NONE
1242 
1243  REAL*8, INTENT(INOUT) :: qp(n_vars)
1244 
1245  !Pressure phase 2
1246  qp(idx_p2) = qp(idx_p2) + 1.d-2
1247 
1248  !Velocity phase 1
1249  IF ( qp(idx_u1) - qp(idx_u2) .GT. -1d-7 ) THEN
1250 
1251  qp(idx_u2) = qp(idx_u2) * ( 1.0001d0 )
1252 
1253  END IF
1254 
1255  END SUBROUTINE perturbe_qp
1256 
1257  !******************************************************************************
1259  !
1268  !******************************************************************************
1269 
1270  SUBROUTINE advance_dz( qp , dz , check_convergence )
1271 
1272  ! external variables
1273  USE parameters, ONLY : alfa_impl
1274 
1275  ! external subroutines
1276  USE constitutive, ONLY : eos
1277  USE equations, ONLY : phys_var_qp
1278 
1279  IMPLICIT none
1280 
1281  REAL*8, INTENT(INOUT) :: qp(n_eqns)
1282  REAL*8, INTENT(IN) :: dz
1283  LOGICAL, INTENT(OUT) :: check_convergence
1284 
1285  REAL*8 :: qp_rel(n_vars)
1286  REAL*8 :: qp_org(n_vars)
1287 
1288  REAL*8 :: left_matrix(n_eqns,n_eqns)
1289  REAL*8 :: right_term(n_eqns)
1290 
1291  REAL*8 :: scal_f , scal_f_old
1292 
1293  INTEGER :: pivot(n_eqns)
1294  INTEGER :: ok
1295 
1296  INTEGER :: idx_max(1)
1297 
1298  REAL*8 :: delta_qp_rel(n_eqns)
1299  REAL*8 :: grad_f(n_eqns)
1300 
1301  REAL*8, DIMENSION(size(qp)) :: qp_NR_old , desc_dir
1302 
1303  REAL*8 :: qp_rel_NR_old(n_vars)
1304 
1305  REAL*8 :: check_NR_error
1306 
1307  REAL*8 :: scal_f_init
1308 
1309  LOGICAL :: opt_search_NL
1310  REAL*8, PARAMETER :: STPMX=100.d0
1311  REAL*8 :: stpmax
1312  LOGICAL :: check
1313  ! LOGICAL :: check_phys_var
1314 
1315  REAL*8, PARAMETER :: tol_rel_NR = 1.d-20 , tol_abs_nr = 1.d-20
1316 
1317 
1318  LOGICAL :: normalize_qp
1319  LOGICAL :: normalize_f
1320 
1321  REAL*8 :: coeff_f(n_eqns)
1322 
1323  REAL*8 :: arg_check(n_eqns)
1324 
1325  INTEGER :: i,j
1326 
1327  check_convergence = .false.
1328 
1329  normalize_qp = .true.
1330  normalize_f = .false.
1331 
1332  IF ( normalize_qp ) THEN
1333 
1334  DO i = 1,n_vars
1335 
1336  IF ( qp(i) .EQ. 0.d0 ) THEN
1337 
1338  qp_org(i) = 1.d0
1339 
1340  ELSE
1341 
1342  qp_org(i) = qp(i)
1343 
1344  END IF
1345 
1346  END DO
1347 
1348  ELSE
1349 
1350  qp_org(1:n_vars) = 1.0d0
1351 
1352  END IF
1353 
1354  qp_rel = qp / qp_org
1355 
1356  IF ( normalize_qp ) THEN
1357 
1358  qp_rel = qp_rel * 1.00001d0
1359 
1360  END IF
1361 
1362 
1363  coeff_f(1:n_eqns) = 1.d0
1364 
1365  IF ( normalize_f ) THEN
1366 
1367  CALL eval_f( qp_rel , qp_org , dz , coeff_f , right_term , scal_f )
1368 
1369  DO i = 1,n_vars
1370 
1371  IF ( abs(right_term(i)) .GE. 1.d-5 ) THEN
1372 
1373  coeff_f(i) = 1.d0 / right_term(i)
1374 
1375  END IF
1376 
1377  END DO
1378 
1379  END IF
1380 
1381  IF ( verbose_level .GE. 3 ) THEN
1382 
1383  CALL eval_f( qp_rel , qp_org , dz , coeff_f , right_term , scal_f )
1384  WRITE(*,*) 'before iter'
1385  WRITE(*,*) 'right_term',right_term
1386  WRITE(*,*) 'qp'
1387  WRITE(*,*) qp
1388 
1389  END IF
1390 
1391  DO nl_iter = 1,max_nl_iter
1392 
1393  CALL eval_jacobian( qp_rel , qp_org , dz , coeff_f , left_matrix )
1394 
1395  CALL eval_f( qp_rel , qp_org , dz , coeff_f , right_term , scal_f )
1396 
1397  IF ( nl_iter .EQ. 1 ) scal_f_init = scal_f
1398 
1399  delta_qp_rel = right_term
1400 
1401  call dgesv(n_eqns, 1, left_matrix , n_eqns, pivot, delta_qp_rel , n_eqns,&
1402  ok)
1403 
1404  IF ( ok .EQ. 0 ) THEN
1405 
1406  qp_rel_nr_old = qp_rel
1407  qp_nr_old = qp_rel * qp_org
1408  scal_f_old = scal_f
1409  desc_dir = - delta_qp_rel
1410 
1411  opt_search_nl = .true.
1412 
1413  IF ( ( opt_search_nl ) .AND. ( nl_iter .GT. 0 ) ) THEN
1414 
1415  stpmax = stpmx * max( dsqrt( dot_product(qp_rel,qp_rel) ) , &
1416  dble(SIZE(qp_rel)) )
1417 
1418  grad_f = matmul( right_term , left_matrix )
1419 
1420  CALL steady_lnsrch( qp_rel_nr_old , qp_org , scal_f_old , grad_f , &
1421  desc_dir , dz , coeff_f , qp_rel , scal_f , right_term , &
1422  stpmax , check , eval_f )
1423 
1424  IF ( check ) THEN
1425 
1426  desc_dir = - 0.5d0 * delta_qp_rel
1427 
1428  qp_rel = qp_rel_nr_old + desc_dir
1429 
1430  CALL eval_f( qp_rel , qp_org , dz , coeff_f , right_term , &
1431  scal_f )
1432 
1433  END IF
1434 
1435  DO i=idx_xd_first,idx_xd_last
1436 
1437  qp_rel(i) = max(qp_rel(i),0.d0)
1438 
1439  END DO
1440 
1441  DO i=idx_beta_first,idx_beta_last
1442 
1443  qp_rel(i) = max(qp_rel(i),0.d0)
1444 
1445  END DO
1446 
1447  ELSE
1448 
1449  qp_rel = qp_rel_nr_old + desc_dir
1450 
1451  DO i=idx_xd_first,idx_xd_last
1452 
1453  qp_rel(i) = max(qp_rel(i),0.d0)
1454 
1455  END DO
1456 
1457  DO i=idx_beta_first,idx_beta_last
1458 
1459  qp_rel(i) = max(qp_rel(i),0.d0)
1460 
1461  END DO
1462 
1463  CALL eval_f( qp_rel , qp_org , dz , coeff_f , right_term , scal_f )
1464 
1465  END IF
1466 
1467  ! Sometimes it happens that the crystal content becomes negative
1468  ! (even if very small, for example -1.e-140). Thus I force the
1469  ! crystal content and the dissolved gas to be non-negative.
1470 
1471 
1472  DO i=idx_xd_first,idx_xd_last
1473 
1474  qp_rel(i) = max(qp_rel(i),0.d0)
1475 
1476  END DO
1477 
1478  DO i=idx_beta_first,idx_beta_last
1479 
1480  qp_rel(i) = max(qp_rel(i),0.d0)
1481 
1482  END DO
1483 
1484  qp = qp_rel * qp_org
1485 
1486  arg_check = abs( qp_rel(1:n_vars) - qp_rel_nr_old(1:n_vars) ) / &
1487  ( tol_abs_nr + tol_rel_nr * max( abs(qp_rel(1:n_vars)) , &
1488  abs(qp_rel_nr_old(1:n_vars)) ) )
1489 
1490  check_nr_error = maxval( arg_check )
1491 
1492  idx_max = maxloc( arg_check )
1493 
1494  IF ( verbose_level .GE. 3 ) THEN
1495 
1496  WRITE(*,*) 'iter = ',nl_iter
1497  WRITE(*,*) 'right_term'
1498  WRITE(*,*) right_term
1499  WRITE(*,*) ''
1500  WRITE(*,*) 'desc_dir'
1501  WRITE(*,*) desc_dir
1502  WRITE(*,*) ''
1503  WRITE(*,*) 'qp_NR_old'
1504  WRITE(*,*) qp_nr_old
1505  WRITE(*,*) ''
1506  WRITE(*,*) 'qp'
1507  WRITE(*,*) qp
1508  WRITE(*,*) ''
1509 
1510  WRITE(*,*) 'scal_f = '
1511  WRITE(*,*) scal_f_init , scal_f , scal_f/scal_f_old , &
1512  scal_f/scal_f_init
1513  WRITE(*,*) ''
1514 
1515  WRITE(*,*) 'check_NR_error = '
1516  WRITE(*,*) check_nr_error, idx_max, &
1517  qp_rel(idx_max) * qp_org(idx_max), &
1518  qp_rel_nr_old(idx_max) * qp_org(idx_max)
1519  WRITE(*,*) ''
1520  WRITE(*,*) 'zeta = '
1521  WRITE(*,*) zeta
1522 
1523 
1524  READ(*,*)
1525 
1526  END IF
1527 
1528  ELSE
1529 
1530  IF ( verbose_level .GE. 3 ) THEN
1531 
1532  WRITE(*,*) 'advance_dz : num_cond too small'
1533  WRITE(*,*) 'left_matrix'
1534 
1535  DO j=1,n_eqns
1536  WRITE(*,*) left_matrix(j,:)
1537  END DO
1538 
1539  READ(*,*)
1540 
1541  END IF
1542 
1543  EXIT
1544 
1545  END IF
1546 
1547 
1548  IF ( qp_rel(1) .LT. 0.0d0 ) THEN
1549 
1550  check_convergence = .false.
1551 
1552  EXIT
1553 
1554  END IF
1555 
1556 
1557  IF ( check_nr_error .LT. 0.10d0 ) THEN
1558 
1559  IF ( scal_f / scal_f_init .LT. 1.d-10 ) check_convergence = .true.
1560 
1561  EXIT
1562 
1563  END IF
1564 
1565  IF ( scal_f / scal_f_init .LT. 1.d-14 ) THEN
1566 
1567  check_convergence = .true.
1568 
1569  EXIT
1570 
1571  END IF
1572 
1573  END DO
1574 
1575  IF ( scal_f / scal_f_init .LT. 1.d-11 ) check_convergence = .true.
1576 
1577  IF ( sum(qp_rel(idx_alfa_first:idx_alfa_last)) .LT. 0.0d0 ) THEN
1578 
1579  check_convergence = .false.
1580 
1581  END IF
1582 
1583  IF ( sum(qp_org(idx_alfa_first:idx_alfa_last) &
1584  * qp_rel(idx_alfa_first:idx_alfa_last)) .GE. 1.0d0 ) THEN
1585 
1586  check_convergence = .false.
1587 
1588  END IF
1589 
1590 
1591  IF ( verbose_level .GE. 1 ) THEN
1592 
1593  WRITE(*,*)''
1594  WRITE(*,*) 'scal_f = ' , nl_iter , idx_max , &
1595  scal_f/scal_f_init,check_nr_error, idx_max, ok
1596  WRITE(*,*)''
1597  WRITE(*,*)'qp_new = '
1598  WRITE(*,*) qp_rel * qp_org
1599 
1600  END IF
1601 
1602  END SUBROUTINE advance_dz
1603 
1604  !******************************************************************************
1606  !
1618  !******************************************************************************
1619 
1620  SUBROUTINE eval_f( qp_rel , qp_org , dz , coeff_f , right_term , scal_f )
1621 
1622  USE equations, ONLY : eval_fluxes_qp
1623  USE equations, ONLY : eval_nonhyperbolic_terms_qp
1624 
1625  USE parameters, ONLY : alfa_impl
1626 
1627  IMPLICIT NONE
1628 
1629  REAL*8, INTENT(IN) :: qp_rel(:)
1630  REAL*8, INTENT(IN) :: qp_org(:)
1631  REAL*8, INTENT(IN) :: dz
1632  REAL*8, INTENT(IN) :: coeff_f(:)
1633  REAL*8, INTENT(OUT) :: right_term(:)
1634  REAL*8, INTENT(OUT) :: scal_f
1635 
1636  REAL*8 :: qp(n_vars)
1637  REAL*8 :: fluxes(n_eqns)
1638  REAL*8 :: nh_terms(n_eqns)
1639 
1640 
1641  qp = qp_rel * qp_org
1642 
1643  CALL eval_fluxes_qp( r_qp=qp , r_flux=fluxes)
1644 
1645  CALL eval_nonhyperbolic_terms_qp( r_qp = qp , r_nh_term_impl = nh_terms )
1646  ! correction for isothermal model
1647 
1648 
1649  !WRITE(*,*) 'fluxes'
1650  !WRITE(*,*) fluxes
1651 
1652  !WRITE(*,*) 'nh_terms'
1653  !WRITE(*,*) nh_terms
1654  !READ(*,*)
1655 
1656  IF ( isothermal ) THEN
1657 
1658  fluxes(idx_mix_engy_eqn) = qp(idx_t)
1659  nh_terms(idx_mix_engy_eqn) = 0.d0
1660  right_term(idx_mix_engy_eqn) = fluxes(idx_mix_engy_eqn) &
1661  - fluxes_old(idx_mix_engy_eqn)
1662 
1663  END IF
1664 
1665  right_term = ( fluxes - fluxes_old ) - dz * ( alfa_impl * nh_terms &
1666  + ( 1.d0 - alfa_impl ) * nh_terms_old )
1667 
1668  right_term = right_term * coeff_f
1669 
1670  IF ( verbose_level .GE. 3 ) THEN
1671 
1672  WRITE(*,*) 'eval_f'
1673  WRITE(*,*) 'fluxes',REAL(fluxes)
1674  WRITE(*,*) 'fluxes_old',REAL(fluxes_old)
1675  WRITE(*,*) 'nh_terms',REAL(nh_terms)
1676  WRITE(*,*) 'nh_terms_old',REAL(nh_terms_old)
1677  READ(*,*)
1678 
1679  END IF
1680 
1681  scal_f = 0.5d0 * dot_product( right_term , right_term )
1682 
1683  END SUBROUTINE eval_f
1684 
1685  !******************************************************************************
1687  !
1699  !******************************************************************************
1700 
1701  SUBROUTINE eval_jacobian( qp_rel , qp_org , dz , coeff_f , left_matrix)
1702 
1703  ! external procedures
1704  USE equations, ONLY : eval_fluxes_qp
1705  USE equations, ONLY : eval_nonhyperbolic_terms_qp
1706 
1707  USE parameters, ONLY : alfa_impl
1708 
1709  IMPLICIT NONE
1710 
1711  REAL*8, INTENT(IN) :: qp_rel(n_vars)
1712  REAL*8, INTENT(IN) :: qp_org(n_vars)
1713  REAL*8, INTENT(IN) :: dz
1714  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1715  REAL*8, INTENT(OUT) :: left_matrix(n_eqns,n_eqns)
1716 
1717  REAL*8 :: h
1718 
1719  REAL*8 :: qp(n_vars)
1720  COMPLEX*16 :: qp_cmplx(n_eqns)
1721  COMPLEX*16 :: fluxes_cmplx(n_eqns)
1722  COMPLEX*16 :: nh_terms_cmplx(n_eqns)
1723 
1724  REAL*8 :: Jacob_fluxes(n_eqns,n_eqns)
1725  REAL*8 :: Jacob_nh(n_eqns,n_eqns)
1726 
1727  INTEGER :: i
1728 
1729  h = n_eqns * epsilon(1.d0)
1730 
1731  qp = qp_rel * qp_org
1732 
1733  qp_cmplx(1:n_eqns) = dcmplx( qp(1:n_eqns) , 0.d0 )
1734 
1735  DO i=1,n_vars
1736 
1737  qp_cmplx(i) = dcmplx( qp_rel(i) , h ) * qp_org(i)
1738 
1739  CALL eval_fluxes_qp( c_qp=qp_cmplx , c_flux=fluxes_cmplx )
1740  CALL eval_nonhyperbolic_terms_qp( c_qp = qp_cmplx , c_nh_term_impl = &
1741  nh_terms_cmplx )
1742 
1743  IF ( isothermal ) THEN
1744 
1745  fluxes_cmplx(idx_mix_engy_eqn) = qp_cmplx(idx_t)
1746  nh_terms_cmplx(idx_mix_engy_eqn) = dcmplx( 0.d0 , 0.d0)
1747 
1748  END IF
1749 
1750  jacob_fluxes(1:n_eqns,i) = dimag(fluxes_cmplx) / h
1751  jacob_nh(1:n_eqns,i) = dimag(nh_terms_cmplx) / h
1752 
1753  jacob_fluxes(1:n_eqns,i) = jacob_fluxes(1:n_eqns,i) * coeff_f
1754  jacob_nh(1:n_eqns,i) = jacob_nh(1:n_eqns,i) * coeff_f
1755 
1756  qp_cmplx(i) = dcmplx( qp(i) , 0.d0 )
1757 
1758  END DO
1759 
1760  left_matrix = jacob_fluxes - dz * alfa_impl * jacob_nh
1761 
1762 
1763  END SUBROUTINE eval_jacobian
1764 
1765  !******************************************************************************
1767  !
1786  !******************************************************************************
1787 
1788  SUBROUTINE steady_lnsrch( x_rel_init , x_org , scal_f_old , grad_f , desc_dir &
1789  , dz , coeff_f , x_rel_new , scal_f , f_nl , stpmax , check , callf )
1790 
1791  IMPLICIT NONE
1792 
1794  REAL*8, DIMENSION(:), INTENT(IN) :: x_rel_init
1795 
1797  REAL*8, DIMENSION(:), INTENT(IN) :: x_org
1798 
1800  REAL*8, DIMENSION(:), INTENT(IN) :: grad_f
1801 
1803  REAL*8, INTENT(IN) :: scal_f_old
1804 
1806  REAL*8, DIMENSION(:), INTENT(INOUT) :: desc_dir
1807 
1808  REAL*8, INTENT(IN) :: dz
1809 
1810  REAL*8, INTENT(IN) :: stpmax
1811 
1812  REAL*8, DIMENSION(:), INTENT(IN) :: coeff_f
1813 
1815  REAL*8, DIMENSION(:), INTENT(OUT) :: x_rel_new
1816 
1818  REAL*8, INTENT(OUT) :: scal_f
1819 
1821  REAL*8, DIMENSION(:), INTENT(OUT) :: f_nl
1822 
1824  LOGICAL, INTENT(OUT) :: check
1825 
1826  INTERFACE
1827  SUBROUTINE callf( x_rel , x_org , dz , coeff_f , f_nl , scal_f )
1828 
1829  IMPLICIT NONE
1830 
1831  REAL*8, INTENT(IN) :: x_rel(:)
1832  REAL*8, INTENT(IN) :: x_org(:)
1833  REAL*8, INTENT(IN) :: dz
1834  REAL*8, INTENT(IN) :: coeff_f(:)
1835  REAL*8, INTENT(OUT) :: f_nl(:)
1836  REAL*8, INTENT(OUT) :: scal_f
1837 
1838  END SUBROUTINE callf
1839  END INTERFACE
1840 
1841  REAL*8, PARAMETER :: TOLX=epsilon(x_rel_init)
1842 
1843  INTEGER, DIMENSION(1) :: ndum
1844  REAL*8 :: ALF , a,alam,alam2,alamin,b,disc
1845  REAL*8 :: scal_f2
1846  REAL*8 :: desc_dir_abs
1847  REAL*8 :: rhs1 , rhs2 , slope, tmplam
1848 
1849  alf = 1.0d-4
1850 
1851  IF ( size(grad_f) == size(desc_dir) .AND. size(grad_f) == size(x_rel_new) &
1852  .AND. size(x_rel_new) == size(x_rel_init) ) THEN
1853 
1854  ndum = size(grad_f)
1855 
1856  ELSE
1857 
1858  WRITE(*,*) 'nrerror: an assert_eq failed with this tag:', 'lnsrch'
1859  stop 'program terminated by assert_eq4'
1860 
1861  END IF
1862 
1863  check = .false.
1864 
1865  desc_dir_abs = dsqrt( dot_product(desc_dir,desc_dir) )
1866 
1867  IF ( desc_dir_abs > stpmax ) desc_dir(:) = desc_dir(:) * stpmax / &
1868  desc_dir_abs
1869 
1870  slope = dot_product(grad_f,desc_dir)
1871 
1872  alamin = tolx / maxval( abs( desc_dir(:)) / max( abs( x_rel_init(:)) , &
1873  1.d0 ) )
1874 
1875  alamin = 1.d-20
1876 
1877  IF ( alamin .EQ. 0.d0) THEN
1878 
1879  x_rel_new(:) = x_rel_init(:)
1880 
1881  WRITE(*,*) 'alam 0'
1882 
1883  RETURN
1884 
1885  END IF
1886 
1887  alam = 1.0d0
1888  alam2 = alam
1889 
1890  optimal_step_search: DO
1891 
1892  IF ( verbose_level .GE. 4 ) THEN
1893 
1894  WRITE(*,*) 'alam',alam
1895 
1896  END IF
1897 
1898  x_rel_new = x_rel_init + alam * desc_dir
1899 
1900  CALL callf( x_rel_new , x_org , dz , coeff_f , f_nl , scal_f )
1901 
1902  scal_f2 = scal_f
1903 
1904  IF ( verbose_level .GE. 4 ) THEN
1905 
1906  WRITE(*,*) 'x',x_rel_new,x_rel_init
1907  WRITE(*,*) 'lnsrch: effe_old,effe',scal_f_old,scal_f
1908  READ(*,*)
1909 
1910  END IF
1911 
1912  IF ( alam < alamin ) THEN
1913  ! convergence on Delta_x
1914 
1915  IF ( verbose_level .GE. 4 ) THEN
1916 
1917  WRITE(*,*) ' convergence on Delta_x',alam,alamin
1918 
1919  END IF
1920 
1921  x_rel_new(:) = x_rel_init(:)
1922  scal_f = scal_f_old
1923  check = .true.
1924 
1925  EXIT optimal_step_search
1926 
1927  ELSE IF ( scal_f .LE. scal_f_old + alf * alam * slope ) THEN
1928 
1929  EXIT optimal_step_search
1930 
1931  ELSE
1932 
1933  IF ( alam .EQ. 1.d0 ) THEN
1934 
1935  tmplam = - slope / ( 2.0d0 * ( scal_f - scal_f_old - slope ) )
1936 
1937  ELSE
1938 
1939  rhs1 = scal_f - scal_f_old - alam*slope
1940  rhs2 = scal_f2 - scal_f_old - alam2*slope
1941 
1942  a = ( rhs1/alam**2.d0 - rhs2/alam2**2.d0 ) / ( alam - alam2 )
1943  b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )
1944 
1945  IF ( a .EQ. 0.d0 ) THEN
1946 
1947  tmplam = - slope / ( 2.0d0 * b )
1948 
1949  ELSE
1950 
1951  disc = b*b - 3.0d0*a*slope
1952 
1953  IF ( disc .LT. 0.d0 ) THEN
1954 
1955  tmplam = 0.5d0 * alam
1956 
1957  ELSE IF ( b .LE. 0.d0 ) THEN
1958 
1959  tmplam = ( - b + dsqrt(disc) ) / ( 3.d0 * a )
1960 
1961  ELSE
1962 
1963  tmplam = - slope / ( b + dsqrt(disc) )
1964 
1965  ENDIF
1966 
1967  END IF
1968 
1969  IF ( tmplam .GT. 0.5d0*alam ) tmplam = 0.5d0 * alam
1970 
1971  END IF
1972 
1973  END IF
1974 
1975  alam2 = alam
1976  scal_f2 = scal_f
1977  alam = max( tmplam , 0.1d0*alam )
1978 
1979  END DO optimal_step_search
1980 
1981  END SUBROUTINE steady_lnsrch
1982 
1983  !******************************************************************************
1985  !
2000  !******************************************************************************
2001 
2002  SUBROUTINE linear_extrapolation(zeta_old,qp_old,zeta,qp,extrap_z,extrap_z_p, &
2003  extrap_z_mach,extrap_flag)
2004 
2005  ! external procedures
2006  USE constitutive, ONLY : eos , alfa_2
2007  USE constitutive, ONLY : sound_speeds
2008  USE equations, ONLY : phys_var_qp
2009 
2010  ! external variables
2011  USE constitutive, ONLY : x_d_md, x_ex_dis_in
2012  USE init, ONLY : p1_in , p2_in
2013 
2014  IMPLICIT NONE
2015 
2016  REAL*8, INTENT(IN) :: zeta_old
2017  REAL*8, INTENT(IN) :: qp_old(n_eqns)
2018  REAL*8, INTENT(IN) :: zeta
2019  REAL*8, INTENT(IN) :: qp(n_eqns)
2020 
2021  REAL*8, INTENT(OUT) :: extrap_z
2022  REAL*8, INTENT(OUT) :: extrap_z_p
2023  REAL*8, INTENT(OUT) :: extrap_z_mach
2024  INTEGER, INTENT(OUT) :: extrap_flag
2025 
2026  REAL*8 :: r_p_2 , r_p_2_old
2027  REAL*8 :: r_p_1 , r_p_1_old
2028  REAL*8 :: mach , mach_old
2029 
2030  REAL*8 :: C_mix
2031 
2032  REAL*8 :: p_check
2033 
2034  REAL*8 :: extrap_z_2
2035  REAL*8 :: extrap_z_1
2036 
2037  REAL*8 :: coeff_p_1 , coeff_p_2
2038 
2039  REAL*8 :: p_out_1 , p_out_2
2040 
2041  REAL*8 :: eps_extrap
2042 
2043  REAL*8 :: r_alfa_2
2044 
2045  eps_extrap = 1.d-10
2046 
2047  extrap_flag = 0
2048 
2049  !--- check on small values of the phisical variables ------------------------
2050 
2051 
2052  r_p_1 = qp(idx_p1)
2053  r_p_2 = qp(idx_p2)
2054 
2055  IF ( qp(n_vars) .GT. 0.d0 ) THEN
2056 
2057  p_out_1 = p_out
2058  p_out_2 = p_out
2059 
2060  ELSE
2061 
2062  p_out_1 = p_out
2063  p_out_2 = p_out
2064 
2065  END IF
2066 
2067  CALL phys_var_qp( r_qp = qp )
2068  CALL eos
2069  CALL sound_speeds(c_mix,mach)
2070 
2071 
2072 
2073  IF ( verbose_level .GT. 1 ) THEN
2074 
2075  WRITE(*,*) 'p1,p2,mach',r_p_1,r_p_2,mach
2076 
2077  END IF
2078 
2079  ! ---- Check if the solution at zeta has already reached a boundary condition
2080 
2081  p_check = min(r_p_2,r_p_1)
2082 
2083  r_alfa_2 = REAL(alfa_2)
2084 
2085  IF ( ( r_p_1 .LT. p_out_1 ) .OR. ( r_p_2 .LT. p_out_2 ) .OR. &
2086  ( mach .GT. 1.d0 ) ) THEN
2087 
2088  IF ( verbose_level .GE. -1 ) THEN
2089 
2090  WRITE(*,*) 'Boundary conditions before the top'
2091  WRITE(*,*) 'z = ',zeta,'Pressures = ',r_p_1,r_p_2
2092  WRITE(*,*) 'Mach =',mach
2093 
2094  END IF
2095 
2096  extrap_z = zeta
2097  extrap_flag = 4
2098 
2099  RETURN
2100 
2101  END IF
2102 
2103  IF ( r_alfa_2 .GE. 1.d0 ) THEN
2104 
2105  IF ( verbose_level .GE. -1 ) THEN
2106 
2107  WRITE(*,*) 'Boundary conditions before the top'
2108  WRITE(*,*) 'z = ',zeta,'alfa gas = ',r_alfa_2
2109 
2110  END IF
2111 
2112  extrap_z = zeta
2113  extrap_flag = 4
2114 
2115  RETURN
2116 
2117  END IF
2118 
2119  IF ( (abs(zeta - zeta_old) .LT. 1e-11) .AND. &
2120  (r_p_1 .LT. 1.0e+6) ) THEN
2121 
2122  IF ( verbose_level .GE. -1 ) THEN
2123 
2124  WRITE(*,*) 'Pressure Conditions reached before the exit'
2125  !READ(*,*)
2126 
2127  END IF
2128 
2129  extrap_z = zeta
2130  extrap_flag = 5
2131 
2132  RETURN
2133 
2134  END IF
2135 
2136  IF ( (abs(zeta - zeta_old) .LT. 1e-6) .AND. &
2137  ( sum(REAL(x_d_md)) .LT. sum(x_ex_dis_in) ) ) then
2138 
2139  IF ( verbose_level .GE. -1 ) THEN
2140 
2141  WRITE(*,*) 'Dissolved gas Conditions reached before the exit'
2142  !READ(*,*)
2143 
2144  END IF
2145 
2146  extrap_z = zeta
2147  extrap_flag = 6
2148 
2149  RETURN
2150 
2151  END IF
2152 
2153 
2154 
2155  ! --------------------------------------------------------------------------
2156 
2157  r_p_1_old = qp_old(idx_p1)
2158  r_p_2_old = qp_old(idx_p2)
2159 
2160  CALL phys_var_qp( r_qp = qp_old )
2161  CALL eos
2162  CALL sound_speeds(c_mix,mach_old)
2163 
2164  ! ---- gas pressure extrapolation -------------------------------------------
2165  IF ( r_p_2 - r_p_2_old .LT. 0.d0 ) THEN
2166 
2167  extrap_z_2 = zeta_old - ( zeta - zeta_old ) / ( r_p_2 - r_p_2_old) * &
2168  ( r_p_2_old )
2169 
2170  IF ( verbose_level .GE. 1 ) THEN
2171 
2172  WRITE(*,*) 'gas pressure extrapolation'
2173  WRITE(*,*) 'extrap_z_2 =',extrap_z_2,r_p_2
2174 
2175  END IF
2176 
2177  ELSE
2178 
2179  extrap_z_2 = zn
2180 
2181  END IF
2182 
2183  ! --- liquid pressure extrapolation -----------------------------------------
2184  IF ( r_p_1 - r_p_1_old .LT. 0.d0 ) THEN
2185 
2186  extrap_z_1 = zeta_old + ( zeta - zeta_old ) / ( r_p_1 - r_p_1_old) * &
2187  ( p_out_1 - r_p_1_old )
2188 
2189  IF ( verbose_level .GE. 1 ) THEN
2190 
2191  WRITE(*,*) 'liquid pressure extrapolation'
2192  WRITE(*,*) 'extrap_z_1 =',extrap_z_1,r_p_1
2193 
2194  END IF
2195 
2196  ELSE
2197 
2198  extrap_z_1 = zn
2199 
2200  END IF
2201 
2202  extrap_z_p = min( extrap_z_1 , extrap_z_2 )
2203 
2204  ! --- Mach number extrapolation ---------------------------------------------
2205  IF ( mach - mach_old .GT. 0.d0 ) THEN
2206 
2207  extrap_z_mach = zeta_old + ( zeta - zeta_old ) / ( mach - mach_old ) * &
2208  ( 1.d0 - mach_old )
2209 
2210  IF ( verbose_level .GE. 1 ) THEN
2211 
2212  WRITE(*,*) 'Mach number extrapolation'
2213  WRITE(*,*) 'extrap_z_mach',extrap_z_mach,mach
2214 
2215  END IF
2216 
2217  ELSE
2218 
2219  extrap_z_mach = zn
2220 
2221  END IF
2222 
2223  ! ---------------------------------------------------------------------------
2224  ! If the top is reached without fulfilling the boundary conditions then save
2225  ! the extrapolated point at which the condition is satisfied
2226 
2227  IF ( zeta .EQ. zn ) THEN
2228 
2229  IF ( extrap_z_1 .GT. zn ) extrap_z_p = extrap_z_1
2230 
2231  IF ( extrap_z_2 .GT. zn ) extrap_z_p = min( extrap_z_p , extrap_z_2 )
2232 
2233  extrap_z = extrap_z_p
2234 
2235  IF ( extrap_z_mach .GT. zn ) THEN
2236 
2237  extrap_z = min( extrap_z_p , extrap_z_mach )
2238  extrap_flag = -3
2239 
2240  END IF
2241 
2242  RETURN
2243 
2244  END IF
2245 
2246  ! --- Check on the extrapolated values --------------------------------------
2247 
2248  IF ( ( extrap_z_mach - zeta ) * ( 1.d0 - mach ) .LT. eps_extrap ) THEN
2249 
2250  IF ( extrap_z_mach .LT. zn ) THEN
2251 
2252  extrap_z = extrap_z_mach
2253  extrap_flag = 3
2254 
2255  RETURN
2256 
2257  END IF
2258 
2259  END IF
2260 
2261  coeff_p_1 = ( p1_in - r_p_1 ) / ( p1_in - p_out_1 )
2262 
2263  IF ( ( extrap_z_1 - zeta ) * ( 1.d0 - coeff_p_1 ) .LT. eps_extrap ) THEN
2264 
2265  IF ( extrap_z_1 .LT. zn ) THEN
2266 
2267  extrap_z = extrap_z_1
2268  extrap_flag = 1
2269 
2270  RETURN
2271 
2272  END IF
2273 
2274  END IF
2275 
2276  coeff_p_2 = ( r_p_2 - p2_in ) / ( p_out_2 - p2_in )
2277 
2278  IF ( ( extrap_z_2 - zeta ) * ( 1.d0 - coeff_p_2 ) .LT. eps_extrap ) THEN
2279 
2280  IF ( extrap_z_2 .LT. zn ) THEN
2281 
2282  extrap_z = extrap_z_2
2283  extrap_flag = 2
2284 
2285  RETURN
2286 
2287  END IF
2288 
2289  END IF
2290 
2291  END SUBROUTINE linear_extrapolation
2292 
2293 END MODULE steady_solver
constitutive::x_ex_dis_in
real *8, dimension(:), allocatable x_ex_dis_in
Definition: constitutive.f90:300
constitutive::u_1
complex *16 u_1
melt-crystals phase local velocity
Definition: constitutive.f90:148
equations::isothermal
logical isothermal
Flag for isothermal runs: .
Definition: equations.f90:33
parameters::idx_p2
integer idx_p2
Index of p2 in the qp array.
Definition: parameters.f90:38
parameters
Parameters.
Definition: parameters.f90:10
parameters::idx_cry_eqn_last
integer idx_cry_eqn_last
Definition: parameters.f90:59
constitutive::s_1
complex *16 s_1
local specific entropy of the melt-crystals phase
Definition: constitutive.f90:104
parameters::idx_ex_gas_eqn_last
integer idx_ex_gas_eqn_last
Definition: parameters.f90:57
parameters::idx_mix_mom_eqn
integer idx_mix_mom_eqn
Definition: parameters.f90:51
constitutive::rho_g
complex *16, dimension(:), allocatable rho_g
exsolved gas local density
Definition: constitutive.f90:123
constitutive::rho_c
complex *16, dimension(:), allocatable rho_c
crystals local density
Definition: constitutive.f90:121
parameters::alfa_impl
real *8, parameter alfa_impl
Parameter for numerical scheme: .
Definition: parameters.f90:91
steady_solver::zeta
real *8 zeta
Definition: steady_solver.f90:40
constitutive::mu_2
complex *16 mu_2
free Gibbs energy of the exsolved gas
Definition: constitutive.f90:112
parameters::idx_u2
integer idx_u2
Index of u2 in the qp array.
Definition: parameters.f90:40
constitutive::rho_m
complex *16 rho_m
melt local density
Definition: constitutive.f90:120
parameters::n_gas
integer n_gas
Numbeer of crystal phases.
Definition: parameters.f90:30
parameters::eps_conv
real *8 eps_conv
Residual for the convergence of the shooting method. The solution is accepted if one of these conditi...
Definition: parameters.f90:122
parameters::tol_abs
real *8, parameter tol_abs
Definition: parameters.f90:98
constitutive::explosive
logical explosive
Flag to choose the eruptive style: .
Definition: constitutive.f90:380
constitutive::e_m
complex *16 e_m
local specific internal energy of the melt
Definition: constitutive.f90:96
parameters::shooting
logical shooting
Flag for the shooting technique: .
Definition: parameters.f90:107
steady_solver::steady_shooting
subroutine steady_shooting
Shooting Method.
Definition: steady_solver.f90:65
constitutive::alfa_2
complex *16 alfa_2
total exsolved gas local volume fraction
Definition: constitutive.f90:130
constitutive::u_2
complex *16 u_2
exsolved gas local velocity
Definition: constitutive.f90:149
parameters::idx_ex_gas_eqn_first
integer idx_ex_gas_eqn_first
Definition: parameters.f90:56
constitutive::f_alfa3
subroutine f_alfa3(r_p_2, xtot, r_beta, r_rho_md, r_rho_g, r_alfa_g)
Bottom exsolved gas.
Definition: constitutive.f90:2103
constitutive::rho_mix
complex *16 rho_mix
mixture local density
Definition: constitutive.f90:169
init::p_out
real *8 p_out
Outlet pressure.
Definition: init.f90:38
geometry::comp_cells
integer comp_cells
Number of control volumes in the computational domain.
Definition: geometry.f90:44
equations::eval_fluxes_qp
subroutine eval_fluxes_qp(c_qp, r_qp, c_flux, r_flux)
Fluxes.
Definition: equations.f90:293
constitutive::e_g
complex *16, dimension(:), allocatable e_g
local specific internal energy of the exsolved gas
Definition: constitutive.f90:99
steady_solver::integrate_equations
subroutine integrate_equations(qp, flag_output, extrap_z, extrap_z_p, extrap_z_mach, extrap_flag, r_p_1, r_p_2, mach)
Steady state integration.
Definition: steady_solver.f90:625
constitutive::frag_thr
real *8 frag_thr
Threshold for the fragmentation.
Definition: constitutive.f90:268
geometry::update_radius
subroutine update_radius(zeta)
Definition: geometry.f90:166
steady_solver::fluxes_old
real *8, dimension(:), allocatable fluxes_old
Definition: steady_solver.f90:41
equations::lateral_degassing
logical lateral_degassing
Flag for lateral degassing: .
Definition: equations.f90:48
inpout
Input/Output module.
Definition: inpout.f90:11
equations::eval_nonhyperbolic_terms_qp
subroutine eval_nonhyperbolic_terms_qp(c_qp, c_nh_term_impl, r_qp, r_nh_term_impl)
Non-Hyperbolic terms.
Definition: equations.f90:491
constitutive::p_2
complex *16 p_2
local pressure of the exsolved gas
Definition: constitutive.f90:102
steady_solver::eval_f
subroutine eval_f(qp_rel, qp_org, dz, coeff_f, right_term, scal_f)
Nonlinear function evaluation.
Definition: steady_solver.f90:1621
constitutive::e_2
complex *16 e_2
local specific internal energy of the exsolved gas
Definition: constitutive.f90:95
parameters::idx_alfa_first
integer idx_alfa_first
First index of alfa in the qp array.
Definition: parameters.f90:44
steady_solver::steady_lnsrch
subroutine steady_lnsrch(x_rel_init, x_org, scal_f_old, grad_f, desc_dir, dz, coeff_f, x_rel_new, scal_f, f_nl, stpmax, check, callf)
Search the descent stepsize.
Definition: steady_solver.f90:1790
equations::lateral_degassing_flag
logical lateral_degassing_flag
Input flag for lateral degassing: .
Definition: equations.f90:42
parameters::idx_mix_engy_eqn
integer idx_mix_engy_eqn
Definition: parameters.f90:53
constitutive::s_g
complex *16, dimension(:), allocatable s_g
local specific entropy of the exsolved gas
Definition: constitutive.f90:109
parameters::tol_rel
real *8, parameter tol_rel
Definition: parameters.f90:99
constitutive::rho_2
complex *16 rho_2
exsolved gas local density
Definition: constitutive.f90:119
geometry::z0
real *8 z0
Left (bottom) of the physical domain.
Definition: geometry.f90:20
parameters::max_nl_iter
integer, parameter max_nl_iter
Maximum iterations of the Newthon-Raphson solver.
Definition: parameters.f90:94
constitutive::f_alfa
subroutine f_alfa(xtot, xmax, r_beta, r_rho_md, r_rho_2, r_alfa_2)
Bottom exsolved gas.
Definition: constitutive.f90:2060
constitutive::sound_speeds
subroutine sound_speeds(C_mix, mach)
Local sound speeds.
Definition: constitutive.f90:647
steady_solver::perturbe_qp
subroutine perturbe_qp(qp)
Solution preturbation.
Definition: steady_solver.f90:1240
constitutive::mu_m
complex *16 mu_m
free Gibbs energy of the melt
Definition: constitutive.f90:113
steady_solver::nl_iter
integer nl_iter
Definition: steady_solver.f90:48
constitutive::u_mix
complex *16 u_mix
mixture velocity
Definition: constitutive.f90:170
constitutive::s_2
complex *16 s_2
local specific entropy of the exsolved gas
Definition: constitutive.f90:105
parameters::idx_mix_mass_eqn
integer idx_mix_mass_eqn
Definition: parameters.f90:49
parameters::idx_alfa_last
integer idx_alfa_last
Last index of alfa in the qp array.
Definition: parameters.f90:45
geometry::radius
real *8 radius
Effective radius.
Definition: geometry.f90:27
steady_solver::eval_jacobian
subroutine eval_jacobian(qp_rel, qp_org, dz, coeff_f, left_matrix)
Jacobian evaluation.
Definition: steady_solver.f90:1702
constitutive::s_m
complex *16 s_m
local specific entropy of the melt
Definition: constitutive.f90:106
parameters::idx_rel_vel_eqn
integer idx_rel_vel_eqn
Definition: parameters.f90:52
equations
Governing equations.
Definition: equations.f90:9
constitutive::zeta_lith
real *8 zeta_lith
Elevation above the bottom for the evaluation of the lithostatic pressure.
Definition: constitutive.f90:358
parameters::idx_beta_first
integer idx_beta_first
First index of beta in the qp array.
Definition: parameters.f90:46
constitutive::s_c
complex *16, dimension(:), allocatable s_c
local specific entropy of the crystals
Definition: constitutive.f90:107
geometry::zeta_exit
real *8 zeta_exit
Right (top) of the physical domain.
Definition: geometry.f90:22
steady_solver::linear_extrapolation
subroutine linear_extrapolation(zeta_old, qp_old, zeta, qp, extrap_z, extrap_z_p, extrap_z_mach, extrap_flag)
Linear extrapolation of the solution.
Definition: steady_solver.f90:2004
init
Initial solution.
Definition: init.f90:11
init::u1_in
real *8 u1_in
Inlet first phase velocity.
Definition: init.f90:31
constitutive::frag_eff
real *8 frag_eff
index of fragmentation in the interval [0;1]
Definition: constitutive.f90:258
parameters::idx_dis_gas_eqn_first
integer idx_dis_gas_eqn_first
Definition: parameters.f90:54
steady_solver::nh_terms_old
real *8, dimension(:), allocatable nh_terms_old
Definition: steady_solver.f90:42
geometry::z_comp
real *8, dimension(:), allocatable z_comp
Location of the centers of the control volume of the domain.
Definition: geometry.f90:12
parameters::idx_xd_first
integer idx_xd_first
First index of xd in the qp array.
Definition: parameters.f90:42
parameters::idx_p1
integer idx_p1
Index of p1 in the qp array.
Definition: parameters.f90:37
geometry::zn
real *8 zn
Right (top) of the physical domain.
Definition: geometry.f90:21
parameters::idx_u1
integer idx_u1
Index of u1 in the qp array.
Definition: parameters.f90:39
constitutive
Constitutive equations.
Definition: constitutive.f90:9
parameters::idx_cry_eqn_first
integer idx_cry_eqn_first
Definition: parameters.f90:58
parameters::idx_beta_last
integer idx_beta_last
Last index of beta in the qp array.
Definition: parameters.f90:47
constitutive::e_c
complex *16, dimension(:), allocatable e_c
local specific internal energy of the crystals
Definition: constitutive.f90:97
constitutive::e_1
complex *16 e_1
local specific internal energy of the melt-crystals phase
Definition: constitutive.f90:94
parameters::idx_vol1_eqn
integer idx_vol1_eqn
Definition: parameters.f90:50
init::p1_in
real *8 p1_in
Inlet first phase pressure.
Definition: init.f90:28
parameters::idx_dis_gas_eqn_last
integer idx_dis_gas_eqn_last
Definition: parameters.f90:55
constitutive::rho_1
complex *16 rho_1
dis_gas+melt+crystals phase local density
Definition: constitutive.f90:118
steady_solver::zeta_old
real *8 zeta_old
Definition: steady_solver.f90:44
steady_solver
Steady state equations integration.
Definition: steady_solver.f90:11
parameters::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters.f90:32
init::p2_in
real *8 p2_in
Inlet second phase pressure.
Definition: init.f90:29
equations::alfa2_lat_thr
real *8 alfa2_lat_thr
Exsolved gas volume fraction threshold for lateral degassing.
Definition: equations.f90:52
parameters::verbose_level
integer verbose_level
Definition: parameters.f90:101
steady_solver::advance_dz
subroutine advance_dz(qp, dz, check_convergence)
Solution advance in space.
Definition: steady_solver.f90:1271
constitutive::mu_1
complex *16 mu_1
free Gibbs energy of the melt-crystals phase
Definition: constitutive.f90:111
constitutive::p_1
complex *16 p_1
local pressure of the melt-crystals phase
Definition: constitutive.f90:101
inpout::output_steady
subroutine output_steady(zeta, qp, radius)
Write the steady solution on the output unit.
Definition: inpout.f90:1204
constitutive::eval_densities
subroutine eval_densities
Phases densities.
Definition: constitutive.f90:729
constitutive::t
complex *16 t
mixture local temperature
Definition: constitutive.f90:168
equations::phys_var_qp
subroutine phys_var_qp(r_qp, c_qp)
Physical variables.
Definition: equations.f90:117
constitutive::x_d_md
complex *16, dimension(:), allocatable x_d_md
dissolved gas mass fraction in the melt+dis.gas phase
Definition: constitutive.f90:158
geometry
Grid module.
Definition: geometry.f90:7
constitutive::eos
subroutine eos
Equation of state.
Definition: constitutive.f90:573
constitutive::mu_g
complex *16, dimension(:), allocatable mu_g
free Gibbs energy of the exsolved gas
Definition: constitutive.f90:116
parameters::idx_t
integer idx_t
Index of T in the qp array.
Definition: parameters.f90:41
steady_solver::increase_flow_rate
logical increase_flow_rate
Definition: steady_solver.f90:46
parameters::idx_xd_last
integer idx_xd_last
Last index of xd in the qp array.
Definition: parameters.f90:43
constitutive::mu_c
complex *16, dimension(:), allocatable mu_c
free Gibbs energy of the crystals
Definition: constitutive.f90:114
init::init_steady
subroutine init_steady(u_0, qp)
Steady problem initialization.
Definition: init.f90:53
steady_solver::dz_max
real *8 dz_max
Definition: steady_solver.f90:44
parameters::n_eqns
integer n_eqns
Number of equations.
Definition: parameters.f90:33
Generated on Tue Mar 7 2017 09:55:58 for MAMMA by   doxygen 1.8.11