Next: H.2.6 Evaluate_Gradient_Cells_Ortho_Diffusion Procedure Up: H.2 Ortho_Diffusion Class Code Previous: H.2.4 Initialized_Ortho_Diffusion Procedure
The main documentation of the Add_to_Matrix_Equation_Ortho_Diffusion Procedure contains additional explanation of this code listing.
subroutine Add_to_Matrix_Equation_Ortho_D (Diff_Term, ELLM, RHS_MV)
! Input variables.
! Diff_Term to be added.
type(Ortho_Diffusion_type), intent(inout) :: Diff_Term
! Input/Output variables.
! ELL_Matrix to be incremented.
type(ELL_Matrix_type), intent(inout) :: ELLM
! RHS mathematic vector.
type(Mathematic_Vector_type), intent(inout) :: RHS_MV
! Internal variables.
! Arrays to manipulate values of the matrix.
type(real,2) :: Matrix_Values
type(real,1) :: RHS_Values
type(integer,1) :: Matrix_Rows
type(integer,2) :: Matrix_Columns
type(integer) :: face, i, other_cell ! Loop parameters.
type(integer) :: shift ! Index shift.
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.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Verify requirements.
VERIFY(Valid_State(Diff_Term),5) ! Diff_Term is valid.
VERIFY(Valid_State(ELLM),5) ! ELLM is valid.
VERIFY(Valid_State(RHS_MV),5) ! RHS_MV is valid.
! 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.
! Create some working vectors.
call Initialize (Matrix_Values, NCells_PE(Diff_Term%Mesh), &
1+Faces_per_Cell)
call Initialize (RHS_Values, NCells_PE(Diff_Term%Mesh))
call Initialize (Matrix_Rows, NCells_PE(Diff_Term%Mesh))
call Initialize (Matrix_Columns, NCells_PE(Diff_Term%Mesh), &
1+Faces_per_Cell)
! 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)
!call Get_VecDeltaR21_Cells_of_Cells (DeltaR21_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)
! Must be dimensioned as follows:
! call Initialize (Diff_Term%Boundary_Condition, &
! NCells_PE(Diff_Term%Mesh), &
! Faces_per_Cell)
! 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
! 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
! Set Matrix_Rows, Matrix_Columns.
!
! To spread out the rows and columns from this equation into a multiple
! equation matrix, translate local numbers via:
! i --> NEquations*(i - 1) + Equation
! and global numbers via:
! i + shift --> NEquations*(i + shift - 1) + Equation
Matrix_Columns(:,2:) = Diff_Term%Mesh%Cells_of_Cells_Index
Matrix_Columns(:,2:) = Diff_Term%NEquations * (Matrix_Columns(:,2:) - 1) &
+ Diff_Term%Equation
shift = First_PE(Diff_Term%Structure) - 1
do i = 1, Length_PE(Diff_Term%Structure)
Matrix_Rows(i) = Diff_Term%NEquations*(i - 1) + Diff_Term%Equation
Matrix_Columns(i,1) = Diff_Term%NEquations*(i + shift - 1) + &
Diff_Term%Equation
end do
! Calculate matrix coefficients for F . A .
! f f
! Matrix structure:
! Matrix_Values(:,1) = Diagonal, the cell center value.
! Matrix_Values(:,2:Faces_per_Cell+1) = The values for the cells on
! the other side of each face, shifted by one from the face number
! to allow for the diagonal entry.
Matrix_Values = zero
do face = 1, Faces_per_Cell
other_cell = face ! Another name for clearer semantics.
where (Diff_Term%Boundary_Condition(:,face) == Internal_or_Periodic_BC)
! For the internal faces and periodic boundary conditions:
!
! ( Phi - Phi )
! h ( 1 2 )
! F . A = A D ---------------
! f f f 12 | delta r |
! | 12 |
!
! The other_cell (Phi_2) coefficient is calculated and stored
! here. The cell (diagonal entry, Phi_1) coefficient is calculated
! below by summing and subtracting all the other_cell coefficients.
Matrix_Values(:,other_cell+1) = &
- Pseudo_Ortho_Area_F_of_C(:,face) &
* Harmonic_Diffusion_Coef_F_of_C(:,face) &
/ DeltaR21_Cells_of_Cells(:,face)
elsewhere
! For the generalized Robin Boundary Conditions:
!
! A ( Beta Phi - Beta Phi )
! f ( 1 1 3 BC )
! F . A = -------------------------------------
! f f ( Beta + Beta | delta r | / D )
! ( 2 1 | 1f | 1 )
!
! The cell (diagonal entry, Phi_1) coefficient and the constant, RHS,
! value are both calculated here. Note the sign flip on the RHS value
! to move it to the right-hand side.
!
! Since a Fick's law extrapolation is used to determine this
! term, if D_1 is set to zero the flow (F.A) is set to zero,
! regardless of the Beta values.
where (Diff_Term%Coefficient(:) /= zero)
Matrix_Values(:,1) = Matrix_Values(:,1) + &
Pseudo_Ortho_Area_F_of_C(:,face) &
* Beta(1,:,face) &
/ (Beta(2,:,face) &
+ Beta(1,:,face) * DeltaR1f_Cells_of_Cells(:,face) &
/ Diff_Term%Coefficient(:) )
RHS_Values(:) = RHS_Values(:) + &
Pseudo_Ortho_Area_F_of_C(:,face) &
* 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 where
end do
! Finish the cell (diagonal entry, Phi_1) coefficient calculation by
! by summing and subtracting all the other_cell coefficients.
Matrix_Values(:,1) = Matrix_Values(:,1) - SUM(Matrix_Values(:,2:),2)
! Add the values to the ELLM matrix and RHS vector.
call Add_Values (ELLM, Matrix_Values, Matrix_Rows, Matrix_Columns, &
Global=.false.)
call Add_Values (RHS_MV, RHS_Values, Matrix_Rows, Global=.false.)
! 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 (Pseudo_Ortho_Area_F_of_C)
call Finalize (Beta)
call Finalize (Area_Faces_of_Cells)
call Finalize (Matrix_Values)
call Finalize (RHS_Values)
call Finalize (Matrix_Rows)
call Finalize (Matrix_Columns)
! Verify guarantees.
VERIFY(Valid_State(ELLM),5) ! ELLM is valid.
VERIFY(Valid_State(RHS_MV),5) ! RHS_MV is valid.
return
end subroutine Add_to_Matrix_Equation_Ortho_D