H.2.8 Get_Harmonic_Diffusion_Coef_Ortho_Diffusion Procedure

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

  subroutine Get_Harmonic_Diffusion_Coef (Harmonic_Diffusion_Coef_F_of_C, &
                                          Diff_Term)

    ! Input variable.

    type(Ortho_Diffusion_type), intent(inout) :: Diff_Term ! Diffusion term.

    ! Output variable.

    ! Harmonic diffusion coefficient at each face of each cell.
    type(real,2) :: Harmonic_Diffusion_Coef_F_of_C

    ! Internal variables.

    type(integer) :: face, other_cell      ! Loop parameters.
    ! Diffusion coefficient variables for cells and cells of cells.
    type(Collected_Array_type) :: Coefficient_Cells_of_Cells_CA
    type(Distributed_Vector_type) :: Coefficient_Cells_DV
    type(real,2) :: Coefficient_Cells_of_Cells
    ! Absolute distances between the cell centers and the faces.
    type(real,2) :: DeltaR21_Cells_of_Cells, DeltaR1f_Cells_of_Cells, &
                    DeltaR2f_Cells_of_Cells
    ! A toggle for different harmonic averages.
    type(character,name_length) :: harmonic_diffusion_option
    type(integer) :: Faces_per_Cell  ! Number of faces per cell.
    type(integer) :: NDimensions     ! Number of dimensions.

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

    ! Verify requirements.

    VERIFY(Valid_State(Diff_Term),5)            ! Diff_Term is valid.
    ! Harmonic_Diffusion_Coef_F_of_C is valid.
    VERIFY(Valid_State(Harmonic_Diffusion_Coef_F_of_C),5) 
    ! Harmonic_Diffusion_Coef_F_of_C is dimensioned correctly.
    VERIFY(SIZE(Harmonic_Diffusion_Coef_F_of_C,1) == dnl
           NCells_PE(Diff_Term%Mesh),5)  
    VERIFY(SIZE(Harmonic_Diffusion_Coef_F_of_C,2) == dnl
           Get_Faces_per_Cell(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.
    !   2. Other pseudo-ortho methods.

    ! 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 Initialize (DeltaR2f_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)
    call Get_DeltaR2f_Cells_of_Cells (DeltaR2f_Cells_of_Cells, Diff_Term%Mesh)

    ! Get the Diffusion Coefficient for the neighboring cells.

    call Initialize (Coefficient_Cells_DV, Diff_Term%Structure, 1)
    call Initialize (Coefficient_Cells_of_Cells_CA, &
                     Diff_Term%Mesh%Cells_of_Cells_Index, 1)
    call Initialize (Coefficient_Cells_of_Cells, NCells_PE(Diff_Term%Mesh), &
                     Faces_per_Cell)
    Coefficient_Cells_DV = Diff_Term%Coefficient
    Coefficient_Cells_of_Cells_CA = Coefficient_Cells_DV
    Coefficient_Cells_of_Cells = Coefficient_Cells_of_Cells_CA
    call Finalize (Coefficient_Cells_of_Cells_CA)
    call Finalize (Coefficient_Cells_DV)

    ! Define Harmonic Average Diffusion Coefficient. Four options:
    !
    !   alpha - DeltaR21 (DeltaR1f/D1 + DeltaR2F/D2)^-1
    !   beta  - (DeltaR1f + DeltaR2F) (DeltaR1f/D1 + DeltaR2F/D2)^-1
    !   gamma - (V1 + V2) (V2/D1 + V1/D2)^-1
    !   delta - (V1 + V2) (V1/D1 + V2/D2)^-1
    !
    ! Options alpha, beta, and delta reduce to the orthogonal method on 
    ! orthogonal meshes. Options beta, gamma, and delta reduce to the 
    ! constant D value if D is constant, regardless of the geometry. Option 
    ! gamma reproduces a bug in Telluride, and is here for comparison 
    ! purposes only (don't use it).

    harmonic_diffusion_option = 'beta'   ! Hardwired to beta for now.
    do face = 1, Faces_per_Cell
      other_cell = face  ! Another name for clearer semantics.

      if (harmonic_diffusion_option == 'alpha') then
        where (Diff_Term%Coefficient(:) == zero .OR. &
               Coefficient_Cells_of_Cells(:,other_cell) == zero) 
          Harmonic_Diffusion_Coef_F_of_C(:,face) = zero
        elsewhere
          Harmonic_Diffusion_Coef_F_of_C(:,face) = &
            DeltaR21_Cells_of_Cells(:,face) * &
            (DeltaR1f_Cells_of_Cells(:,face) / &
             Diff_Term%Coefficient(:) + &
             DeltaR2f_Cells_of_Cells(:,face) / &
             Coefficient_Cells_of_Cells(:,other_cell))**(-1)
        end where
      else if (harmonic_diffusion_option == 'beta') then
        where (Diff_Term%Coefficient(:) == zero .OR. &
               Coefficient_Cells_of_Cells(:,other_cell) == zero) 
          Harmonic_Diffusion_Coef_F_of_C(:,face) = zero
        elsewhere
          Harmonic_Diffusion_Coef_F_of_C(:,face) = &
            (DeltaR1f_Cells_of_Cells(:,face) + &
             DeltaR2f_Cells_of_Cells(:,face)) * &
            (DeltaR1f_Cells_of_Cells(:,face) / &
             Diff_Term%Coefficient(:) + &
             DeltaR2f_Cells_of_Cells(:,face) / &
             Coefficient_Cells_of_Cells(:,other_cell))**(-1)
        end where
      else
        VERIFY(.false.,1) ! Not implemented yet.
      end if
    end do

    ! Finalize working structures.

    call Finalize (DeltaR21_Cells_of_Cells)
    call Finalize (DeltaR1f_Cells_of_Cells)
    call Finalize (DeltaR2f_Cells_of_Cells)
    call Finalize (Coefficient_Cells_of_Cells)

    ! Verify guarantees.

    ! Harmonic_Diffusion_Coef_F_of_C is valid.
    VERIFY(Valid_State(Harmonic_Diffusion_Coef_F_of_C),5) 
    VERIFY(Valid_State(Diff_Term),5)      ! Diff_Term is valid.

    return
  end subroutine Get_Harmonic_Diffusion_Coef