Next: H.2.8 Get_Harmonic_Diffusion_Coef_Ortho_Diffusion Procedure Up: H.2 Ortho_Diffusion Class Code Previous: H.2.6 Evaluate_Gradient_Cells_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