IMEXSfloW  1.0
templategithubproject
 All Classes Files Functions Variables Pages
solver.f90
Go to the documentation of this file.
1 !*******************************************************************************
3 !
6 !
10 !
11 !********************************************************************************
12 MODULE solver
13 
14  ! external variables
15 
16  USE constitutive, ONLY : implicit_flag
17 
18  USE geometry, ONLY : comp_cells
19  USE geometry, ONLY : comp_interfaces
20 
21  USE geometry, ONLY : b_cent , b_prime , b_stag
22 
23  USE parameters, ONLY : n_eqns , n_vars , n_nh
24  USE parameters, ONLY : n_rk
25  USE parameters, ONLY : reconstr_variables
26  USE parameters, ONLY : verbose_level
27 
28  USE parameters, ONLY : bcl , bcr
29 
30  IMPLICIT none
31 
33  REAL*8, ALLOCATABLE :: q(:,:)
34 
35  REAL*8, ALLOCATABLE :: q0(:,:)
36 
37  REAL*8, ALLOCATABLE :: q_fv(:,:)
38 
39  REAL*8, ALLOCATABLE :: q_interfaceL(:,:)
40 
41  REAL*8, ALLOCATABLE :: q_interfaceR(:,:)
42 
43  REAL*8, ALLOCATABLE :: a_interfaceL(:,:)
44 
45  REAL*8, ALLOCATABLE :: a_interfaceR(:,:)
46 
47  REAL*8, ALLOCATABLE :: H_interface(:,:)
48 
49  REAL*8, ALLOCATABLE :: qp(:,:)
50 
51 
53  REAL*8 :: dt
54 
55  LOGICAL, ALLOCATABLE :: mask22(:,:) , mask21(:,:) , mask11(:,:) , mask12(:,:)
56 
58  REAL*8, ALLOCATABLE :: a_tilde_ij(:,:)
59 
60  REAL*8, ALLOCATABLE :: a_dirk_ij(:,:)
61 
63  REAL*8, ALLOCATABLE :: omega_tilde(:)
64 
65  REAL*8, ALLOCATABLE :: omega(:)
66 
68  REAL*8, ALLOCATABLE :: a_tilde(:)
69 
70  REAL*8, ALLOCATABLE :: a_dirk(:)
71 
72  REAL*8 :: a_diag
73 
75  REAL*8, ALLOCATABLE :: q_rk(:,:,:)
76 
77  REAL*8, ALLOCATABLE :: F_x(:,:,:)
78 
79  REAL*8, ALLOCATABLE :: NH(:,:,:)
80 
82  REAL*8, ALLOCATABLE :: expl_terms(:,:,:)
83 
85  REAL*8, ALLOCATABLE :: Fxj(:,:)
86 
87  REAL*8, ALLOCATABLE :: NHj(:,:)
88 
89  REAL*8, ALLOCATABLE :: expl_terms_j(:,:)
90 
91  LOGICAL :: normalize_q
92  LOGICAL :: normalize_f
93  LOGICAL :: opt_search_NL
94 
95  REAL*8, ALLOCATABLE :: residual_term(:,:)
96 
97 
98 CONTAINS
99 
100  !*****************************************************************************
102  !
105  !
109  !
110  !*****************************************************************************
111 
112  SUBROUTINE allocate_solver_variables
113 
114  IMPLICIT NONE
115 
116  REAL*8 :: gamma
117 
118  INTEGER :: i,j
119 
120  ALLOCATE( q( n_vars , comp_cells ) , q0( n_vars , comp_cells ) )
121 
122  ALLOCATE( q_fv( n_vars , comp_cells ) )
123 
124  ALLOCATE( q_interfacel( n_vars , comp_interfaces ) )
125  ALLOCATE( q_interfacer( n_vars , comp_interfaces ) )
126 
127  ALLOCATE( a_interfacel( n_eqns , comp_interfaces ) )
128  ALLOCATE( a_interfacer( n_eqns , comp_interfaces ) )
129 
130  ALLOCATE( h_interface( n_eqns , comp_interfaces ) )
131 
132  ALLOCATE( qp( n_vars , comp_cells ) )
133 
134  ALLOCATE( a_tilde_ij(n_rk,n_rk) )
135  ALLOCATE( a_dirk_ij(n_rk,n_rk) )
136  ALLOCATE( omega_tilde(n_rk) )
137  ALLOCATE( omega(n_rk) )
138 
139  ALLOCATE( mask22(n_eqns,n_eqns) )
140  ALLOCATE( mask21(n_eqns,n_eqns) )
141  ALLOCATE( mask11(n_eqns,n_eqns) )
142  ALLOCATE( mask12(n_eqns,n_eqns) )
143 
144  mask11(1:n_eqns,1:n_eqns) = .false.
145  mask12(1:n_eqns,1:n_eqns) = .false.
146  mask22(1:n_eqns,1:n_eqns) = .false.
147  mask21(1:n_eqns,1:n_eqns) = .false.
148 
149  DO i = 1,n_eqns
150 
151  DO j = 1,n_eqns
152 
153  IF ( .NOT.implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
154  mask11(j,i) = .true.
155  IF ( implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
156  mask12(j,i) = .true.
157  IF ( implicit_flag(i) .AND. implicit_flag(j) ) &
158  mask22(j,i) = .true.
159  IF ( .NOT.implicit_flag(i) .AND. implicit_flag(j) ) &
160  mask21(j,i) = .true.
161 
162  END DO
163 
164  END DO
165 
166 
167  a_tilde_ij = 0.d0
168  a_dirk_ij = 0.d0
169  omega_tilde = 0.d0
170  omega = 0.d0
171 
172  gamma = 1.d0 - 1.d0 / sqrt(2.d0)
173 
174  IF ( n_rk .EQ. 1 ) THEN
175 
176  a_tilde_ij(1,1) = 1.d0
177 
178  omega_tilde(1) = 1.d0
179 
180  a_dirk_ij(1,1) = 0.d0
181 
182  omega(1) = 0.d0
183 
184  ELSEIF ( n_rk .EQ. 2 ) THEN
185 
186  a_tilde_ij(2,1) = 1.0d0
187 
188  omega_tilde(1) = 1.0d0
189  omega_tilde(2) = 0.0d0
190 
191  a_dirk_ij(2,2) = 1.0d0
192 
193  omega(1) = 0.d0
194  omega(2) = 1.d0
195 
196  ELSEIF ( n_rk .EQ. 3 ) THEN
197 
198  a_tilde_ij(2,1) = 0.5d0
199  a_tilde_ij(3,1) = 0.5d0
200  a_tilde_ij(3,2) = 0.5d0
201 
202  omega_tilde(1) = 1.0d0 / 3.0d0
203  omega_tilde(2) = 1.0d0 / 3.0d0
204  omega_tilde(3) = 1.0d0 / 3.0d0
205 
206  a_dirk_ij(1,1) = 0.25d0
207  a_dirk_ij(2,2) = 0.25d0
208  a_dirk_ij(3,1) = 1.0d0 / 3.0d0
209  a_dirk_ij(3,2) = 1.0d0 / 3.0d0
210  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
211 
212  omega(1) = 1.0d0 / 3.0d0
213  omega(2) = 1.0d0 / 3.0d0
214  omega(3) = 1.0d0 / 3.0d0
215 
216  ELSEIF ( n_rk .EQ. 4 ) THEN
217 
218  ! LRR(3,2,2)
219 
220  a_tilde_ij(2,1) = 0.5d0
221  a_tilde_ij(3,1) = 1.d0 / 3.d0
222  a_tilde_ij(4,2) = 1.0d0
223 
224  omega_tilde(1) = 0.d0
225  omega_tilde(2) = 1.0d0
226  omega_tilde(3) = 0.0d0
227  omega_tilde(4) = 0.d0
228 
229  a_dirk_ij(2,2) = 0.5d0
230  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
231  a_dirk_ij(4,3) = 0.75d0
232  a_dirk_ij(4,4) = 0.25d0
233 
234  omega(1) = 0.d0
235  omega(2) = 0.d0
236  omega(3) = 0.75d0
237  omega(4) = 0.25d0
238 
239  END IF
240 
241  ALLOCATE( a_tilde(n_rk) )
242  ALLOCATE( a_dirk(n_rk) )
243 
244  ALLOCATE( q_rk( n_vars , comp_cells , n_rk ) )
245  ALLOCATE( f_x( n_eqns , comp_cells , n_rk ) )
246  ALLOCATE( nh( n_eqns , comp_cells , n_rk ) )
247 
248  ALLOCATE( expl_terms( n_eqns , comp_cells , n_rk ) )
249 
250  ALLOCATE( fxj(n_eqns,n_rk) )
251  ALLOCATE( nhj(n_eqns,n_rk) )
252  ALLOCATE( expl_terms_j(n_eqns,n_rk) )
253 
254  ALLOCATE( residual_term( n_vars , comp_cells ) )
255 
256  END SUBROUTINE allocate_solver_variables
257 
258  !******************************************************************************
260  !
263  !
267  !
268  !******************************************************************************
269 
270  SUBROUTINE deallocate_solver_variables
271 
272  DEALLOCATE( q , q0 )
273 
274  DEALLOCATE( q_fv )
275 
276  DEALLOCATE( q_interfacel )
277  DEALLOCATE( q_interfacer )
278 
279  DEALLOCATE( a_interfacel )
280  DEALLOCATE( a_interfacer )
281 
282  DEALLOCATE( h_interface )
283 
284  Deallocate( qp )
285 
286  DEALLOCATE( a_tilde_ij )
287  DEALLOCATE( a_dirk_ij )
288  DEALLOCATE( omega_tilde )
289  DEALLOCATE( omega )
290 
291  DEALLOCATE( implicit_flag )
292 
293  DEALLOCATE( a_tilde )
294  DEALLOCATE( a_dirk )
295 
296  DEALLOCATE( q_rk )
297  DEALLOCATE( f_x )
298  DEALLOCATE( nh )
299 
300  DEALLOCATE( fxj )
301  DEALLOCATE( nhj )
302 
303  DEALLOCATE( mask22 , mask21 , mask11 , mask12 )
304 
305  DEALLOCATE( residual_term )
306 
307  END SUBROUTINE deallocate_solver_variables
308 
309  !*****************************************************************************
311  !
315  !
319  !
320  !*****************************************************************************
321 
322  SUBROUTINE timestep
323 
324  ! External variables
325  USE geometry, ONLY : dx
326  USE parameters, ONLY : max_dt , cfl
327 
328  ! External procedures
329  USE constitutive, ONLY : eval_local_speeds
330 
331  IMPLICIT none
332 
333  REAL*8 :: vel_max(n_vars)
334  REAL*8 :: vel_min(n_vars)
335  REAL*8 :: vel_j
336  REAL*8 :: dt_cfl
337  REAL*8 :: qj(n_vars)
338 
339  REAL*8 :: a_interface_max(n_eqns,comp_interfaces)
340  REAL*8 :: dt_interface
341 
342  INTEGER :: i, j
343 
344  dt = max_dt
345 
346  IF ( cfl .NE. -1.d0 ) THEN
347 
348  IF ( reconstr_variables .EQ. 0 ) THEN
349 
350  ! Linear reconstruction of the conservative variables at the interfaces
351  CALL reconstruction_cons
352 
353  ELSEIF ( reconstr_variables .EQ. 0 ) THEN
354 
355  ! Linear reconstruction of the physical var. (h+B,u) at the interfaces
356  CALL reconstruction_phys1
357 
358  ELSE
359 
360  ! Linear reconstruction of the physical var. (h,u) at the interfaces
361  CALL reconstruction_phys2
362 
363  END IF
364 
365  ! Evaluation of the maximum local speeds at the interfaces
366  CALL eval_speeds
367 
368  DO i=1,n_vars
369 
370  a_interface_max(i,:) = max( a_interfacer(i,:) , -a_interfacel(i,:) )
371 
372  END DO
373 
374  DO j = 1,comp_cells
375 
376  ! qj = q( 1:n_vars , j )
377 
378  ! CALL eval_local_speeds(qj,B_cent(j),vel_min,vel_max)
379 
380  ! vel_j = MAX( MAXVAL(ABS(vel_min)) , MAXVAL(ABS(vel_max)) )
381 
382  ! dt_cfl = cfl * dx / vel_j
383 
384  ! dt = MIN( dt , dt_cfl )
385 
386  dt_interface = cfl * dx / max( maxval(a_interface_max(:,j)) , &
387  maxval(a_interface_max(:,j+1)) )
388 
389  dt = min( dt , dt_interface )
390 
391  END DO
392 
393  END IF
394 
395 
396 
397 
398  END SUBROUTINE timestep
399 
400  !*****************************************************************************
402  !
407  !
411  !
412  !******************************************************************************
413 
414  SUBROUTINE imex_rk_solver
415 
416  USE constitutive, ONLY : eval_nonhyperbolic_terms
417 
418  USE constitutive, ONLY : qc_to_qp
419 
420  IMPLICIT NONE
421 
422  REAL*8 :: q_guess(n_vars)
423  INTEGER :: i_rk
424  INTEGER :: j
425 
426  REAL*8 :: nh_term(n_eqns)
427 
428  INTEGER :: i
429 
430  q0( 1:n_vars , 1:comp_cells ) = q( 1:n_vars , 1:comp_cells )
431 
432  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: beginning'
433 
434  ! Initialization of the variables for the Runge-Kutta scheme
435  q_rk(1:n_vars,1:comp_cells,1:n_rk) = 0.d0
436 
437  f_x(1:n_eqns,1:comp_cells,1:n_rk) = 0.d0
438 
439  nh(1:n_eqns,1:comp_cells,1:n_rk) = 0.d0
440 
441  expl_terms(1:n_eqns,1:comp_cells,1:n_rk) = 0.d0
442 
443  DO i_rk = 1,n_rk
444 
445  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: i_RK',i_rk
446 
447  ! define the explicits coefficients for the i-th step of the Runge-Kutta
448  a_tilde = 0.d0
449  a_dirk = 0.d0
450 
451  a_tilde(1:i_rk-1) = a_tilde_ij(i_rk,1:i_rk-1)
452  a_dirk(1:i_rk-1) = a_dirk_ij(i_rk,1:i_rk-1)
453 
454  ! define the implicit coefficient for the i-th step of the Runge-Kutta
455  a_diag = a_dirk_ij(i_rk,i_rk)
456 
457  DO j = 1,comp_cells
458 
459  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: j',j
460 
461  IF ( i_rk .EQ. 1 ) THEN
462 
463  q_guess(1:n_vars) = q0( 1:n_vars , j )
464 
465  ELSE
466 
467  q_guess(1:n_vars) = q_rk( 1:n_vars , j , max(1,i_rk-1) )
468 
469  END IF
470 
471  fxj(1:n_eqns,1:n_rk) = f_x( 1:n_eqns , j , 1:n_rk )
472 
473  nhj(1:n_eqns,1:n_rk) = nh( 1:n_eqns , j , 1:n_rk )
474 
475  expl_terms_j(1:n_eqns,1:n_rk) = expl_terms( 1:n_eqns , j , 1:n_rk )
476 
477  IF ( verbose_level .GE. 2 ) THEN
478 
479  WRITE(*,*) 'q_guess',q_guess
480  CALL qc_to_qp( q_guess , b_cent(j) , qp(1:n_vars,j) )
481  WRITE(*,*) 'q_guess: qp',qp(1:n_vars,j)
482 
483  END IF
484 
485  IF ( a_diag .NE. 0.d0 ) THEN
486 
487  CALL solve_rk_step( b_cent(j) , b_prime(j) , q_guess , &
488  q0(1:n_vars,j ) , a_tilde , a_dirk , a_diag )
489 
490  END IF
491 
492  q_rk( 1:n_vars , j , i_rk ) = q_guess
493 
494 
495  ! store the non-hyperbolic term for the explicit computations
496  IF ( a_diag .EQ. 0.d0 ) THEN
497 
498  CALL eval_nonhyperbolic_terms( b_cent(j) , b_prime(j) , &
499  r_qj = q_guess , r_nh_term_impl = nh(1:n_eqns,j,i_rk) )
500 
501  ELSE
502 
503  nh( 1:n_eqns , j , i_rk ) = 1.d0 / a_diag * ( ( q_guess - &
504  q0( 1:n_vars , j ) ) / dt + &
505  ( matmul(fxj,a_tilde) - matmul(nhj,a_dirk) ) )
506 
507  END IF
508 
509  IF ( verbose_level .GE. 2 ) THEN
510 
511  WRITE(*,*) 'imex_RK_solver: qc',q_guess
512  CALL qc_to_qp( q_guess, b_cent(j) , qp(1:n_vars,j) )
513  WRITE(*,*) 'imex_RK_solver: qp',qp(1:n_vars,j)
514  READ(*,*)
515 
516  END IF
517 
518  END DO
519 
520  ! Eval and save the explicit hyperbolic (fluxes) terms
521  CALL eval_hyperbolic_terms( q_rk(1:n_vars,1:comp_cells,i_rk) , &
522  f_x(1:n_eqns,1:comp_cells,i_rk) )
523 
524  ! Eval and save the other explicit terms (e.g. gravity or viscous forces)
525  CALL eval_explicit_terms( q_rk(1:n_vars,:,i_rk) , &
526  expl_terms(1:n_eqns,:,i_rk) )
527 
528  IF ( verbose_level .GE. 1 ) THEN
529 
530  WRITE(*,*) 'Fx, expl_terms'
531 
532  DO j = 1,comp_cells
533 
534  WRITE(*,*) f_x(1:n_vars,j,i_rk),expl_terms(2,j,i_rk)
535 
536  END DO
537 
538  READ(*,*)
539 
540  END IF
541 
542  END DO
543 
544  DO j = 1,comp_cells
545 
546  residual_term(1:n_vars,j) = matmul( f_x(1:n_eqns,j,1:n_rk) &
547  + expl_terms(1:n_eqns,j,1:n_rk) , omega_tilde ) &
548  - matmul( nh(1:n_eqns,j,1:n_rk) , omega )
549 
550  END DO
551 
552 
553 
554  DO j = 1,comp_cells
555 
556  IF ( verbose_level .GE. 1 ) THEN
557 
558  WRITE(*,*) 'cell j =',j
559  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j)
560  CALL qc_to_qp(q0(1:n_vars,j) , b_cent(j) , qp(1:n_vars,j))
561  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j)
562 
563  END IF
564 
565  q(1:n_vars,j) = q0(1:n_vars,j) - dt * residual_term(1:n_vars,j)
566 
567  IF ( verbose_level .GE. 1 ) THEN
568 
569  CALL qc_to_qp(q(1:n_vars,j) , b_cent(j) , qp(1:n_vars,j))
570 
571  WRITE(*,*) 'after imex_RK_solver: qc',q(1:n_vars,j)
572  WRITE(*,*) 'after imex_RK_solver: qp',qp(1:n_vars,j)
573  READ(*,*)
574 
575  END IF
576 
577  END DO
578 
579  END SUBROUTINE imex_rk_solver
580 
581  !******************************************************************************
583  !
590  !
598  !
602  !
603  !******************************************************************************
604 
605  SUBROUTINE solve_rk_step( Bj , Bprimej, qj, qj_old, a_tilde , a_dirk , a_diag )
606 
607  USE parameters, ONLY : max_nl_iter , tol_rel , tol_abs
608 
609  USE constitutive, ONLY : qc_to_qp
610 
611  IMPLICIT NONE
612 
613  REAL*8, INTENT(IN) :: bj
614  REAL*8, INTENT(IN) :: bprimej
615  REAL*8, INTENT(INOUT) :: qj(n_vars)
616  REAL*8, INTENT(IN) :: qj_old(n_vars)
617  REAL*8, INTENT(IN) :: a_tilde(n_rk)
618  REAL*8, INTENT(IN) :: a_dirk(n_rk)
619  REAL*8, INTENT(IN) :: a_diag
620 
621  REAL*8 :: qj_org(n_vars) , qj_rel(n_vars)
622 
623  REAL*8 :: left_matrix(n_eqns,n_vars)
624  REAL*8 :: right_term(n_eqns)
625 
626  REAL*8 :: scal_f
627 
628  REAL*8 :: coeff_f(n_eqns)
629 
630  REAL*8 :: qj_rel_nr_old(n_vars)
631  REAL*8 :: scal_f_old
632  REAL*8 :: desc_dir(n_vars)
633  REAL*8 :: grad_f(n_vars)
634 
635  INTEGER :: pivot(n_vars)
636 
637  REAL*8 :: desc_dir_small(n_vars)
638 
639  REAL*8 :: left_matrix_small22(n_nh,n_nh)
640  REAL*8 :: left_matrix_small21(n_eqns-n_nh,n_nh)
641  REAL*8 :: left_matrix_small11(n_eqns-n_nh,n_vars-n_nh)
642  REAL*8 :: left_matrix_small12(n_nh,n_vars-n_nh)
643 
644  REAL*8 :: right_term_small2(n_nh)
645  REAL*8 :: desc_dir_small2(n_nh)
646  INTEGER :: pivot_small2(n_nh)
647 
648  REAL*8 :: right_term_small1(n_eqns-n_nh)
649  REAL*8 :: desc_dir_small1(n_vars-n_nh)
650  INTEGER :: pivot_small1(n_vars-n_nh)
651 
652  INTEGER :: ok
653 
654  INTEGER :: i
655  INTEGER :: nl_iter
656 
657  REAL*8, PARAMETER :: stpmx=100.d0
658  REAL*8 :: stpmax
659  LOGICAL :: check
660 
661  REAL*8, PARAMETER :: tolf=1.d-10 , tolmin=1.d-6
662  REAL*8 :: tolx
663 
664  REAL*8 :: scal_f_init
665 
666  REAL*8 :: qpj(n_vars)
667 
668  REAL*8 :: desc_dir2(n_vars)
669 
670  REAL*8 :: desc_dir_temp(n_vars)
671 
672  normalize_q = .true.
673  normalize_f = .false.
674  opt_search_nl = .true.
675 
676  coeff_f(1:n_eqns) = 1.d0
677 
678  ! normalize the functions of the nonlinear system
679 
680  IF ( normalize_f ) THEN
681 
682  qj = qj_old - dt * ( matmul(fxj,a_tilde) - matmul(nhj,a_dirk) )
683 
684  CALL eval_f( bj , bprimej , qj , qj_old , a_tilde , a_dirk , a_diag , &
685  coeff_f , right_term , scal_f )
686 
687  IF ( verbose_level .GE. 3 ) THEN
688 
689  WRITE(*,*) 'solve_rk_step: non-normalized right_term'
690  WRITE(*,*) right_term
691  WRITE(*,*) 'scal_f',scal_f
692 
693  END IF
694 
695  DO i=1,n_eqns
696 
697  IF ( abs(right_term(i)) .GE. 1.d0 ) coeff_f(i) = 1.d0 / right_term(i)
698 
699  END DO
700 
701  right_term = coeff_f * right_term
702 
703  scal_f = 0.5d0 * dot_product( right_term , right_term )
704 
705  IF ( verbose_level .GE. 3 ) THEN
706  WRITE(*,*) 'solve_rk_step: after normalization',scal_f
707  END IF
708 
709  END IF
710 
711  ! set the initial guess for the NR iterative solver
712 
713 !!$ qj = qj_old - dt * ( MATMUL(Fxj,a_tilde) - MATMUL(NHj,a_dirk) )
714 !!$
715 !!$ DO i=1,n_eqns
716 !!$
717 !!$ IF ( implicit_flag(i) ) qj(i) = qj_old(i)
718 !!$
719 !!$ END DO
720 
721 ! qj = qj_old
722 
723  !---- normalize the conservative variables ------
724 
725  IF ( normalize_q ) THEN
726 
727  qj_org = qj
728 
729  qj_org = max( abs(qj_org) , 1.d-3 )
730 
731  ELSE
732 
733  qj_org(1:n_vars) = 1.d0
734 
735  END IF
736 
737  qj_rel = qj / qj_org
738 
739  ! -----------------------------------------------
740 
741  DO nl_iter=1,max_nl_iter
742 
743  tolx = epsilon(qj_rel)
744 
745  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solve_rk_step: nl_iter',nl_iter
746 
747  CALL eval_f( bj , bprimej , qj , qj_old , a_tilde , a_dirk , a_diag , &
748  coeff_f , right_term , scal_f )
749 
750  IF ( verbose_level .GE. 2 ) THEN
751 
752  WRITE(*,*) 'solve_rk_step: right_term',right_term
753 
754  END IF
755 
756  IF ( verbose_level .GE. 2 ) THEN
757 
758  WRITE(*,*) 'before_lnsrch: scal_f',scal_f
759 
760  END IF
761 
762  ! check the residual of the system
763 
764  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
765 
766  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
767  RETURN
768 
769  END IF
770 
771  IF ( ( normalize_f ) .AND. ( scal_f < 1.d-6 ) ) THEN
772 
773  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check scal_f',check
774  RETURN
775 
776  END IF
777 
778  ! ---- evaluate the descent direction ------------------------------------
779 
780  IF ( count( implicit_flag ) .EQ. n_eqns ) THEN
781 
782  CALL eval_jacobian(bj,bprimej,qj_rel,qj_org,coeff_f,left_matrix)
783 
784  desc_dir_temp = - right_term
785 
786  CALL dgesv(n_eqns,1, left_matrix , n_eqns, pivot, desc_dir_temp , &
787  n_eqns, ok)
788 
789  desc_dir = desc_dir_temp
790 
791  ELSE
792 
793  CALL eval_jacobian(bj,bprimej,qj_rel,qj_org,coeff_f,left_matrix)
794 
795  left_matrix_small11 = reshape(pack(left_matrix, mask11), &
796  [n_eqns-n_nh,n_eqns-n_nh])
797 
798  left_matrix_small12 = reshape(pack(left_matrix, mask12), &
799  [n_nh,n_eqns-n_nh])
800 
801  left_matrix_small22 = reshape(pack(left_matrix, mask22), &
802  [n_nh,n_nh])
803 
804  left_matrix_small21 = reshape(pack(left_matrix, mask21), &
805  [n_eqns-n_nh,n_nh])
806 
807 
808  desc_dir_small1 = pack( right_term, .NOT.implicit_flag )
809  desc_dir_small2 = pack( right_term , implicit_flag )
810 
811  DO i=1,n_vars-n_nh
812 
813  desc_dir_small1(i) = desc_dir_small1(i) / left_matrix_small11(i,i)
814 
815  END DO
816 
817  desc_dir_small2 = desc_dir_small2 - &
818  matmul( desc_dir_small1 , left_matrix_small21 )
819 
820  CALL dgesv(n_nh,1, left_matrix_small22 , n_nh , pivot_small2 , &
821  desc_dir_small2 , n_nh, ok)
822 
823  desc_dir = unpack( - desc_dir_small2 , implicit_flag , 0.0d0 ) &
824  + unpack( - desc_dir_small1 , .NOT.implicit_flag , 0.0d0 )
825 
826  END IF
827 
828  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'desc_dir',desc_dir
829 
830  qj_rel_nr_old = qj_rel
831  scal_f_old = scal_f
832 
833  IF ( ( opt_search_nl ) .AND. ( nl_iter .GT. 1 ) ) THEN
834  ! Search for the step lambda giving a sufficient decrease in the solution
835 
836  stpmax = stpmx * max( sqrt( dot_product(qj_rel,qj_rel) ) , &
837  dble( SIZE(qj_rel) ) )
838 
839  grad_f = matmul( right_term , left_matrix )
840 
841  desc_dir2 = desc_dir
842 
843  CALL lnsrch( bj , bprimej , qj_rel_nr_old , qj_org , qj_old , &
844  scal_f_old , grad_f , desc_dir , coeff_f , qj_rel , scal_f , &
845  right_term , stpmax , check )
846 
847  ELSE
848 
849  qj_rel = qj_rel_nr_old + desc_dir
850 
851  qj = qj_rel * qj_org
852 
853  CALL eval_f( bj , bprimej , qj , qj_old , a_tilde , a_dirk , a_diag , &
854  coeff_f , right_term , scal_f )
855 
856  END IF
857 
858  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'after_lnsrch: scal_f',scal_f
859 
860  qj = qj_rel * qj_org
861 
862  IF ( verbose_level .GE. 3 ) THEN
863 
864  WRITE(*,*) 'qj',qj
865  CALL qc_to_qp( qj , bj , qpj)
866  WRITE(*,*) 'qp',qpj
867 
868  END IF
869 
870  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
871 
872  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
873  check= .false.
874  RETURN
875 
876  END IF
877 
878  IF (check) THEN
879 
880  check = ( maxval( abs(grad_f(:)) * max( abs( qj_rel(:) ),1.d0 ) / &
881  max( scal_f , 0.5d0 * SIZE(qj_rel) ) ) < tolmin )
882 
883  IF ( verbose_level .GE. 3 ) WRITE(*,*) '2: check',check
884 ! RETURN
885 
886  END IF
887 
888  IF ( maxval( abs( qj_rel(:) - qj_rel_nr_old(:) ) / max( abs( qj_rel(:)) ,&
889  1.d0 ) ) < tolx ) THEN
890 
891  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check',check
892  RETURN
893 
894  END IF
895 
896  END DO
897 
898  END SUBROUTINE solve_rk_step
899 
900  !******************************************************************************
902  !
922  !******************************************************************************
923 
924  SUBROUTINE lnsrch( Bj , Bprimej , qj_rel_NR_old , qj_org , qj_old , &
925  scal_f_old , grad_f , desc_dir , coeff_f , qj_rel , scal_f , right_term ,&
926  stpmax , check )
927 
928  IMPLICIT NONE
929 
930  REAL*8, INTENT(IN) :: bj
931 
932  REAL*8, INTENT(IN) :: bprimej
933 
935  REAL*8, DIMENSION(:), INTENT(IN) :: qj_rel_nr_old
936 
938  REAL*8, DIMENSION(:), INTENT(IN) :: qj_org
939 
941  REAL*8, DIMENSION(:), INTENT(IN) :: qj_old
942 
944  REAL*8, DIMENSION(:), INTENT(IN) :: grad_f
945 
947  REAL*8, INTENT(IN) :: scal_f_old
948 
950  REAL*8, DIMENSION(:), INTENT(INOUT) :: desc_dir
951 
952  REAL*8, INTENT(IN) :: stpmax
953 
955  REAL*8, DIMENSION(:), INTENT(IN) :: coeff_f
956 
958  REAL*8, DIMENSION(:), INTENT(OUT) :: qj_rel
959 
961  REAL*8, INTENT(OUT) :: scal_f
962 
964  REAL*8, INTENT(OUT) :: right_term(n_eqns)
965 
967  LOGICAL, INTENT(OUT) :: check
968 
969  REAL*8, PARAMETER :: tolx=epsilon(qj_rel)
970 
971  INTEGER, DIMENSION(1) :: ndum
972  REAL*8 :: alf , a,alam,alam2,alamin,b,disc
973  REAL*8 :: scal_f2
974  REAL*8 :: desc_dir_abs
975  REAL*8 :: rhs1 , rhs2 , slope, tmplam
976 
977  REAL*8 :: scal_f_min , alam_min
978 
979  REAL*8 :: qj(n_vars)
980 
981  alf = 1.0d-4
982 
983  IF ( size(grad_f) == size(desc_dir) .AND. size(grad_f) == size(qj_rel) .AND. &
984  size(qj_rel) == size(qj_rel_nr_old) ) THEN
985 
986  ndum = size(grad_f)
987 
988  ELSE
989 
990  WRITE(*,*) 'nrerror: an assert_eq failed with this tag:', 'lnsrch'
991  stop 'program terminated by assert_eq4'
992 
993  END IF
994 
995  check = .false.
996 
997  desc_dir_abs = sqrt( dot_product(desc_dir,desc_dir) )
998 
999  IF ( desc_dir_abs > stpmax ) desc_dir(:) = desc_dir(:) * stpmax / desc_dir_abs
1000 
1001  slope = dot_product(grad_f,desc_dir)
1002 
1003  alamin = tolx / maxval( abs( desc_dir(:)) / max( abs(qj_rel_nr_old(:)),1.d0 ) )
1004 
1005  IF ( alamin .EQ. 0.d0) THEN
1006 
1007  qj_rel(:) = qj_rel_nr_old(:)
1008 
1009  RETURN
1010 
1011  END IF
1012 
1013  alam = 1.0d0
1014 
1015  scal_f_min = scal_f_old
1016 
1017  optimal_step_search: DO
1018 
1019  IF ( verbose_level .GE. 4 ) THEN
1020 
1021  WRITE(*,*) 'alam',alam
1022 
1023  END IF
1024 
1025  qj_rel = qj_rel_nr_old + alam * desc_dir
1026 
1027  qj = qj_rel * qj_org
1028 
1029  CALL eval_f( bj , bprimej , qj , qj_old , a_tilde , a_dirk , a_diag , &
1030  coeff_f , right_term , scal_f )
1031 
1032  IF ( verbose_level .GE. 4 ) THEN
1033 
1034  WRITE(*,*) 'lnsrch: effe_old,effe',scal_f_old,scal_f
1035  READ(*,*)
1036 
1037  END IF
1038 
1039  IF ( scal_f .LT. scal_f_min ) THEN
1040 
1041  scal_f_min = scal_f
1042  alam_min = alam
1043 
1044  END IF
1045 
1046  IF ( scal_f .LE. 0.9 * scal_f_old ) THEN
1047  ! sufficient function decrease
1048 
1049  IF ( verbose_level .GE. 4 ) THEN
1050 
1051  WRITE(*,*) 'sufficient function decrease'
1052 
1053  END IF
1054 
1055  EXIT optimal_step_search
1056 
1057  ELSE IF ( alam < alamin ) THEN
1058  ! convergence on Delta_x
1059 
1060  IF ( verbose_level .GE. 4 ) THEN
1061 
1062  WRITE(*,*) ' convergence on Delta_x',alam,alamin
1063 
1064  END IF
1065 
1066  qj_rel(:) = qj_rel_nr_old(:)
1067  scal_f = scal_f_old
1068  check = .true.
1069 
1070  EXIT optimal_step_search
1071 
1072 ! ELSE IF ( scal_f .LE. scal_f_old + ALF * alam * slope ) THEN
1073  ELSE
1074 
1075  IF ( alam .EQ. 1.d0 ) THEN
1076 
1077  tmplam = - slope / ( 2.0d0 * ( scal_f - scal_f_old - slope ) )
1078 
1079  ELSE
1080 
1081  rhs1 = scal_f - scal_f_old - alam*slope
1082  rhs2 = scal_f2 - scal_f_old - alam2*slope
1083 
1084  a = ( rhs1/alam**2.d0 - rhs2/alam2**2.d0 ) / ( alam - alam2 )
1085  b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )
1086 
1087  IF ( a .EQ. 0.d0 ) THEN
1088 
1089  tmplam = - slope / ( 2.0d0 * b )
1090 
1091  ELSE
1092 
1093  disc = b*b - 3.0d0*a*slope
1094 
1095  IF ( disc .LT. 0.d0 ) THEN
1096 
1097  tmplam = 0.5d0 * alam
1098 
1099  ELSE IF ( b .LE. 0.d0 ) THEN
1100 
1101  tmplam = ( - b + sqrt(disc) ) / ( 3.d0 * a )
1102 
1103  ELSE
1104 
1105  tmplam = - slope / ( b + sqrt(disc) )
1106 
1107  ENDIF
1108 
1109  END IF
1110 
1111  IF ( tmplam .GT. 0.5d0*alam ) tmplam = 0.5d0 * alam
1112 
1113  END IF
1114 
1115  END IF
1116 
1117  alam2 = alam
1118  scal_f2 = scal_f
1119  alam = max( tmplam , 0.5d0*alam )
1120 
1121  END DO optimal_step_search
1122 
1123  END SUBROUTINE lnsrch
1124 
1125  !******************************************************************************
1127  !
1143  !******************************************************************************
1144 
1145  SUBROUTINE eval_f( Bj , Bprimej , qj , qj_old , a_tilde , a_dirk , a_diag , &
1146  coeff_f , f_nl , scal_f )
1147 
1148  USE constitutive, ONLY : eval_nonhyperbolic_terms
1149 
1150  IMPLICIT NONE
1151 
1152  REAL*8, INTENT(IN) :: bj
1153  REAL*8, INTENT(IN) :: bprimej
1154  REAL*8, INTENT(IN) :: qj(n_vars)
1155  REAL*8, INTENT(IN) :: qj_old(n_vars)
1156  REAL*8, INTENT(IN) :: a_tilde(n_rk)
1157  REAL*8, INTENT(IN) :: a_dirk(n_rk)
1158  REAL*8, INTENT(IN) :: a_diag
1159  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1160  REAL*8, INTENT(OUT) :: f_nl(n_eqns)
1161  REAL*8, INTENT(OUT) :: scal_f
1162 
1163  REAL*8 :: nh_term_impl(n_eqns)
1164  REAL*8 :: rj(n_eqns)
1165 
1166  CALL eval_nonhyperbolic_terms( bj , bprimej , r_qj = qj , &
1167  r_nh_term_impl = nh_term_impl )
1168 
1169  rj = ( matmul(fxj,a_tilde) - matmul(nhj,a_dirk) ) - a_diag * nh_term_impl
1170 
1171  f_nl = qj - qj_old + dt * rj
1172 
1173  f_nl = coeff_f * f_nl
1174 
1175  scal_f = 0.5d0 * dot_product( f_nl , f_nl )
1176 
1177  END SUBROUTINE eval_f
1178 
1179  !******************************************************************************
1181  !
1184  !
1191  !
1195  !******************************************************************************
1196 
1197  SUBROUTINE eval_jacobian(Bj , Bprimej , qj_rel , qj_org , coeff_f, left_matrix)
1198 
1199  USE constitutive, ONLY : eval_nonhyperbolic_terms
1200 
1201  IMPLICIT NONE
1202 
1203  REAL*8, INTENT(IN) :: bj
1204  REAL*8, INTENT(IN) :: bprimej
1205  REAL*8, INTENT(IN) :: qj_rel(n_vars)
1206  REAL*8, INTENT(IN) :: qj_org(n_vars)
1207  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1208  REAL*8, INTENT(OUT) :: left_matrix(n_eqns,n_vars)
1209 
1210  REAL*8 :: jacob_relax(n_eqns,n_vars)
1211  COMPLEX*16 :: nh_terms_cmplx_impl(n_eqns)
1212  COMPLEX*16 :: qj_cmplx(n_vars) , qj_rel_cmplx(n_vars)
1213 
1214  REAL*8 :: h
1215 
1216  INTEGER :: i
1217 
1218  h = n_vars * epsilon(1.d0)
1219 
1220  ! initialize the matrix of the linearized system and the Jacobian
1221 
1222  left_matrix(1:n_eqns,1:n_vars) = 0.d0
1223  jacob_relax(1:n_eqns,1:n_vars) = 0.d0
1224 
1225  ! evaluate the jacobian of the non-hyperbolic terms
1226 
1227  DO i=1,n_vars
1228 
1229  left_matrix(i,i) = coeff_f(i) * qj_org(i)
1230 
1231  IF ( implicit_flag(i) ) THEN
1232 
1233  qj_rel_cmplx(1:n_vars) = qj_rel(1:n_vars)
1234  qj_rel_cmplx(i) = dcmplx(qj_rel(i), h)
1235 
1236  qj_cmplx = qj_rel_cmplx * qj_org
1237 
1238  CALL eval_nonhyperbolic_terms( bj , bprimej , c_qj = qj_cmplx , &
1239  c_nh_term_impl = nh_terms_cmplx_impl )
1240 
1241  jacob_relax(1:n_eqns,i) = coeff_f(i) * &
1242  aimag(nh_terms_cmplx_impl) / h
1243 
1244  left_matrix(1:n_eqns,i) = left_matrix(1:n_eqns,i) - dt * a_diag &
1245  * jacob_relax(1:n_eqns,i)
1246 
1247  END IF
1248 
1249  END DO
1250 
1251  END SUBROUTINE eval_jacobian
1252 
1253  !******************************************************************************
1255  !
1263  !******************************************************************************
1264 
1265  SUBROUTINE eval_explicit_terms( q_expl , expl_terms )
1266 
1267  USE constitutive, ONLY : eval_explicit_forces
1268 
1269  IMPLICIT NONE
1270 
1271  REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells)
1272  REAL*8, INTENT(OUT) :: expl_terms(n_eqns,comp_cells)
1273 
1274  REAL*8 :: qc(n_vars)
1275  REAL*8 :: expl_forces_term(n_eqns)
1276 
1277  INTEGER :: j
1278 
1279  DO j = 1,comp_cells
1280 
1281  qc = q_expl(1:n_vars,j)
1282 
1283  CALL eval_explicit_forces( b_cent(j) , b_prime(j) , qc , expl_forces_term )
1284 
1285  expl_terms(1:n_eqns,j) = expl_forces_term
1286 
1287  END DO
1288 
1289  END SUBROUTINE eval_explicit_terms
1290 
1291  !******************************************************************************
1293  !
1301  !******************************************************************************
1302 
1303  SUBROUTINE eval_hyperbolic_terms( q_expl , F_x )
1304 
1305  ! External variables
1306  USE geometry, ONLY : dx
1307  USE parameters, ONLY : solver_scheme
1308 
1309  IMPLICIT NONE
1310 
1311  REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells)
1312  REAL*8, INTENT(OUT) :: f_x(n_eqns,comp_cells)
1313 
1314  REAL*8 :: q_old(n_vars,comp_cells)
1315 
1316  INTEGER :: i, j
1317 
1318  REAL*8 :: h_new
1319 
1320  q_old = q
1321 
1322  q = q_expl
1323 
1324  IF ( reconstr_variables .EQ. 0 ) THEN
1325 
1326  ! Linear reconstruction of the conservative variables at the interfaces
1327  CALL reconstruction_cons
1328 
1329  ELSEIF ( reconstr_variables .EQ. 0 ) THEN
1330 
1331  ! Linear reconstruction of the physical var. (h+B,u) at the interfaces
1332  CALL reconstruction_phys1
1333 
1334  ELSE
1335 
1336  ! Linear reconstruction of the physical var. (h,u) at the interfaces
1337  CALL reconstruction_phys2
1338 
1339  END IF
1340 
1341  ! Evaluation of the maximum local speeds at the interfaces
1342  CALL eval_speeds
1343 
1344  ! Evaluation of the numerical fluxes
1345  SELECT CASE ( solver_scheme )
1346 
1347  CASE ("LxF")
1348 
1349  CALL eval_flux_lxf
1350 
1351  CASE ("GFORCE")
1352 
1353  CALL eval_flux_gforce
1354 
1355  CASE ("KT")
1356 
1357  CALL eval_flux_kt
1358 
1359  END SELECT
1360 
1361 
1362  ! Advance in time the solution
1363  DO j = 1,comp_cells
1364 
1365  DO i=1,n_eqns
1366 
1367  f_x(i,j) = ( h_interface(i,j+1) - h_interface(i,j) ) / dx
1368 
1369  END DO
1370 
1371  h_new = q_expl(1,j) - dt * f_x(1,j) - b_cent(j)
1372 
1373  IF ( h_new .LT. 0.d0 ) THEN
1374 
1375  WRITE(*,*) 'j,h',j,h_new
1376  WRITE(*,*) 'dt',dt
1377  DO i=-2,2
1378 
1379  WRITE(*,*) b_cent(j+i) , b_prime(j+i) , q_expl(1,j+i) , q_expl(2,j+i)
1380 
1381  END DO
1382 
1383  WRITE(*,*) 'w_interface(j) ',q_interfacel(1,j) , q_interfacer(1,j)
1384  WRITE(*,*) 'w_interface(j+1) ',q_interfacel(1,j+1) , q_interfacer(1,j+1)
1385 
1386  WRITE(*,*) 'uh_interface(j) ',q_interfacel(2,j) , q_interfacer(2,j)
1387  WRITE(*,*) 'uh_interface(j+1) ',q_interfacel(2,j+1) , q_interfacer(2,j+1)
1388 
1389  WRITE(*,*) 'a(j) ',a_interfacel(1,j) , a_interfacer(1,j)
1390  WRITE(*,*) 'a(j+1) ',a_interfacel(1,j+1) , a_interfacer(1,j+1)
1391 
1392  WRITE(*,*) 'H_interface ',h_interface(i,j) , h_interface(i,j+1)
1393  WRITE(*,*)
1394  READ(*,*)
1395 
1396  END IF
1397 
1398  IF ( verbose_level .GE. 1 ) THEN
1399 
1400  WRITE(*,*) 'F_x',j,f_x(:,j), h_interface(1,j+1) , h_interface(1,j)
1401  READ(*,*)
1402 
1403  END IF
1404 
1405  END DO
1406 
1407  !WRITE(*,*) 'flux left', H_interface(1,1)
1408  !WRITE(*,*) 'flux right', H_interface(1,comp_cells+1)
1409  !WRITE(*,*) 'sum fluxes',SUM(F_x(1,:))
1410  !READ(*,*)
1411 
1412 
1413  q = q_old
1414 
1415  END SUBROUTINE eval_hyperbolic_terms
1416 
1417 
1418  !******************************************************************************
1420  !
1426  !******************************************************************************
1427 
1428  SUBROUTINE eval_flux_kt
1429 
1430  ! External procedures
1431  USE constitutive, ONLY : eval_fluxes
1432 
1433  ! External variables
1434  USE geometry, ONLY : dx
1435 
1436  IMPLICIT NONE
1437 
1438  REAL*8 :: fluxl(n_eqns)
1439  REAL*8 :: fluxr(n_eqns)
1440 
1441  REAL*8 :: flux_avg(n_eqns)
1442 
1443  REAL*8 :: nc_fluxl(n_eqns)
1444  REAL*8 :: nc_fluxr(n_eqns)
1445 
1446  REAL*8 :: nc_flux_avg(n_eqns)
1447 
1448  REAL*8 :: nc_terml(n_eqns)
1449  REAL*8 :: nc_termr(n_eqns)
1450 
1451  REAL*8 :: nc_term_avg(n_eqns)
1452 
1453  INTEGER :: j
1454  INTEGER :: i
1455 
1456  REAL*8 :: dt_loc , vel_loc
1457 
1458  DO j = 1 , comp_interfaces
1459 
1460  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacer(1:n_vars,j) , &
1461  r_flux = fluxr )
1462 
1463  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacel(1:n_vars,j) , &
1464  r_flux = fluxl )
1465 
1466  CALL average_kt( a_interfacel(:,j) , a_interfacer(:,j) , fluxl , fluxr , &
1467  flux_avg )
1468 
1469  DO i=1,n_eqns
1470 
1471  IF ( a_interfacel(i,j) .EQ. a_interfacer(i,j) ) THEN
1472 
1473  h_interface(i,j) = 0.d0
1474 
1475  ELSE
1476 
1477  h_interface(i,j) = flux_avg(i) &
1478  + ( a_interfacer(i,j) * a_interfacel(i,j) ) &
1479  / ( a_interfacer(i,j) - a_interfacel(i,j) ) &
1480  * ( q_interfacer(i,j) - q_interfacel(i,j) )
1481 
1482  END IF
1483 
1484  END DO
1485 
1486 ! IF ( j .EQ. comp_interfaces ) THEN
1487 !
1488 ! WRITE(*,*) 'flux', H_interface(1,j) , flux_avg(1), &
1489 ! ( a_interfaceR(1,j) * a_interfaceL(1,j) ) &
1490 ! / ( a_interfaceR(1,j) - a_interfaceL(1,j) ) &
1491 ! * ( q_interfaceR(1,j) - q_interfaceL(1,j) )
1492 !
1493 ! END IF
1494 
1495  END DO
1496 
1497  END SUBROUTINE eval_flux_kt
1498 
1499  SUBROUTINE average_kt( aL , aR , wL , wR , w_avg )
1500 
1501  IMPLICIT NONE
1502 
1503  REAL*8, INTENT(IN) :: al(:) , ar(:)
1504  REAL*8, INTENT(IN) :: wl(:) , wr(:)
1505  REAL*8, INTENT(OUT) :: w_avg(:)
1506 
1507  INTEGER :: n
1508  INTEGER :: i
1509 
1510  n = SIZE( al )
1511 
1512  DO i=1,n
1513 
1514  IF ( al(i) .EQ. ar(i) ) THEN
1515 
1516  ! w_avg(i) = 0.5D0 * ( wL(i) + wR(i) )
1517  w_avg(i) = 0.0d0
1518 
1519  ELSE
1520 
1521  w_avg(i) = ( ar(i) * wl(i) - al(i) * wr(i) ) / ( ar(i) - al(i) )
1522 
1523  END IF
1524 
1525  END DO
1526 
1527  END SUBROUTINE average_kt
1528 
1529  !******************************************************************************
1534  !******************************************************************************
1535 
1536  SUBROUTINE eval_flux_gforce
1537 
1538  ! External procedures
1539  USE constitutive, ONLY : eval_fluxes
1540 
1541  ! External variables
1542  USE geometry, ONLY : dx
1543 
1544  IMPLICIT NONE
1545 
1546  REAL*8 :: fluxl(n_eqns)
1547  REAL*8 :: fluxr(n_eqns)
1548  REAL*8 :: flux_lf(n_eqns)
1549  REAL*8 :: flux_lw(n_eqns)
1550  REAL*8 :: q_lw(n_vars)
1551 
1552  INTEGER :: j
1553  INTEGER :: i
1554 
1555  REAL*8 :: cfl_loc
1556  REAL*8 :: vel_loc
1557  REAL*8 :: dt_loc
1558  REAL*8 :: omega
1559 
1560  cfl_loc = 0.9
1561 
1562  omega = 1.d0 / ( 1.d0 + cfl_loc )
1563 
1564  DO j = 1,comp_interfaces
1565 
1566  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacel(1:n_vars,j) , &
1567  r_flux = fluxl )
1568 
1569  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacer(1:n_vars,j) , &
1570  r_flux = fluxr )
1571 
1572  vel_loc = max( abs( a_interfacel(1,j) ) , abs( a_interfacer(1,j) ) )
1573 
1574  dt_loc = cfl_loc * dx / vel_loc
1575 
1576  DO i=1,n_vars
1577 
1578  flux_lf(i) = 0.5d0 * ( fluxl(i) + fluxr(i) ) - 0.5d0 * dx / &
1579  dt_loc * ( q_interfacer(i,j) - q_interfacel(i,j) )
1580 
1581  q_lw(i) = 0.5d0 * ( q_interfacer(i,j) + q_interfacel(i,j) ) &
1582  - 0.5d0 * dt_loc / dx * ( fluxr(i) - fluxl(i) )
1583 
1584  END DO
1585 
1586  CALL eval_fluxes( b_stag(j) , r_qj = q_lw, r_flux = flux_lw )
1587 
1588  DO i=1,n_eqns
1589 
1590  h_interface(i,j) = ( 1.d0 - omega ) * flux_lf(i) + omega * flux_lw(i)
1591 
1592  END DO
1593 
1594  END DO
1595 
1596  END SUBROUTINE eval_flux_gforce
1597 
1598  !*****************************************************************************
1603  !*****************************************************************************
1604 
1605  SUBROUTINE eval_flux_lxf
1606 
1607  ! External procedures
1608  USE constitutive, ONLY : eval_fluxes
1609 
1610  ! External variables
1611  USE geometry, ONLY : dx
1612 
1613  IMPLICIT NONE
1614 
1615  REAL*8 :: fluxl(n_eqns)
1616  REAL*8 :: fluxr(n_eqns)
1617  REAL*8 :: flux_lf(n_eqns)
1618 
1619  INTEGER :: j
1620  INTEGER :: i
1621  REAL*8 :: vel_loc
1622  REAL*8 :: dt_loc !
1623 
1624  DO j = 1,comp_interfaces
1625 
1626  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacel(1:n_eqns,j) , &
1627  r_flux = fluxl )
1628 
1629  CALL eval_fluxes( b_stag(j) , r_qj = q_interfacer(1:n_eqns,j) , &
1630  r_flux = fluxr )
1631 
1632  vel_loc = max( abs( a_interfacel(1,j) ) , abs( a_interfacer(1,j) ) )
1633 
1634  dt_loc = 0.9d0 * dx / vel_loc
1635 
1636  DO i = 1,n_eqns
1637 
1638  flux_lf(i) = 0.5d0 * ( fluxr(i) + fluxl(i) ) - 0.5d0 * dx / &
1639  dt_loc *( q_interfacer(i,j) - q_interfacel(i,j) )
1640 
1641  END DO
1642 
1643  DO i=1,n_eqns
1644 
1645  h_interface(i,j) = flux_lf(i)
1646 
1647  END DO
1648 
1649  END DO
1650 
1651  END SUBROUTINE eval_flux_lxf
1652 
1653  !******************************************************************************
1655  !
1662  !******************************************************************************
1663 
1664  SUBROUTINE reconstruction_cons
1665 
1666  ! External procedures
1667  USE constitutive, ONLY : qc_to_qp , qp_to_qc
1668  USE parameters, ONLY : limiter
1669 
1670  ! External variables
1671  USE geometry, ONLY : x_comp , x_stag
1672  USE parameters, ONLY : reconstr_coeff
1673 
1674  IMPLICIT NONE
1675 
1676  REAL*8 :: qc(n_vars)
1677  REAL*8 :: qcl(n_vars)
1678  REAL*8 :: qcr(n_vars)
1679 
1680  REAL*8 :: qc_stencil(3)
1681  REAL*8 :: x_stencil(3)
1682  REAL*8 :: qc_prime
1683 
1684  REAL*8 :: qp_check(n_vars)
1685 
1686  REAL*8 :: dxl , dxr
1687 
1688  INTEGER :: j
1689  INTEGER :: i
1690 
1691  DO i=1,n_vars
1692 
1693  ! Slope in the first cell
1694  IF ( bcl(i)%flag .EQ. 0 ) THEN
1695  ! Dirichelet boundary condition
1696 
1697  qc_stencil(1) = bcl(i)%value
1698  qc_stencil(2:3) = q(i,1:2)
1699 
1700  x_stencil(1) = x_stag(1)
1701  x_stencil(2:3) = x_comp(1:2)
1702 
1703  CALL limit( qc_stencil , x_stencil , limiter(i) , qc_prime )
1704 
1705  ELSEIF ( bcl(i)%flag .EQ. 1 ) THEN
1706  ! Neumann boundary condition (fixed slope)
1707 
1708  qc_prime = bcl(i)%value
1709 
1710  ELSEIF ( bcl(i)%flag .EQ. 2 ) THEN
1711  ! Linear extrapolation from internal values
1712 
1713  qc_prime = ( q(i,2) - q(i,1) ) / ( x_comp(2) - x_comp(1) )
1714 
1715  END IF
1716 
1717  ! Linear reconstruction of the conservative variables at the boundaries
1718  ! of for the first cell
1719  dxl = x_comp(1) - x_stag(1)
1720  qcl(i) = q(i,1) - reconstr_coeff * dxl * qc_prime
1721 
1722  dxr = x_stag(2) - x_comp(1)
1723  qcr(i) = q(i,1) + reconstr_coeff * dxr * qc_prime
1724 
1725  ! Positivity preserving reconstruction for h
1726  IF (i.eq.1) THEN
1727 
1728  IF ( qcr(i) .LT. b_stag(2) ) THEN
1729 
1730  qc_prime = ( b_stag(2) - q(i,1) ) / dxr
1731 
1732  qcl(i) = q(i,1) - reconstr_coeff * dxl * qc_prime
1733 
1734  qcr(i) = q(i,1) + reconstr_coeff * dxr * qc_prime
1735 
1736  ENDIF
1737 
1738  IF ( qcl(i) .LT. b_stag(1) ) THEN
1739 
1740  qc_prime = ( q(i,1) - b_stag(1) ) / dxl
1741 
1742  qcl(i) = q(i,1) - reconstr_coeff * dxl * qc_prime
1743 
1744  qcr(i) = q(i,1) + reconstr_coeff * dxr * qc_prime
1745 
1746  ENDIF
1747 
1748  ENDIF
1749 
1750  END DO
1751 
1752  ! Correction for small values (Eqs. 2.17 and 2.21 K&P)
1753  CALL qc_to_qp( qcl , b_stag(1) , qp_check )
1754  CALL qp_to_qc( qp_check , b_stag(1) , qcl )
1755 
1756  CALL qc_to_qp( qcr , b_stag(2) , qp_check )
1757  CALL qp_to_qc( qp_check , b_stag(2) , qcr )
1758 
1759 
1760  ! Value at the right of the first interface
1761  q_interfacer(1:n_vars,1) = qcl
1762  ! Value at the left of the second interface
1763  q_interfacel(1:n_vars,2) = qcr
1764 
1765  ! Values at the left of the first interface (out of the domain)
1766  DO i=1,n_vars
1767 
1768  IF ( bcl(i)%flag .EQ. 0 ) THEN
1769 
1770  q_interfacer(i,1) = bcl(i)%value
1771  q_interfacel(i,1) = bcl(i)%value
1772 
1773  ELSEIF ( bcl(i)%flag .EQ. 1 ) THEN
1774 
1775  q_interfacel(i,1) = q_interfacer(i,1)
1776 
1777  ELSE
1778 
1779  q_interfacel(i,1) = qcl(i)
1780 
1781  END IF
1782 
1783  END DO
1784 
1785  ! Linear reconstruction of the physical variables for the internal cells
1786  DO j = 2,comp_cells-1
1787 
1788  x_stencil(1:3) = x_comp(j-1:j+1)
1789 
1790  DO i=1,n_vars
1791 
1792  qc_stencil = q(i,j-1:j+1)
1793 
1794  CALL limit( qc_stencil , x_stencil , limiter(i) , qc_prime )
1795 
1796  ! Linear reconstruction at the left bdry of the cell
1797  dxl = x_comp(j) - x_stag(j)
1798  qcl(i) = q(i,j) - reconstr_coeff * dxl * qc_prime
1799 
1800  ! Linear reconstruction at the right bdry of the cell
1801  dxr = x_stag(j+1) - x_comp(j)
1802  qcr(i) = q(i,j) + reconstr_coeff * dxr * qc_prime
1803 
1804  ! Positivity preserving reconstruction for h
1805  IF ( i .EQ. 1 ) THEN
1806 
1807  IF ( qcr(i) .LT. b_stag(j+1) ) THEN
1808 
1809  qc_prime = ( b_stag(j+1) - q(i,j) ) / dxr
1810 
1811  qcl(i) = q(i,j) - reconstr_coeff * dxl * qc_prime
1812 
1813  qcr(i) = q(i,j) + reconstr_coeff * dxr * qc_prime
1814 
1815  ENDIF
1816 
1817  IF ( qcl(i) .LT. b_stag(j) ) THEN
1818 
1819  qc_prime = ( q(i,j) - b_stag(j) ) / dxl
1820 
1821  qcl(i) = q(i,j) - reconstr_coeff * dxl * qc_prime
1822 
1823  qcr(i) = q(i,j) + reconstr_coeff * dxr * qc_prime
1824 
1825  ENDIF
1826 
1827  ENDIF
1828 
1829  END DO
1830 
1831 
1832  ! Correction for small values (Eqs. 2.17 and 2.21 K&P)
1833  CALL qc_to_qp( qcl , b_stag(j) , qp_check )
1834  CALL qp_to_qc( qp_check , b_stag(j) , qcl )
1835 
1836 
1837  CALL qc_to_qp( qcr , b_stag(j+1) , qp_check )
1838  CALL qp_to_qc( qp_check , b_stag(j+1) , qcr )
1839 
1840  ! Value at the right of the j-th interface
1841  q_interfacer(:,j) = qcl
1842  ! Value at the left of the j+1-th interface
1843  q_interfacel(:,j+1) = qcr
1844 
1845  END DO
1846 
1847  ! Linear reconstruction of the physical variables at the boundaries
1848  ! of the last cell (j=comp_cells)
1849 
1850  DO i=1,n_vars
1851 
1852  IF ( bcr(i)%flag .EQ. 0 ) THEN
1853 
1854  qc_stencil(1:2) = q(i,comp_cells-1:comp_cells)
1855  qc_stencil(3) = bcr(i)%value
1856 
1857  x_stencil(1:2) = x_comp(comp_cells-1:comp_cells)
1858  x_stencil(3) = x_stag(comp_interfaces)
1859 
1860  CALL limit( qc_stencil , x_stencil , limiter(i) , qc_prime )
1861 
1862  ELSEIF ( bcr(i)%flag .EQ. 1 ) THEN
1863 
1864  qc_prime = bcr(i)%value
1865 
1866  ELSEIF ( bcr(i)%flag .EQ. 2 ) THEN
1867 
1868  qc_prime = ( q(i,comp_cells) - q(i,comp_cells-1) ) / &
1869  ( x_comp(comp_cells) - x_comp(comp_cells-1) )
1870 
1871  END IF
1872 
1873  dxl = x_comp(comp_cells) - x_stag(comp_interfaces-1)
1874  qcl(i) = q(i,comp_cells) - reconstr_coeff * dxl * qc_prime
1875 
1876  dxr = x_stag(comp_interfaces) - x_comp(comp_cells)
1877  qcr(i) = q(i,comp_cells) + reconstr_coeff * dxr * qc_prime
1878 
1879  ! Positivity preserving reconstruction for h
1880  IF (i.eq.1) THEN
1881 
1882  IF ( qcr(i) .LT. b_stag(comp_interfaces) ) THEN
1883 
1884  qc_prime = ( b_stag(comp_interfaces) - q(i,comp_cells) ) / dxr
1885 
1886  qcl(i) = q(i,comp_cells) - reconstr_coeff * dxl * qc_prime
1887 
1888  qcr(i) = q(i,comp_cells) + reconstr_coeff * dxr * qc_prime
1889 
1890  ENDIF
1891 
1892  IF ( qcl(i) .LT. b_stag(comp_interfaces-1) ) THEN
1893 
1894  qc_prime = ( q(i,comp_cells) - b_stag(comp_cells) ) / dxl
1895 
1896  qcl(i) = q(i,comp_cells) - reconstr_coeff * dxl * qc_prime
1897 
1898  qcr(i) = q(i,comp_cells) + reconstr_coeff * dxr * qc_prime
1899 
1900  ENDIF
1901 
1902  ENDIF
1903 
1904  END DO
1905 
1906  ! Correction for small values (Eqs. 2.17 and 2.21 K&P)
1907  CALL qc_to_qp( qcl , b_stag(comp_cells) , qp_check )
1908  CALL qp_to_qc( qp_check , b_stag(comp_cells) , qcl )
1909 
1910  CALL qc_to_qp( qcr , b_stag(comp_cells+1) , qp_check )
1911  CALL qp_to_qc( qp_check , b_stag(comp_cells+1) , qcr )
1912 
1913  q_interfacer(:,comp_interfaces-1) = qcl
1914  q_interfacel(:,comp_interfaces) = qcr
1915 
1916  ! Values at the right of the last interface (out of the domain)
1917  DO i=1,n_vars
1918 
1919  IF ( bcr(i)%flag .EQ. 0 ) THEN
1920 
1921  q_interfacel(i,comp_interfaces) = bcr(i)%value
1922  q_interfacer(i,comp_interfaces) = bcr(i)%value
1923 
1924  ELSEIF ( bcr(i)%flag .EQ. 1 ) THEN
1925 
1926  q_interfacer(i,comp_interfaces) = q_interfacel(i,comp_interfaces)
1927 
1928  ELSE
1929 
1930  q_interfacer(i,comp_interfaces) = qcr(i)
1931 
1932  END IF
1933 
1934  END DO
1935 
1936  END SUBROUTINE reconstruction_cons
1937 
1938  !******************************************************************************
1940  !
1947  !******************************************************************************
1948 
1949 
1950  SUBROUTINE reconstruction_phys1
1951 
1952  ! External procedures
1953  USE constitutive, ONLY : qc_to_qp , qp_to_qc
1954  USE constitutive, ONLY : qc_to_qp2 , qp2_to_qc
1955  USE parameters, ONLY : limiter
1956 
1957  ! External variables
1958  USE geometry, ONLY : x_comp , x_stag
1959  USE parameters, ONLY : reconstr_coeff
1960 
1961  IMPLICIT NONE
1962 
1963  REAL*8 :: qc(n_vars)
1964  REAL*8 :: qpl(n_vars)
1965  REAL*8 :: qpr(n_vars)
1966  REAL*8 :: qp_bdry(n_vars)
1967 
1968  REAL*8 :: qp_stencil(3)
1969  REAL*8 :: x_stencil(3)
1970  REAL*8 :: qp_prime
1971 
1972  REAL*8 :: dxl , dxr
1973 
1974  INTEGER :: j
1975  INTEGER :: i
1976 
1977 
1978  ! Convert the conservative variables to the physical variables
1979  DO j = 1,comp_cells
1980 
1981  qc = q(1:n_vars,j)
1982 
1983  CALL qc_to_qp( qc , b_cent(j) , qp(1:n_vars,j) )
1984 
1985  END DO
1986 
1987  DO i=1,n_vars
1988 
1989  ! Slope in the first cell
1990  IF ( bcl(i)%flag .EQ. 0 ) THEN
1991 
1992  qp_stencil(1) = bcl(i)%value
1993  qp_stencil(2:3) = qp(i,1:2)
1994 
1995  x_stencil(1) = x_stag(1)
1996  x_stencil(2:3) = x_comp(1:2)
1997 
1998  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
1999 
2000  ELSEIF ( bcl(i)%flag .EQ. 1 ) THEN
2001 
2002  qp_prime = bcl(i)%value
2003 
2004  ELSEIF ( bcl(i)%flag .EQ. 2 ) THEN
2005 
2006  qp_prime = ( qp(i,2) - qp(i,1) ) / ( x_comp(2) - x_comp(1) )
2007 
2008  END IF
2009 
2010  ! Linear reconstruction of the physical variables at the boundaries
2011  ! of for the first cell
2012  dxl = x_comp(1) - x_stag(1)
2013  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2014 
2015  dxr = x_stag(2) - x_comp(1)
2016  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2017 
2018  ! positivity preserving reconstruction for h
2019  IF(i.eq.1)THEN
2020 
2021  IF(qpr(i).LT.b_stag(2))THEN
2022 
2023  qp_prime=(b_stag(2)-qp(i,1))/dxr
2024 
2025  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2026 
2027  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2028 
2029  ENDIF
2030 
2031  IF(qpl(i).LT.b_stag(1))THEN
2032 
2033  qp_prime=(qp(i,1)-b_stag(1))/dxl
2034 
2035  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2036 
2037  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2038 
2039  ENDIF
2040 
2041  ENDIF
2042 
2043  END DO
2044 
2045  DO i=1,n_vars
2046 
2047  IF ( bcl(i)%flag .EQ. 0 ) THEN
2048 
2049  qpl(i) = bcl(i)%value
2050 
2051  END IF
2052 
2053  END DO
2054 
2055 
2056  ! Convert back from physical to conservative variables
2057  CALL qp_to_qc( qpl , b_stag(1) , q_interfacer(1:n_vars,1) )
2058  CALL qp_to_qc( qpr , b_stag(2) , q_interfacel(1:n_vars,2) )
2059 
2060 
2061  DO i=1,n_vars
2062 
2063  IF ( bcl(i)%flag .EQ. 0 ) THEN
2064 
2065  qp_bdry(i) = bcl(i)%value
2066 
2067  ELSE
2068 
2069  !qp_bdry(i) = qpL(i)
2070 
2071  ! fixed wall
2072  IF(i.eq.2)THEN
2073 
2074  qp_bdry(i) = -qpl(i)
2075 
2076  ELSE
2077 
2078  qp_bdry(i) = qpl(i)
2079 
2080  ENDIF
2081 
2082  END IF
2083 
2084  END DO
2085 
2086  CALL qp_to_qc( qp_bdry , b_stag(1) , q_interfacel(:,1) )
2087 
2088  ! Linear reconstruction of the physical variables for the internal cells
2089  DO j = 2,comp_cells-1
2090 
2091  x_stencil(1:3) = x_comp(j-1:j+1)
2092 
2093  DO i=1,n_vars
2094 
2095  qp_stencil = qp(i,j-1:j+1)
2096 
2097  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
2098 
2099  dxl = x_comp(j) - x_stag(j)
2100  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2101 
2102  dxr = x_stag(j+1) - x_comp(j)
2103  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2104 
2105  ! positivity preserving reconstruction for h
2106  IF(i.eq.1)THEN
2107 
2108  IF(qpr(i).LT.b_stag(j+1))THEN
2109 
2110  qp_prime=(b_stag(j+1)-qp(i,j))/dxr
2111 
2112  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2113 
2114  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2115 
2116  ENDIF
2117 
2118  IF(qpl(i).LT.b_stag(j))THEN
2119 
2120  qp_prime=(qp(i,j)-b_stag(j))/dxl
2121 
2122  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2123 
2124  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2125 
2126  ENDIF
2127 
2128  ENDIF
2129 
2130  END DO
2131 
2132  ! Convert back from physical to conservative variables
2133  CALL qp_to_qc( qpl , b_stag(j) , q_interfacer(:,j) )
2134  CALL qp_to_qc( qpr , b_stag(j+1) , q_interfacel(:,j+1) )
2135 
2136  END DO
2137 
2138  ! Linear reconstruction of the physical variables at the boundaries
2139  ! of the last cell (j=comp_cells)
2140 
2141  DO i=1,n_vars
2142 
2143  IF ( bcr(i)%flag .EQ. 0 ) THEN
2144 
2145  qp_stencil(1:2) = qp(i,comp_cells-1:comp_cells)
2146  qp_stencil(3) = bcr(i)%value
2147 
2148  x_stencil(1:2) = x_comp(comp_cells-1:comp_cells)
2149  x_stencil(3) = x_stag(comp_interfaces)
2150 
2151  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
2152 
2153  ELSEIF ( bcr(i)%flag .EQ. 1 ) THEN
2154 
2155  qp_prime = bcr(i)%value
2156 
2157  ELSEIF ( bcr(i)%flag .EQ. 2 ) THEN
2158 
2159  qp_prime = ( qp(i,comp_cells) - qp(i,comp_cells-1) ) / &
2160  ( x_comp(comp_cells) - x_comp(comp_cells-1) )
2161 
2162  END IF
2163 
2164  dxl = x_comp(comp_cells) - x_stag(comp_interfaces-1)
2165  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2166 
2167  dxr = x_stag(comp_interfaces) - x_comp(comp_cells)
2168  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2169 
2170  ! positivity preserving reconstruction for h
2171  IF(i.eq.1)THEN
2172 
2173  IF(qpr(i).LT.b_stag(comp_interfaces))THEN
2174 
2175  qp_prime=(b_stag(comp_interfaces)-qp(i,comp_cells))/dxr
2176 
2177  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2178 
2179  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2180 
2181  ENDIF
2182 
2183  IF(qpl(i).LT.b_stag(comp_interfaces-1))THEN
2184 
2185  qp_prime=(qp(i,comp_cells)-b_stag(comp_cells))/dxl
2186 
2187  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2188 
2189  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2190 
2191  ENDIF
2192 
2193  ENDIF
2194 
2195  END DO
2196 
2197  DO i=1,n_vars
2198 
2199  IF ( bcr(i)%flag .EQ. 0 ) THEN
2200 
2201  qpr(i) = bcr(i)%value
2202 
2203  END IF
2204 
2205  END DO
2206 
2207  ! Convert back from physical to conservative variables
2208  CALL qp_to_qc( qpl , b_stag(comp_interfaces-1) , q_interfacer(:,comp_interfaces-1) )
2209  CALL qp_to_qc( qpr , b_stag(comp_interfaces) , q_interfacel(:,comp_interfaces) )
2210 
2211  DO i=1,n_vars
2212 
2213  IF ( bcr(i)%flag .EQ. 0 ) THEN
2214 
2215  qp_bdry(i) = bcr(i)%value
2216 
2217  ELSE
2218 
2219  !qp_bdry(i) = qpR(i)
2220 
2221  ! fixed wall
2222  IF(i.eq.2)THEN
2223 
2224  qp_bdry(i) = -qpr(i)
2225 
2226  ELSE
2227 
2228  qp_bdry(i) = qpr(i)
2229 
2230  ENDIF
2231 
2232  END IF
2233 
2234  END DO
2235 
2236  CALL qp_to_qc( qp_bdry , b_stag(comp_interfaces) , &
2237  q_interfacer(:,comp_interfaces) )
2238 
2239  END SUBROUTINE reconstruction_phys1
2240 
2241  !******************************************************************************
2243  !
2250  !******************************************************************************
2251 
2252  SUBROUTINE reconstruction_phys2
2253 
2254  ! External procedures
2255  USE constitutive, ONLY : qc_to_qp , qp_to_qc
2256  USE constitutive, ONLY : qc_to_qp2 , qp2_to_qc
2257  USE parameters, ONLY : limiter
2258 
2259  ! External variables
2260  USE geometry, ONLY : x_comp , x_stag
2261  USE parameters, ONLY : reconstr_coeff
2262 
2263  IMPLICIT NONE
2264 
2265  REAL*8 :: qc(n_vars)
2266  REAL*8 :: qpl(n_vars)
2267  REAL*8 :: qpr(n_vars)
2268  REAL*8 :: qp_bdry(n_vars)
2269 
2270  REAL*8 :: qp_stencil(3)
2271  REAL*8 :: x_stencil(3)
2272  REAL*8 :: qp_prime
2273 
2274  REAL*8 :: dxl , dxr
2275 
2276  INTEGER :: j
2277  INTEGER :: i
2278 
2279 
2280  ! Convert the conservative variables to the physical variables
2281  DO j = 1,comp_cells
2282 
2283  qc = q(1:n_vars,j)
2284 
2285  CALL qc_to_qp2( qc , b_cent(j) , qp(1:n_vars,j) )
2286 
2287  END DO
2288 
2289  DO i=1,n_vars
2290 
2291  ! Slope in the first cell
2292  IF ( bcl(i)%flag .EQ. 0 ) THEN
2293 
2294  qp_stencil(1) = bcl(i)%value
2295  qp_stencil(2:3) = qp(i,1:2)
2296 
2297  x_stencil(1) = x_stag(1)
2298  x_stencil(2:3) = x_comp(1:2)
2299 
2300  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
2301 
2302  ELSEIF ( bcl(i)%flag .EQ. 1 ) THEN
2303 
2304  qp_prime = bcl(i)%value
2305 
2306  ELSEIF ( bcl(i)%flag .EQ. 2 ) THEN
2307 
2308  qp_prime = ( qp(i,2) - qp(i,1) ) / ( x_comp(2) - x_comp(1) )
2309 
2310  END IF
2311 
2312  ! Linear reconstruction of the physical variables at the boundaries
2313  ! of for the first cell
2314  dxl = x_comp(1) - x_stag(1)
2315  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2316 
2317  dxr = x_stag(2) - x_comp(1)
2318  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2319 
2320  ! positivity preserving reconstruction for h
2321  IF ( i .EQ. 1 ) THEN
2322 
2323  IF ( qpr(i) .LT. 0.d0 ) THEN
2324 
2325  qp_prime = - qp(i,1) / dxr
2326 
2327  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2328 
2329  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2330 
2331  ENDIF
2332 
2333  IF ( qpl(i) .LT. 0.d0 ) THEN
2334 
2335  qp_prime = qp(i,1) / dxl
2336 
2337  qpl(i) = qp(i,1) - reconstr_coeff * dxl * qp_prime
2338 
2339  qpr(i) = qp(i,1) + reconstr_coeff * dxr * qp_prime
2340 
2341  ENDIF
2342 
2343  ENDIF
2344 
2345  END DO
2346 
2347  DO i=1,n_vars
2348 
2349  IF ( bcl(i)%flag .EQ. 0 ) THEN
2350 
2351  qpl(i) = bcl(i)%value
2352 
2353  END IF
2354 
2355  END DO
2356 
2357 
2358  ! Convert back from physical to conservative variables
2359  CALL qp2_to_qc( qpl , b_stag(1) , q_interfacer(1:n_vars,1) )
2360  CALL qp2_to_qc( qpr , b_stag(2) , q_interfacel(1:n_vars,2) )
2361 
2362 
2363  DO i=1,n_vars
2364 
2365  IF ( bcl(i)%flag .EQ. 0 ) THEN
2366 
2367  qp_bdry(i) = bcl(i)%value
2368 
2369  ELSE
2370 
2371  !qp_bdry(i) = qpL(i)
2372 
2373  ! fixed wall
2374  IF(i.eq.2)THEN
2375 
2376  qp_bdry(i) = -qpl(i)
2377 
2378  ELSE
2379 
2380  qp_bdry(i) = qpl(i)
2381 
2382  ENDIF
2383 
2384  END IF
2385 
2386  END DO
2387 
2388  CALL qp2_to_qc( qp_bdry , b_stag(1) , q_interfacel(:,1) )
2389 
2390  ! Linear reconstruction of the physical variables for the internal cells
2391  DO j = 2,comp_cells-1
2392 
2393  x_stencil(1:3) = x_comp(j-1:j+1)
2394 
2395  DO i=1,n_vars
2396 
2397  qp_stencil = qp(i,j-1:j+1)
2398 
2399  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
2400 
2401  dxl = x_comp(j) - x_stag(j)
2402  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2403 
2404  dxr = x_stag(j+1) - x_comp(j)
2405  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2406 
2407  ! positivity preserving reconstruction for h
2408  IF ( i .EQ. 1 ) THEN
2409 
2410  IF ( qpr(i) .LT. 0.d0 ) THEN
2411 
2412  qp_prime = - qp(i,j) / dxr
2413 
2414  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2415 
2416  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2417 
2418  ENDIF
2419 
2420  IF ( qpl(i) .LT. 0.d0 ) THEN
2421 
2422  qp_prime = qp(i,j) / dxl
2423 
2424  qpl(i) = qp(i,j) - reconstr_coeff * dxl * qp_prime
2425 
2426  qpr(i) = qp(i,j) + reconstr_coeff * dxr * qp_prime
2427 
2428  ENDIF
2429 
2430  ENDIF
2431 
2432  END DO
2433 
2434  ! Convert back from physical to conservative variables
2435  CALL qp2_to_qc( qpl , b_stag(j) , q_interfacer(:,j) )
2436  CALL qp2_to_qc( qpr , b_stag(j+1) , q_interfacel(:,j+1) )
2437 
2438  END DO
2439 
2440  ! Linear reconstruction of the physical variables at the boundaries
2441  ! of the last cell (j=comp_cells)
2442 
2443  DO i=1,n_vars
2444 
2445  IF ( bcr(i)%flag .EQ. 0 ) THEN
2446 
2447  qp_stencil(1:2) = qp(i,comp_cells-1:comp_cells)
2448  qp_stencil(3) = bcr(i)%value
2449 
2450  x_stencil(1:2) = x_comp(comp_cells-1:comp_cells)
2451  x_stencil(3) = x_stag(comp_interfaces)
2452 
2453  CALL limit( qp_stencil , x_stencil , limiter(i) , qp_prime )
2454 
2455  ELSEIF ( bcr(i)%flag .EQ. 1 ) THEN
2456 
2457  qp_prime = bcr(i)%value
2458 
2459  ELSEIF ( bcr(i)%flag .EQ. 2 ) THEN
2460 
2461  qp_prime = ( qp(i,comp_cells) - qp(i,comp_cells-1) ) / &
2462  ( x_comp(comp_cells) - x_comp(comp_cells-1) )
2463 
2464  END IF
2465 
2466  dxl = x_comp(comp_cells) - x_stag(comp_interfaces-1)
2467  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2468 
2469  dxr = x_stag(comp_interfaces) - x_comp(comp_cells)
2470  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2471 
2472  ! positivity preserving reconstruction for h
2473  IF ( i .EQ. 1 ) THEN
2474 
2475  IF ( qpr(i) .LT. 0.d0 ) THEN
2476 
2477  qp_prime = - qp(i,comp_cells) / dxr
2478 
2479  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2480 
2481  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2482 
2483  ENDIF
2484 
2485  IF ( qpl(i) .LT. 0.d0 ) THEN
2486 
2487  qp_prime = qp(i,comp_cells) / dxl
2488 
2489  qpl(i) = qp(i,comp_cells) - reconstr_coeff * dxl * qp_prime
2490 
2491  qpr(i) = qp(i,comp_cells) + reconstr_coeff * dxr * qp_prime
2492 
2493  ENDIF
2494 
2495  ENDIF
2496 
2497  END DO
2498 
2499  DO i=1,n_vars
2500 
2501  IF ( bcr(i)%flag .EQ. 0 ) THEN
2502 
2503  qpr(i) = bcr(i)%value
2504 
2505  END IF
2506 
2507  END DO
2508 
2509  ! Convert back from physical to conservative variables
2510  CALL qp2_to_qc( qpl , b_stag(comp_interfaces-1) , q_interfacer(:,comp_interfaces-1) )
2511  CALL qp2_to_qc( qpr , b_stag(comp_interfaces) , q_interfacel(:,comp_interfaces) )
2512 
2513  DO i=1,n_vars
2514 
2515  IF ( bcr(i)%flag .EQ. 0 ) THEN
2516 
2517  qp_bdry(i) = bcr(i)%value
2518 
2519  ELSE
2520 
2521  !qp_bdry(i) = qpR(i)
2522 
2523  ! fixed wall
2524  IF(i.eq.2)THEN
2525 
2526  qp_bdry(i) = -qpr(i)
2527 
2528  ELSE
2529 
2530  qp_bdry(i) = qpr(i)
2531 
2532  ENDIF
2533 
2534  END IF
2535 
2536  END DO
2537 
2538  CALL qp2_to_qc( qp_bdry , b_stag(comp_interfaces) , &
2539  q_interfacer(:,comp_interfaces) )
2540 
2541  END SUBROUTINE reconstruction_phys2
2542 
2543  !******************************************************************************
2545  !
2551  !******************************************************************************
2552 
2553  SUBROUTINE eval_speeds
2554 
2555  ! External procedures
2556  USE constitutive, ONLY : eval_local_speeds2
2557 
2558  IMPLICIT NONE
2559 
2560  REAL*8 :: abslambdal_min(n_vars) , abslambdal_max(n_vars)
2561  REAL*8 :: abslambdar_min(n_vars) , abslambdar_max(n_vars)
2562  REAL*8 :: min_r(n_vars) , max_r(n_vars)
2563 
2564  INTEGER :: j
2565 
2566  DO j = 1 , comp_interfaces
2567 
2568  CALL eval_local_speeds2( q_interfacel(:,j) , b_stag(j) , abslambdal_min , &
2569  abslambdal_max )
2570  CALL eval_local_speeds2( q_interfacer(:,j) , b_stag(j) , abslambdar_min , &
2571  abslambdar_max )
2572 
2573  min_r = min(abslambdal_min , abslambdar_min , 0.0d0)
2574  max_r = max(abslambdal_max , abslambdar_max , 0.0d0)
2575 
2576  a_interfacel(:,j) = min_r
2577  a_interfacer(:,j) = max_r
2578 
2579  END DO
2580 
2581  END SUBROUTINE eval_speeds
2582 
2583 
2584 
2585  !******************************************************************************
2587  !
2601  !******************************************************************************
2602 
2603  SUBROUTINE limit( v , z , limiter , slope_lim )
2604 
2605  USE parameters, ONLY : theta
2606 
2607  IMPLICIT none
2608 
2609  REAL*8, INTENT(IN) :: v(3)
2610  REAL*8, INTENT(IN) :: z(3)
2611  INTEGER, INTENT(IN) :: limiter
2612 
2613  REAL*8, INTENT(OUT) :: slope_lim
2614 
2615  REAL*8 :: a , b , c
2616 
2617  REAL*8 :: sigma1 , sigma2
2618 
2619  a = ( v(3) - v(2) ) / ( z(3) - z(2) )
2620  b = ( v(2) - v(1) ) / ( z(2) - z(1) )
2621  c = ( v(3) - v(1) ) / ( z(3) - z(1) )
2622 
2623  SELECT CASE (limiter)
2624 
2625  CASE ( 0 )
2626 
2627  slope_lim = 0.d0
2628 
2629  CASE ( 1 )
2630 
2631  ! minmod
2632 
2633  slope_lim = minmod(a,b)
2634 
2635  CASE ( 2 )
2636 
2637  ! superbee
2638 
2639  sigma1 = minmod( a , 2.d0*b )
2640  sigma2 = minmod( 2.d0*a , b )
2641  slope_lim = maxmod( sigma1 , sigma2 )
2642 
2643  CASE ( 3 )
2644 
2645  ! van_leer
2646 
2647  slope_lim = minmod( c , theta * minmod( a , b ) )
2648 
2649  END SELECT
2650 
2651  END SUBROUTINE limit
2652 
2653 
2654  REAL*8 FUNCTION minmod(a,b)
2655 
2656  IMPLICIT none
2657 
2658  REAL*8 :: a , b , sa , sb
2659 
2660  IF ( a*b .LE. 0.d0 ) THEN
2661 
2662  minmod = 0.d0
2663 
2664  ELSE
2665 
2666  sa = a / abs(a)
2667  sb = b / abs(b)
2668 
2669  minmod = 0.5 * ( sa+sb ) * min( abs(a) , abs(b) )
2670 
2671  END IF
2672 
2673  END FUNCTION minmod
2674 
2675  REAL*8 function maxmod(a,b)
2676 
2677  IMPLICIT none
2678 
2679  REAL*8 :: a , b , sa , sb
2680 
2681  IF ( a*b .EQ. 0.d0 ) THEN
2682 
2683  maxmod = 0.d0
2684 
2685  ELSE
2686 
2687  sa = a / abs(a)
2688  sb = b / abs(b)
2689 
2690  maxmod = 0.5 * ( sa+sb ) * max( abs(a) , abs(b) )
2691 
2692  END IF
2693 
2694  END function maxmod
2695 
2696 END MODULE solver
solver::imex_rk_solver
subroutine imex_rk_solver
Runge-Kutta integration.
Definition: solver.f90:414
constitutive::qc_to_qp2
subroutine qc_to_qp2(qc, B, qp)
Conservative to 2nd set of physical variables.
Definition: constitutive.f90:241
solver::timestep
subroutine timestep
Time-step computation.
Definition: solver.f90:322
solver
Numerical solver.
Definition: solver.f90:12
constitutive::eval_explicit_forces
subroutine eval_explicit_forces(Bj, Bprimej, qj, expl_forces_term)
Explicit Forces term.
Definition: constitutive.f90:453
solver::eval_f
subroutine eval_f(Bj, Bprimej, qj, qj_old, a_tilde, a_dirk, a_diag, coeff_f, f_nl, scal_f)
Evaluate the nonlinear system.
Definition: solver.f90:1145
solver::eval_flux_kt
subroutine eval_flux_kt
Semidiscrete numerical fluxes.
Definition: solver.f90:1428
solver::reconstruction_phys2
subroutine reconstruction_phys2
Linear reconstruction.
Definition: solver.f90:2252
solver::eval_hyperbolic_terms
subroutine eval_hyperbolic_terms(q_expl, F_x)
Semidiscrete finite volume central scheme.
Definition: solver.f90:1303
solver::solve_rk_step
subroutine solve_rk_step(Bj, Bprimej, qj, qj_old, a_tilde, a_dirk, a_diag)
Runge-Kutta single step integration.
Definition: solver.f90:605
constitutive::qp2_to_qc
subroutine qp2_to_qc(qp, B, qc)
From 2nd set of physical to conservative variables.
Definition: constitutive.f90:271
solver::eval_flux_lxf
subroutine eval_flux_lxf
Numerical fluxes Lax-Friedrichs.
Definition: solver.f90:1605
solver::deallocate_solver_variables
subroutine deallocate_solver_variables
Memory deallocation.
Definition: solver.f90:270
constitutive::qp_to_qc
subroutine qp_to_qc(qp, B, qc)
Physical to conservative variables.
Definition: constitutive.f90:206
geometry
Grid module.
Definition: geometry.f90:7
constitutive::eval_fluxes
subroutine eval_fluxes(Bj, c_qj, r_qj, c_flux, r_flux)
Hyperbolic Fluxes.
Definition: constitutive.f90:304
parameters
Parameters.
Definition: parameters.f90:7
constitutive
Constitutive equations.
Definition: constitutive.f90:4
solver::limit
subroutine limit(v, z, limiter, slope_lim)
Slope limiter.
Definition: solver.f90:2603
solver::allocate_solver_variables
subroutine allocate_solver_variables
Memory allocation.
Definition: solver.f90:112
solver::lnsrch
subroutine lnsrch(Bj, Bprimej, qj_rel_NR_old, qj_org, qj_old, scal_f_old, grad_f, desc_dir, coeff_f, qj_rel, scal_f, right_term, stpmax, check)
Search the descent stepsize.
Definition: solver.f90:924
constitutive::qc_to_qp
subroutine qc_to_qp(qc, B, qp)
Conservative to physical variables.
Definition: constitutive.f90:177
solver::reconstruction_phys1
subroutine reconstruction_phys1
Linear reconstruction.
Definition: solver.f90:1950
solver::reconstruction_cons
subroutine reconstruction_cons
Linear reconstruction.
Definition: solver.f90:1664
solver::maxmod
real *8 function maxmod(a, b)
Definition: solver.f90:2675
solver::eval_jacobian
subroutine eval_jacobian(Bj, Bprimej, qj_rel, qj_org, coeff_f, left_matrix)
Evaluate the jacobian.
Definition: solver.f90:1197
solver::eval_explicit_terms
subroutine eval_explicit_terms(q_expl, expl_terms)
Evaluate the explicit terms.
Definition: solver.f90:1265
solver::average_kt
subroutine average_kt(aL, aR, wL, wR, w_avg)
Definition: solver.f90:1499
solver::eval_speeds
subroutine eval_speeds
Characteristic speeds.
Definition: solver.f90:2553
solver::eval_flux_gforce
subroutine eval_flux_gforce
Numerical fluxes GFORCE.
Definition: solver.f90:1536
constitutive::eval_local_speeds2
subroutine eval_local_speeds2(qj, Bj, vel_min, vel_max)
Local Characteristic speeds.
Definition: constitutive.f90:129
solver::minmod
real *8 function minmod(a, b)
Definition: solver.f90:2654
constitutive::eval_nonhyperbolic_terms
subroutine eval_nonhyperbolic_terms(Bj, Bprimej, c_qj, c_nh_term_impl, r_qj, r_nh_term_impl)
Non-Hyperbolic terms.
Definition: constitutive.f90:378
constitutive::eval_local_speeds
subroutine eval_local_speeds(qj, Bj, vel_min, vel_max)
Local Characteristic speeds.
Definition: constitutive.f90:105
Generated on Tue Nov 22 2016 16:23:06 for IMEXSfloW by   doxygen 1.8.7