init_2d.f90

Go to the documentation of this file.

!********************************************************************************
!
!********************************************************************************

MODULE init_2d

  USE parameters_2d, ONLY : verbose_level
  USE parameters_2d, ONLY : n_solid

  IMPLICIT none

  REAL*8, ALLOCATABLE :: q_init(:,:,:)

  REAL*8, ALLOCATABLE :: thickness_init(:,:)
  
  REAL*8 :: riemann_interface  

  REAL*8 :: hb_w
  REAL*8 :: u_w
  REAL*8 :: v_w
  REAL*8,ALLOCATABLE :: alphas_w(:)
  REAL*8 :: t_w

  REAL*8 :: hb_e
  REAL*8 :: u_e
  REAL*8 :: v_e
  REAL*8,ALLOCATABLE :: alphas_e(:)
  REAL*8 :: t_e


CONTAINS

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE riemann_problem

    USE constitutive_2d, ONLY : qp_to_qc

    USE geometry_2d, ONLY : x_comp , comp_cells_x , comp_cells_y , b_cent

    ! USE geometry_2d, ONLY : x0 , xN , y0 , yN

    USE parameters_2d, ONLY : n_vars , verbose_level , n_solid

    USE solver_2d, ONLY : q

    IMPLICIT none

    ! REAL*8 :: hB            !< height + topography
    ! REAL*8 :: u             !< velocity
    ! REAL*8 :: v             !< velocity

    REAL*8 :: qp(n_vars,comp_cells_x,comp_cells_y) , qj(n_vars)

    INTEGER :: j,k
    INTEGER :: i1

    INTEGER :: i_solid

    REAL*8 :: eps

    IF ( verbose_level .GE. 1 ) THEN

       WRITE(*,*) 'Riemann problem initialization'
       WRITE(*,*) 'x_comp(1)',x_comp(1)
       WRITE(*,*) 'x_comp(comp_cells_x)',x_comp(comp_cells_x)
       WRITE(*,*) 'riemann_interface',riemann_interface

    END IF


    i1 = 0
    
    riemann_int_search:DO j = 1,comp_cells_x
       
       IF ( x_comp(j) .LT. riemann_interface ) THEN

          i1 = j

       ELSE

          EXIT riemann_int_search

       END IF

    END DO riemann_int_search
    
    eps = 1.d-10

    ! Left initial state
    qp(1,1:i1,:) = hb_w
    qp(2,1:i1,:) = u_w
    qp(3,1:i1,:) = v_w
    qp(4,1:i1,:) = t_w

    ALLOCATE( alphas_w(n_solid) )

    DO i_solid=1,n_solid

       qp(4+i_solid,1:i1,:) = alphas_w(i_solid)

    END DO

    qp(n_vars+1,1:i1,:) = 0.d0
    qp(n_vars+2,1:i1,:) = 0.d0

    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'Left state'

    DO j = 1,i1

       DO k = 1,comp_cells_y

         ! evaluate the vector of conservative variables
         CALL qp_to_qc( qp(:,j,k) , b_cent(j,k) , qj )

         q(1:n_vars,j,k) = qj

         IF ( verbose_level .GE. 1 ) THEN 
            
            WRITE(*,*) j,k,b_cent(j,k)
            WRITE(*,*) qp(:,j,k)
            WRITE(*,*) q(1:n_vars,j,k)

         END IF

       ENDDO

    END DO

    IF ( verbose_level .GE. 1 ) READ(*,*)

    ! Right initial state
    qp(1,i1+1:comp_cells_x,:) = hb_e
    qp(2,i1+1:comp_cells_x,:) = u_e
    qp(3,i1+1:comp_cells_x,:) = v_e
    qp(4,i1+1:comp_cells_x,:) = t_e

    ALLOCATE( alphas_e(n_solid) )

    DO i_solid=1,n_solid

       qp(4+i_solid,1:i1,:) = alphas_e(i_solid)

    END DO

    qp(n_vars+1,i1+1:comp_cells_x,:) = 0.d0
    qp(n_vars+2,i1+1:comp_cells_x,:) = 0.d0


    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'Right state'

    DO j = i1+1,comp_cells_x

       DO k = 1,comp_cells_y

         ! evaluate the vector of conservative variables
         CALL qp_to_qc( qp(:,j,k) , b_cent(j,k) , qj )

         q(1:n_vars,j,k) = qj

         IF ( verbose_level .GE. 1 ) THEN 
            
            WRITE(*,*) j,k,b_cent(j,k)
            WRITE(*,*) qp(:,j,k)
            WRITE(*,*) q(1:n_vars,j,k)

         END IF
    
      END DO

    ENDDO

    IF ( verbose_level .GE. 1 ) READ(*,*)

    RETURN

  END SUBROUTINE riemann_problem

  !******************************************************************************
  !
  !
  !
  !******************************************************************************

  SUBROUTINE collapsing_volume

    USE constitutive_2d, ONLY : qp_to_qc

    USE geometry_2d, ONLY : compute_cell_fract

    USE geometry_2d, ONLY : comp_cells_x , comp_cells_y , b_cent

    USE parameters_2d, ONLY : n_vars

    USE parameters_2d, ONLY : x_collapse , y_collapse , r_collapse , t_collapse , &
       h_collapse , alphas_collapse

    USE solver_2d, ONLY : q

    IMPLICIT NONE

    INTEGER :: j,k

    REAL*8 :: qp_init(n_vars+2) ,  qp0_init(n_vars+2)

    REAL*8 :: cell_fract(comp_cells_x,comp_cells_y)

    CALL compute_cell_fract(x_collapse,y_collapse,r_collapse,cell_fract)

    ! values outside the collapsing volume
    qp0_init(1) = 0.d0                  ! h
    qp0_init(2) = 0.d0                  ! hu
    qp0_init(3) = 0.d0                  ! hv
    qp0_init(4) = t_collapse            ! T
    qp0_init(5:4+n_solid) = 0.d0        ! alphas
    qp0_init(n_vars+1:n_vars+2) = 0.d0  ! u,v

    ! values within the collapsing volume
    qp_init(1) = h_collapse
    qp_init(2) = 0.d0
    qp_init(3) = 0.d0
    qp_init(4) = t_collapse
    qp_init(5:4+n_solid) = alphas_collapse(1:n_solid)
    qp_init(n_vars+1:n_vars+2) = 0.d0
    
    DO j = 1,comp_cells_x
       
       DO k = 1,comp_cells_y
          
          IF ( cell_fract(j,k) .GT. 0.d0 ) THEN
             
             qp_init(1) = cell_fract(j,k) * h_collapse

             CALL qp_to_qc( qp_init(1:n_vars+2) , b_cent(j,k) , q(1:n_vars,j,k) )
             
          ELSE
             
             CALL qp_to_qc( qp0_init(1:n_vars+2) , b_cent(j,k) , q(1:n_vars,j,k) )
             
          END IF
          
       END DO
       
    END DO

    RETURN

  END SUBROUTINE collapsing_volume

END MODULE init_2d