compute_face_states Subroutine

  1. weno.f90
  2. weno
  3. compute_face_states

private subroutine compute_face_states(qp, f_qp_left, f_qp_right, flags, dims)

Subroutine to calculate state at the face, generalized for all direction : I,J, and K.

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var):: qp

Store primitive variable at cell center
real(kind=wp), intent(inout), dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var):: f_qp_left

primitive state variable at faces
real(kind=wp), intent(inout), dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var):: f_qp_right

primitive state variable at faces
integer, intent(in), dimension(3):: flags

flags for direction switch
type(extent), intent(in) :: dims

Extent of the domain:imx,jmx,kmx

Called by

proc~~compute_face_states~~CalledByGraph proc~compute_face_states compute_face_states proc~compute_weno_states compute_weno_states proc~compute_weno_states->proc~compute_face_states proc~compute_face_interpolant compute_face_interpolant proc~compute_face_interpolant->proc~compute_weno_states proc~get_total_conservative_residue get_total_conservative_Residue proc~get_total_conservative_residue->proc~compute_face_interpolant proc~get_next_solution get_next_solution proc~get_next_solution->proc~get_total_conservative_residue proc~iterate_one_more_time_step iterate_one_more_time_step proc~iterate_one_more_time_step->proc~get_next_solution program~main main program~main->proc~iterate_one_more_time_step

Contents

Source Code

compute_face_states

Source Code

        subroutine compute_face_states(qp, f_qp_left, f_qp_right, flags, dims)
          !< Subroutine to calculate state at the face, generalized for
          !< all direction : I,J, and K.
            implicit none

            type(extent), intent(in) :: dims
            !< Extent of the domain:imx,jmx,kmx
            integer, dimension(3), intent(in) :: flags
            !< flags for direction switch
            real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(in) :: qp
            !< Store primitive variable at cell center
            real(wp), dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var), intent(inout) :: f_qp_left, f_qp_right
            !< primitive state variable at faces
            integer :: i, j, k, l
            integer :: i_f=0, j_f=0, k_f=0
            real(wp), dimension(3) :: P !< polynomial approximation
            real(wp), dimension(3) :: B !< smoothness factor
            real(wp), dimension(3) :: w !< wieght
            real(wp), dimension(3) :: g !< linear wieght
            real(wp), dimension(-2:2) :: u !< state_variable
            real(wp)               :: eps=1e-6

            g(1) = 1.0/10.0
            g(2) = 6.0/10.0
            g(3) = 3.0/10.0


            i_f = flags(1)
            j_f = flags(2)
            k_f = flags(3)

            do l = 1, dims%n_var
             do k = 1-k_f, dims%kmx-1+k_f
              do j = 1-j_f, dims%jmx-1+j_f
               do i = 1-i_f, dims%imx-1+i_f
                 U(-2) = qp(i-2*i_f,j-2*j_f,k-2*k_f,l)  !u_{i-2}
                 U(-1) = qp(i-1*i_f,j-1*j_f,k-1*k_f,l)  !u_{i-1}
                 u( 0) = qp(i      ,j      ,k      ,l)  !u_{i}
                 U( 1) = qp(i+1*i_f,j+1*j_f,k+1*k_f,l)  !u_{i+1}
                 U( 2) = qp(i+2*i_f,j+2*j_f,k+2*k_f,l)  !u_{i+2}

                 P(1) = ( 2.0*U(-2) -  7.0*U(-1) + 11.0*U(0))/6.0
                 P(2) = (-1.0*U(-1) +  5.0*U( 0) +  2.0*U(1))/6.0
                 P(3) = ( 2.0*U( 0) +  5.0*U( 1) -  1.0*U(2))/6.0

                 B(1) = (13.0/12.0)*(U(-2)-2.0*U(-1)+U(0))**2 + (1.0/4.0)*(    U(-2)-4.0*U(-1)+3.0*U(0))**2
                 B(2) = (13.0/12.0)*(U(-1)-2.0*U( 0)+U(1))**2 + (1.0/4.0)*(    U(-1)-              U(1))**2
                 B(3) = (13.0/12.0)*(U( 0)-2.0*U( 1)+U(2))**2 + (1.0/4.0)*(3.0*U( 0)-4.0*U( 1)+    U(2))**2

                 w(:) = g(:)/(eps + B(:))**2 
                 f_qp_left(i+i_f,j+j_f,k+k_f,l)  = SUM(w*P)/SUM(w)

                 P(1) = ( 2.0*U(2) -  7.0*U( 1) + 11.0*U( 0))/6.0 
                 P(2) = (-1.0*U(1) +  5.0*U( 0) +  2.0*U(-1))/6.0
                 P(3) = ( 2.0*U(0) +  5.0*U(-1) -  1.0*U(-2))/6.0

                 !B(1) = (13.0/12.0)*(U( 2)-2.0*U( 1)+U( 0))**2 + (1.0/4.0)*(    U(2)-4.0*U( 1)+3.0*U( 0))**2
                 !B(2) = (13.0/12.0)*(U( 1)-2.0*U( 0)+U(-1))**2 + (1.0/4.0)*(    U(1)-              U(-1))**2
                 !B(3) = (13.0/12.0)*(U( 0)-2.0*U(-1)+U(-2))**2 + (1.0/4.0)*(3.0*U(0)-4.0*U(-1)+    U(-2))**2

                 w(1) = g(1)/(eps + B(3))**2 
                 w(2) = g(2)/(eps + B(2))**2 
                 w(3) = g(3)/(eps + B(1))**2 
                 f_qp_right(i,j,k,l) = SUM(w*P)/SUM(w)
               end do
              end do
             end do
            end do

        end subroutine compute_face_states