H.2.7 Evaluate_Ortho_Diffusion Procedure

The main documentation of the Evaluate_Ortho_Diffusion Procedure contains additional explanation of this code listing.

  subroutine Evaluate_Ortho_Diffusion (Diff_Term, Phi_MV, Div_D_Grad_Phi_Cells)

    ! Input variables.

    type(Ortho_Diffusion_type), intent(inout) :: Diff_Term ! Diffusion term.
    type(Mathematic_Vector_type), intent(inout) :: Phi_MV  ! Phi math vector.

    ! Output variable.

    type(real,1) :: Div_D_Grad_Phi_Cells  ! Div D Grad at each cell center.

    ! Internal variables.

    type(integer) :: face, other_cell               ! Loop parameters.
    type(real,3) :: Area_Faces_of_Cells, Beta
    type(real,2) :: Pseudo_Ortho_Area_F_of_C
    type(real,2) :: Harmonic_Diffusion_Coef_F_of_C  ! Harmonic D at faces.
    type(real,2) :: DeltaR21_Cells_of_Cells, DeltaR1f_Cells_of_Cells
    type(integer) :: Faces_per_Cell                 ! Number of faces per cell.
    type(integer) :: NDimensions                    ! Number of dimensions.
    type(real,1) :: Volume_Cells                    ! Volume of the cells.

    ! Phi manipulation variables.
    type(real,1) :: Phi_Cells
    type(Distributed_Vector_type) :: Phi_Cells_DV
    type(Collected_Array_type) :: Phi_Cells_of_Cells_CA
    type(real,2) :: Phi_Cells_of_Cells

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Verify requirements.

    VERIFY(Valid_State(Diff_Term),5)            ! Diff_Term is valid.
    VERIFY(Valid_State(Phi_MV),5)               ! Phi_MV is valid.
    VERIFY(Valid_State(Div_D_Grad_Phi_Cells),5) ! Div_D_Grad_Phi_Cells is valid.
    ! Div_D_Grad_Phi_Cells is dimensioned correctly.
    VERIFY(SIZE(Div_D_Grad_Phi_Cells,1) == NCells_PE(Diff_Term%Mesh),5) 

    ! Query the mesh.

    Faces_per_Cell = Get_Faces_per_Cell (Diff_Term%Mesh)
    NDimensions = Get_NDimensions (Diff_Term%Mesh)

    ! Future mods: 
    !   1. Access routine for Mesh%Cells_of_Cells_Index.

    ! Define Delta R's.

    call Initialize (DeltaR21_Cells_of_Cells, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    call Initialize (DeltaR1f_Cells_of_Cells, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    call Get_DeltaR21_Cells_of_Cells (DeltaR21_Cells_of_Cells, Diff_Term%Mesh)
    call Get_DeltaR1f_Cells_of_Cells (DeltaR1f_Cells_of_Cells, Diff_Term%Mesh)

    ! Define Harmonic Average Diffusion Coefficient.

    call Initialize (Harmonic_Diffusion_Coef_F_of_C, &
                     NCells_PE(Diff_Term%Mesh), Faces_per_Cell)
    call Get_Harmonic_Diffusion_Coef (Harmonic_Diffusion_Coef_F_of_C, &
                                      Diff_Term)

    ! Define Pseudo-Ortho Area. Options are:
    !   1. Scalar area.
    !   2. Vector area dotted with unit vector joining the cell centers.
    !   3. Middle area of face-center box.

    call Initialize (Area_Faces_of_Cells, NDimensions, &
                     NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    call Initialize (Pseudo_Ortho_Area_F_of_C, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    call Get_Area_Faces_of_Cells (Area_Faces_of_Cells, Diff_Term%Mesh)
    if (.true.) then 
      ! Option 1, Scalar Area.
      Pseudo_Ortho_Area_F_of_C = ABS(SUM(Area_Faces_of_Cells(:,:,:),1))
    end if

    ! Get cell volumes.

    call Initialize (Volume_Cells, NCells_PE(Diff_Term%Mesh))
    call Get_Volume_Cells (Volume_Cells, Diff_Term%Mesh)

    ! Get Phi for Cells and Cells of Cells.

    call Initialize (Phi_Cells, NCells_PE(Diff_Term%Mesh))
    call Initialize (Phi_Cells_DV, Diff_Term%Structure, 1)
    call Initialize (Phi_Cells_of_Cells_CA, &
                     Diff_Term%Mesh%Cells_of_Cells_Index, 1)
    call Initialize (Phi_Cells_of_Cells, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    Phi_Cells = Phi_MV
    Phi_Cells_DV = Phi_Cells
    Phi_Cells_of_Cells_CA = Phi_Cells_DV
    Phi_Cells_of_Cells = Phi_Cells_of_Cells_CA
    call Finalize (Phi_Cells_DV)
    call Finalize (Phi_Cells_of_Cells_CA)

    ! Define Beta values for Boundary Conditions.

    call Initialize (Beta, 3, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    where (Diff_Term%Boundary_Condition(:,:) == Dirichlet_BC)
      Beta(1,:,:) = one
      Beta(2,:,:) = zero
      Beta(3,:,:) = one
    elsewhere (Diff_Term%Boundary_Condition(:,:) == Homogeneous_BC)
      Beta(1,:,:) = one
      Beta(2,:,:) = zero
      Beta(3,:,:) = zero
    elsewhere (Diff_Term%Boundary_Condition(:,:) == Neumann_BC)
      Beta(1,:,:) = zero
      Beta(2,:,:) = -one
      Beta(3,:,:) = one
    elsewhere (Diff_Term%Boundary_Condition(:,:) == Reflective_BC)
      Beta(1,:,:) = zero
      Beta(2,:,:) = -one
      Beta(3,:,:) = zero
    elsewhere (Diff_Term%Boundary_Condition(:,:) == Source_BC)
      Beta(1,:,:) = Diff_Term%Extrapolation
      Beta(2,:,:) = one
      Beta(3,:,:) = Diff_Term%Extrapolation
    elsewhere (Diff_Term%Boundary_Condition(:,:) == Vacuum_BC)
      Beta(1,:,:) = Diff_Term%Extrapolation
      Beta(2,:,:) = one
      Beta(3,:,:) = zero
    end where

    ! Calculate -Div D Grad Phi, by integrating over the cell and  
    ! dividing by the volume:
    !
    !                             ---            
    !                         1   \              
    !   - Div D Grad Phi  =  ---  /    F . A     
    !                   c     V   ---   f   f  
    !                          c faces        
    !

    Div_D_Grad_Phi_Cells = zero
    do face = 1, Faces_per_Cell
      other_cell = face  ! Another name for clearer semantics.

      where (Diff_Term%Coefficient(:) == zero)

        ! No op. Flux is zero.

      elsewhere (Diff_Term%Boundary_Condition(:,face) == &
                 Internal_or_Periodic_BC)

        ! For the internal faces and periodic boundary conditions:
        !
        !                  h 
        !   F . A      A  D     ( Phi  - Phi  )
        !    f   f      f  12   (    1      2 )
        !   ------  =  ------   ---------------
        !      V         V       | delta r   |
        !       c         c      |        12 |

        Div_D_Grad_Phi_Cells(:) = Div_D_Grad_Phi_Cells(:) + &
          Pseudo_Ortho_Area_F_of_C(:,face) / Volume_Cells(:) &
          * Harmonic_Diffusion_Coef_F_of_C(:,face) &
          * (Phi_Cells(:) - Phi_Cells_of_Cells(:,other_cell)) &
          / DeltaR21_Cells_of_Cells(:,face)

      elsewhere

        ! For the generalized Robin Boundary Conditions:
        !                                
        !   F . A        A  ( Beta  Phi  - Beta  Phi   )
        !    f   f        f (     1    1       3    BC )
        !   ------  =  ----------------------------------------
        !      V        V  ( Beta  + Beta  | delta r   | / D  )
        !       c        c (     2       1 |        1f |    1 )

        Div_D_Grad_Phi_Cells(:) = Div_D_Grad_Phi_Cells(:) + &
          Pseudo_Ortho_Area_F_of_C(:,face) / Volume_Cells(:) &
          * (Beta(1,:,face) * Phi_Cells(:) &
          - Beta(3,:,face) * Diff_Term%Phi_BC(:,face) ) &
          / (Beta(2,:,face) &
          + Beta(1,:,face) * DeltaR1f_Cells_of_Cells(:,face) &
          / Diff_Term%Coefficient(:) )

      end where
    end do

    ! Finalize working structures.

    call Finalize (Area_Faces_of_Cells)
    call Finalize (Pseudo_Ortho_Area_F_of_C)
    call Finalize (Volume_Cells)
    call Finalize (DeltaR21_Cells_of_Cells)
    call Finalize (DeltaR1f_Cells_of_Cells)
    call Finalize (Harmonic_Diffusion_Coef_F_of_C)
    call Finalize (Beta)
    call Finalize (Phi_Cells)
    call Finalize (Phi_Cells_of_Cells)

    ! Verify guarantees.

    VERIFY(Valid_State(Div_D_Grad_Phi_Cells),5) ! Div_D_Grad_Phi_Cells is valid.
    VERIFY(Valid_State(Phi_MV),5)               ! Phi_MV is valid.
    VERIFY(Valid_State(Diff_Term),5)            ! Diff_Term is valid.

    return
  end subroutine Evaluate_Ortho_Diffusion



Michael L. Hall