init_2d.f90

Go to the documentation of this file.

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

MODULE init_2d

  USE parameters_2d, ONLY : temperature_flag , verbose_level ,                  &
       solid_transport_flag

  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 :: xs_w
  REAL*8 :: t_w

  REAL*8 :: hb_e
  REAL*8 :: u_e
  REAL*8 :: v_e
  REAL*8 :: xs_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

    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

    REAL*8 :: eps

    IF ( verbose_level .GE. 1 ) WRITE(*,*) 'Riemann problem initialization'

    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

    IF ( solid_transport_flag ) THEN

       qp(4,1:i1,:) = xs_w

       IF ( temperature_flag ) qp(5,1:i1,:) = t_w

    ELSE

       IF ( temperature_flag ) qp(4,1:i1,:) = t_w

    END IF


    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

    IF ( solid_transport_flag ) THEN

       qp(4,i1+1:comp_cells_x,:) = xs_e

       IF ( temperature_flag ) qp(5,i1+1:comp_cells_x,:) = t_e

    ELSE

       IF ( temperature_flag ) qp(4,i1+1:comp_cells_x,:) = t_e

    END IF


    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 initial_conditions

    USE constitutive_2d, ONLY : qp_to_qc

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

    USE geometry_2d, ONLY : dx , dy

    USE parameters_2d, ONLY : n_vars , released_volume , verbose_level ,        &
         source_flag

    USE solver_2d, ONLY : q

    IMPLICIT none

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

    INTEGER :: j,k

    DO j=1,comp_cells_x

      DO k=1,comp_cells_y

         qp(1,j,k) = thickness_function(x_comp(j),y_comp(k),b_cent(j,k))

         qp(2,j,k) = velocity_u_function(x_comp(j),y_comp(k))

         qp(3,j,k) = velocity_v_function(x_comp(j),y_comp(k))

         IF ( solid_transport_flag ) THEN

            qp(4,j,k) = sediment_function(x_comp(j),y_comp(k))

            IF ( temperature_flag ) THEN
               
               qp(5,j,k) = temperature_function(x_comp(j),y_comp(k))
               
            END IF

         ELSE

            IF ( temperature_flag ) THEN
               
               qp(4,j,k) = temperature_function(x_comp(j),y_comp(k))
               
            END IF

         END IF
            
      ENDDO

    ENDDO

    IF ( released_volume .GT. ( dx * dy * sum( qp(1,:,:)-b_cent(:,:) ) ) ) THEN

       ! Correction for the released volume
       qp(1,:,:) = b_cent(:,:) + ( qp(1,:,:)-b_cent(:,:) ) * released_volume    &
            / ( dx * dy * sum( qp(1,:,:)-b_cent(:,:) ) )
       
    END IF

    WRITE(*,*) 'Initial volume =',dx * dy * sum( qp(1,:,:)-b_cent(:,:) ) 


    DO j = 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 ) WRITE(*,*) j,k,b_cent(j,k),qp(:,j,k)
    
      END DO

    ENDDO

    IF ( verbose_level .GE. 1 ) READ(*,*)
    
    IF ( source_flag ) CALL init_source
    
    RETURN

  END SUBROUTINE initial_conditions

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

  SUBROUTINE init_source

    USE geometry_2d, ONLY : x_comp , comp_cells_x, y_comp , comp_cells_y
    USE geometry_2d, ONLY : dx , dy
    
    USE parameters_2d, ONLY : x_source , y_source , r_source
    USE parameters_2d, ONLY : vfr_source , vel_source

    USE solver_2d, ONLY : source_xy
    
    IMPLICIT NONE

    INTEGER :: h,j,k
    
    REAL*8, ALLOCATABLE :: x_subgrid(:) , y_subgrid(:)
    INTEGER, ALLOCATABLE :: check_subgrid(:)

    INTEGER n_points , n_points2
    REAL*8 :: source_area
    
    n_points = 200
    n_points2 = n_points**2

    ALLOCATE( x_subgrid(n_points2) )
    ALLOCATE( y_subgrid(n_points2) )
    ALLOCATE( check_subgrid(n_points2) )

    x_subgrid = 0.d0
    y_subgrid = 0.d0
    
    DO h = 1,n_points

       x_subgrid(h:n_points2:n_points) = h
       y_subgrid((h-1)*n_points+1:h*n_points) = h

    END DO

    x_subgrid = x_subgrid / ( n_points +1 )
    y_subgrid = y_subgrid / ( n_points +1 )
    
    x_subgrid = ( x_subgrid - 0.5d0 ) * dx
    y_subgrid = ( y_subgrid - 0.5d0 ) * dy
    
    source_xy(:,:) = 0.d0
    
    DO j=1,comp_cells_x

      DO k=1,comp_cells_y

         check_subgrid = 0
         
         WHERE ( ( x_comp(j) + x_subgrid - x_source )**2                        &
              + ( y_comp(k) + y_subgrid - y_source )**2 < r_source**2 )
         
            check_subgrid = 1

         END WHERE

         ! WRITE(*,*) x_comp(j),y_comp(k),REAL(SUM(check_subgrid))/10201.D0

         source_xy(j,k) = REAL(sum(check_subgrid))/n_points2
         
      ENDDO

    ENDDO

    source_area = dx*dy*sum(source_xy)
    
    WRITE(*,*) 'Source area =',source_area,' Error =',abs( 1.d0 -      &
         dx*dy*sum(source_xy) / ( 4.d0*atan(1.d0)*r_source**2 ) )

    vel_source = vfr_source / source_area

    WRITE(*,*) 'Vel source =',vel_source
        
  END SUBROUTINE init_source
  
  !******************************************************************************
  !
  !******************************************************************************

  REAL*8 FUNCTION thickness_function(x,y,Bj)

    USE parameters_2d, ONLY : released_volume , x_release , y_release

    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x,y
    REAL*8, INTENT(IN) :: Bj
    
    REAL*8, PARAMETER :: pig = 4.d0*atan(1.d0)
    REAL*8 :: R

    ! example 1D from Kurganov and Petrova 2007    
    !IF(y.LE.0.d0)THEN
    !   thickness_function=0.4d0+Bj
    !ELSE
    !   thickness_function=Bj
    !ENDIF

    ! example 2D from Kurganov and Petrova 2007    
    ! thickness_function=MAX(0.25,Bj)
    
    r = ( released_volume / pig )**(1.0/3.0)
    
    IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r ) THEN

      thickness_function = bj + r

    ELSE

      thickness_function = bj

    ENDIF

  END FUNCTION thickness_function

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

  REAL*8 FUNCTION velocity_u_function(x,y)
  
    USE parameters_2d, ONLY : released_volume , x_release , y_release
    USE parameters_2d, ONLY : velocity_mod_release, velocity_ang_release

    USE geometry_2d, ONLY : x0 , y0 , dx , dy
    USE geometry_2d, ONLY : x_stag , y_stag , b_ver
    
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x,y
    
    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8 :: R

    INTEGER :: idx1 , idy1
    REAL*8 :: coeff_x , coeff_y

    REAL*8 :: dBdx1 , dBdx2 , dBdx
    REAL*8 :: dBdy1 , dBdy2 , dBdy

    REAL*8 :: max_slope_angle_rad
    REAL*8 :: velocity_ang_release_rad

    
    r = ( released_volume / pig )**(1.d0/3.d0)

    ! indexes of the lower-left vertex of the cell in the computational grid
    idx1 = ceiling( ( x_release - x0 ) / dx )
    idy1 = ceiling( ( y_release - y0 ) / dy )

    ! relative coordinates of (x_release,y_release) within the cell
    coeff_x = ( x_release - x_stag(idx1) ) / dx
    coeff_y = ( y_release - y_stag(idy1) ) / dy


    ! x-derivatives at bottom and top edges of the cell
    dbdx1 = ( b_ver(idx1+1,idy1) - b_ver(idx1,idy1) ) / dx
    dbdx2 = ( b_ver(idx1+1,idy1+1) - b_ver(idx1,idy1+1) ) / dx

    ! x-derivative of terrain at (x_release,y_release)
    dbdx = coeff_y * dbdx2 + ( 1.d0 - coeff_y ) * dbdx1
    
    ! y-derivatives at left and right edges of the cell
    dbdy1 = ( b_ver(idx1,idy1+1) - b_ver(idx1,idy1) ) / dy
    dbdy2 = ( b_ver(idx1+1,idy1+1) - b_ver(idx1+1,idy1) ) / dx

    ! y-derivative of terrain at (x_release,y_release)
    dbdy = coeff_x * dbdy2 + ( 1.d0 - coeff_x ) * dbdy1


    ! direction of maximum slope in radians
    max_slope_angle_rad = datan2(dbdy,dbdx)

    IF ( verbose_level .GE. 1 ) THEN

       WRITE(*,*) x0,x_stag(1),y0,y_stag(1),dx,dy
       WRITE(*,*) '---',x_stag(idx1),x_release,x_stag(idx1+1),coeff_x
       WRITE(*,*) '---',y_stag(idy1),y_release,y_stag(idy1+1),coeff_y
       
       WRITE(*,*) b_ver(idx1,idy1),b_ver(idx1+1,idy1)
       WRITE(*,*) b_ver(idx1,idy1+1),b_ver(idx1+1,idy1+1)
       WRITE(*,*) 'dbdx,dbdy,angle',dbdx,dbdy,max_slope_angle_rad*180.d0/pig
    
       READ(*,*)

    END IF
    
    IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r ) THEN

       velocity_ang_release_rad = velocity_ang_release * pig / 180.d0
       
       IF ( velocity_ang_release .GE. 0.d0 ) THEN 

          ! angle departing from positive x-axis
          velocity_u_function = velocity_mod_release                            &
               * cos( velocity_ang_release * ( 2.d0 * pig / 360.d0 ) )

       ELSE

          ! angle departing from maximum slope direction
          velocity_u_function = velocity_mod_release                            &
               * cos( max_slope_angle_rad + velocity_ang_release_rad )

       END IF
          
    ELSE

      velocity_u_function = 0.d0

    ENDIF

  END FUNCTION velocity_u_function

  !******************************************************************************
  !
  !******************************************************************************
  
  REAL*8 FUNCTION velocity_v_function(x,y)
  
    USE parameters_2d, ONLY : released_volume , x_release , y_release
    USE parameters_2d, ONLY : velocity_mod_release, velocity_ang_release

    USE geometry_2d, ONLY : x0 , y0 , dx , dy
    USE geometry_2d, ONLY : x_stag , y_stag , b_ver

    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x,y
   
    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8 :: R
    
    INTEGER :: idx1 , idy1
    REAL*8 :: coeff_x , coeff_y
    
    REAL*8 :: dBdx1 , dBdx2 , dBdx
    REAL*8 :: dBdy1 , dBdy2 , dBdy

    REAL*8 :: max_slope_angle_rad
    REAL*8 :: velocity_ang_release_rad


    r = ( released_volume / pig )**(1.0/3.0)

    ! indexes of the lower-left vertex of the cell in the computational grid
    idx1 = ceiling( ( x_release - x0 ) / dx )
    idy1 = ceiling( ( y_release - y0 ) / dy )

    ! relative coordinates of (x_release,y_release) within the cell
    coeff_x = ( x_release - x_stag(idx1) ) / dx
    coeff_y = ( y_release - y_stag(idy1) ) / dy


    ! x-derivatives at bottom and top edges of the cell
    dbdx1 = ( b_ver(idx1+1,idy1) - b_ver(idx1,idy1) ) / dx
    dbdx2 = ( b_ver(idx1+1,idy1+1) - b_ver(idx1,idy1+1) ) / dx

    ! x-derivative of terrain at (x_release,y_release)
    dbdx = coeff_y * dbdx2 + ( 1.d0 - coeff_y ) * dbdx1
    
    ! y-derivatives at left and right edges of the cell
    dbdy1 = ( b_ver(idx1,idy1+1) - b_ver(idx1,idy1) ) / dy
    dbdy2 = ( b_ver(idx1+1,idy1+1) - b_ver(idx1+1,idy1) ) / dx

    ! y-derivative of terrain at (x_release,y_release)
    dbdy = coeff_x * dbdy2 + ( 1.d0 - coeff_x ) * dbdy1


    
    IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r ) THEN

       velocity_ang_release_rad = velocity_ang_release * pig / 180.d0
       
       IF ( velocity_ang_release .GE. 0.d0 ) THEN 
          
          ! angle departing from positive x-axis
          velocity_v_function = velocity_mod_release                            &
               * sin( velocity_ang_release * ( 2.d0 * pig / 360.d0 ) )
          
       ELSE
          
          ! angle departing from maximum slope direction
          velocity_v_function = velocity_mod_release                            &
               * sin( max_slope_angle_rad + velocity_ang_release_rad )
          
       END IF

    ELSE

      velocity_v_function = 0.d0

    ENDIF

  END FUNCTION velocity_v_function
  
  !******************************************************************************
  !
  !******************************************************************************

  REAL*8 FUNCTION sediment_function(x,y)
  
    USE parameters_2d, ONLY : released_volume , x_release , y_release
    USE parameters_2d, ONLY : xs_init , xs_ambient
  
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x,y
   
    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8 :: R
   
    r = ( released_volume / pig )**(1.0/3.0)
    
    IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r ) THEN

       sediment_function = xs_init

    ELSE

       sediment_function = xs_ambient

    ENDIF

  END FUNCTION sediment_function
  


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

  REAL*8 FUNCTION temperature_function(x,y)
  
    USE parameters_2d, ONLY : released_volume , x_release , y_release
    USE parameters_2d, ONLY : t_init , t_ambient
  
    IMPLICIT NONE
    
    REAL*8, INTENT(IN) :: x,y
   
    REAL*8, PARAMETER :: pig = 4.0*atan(1.0)
    REAL*8 :: R
   
    r = ( released_volume / pig )**(1.0/3.0)
    
    IF ( dsqrt( (x-x_release)**2 + (y-y_release)**2 ) .LE. r ) THEN

       temperature_function = t_init

    ELSE

       temperature_function = t_ambient

    ENDIF

  END FUNCTION temperature_function
  
END MODULE init_2d