Next: H.2.7 Get_Harmonic_Diffusion_Coef_Ortho_Diffusion Procedure Up: H.2 Ortho_Diffusion Class Code Previous: H.2.5 Add_to_Matrix_Equation_Ortho_Diffusion Procedure
The main documentation of the Evaluate_Gradient_Cells_Ortho_Diffusion Procedure contains additional explanation of this code listing.
subroutine Evaluate_Grad_Cells_Ortho_Diff (Diff_Term, Phi_MV, Grad_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.
! Vector Gradients at each cell center.
type(real,2) :: Grad_Cells
! Internal variables.
type(integer) :: dim, other_cell ! Loop parameters.
! Vector Gradients dotted with the normals at each face of each cell.
type(real,2) :: Grad_dot_N_Faces_of_Cells
type(real,3) :: Beta
type(real,2) :: Harmonic_Diffusion_Coef_F_of_C
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.
! 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(Grad_Cells),5) ! Grad_Cells is valid.
! Grad_Cells is dimensioned correctly.
VERIFY(SIZE(Grad_Cells,1) == Get_NDimensions(Diff_Term%Mesh),5)
VERIFY(SIZE(Grad_Cells,2) == 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)
! 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(:,:) == 1)
! Dirichlet BC.
Beta(1,:,:) = one
Beta(2,:,:) = zero
Beta(3,:,:) = one
elsewhere (Diff_Term%Boundary_Condition(:,:) == -1)
! Homogeneous BC.
Beta(1,:,:) = one
Beta(2,:,:) = zero
Beta(3,:,:) = zero
elsewhere (Diff_Term%Boundary_Condition(:,:) == 2)
! Neumann BC.
Beta(1,:,:) = zero
Beta(2,:,:) = -one
Beta(3,:,:) = one
elsewhere (Diff_Term%Boundary_Condition(:,:) == -2)
! Reflective BC.
Beta(1,:,:) = zero
Beta(2,:,:) = -one
Beta(3,:,:) = zero
elsewhere (Diff_Term%Boundary_Condition(:,:) == 3)
! Source BC.
Beta(1,:,:) = Diff_Term%Extrapolation
Beta(2,:,:) = one
Beta(3,:,:) = Diff_Term%Extrapolation
elsewhere (Diff_Term%Boundary_Condition(:,:) == -3)
! Vacuum BC.
Beta(1,:,:) = Diff_Term%Extrapolation
Beta(2,:,:) = one
Beta(3,:,:) = zero
end where
! Calculate Grad Phi . N.
!
! 1
! Grad Phi . N (outward) = - --- F . N
! f f D f f
! 1
call Initialize (Grad_dot_N_Faces_of_Cells, NCells_PE(Diff_Term%Mesh), &
Faces_per_Cell)
Grad_dot_N_Faces_of_Cells = zero
do other_cell = 1, Faces_per_Cell
where (Diff_Term%Coefficient(:) == zero)
Grad_dot_N_Faces_of_Cells(:,other_cell) = zero
elsewhere (Diff_Term%Boundary_Condition(:,other_cell) == 0)
! For the internal faces and periodic boundary conditions:
!
! h
! D ( Phi - Phi )
! 12 ( 1 2 )
! Grad Phi . N (outward) = - ----- ---------------
! f f D | delta r |
! 1 | 12 |
Grad_dot_N_Faces_of_Cells(:,other_cell) = &
- Harmonic_Diffusion_Coef_F_of_C(:,other_cell) &
/ Diff_Term%Coefficient(:) &
* (Phi_Cells(:) - Phi_Cells_of_Cells(:,other_cell)) &
/ DeltaR21_Cells_of_Cells(:,other_cell)
elsewhere
! For the generalized Robin Boundary Conditions:
!
! ( Beta Phi - Beta Phi )
! ( 1 1 3 BC )
! Grad Phi . N (outward) = - ----------------------------------------
! f f D ( Beta + Beta | delta r | / D )
! 1 ( 2 1 | 1f | 1 )
Grad_dot_N_Faces_of_Cells(:,other_cell) = &
- (Beta(1,:,other_cell) * Phi_Cells(:) &
- Beta(3,:,other_cell) * Diff_Term%Phi_BC(:,other_cell) ) &
/ Diff_Term%Coefficient(:) &
/ (Beta(2,:,other_cell) &
+ Beta(1,:,other_cell) * DeltaR1f_Cells_of_Cells(:,other_cell) &
/ Diff_Term%Coefficient(:) )
end where
end do
! Multiply by N (outward unit normal, simply a plus or minus sign) and
! average face values to get cell center value for each direction.
do dim = 1, NDimensions
Grad_Cells(dim,:) = (- Grad_dot_N_Faces_of_Cells(:,2*dim-1) &
+ Grad_dot_N_Faces_of_Cells(:,2*dim )) &
/ two
end do
! Finalize working structures.
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)
call Finalize (Grad_dot_N_Faces_of_Cells)
! Verify guarantees.
VERIFY(Valid_State(Grad_Cells),5) ! Grad_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_Grad_Cells_Ortho_Diff