Next: H.2.7 Evaluate_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, face, 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(:,:) == 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 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 face = 1, Faces_per_Cell
other_cell = face ! Another name for clearer semantics.
where (Diff_Term%Coefficient(:) == zero)
Grad_dot_N_Faces_of_Cells(:,face) = zero
elsewhere (Diff_Term%Boundary_Condition(:,face) == &
Internal_or_Periodic_BC)
! 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(:,face) = &
- Harmonic_Diffusion_Coef_F_of_C(:,face) &
/ Diff_Term%Coefficient(:) &
* (Phi_Cells(:) - Phi_Cells_of_Cells(:,other_cell)) &
/ DeltaR21_Cells_of_Cells(:,face)
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(:,face) = &
- (Beta(1,:,face) * Phi_Cells(:) &
- Beta(3,:,face) * Diff_Term%Phi_BC(:,face) ) &
/ Diff_Term%Coefficient(:) &
/ (Beta(2,:,face) &
+ Beta(1,:,face) * DeltaR1f_Cells_of_Cells(:,face) &
/ 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
!write(6,*) 'Phi_Cells(5) = ', Phi_Cells(5)
!write(6,*) 'Phi_Cells_of_Cells(5,1) = ', Phi_Cells_of_Cells(5,1)
!write(6,*) 'Phi_Cells_of_Cells(5,2) = ', Phi_Cells_of_Cells(5,2)
!write(6,*) 'Grad_dot_N_Faces_of_Cells(5,1) = ', Grad_dot_N_Faces_of_Cells(5,1)
!write(6,*) 'Grad_dot_N_Faces_of_Cells(5,2) = ', Grad_dot_N_Faces_of_Cells(5,2)
!write(6,*) 'Grad_Cells(1,5) = ', Grad_Cells(1,5)
!write(6,*) 'DeltaR21_Cells_of_Cells(5,1) = ', DeltaR21_Cells_of_Cells(5,1)
!write(6,*) 'DeltaR21_Cells_of_Cells(5,2) = ', DeltaR21_Cells_of_Cells(5,2)
!write(6,*) 'Diff_Term%Coefficient(5) = ', Diff_Term%Coefficient(5)
!write(6,*) 'Harmonic_Diffusion_Coef_F_of_C(5,2) = ', Harmonic_Diffusion_Coef_F_of_C(5,2)
! 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