Next: I. Mesh Module Code Up: H.2 Ortho_Diffusion Class Code Previous: H.2.10 Output_Ortho_Diffusion Procedure
This lightly commented program performs a unit test on the Ortho_Diffusion Class.
program Unit_Test
use Caesar_Linear_Algebra_Module
use Caesar_Multi_Mesh_Class
use Caesar_Ortho_Diffusion_Class
use Caesar_Monomial_Class
use Caesar_Numbers_Module, only: zero, half, one, two, three, four, &
six, eight, twelve, sixteen
implicit none
type(Communication_type) :: Comm
type(Status_type) :: status
type(integer) :: c, NRows
! Mesh variables.
type(real,2) :: Coordinates_Cells ! Coordinates of the cell centers.
type(real,1) :: Lengths ! Physical extent of the
! domain in each direction.
type(Multi_Mesh_type) :: Mesh ! Mesh object.
type(integer) :: NDimensions ! Number of dimensions.
type(integer) :: NCells_X ! Number of cells in X.
type(integer) :: NCells_Y ! Number of cells in Y.
type(integer) :: NCells_Z ! Number of cells in Z.
type(integer) :: Faces_per_Cell ! Number of local faces per cell.
type(integer,2) :: Flag_Faces_of_Cells ! Mesh boundary flags.
type(real,3) :: Coordinates_Faces_of_Cells ! Coordinates of face centers.
type(real) :: X ! Local X coordinate.
! Physics variables.
type(real,1) :: Diffusion_Coefficient ! Diffusion Coefficient.
type(Mathematic_Vector_type) :: Phi_MV ! Diffused Quantity (MV).
type(real,1) :: Phi ! Diffused Quantity (BNV).
type(real,2) :: Gradient_Cells ! Gradient vectors at the cells.
type(real) :: Phi_0, Phi_1 ! Dirichlet BC constants.
type(real,1) :: Sigma_a ! Absorption cross-section.
type(real,1) :: Source ! Source term.
type(integer,2) :: BC_Faces_of_Cells ! Boundary condition flags.
type(real,2) :: Phi_BC_Faces_of_Cells ! Boundary condition constants.
type(real) :: Extrapolation ! Boundary condition factor.
type(real) :: D0, D1, D2, Q0
! Equation variables.
type(Ortho_Diffusion_type) :: Ortho_Diff_Term ! Diffusion term.
type(Monomial_type) :: Linear_Term ! Absorption term.
type(Monomial_type) :: Constant_Term ! Source term.
type(logical) :: Derivative ! Time derivative toggle.
type(real) :: Delta_t ! Time step size.
type(real) :: Exponent ! Exponent used in setting Monomial terms.
! Matrix variables.
type(ELL_Matrix_type) :: Matrix_ELLM ! Diffusion equation matrix.
type(integer) :: MaxNonzeros ! Maximum number of nonzero entries per row.
type(real,1) :: X_BNV ! Solution (BNV).
type(Mathematic_Vector_type) :: X_MV ! Solution (MV).
type(Mathematic_Vector_type) :: RHS_MV ! Right-hand side of diffusion eqn.
type(Mathematic_Vector_type) :: Error_MV ! Error Vector (MV).
type(Mathematic_Vector_type) :: Exact_MV ! Exact Vector (MV).
type(real,1) :: Exact ! Exact Vector (BNV).
type(real) :: Relative_L2_Error_Norm ! Relative L2 Error Norm.
! Base structures for the matrix.
type(Base_Structure_type) :: Row_Structure, Column_Structure
type(Solver_type) :: Solver
type(integer) :: Variable ! Variable number.
type(integer) :: Equation, eq ! Equation number.
type(integer) :: NEquations ! Number of equations.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Initializations.
call Initialize (Comm)
call Output (Comm)
call Initialize (status)
! Create a mesh.
NCells_X = 10
NCells_Y = 2
NCells_Z = 12
!NCells_X = 2
!NCells_Y = 1
!NCells_Z = 1
NDimensions = 3
call Initialize (Lengths, NDimensions)
Lengths = (/ 1.d0, 1.d0, 1.d0 /)
call Initialize (Mesh, NDimensions, Lengths, NCells_X, NCells_Y, NCells_Z, &
'Uniform Mesh', status)
call Initialize (Coordinates_Cells, NDimensions, NCells_PE(Mesh), &
status)
call Get_Coordinates_Cells (Coordinates_Cells, Mesh)
Delta_t = zero
Derivative = .false.
Variable = 3
Equation = 2
NEquations = 3
! Set Diffusion Coefficient for every face of each cell.
! Faces_of_Cells way -- for later:
!
!Faces_per_Cell = two*NDimensions
!call Initialize (Diffusion_Coefficient, NCells_PE(Mesh), Faces_per_Cell, &
! status(4))
!do lf = 1, Faces_per_Cell
! where (Coordinates_Cells(1,:) < 0.5)
! Diffusion_Coefficient(:,lf) = Diff_Coeff(1)
! elsewhere
! Diffusion_Coefficient(:,lf) = Diff_Coeff(2)
! end where
!end do
! Set Boundary Conditions.
! (Should do by face, but doing by faces_of_cells is easier for now.)
! call Initialize (Boundary_Conditions, NFaces_PE(Mesh), status)
Faces_per_Cell = Get_Faces_per_Cell (Mesh)
call Initialize (BC_Faces_of_Cells, NCells_PE(Mesh), Faces_per_Cell, &
status)
call Initialize (Flag_Faces_of_Cells, NCells_PE(Mesh), Faces_per_Cell, &
status)
call Initialize (Phi_BC_Faces_of_Cells, NCells_PE(Mesh), Faces_per_Cell, &
status)
call Get_Flag_Faces_of_Cells (Flag_Faces_of_Cells, Mesh)
! Extrapolation Factor is set to the default value -- cannot leave it out of
! the initialize call because it appears before the status variable.
Extrapolation = half
! Get Cell Face Coordinates.
call Initialize (Coordinates_Faces_of_Cells, NDimensions, &
NCells_PE(Mesh), Faces_per_Cell)
call Get_Coordinates_Faces_of_Cells (Coordinates_Faces_of_Cells, Mesh)
! Physics initializations.
call Initialize (Sigma_a, NCells_PE(Mesh), status)
call Initialize (Source, NCells_PE(Mesh), status)
call Initialize (Diffusion_Coefficient, NCells_PE(Mesh), status)
call Initialize (Phi, NCells_PE(Mesh), status)
call Initialize (Gradient_Cells, NDimensions, NCells_PE(Mesh), status)
call Initialize (Exact, NCells_PE(Mesh), status)
call Initialize (Phi_MV, Cell_Structure(Mesh), 'Phi_Cells', status)
call Initialize (Error_MV, Cell_Structure(Mesh), 'Error = Solution - Exact', &
status)
call Initialize (Exact_MV, Cell_Structure(Mesh), 'Exact', status)
! Set up Matrix Equation.
MaxNonzeros = 2*NDimensions + 1
call Generate_Multiple (Row_Structure, NEquations, Cell_Structure(Mesh), &
'Rows')
call Generate_Multiple (Column_Structure, NEquations, Cell_Structure(Mesh), &
'Columns')
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (X_BNV, Length_PE(Column_Structure), status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.1: 1-D Linear Solution with Dirichlet BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.1: 1-D Linear Solution with Dirichlet BCs'
end if
! Problem parameters.
Phi_0 = 47.d0
D0 = one/3.d0
Q0 = zero
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
Source = -Q0
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Dirichlet_BC
Phi_BC_Faces_of_Cells = Phi_0
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Homogeneous_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
Exact(c) = Phi_0*(one-Coordinates_Cells(1,c))
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
if (ABS(Gradient_Cells(1,c) - (-Phi_0)) > 9.d-8 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 2.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 5.d-8) then
write (6,*) 'Error (Problem X.1): Grad = ', Gradient_Cells(:,c)
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.2: 1-D Linear Solution with Source/Vacuum BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.2: 1-D Linear Solution with Source/Vacuum BCs'
end if
! Problem parameters.
Phi_0 = 47.d0
D0 = one/3.d0
Q0 = zero
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
Source = -Q0
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Source_BC
Phi_BC_Faces_of_Cells = Phi_0
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
Exact(c) = Phi_0 * (one + two*D0 - Coordinates_Cells(1,c)) / &
(one + four*D0)
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
if (ABS(Gradient_Cells(1,c) - (-Phi_0/(one + four*D0))) > 5.d-9 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 4.d-9 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 4.d-9) then
write (6,*) 'Error (Problem X.2): Grad = ', Gradient_Cells(:,c)
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.3: Multi-D Linear Solution with Dirichlet BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.3: Multi-D Linear Solution with Dirichlet BCs'
end if
! Problem parameters.
Phi_0 = 47.d0
D0 = one/3.d0
Q0 = zero
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
Source = -Q0
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Dirichlet_BC
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Dirichlet_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Dirichlet_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Dirichlet_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Dirichlet_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Dirichlet_BC
end where
Phi_BC_Faces_of_Cells(:,:) = Phi_0 * SUM(Coordinates_Faces_of_Cells(:,:,:),1)
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
Exact(c) = Phi_0 * SUM(Coordinates_Cells(:,c),1)
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
if (ABS(Gradient_Cells(1,c) - (Phi_0) ) > 7.d-7 .or. &
ABS(Gradient_Cells(2,c) - (Phi_0) ) > 7.d-7 .or. &
ABS(Gradient_Cells(3,c) - (Phi_0) ) > 7.d-7) then
write (6,*) 'Error (Problem X.3): Grad = ', Gradient_Cells(:,c)
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.4: 1-D Quadratic Solution with Dirichlet BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.4: 1-D Quadratic Solution with Dirichlet BCs'
end if
! Problem parameters.
D0 = one/3.d0
Q0 = 23.d0
Phi_0 = 6.d0
Phi_1 = 19.d0
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
Source = -Q0
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
Phi_BC_Faces_of_Cells = zero
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Dirichlet_BC
Phi_BC_Faces_of_Cells = Phi_0
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Dirichlet_BC
Phi_BC_Faces_of_Cells = Phi_1
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
! Note: I tested problem X.4 and it converged by 2nd Order when
! refined in the X-direction.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
Exact(c) = Q0 / (two*D0) * (X - X**2) + (Phi_1 - Phi_0) * X + Phi_0
!write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (ABS(Gradient_Cells(1,c) - &
(Q0 / (two*D0) * (1 - 2*X) + (Phi_1 - Phi_0)) ) > 1.d-8 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 1.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 1.d-8) then
write (6,*) 'Error (Problem X.4): Grad = ', Gradient_Cells(:,c)
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.5: 1-D Quadratic Solution with Vacuum BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.5: 1-D Quadratic Solution with Vacuum BCs'
end if
! Problem parameters.
D0 = one/3.d0
Q0 = 23.d0
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
Source = -Q0
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
Phi_BC_Faces_of_Cells = zero
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
! Note: I tested problem X.5 and it converged by 2nd Order when
! refined in the X-direction.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
Exact(c) = Q0 / (two*D0) * ( two*D0 + X - X**2)
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (ABS(Gradient_Cells(1,c) - (Q0 / (two*D0) * (1 - 2*X)) ) > 1.d-8 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 1.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 1.d-8) then
write (6,*) 'Error (Problem X.5): Grad = ', Gradient_Cells(:,c)
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.6: 1-D Quartic Solution with Vacuum BCs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.6: 1-D Quartic Solution with Vacuum BCs'
end if
! Problem parameters.
D0 = one/3.d0
Q0 = 23.d0
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
do c = 1, NCells_PE(Mesh)
Source(c) = -Q0 * Coordinates_Cells(1,c)**2
end do
! Set Diffusion Coefficient for each cell.
Diffusion_Coefficient = D0
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
Phi_BC_Faces_of_Cells = zero
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
call Set_0_Diagonal_to_1 (Matrix_ELLM)
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
! Note: I tested problem X.6 and it converged by 2nd Order when
! refined in the X-direction.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
Exact(c) = (Q0 / six) * ((one + eight * D0)/(one + four*D0) * &
(one + X / (two*D0)) - X**4 / (two*D0))
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (ABS(Gradient_Cells(1,c) - &
((Q0 / six) * ((one + eight * D0)/(one + four*D0) * &
(one / (two*D0)) - four*X**3 / (two*D0))) ) > 0.12d0 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 1.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 1.d-8) then
write (6,*) 'Error (Problem X.6): Grad = ', Gradient_Cells(:,c)
write (6,*) 'Error (Problem X.6): Diff = ', ABS(Gradient_Cells(1,c) - &
((Q0 / six) * ((one + eight * D0)/(one + four*D0) * &
(one / (two*D0)) - four*X**3 / (two*D0))))
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Problem X.7: 1-D Quartic Solution with Two Materials
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (this_is_IO_PE) then
write (6,100) 'Problem X.7: 1-D Quartic Solution with Two Materials'
end if
! Problem parameters.
D1 = one/3.d0
D2 = one/3.d6
Q0 = 0.00023d0
! Set Absorption (Removal) cross-section for each cell.
Sigma_a = zero
! Set Source term for each cell.
do c = 1, NCells_PE(Mesh)
Source(c) = -Q0 * Coordinates_Cells(1,c)**2
end do
! Set Diffusion Coefficient for each cell.
where (Coordinates_Cells(1,:) < 0.5)
Diffusion_Coefficient(:) = D1
elsewhere
Diffusion_Coefficient(:) = D2
end where
! Set Boundary Conditions.
where (Flag_Faces_of_Cells == 1) ! Left (-x) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 2) ! Right (+x) Face.
BC_Faces_of_Cells = Vacuum_BC
elsewhere (Flag_Faces_of_Cells == 3) ! Front (-y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 4) ! Back (+y) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 5) ! Bottom (-z) Face.
BC_Faces_of_Cells = Reflective_BC
elsewhere (Flag_Faces_of_Cells == 6) ! Top (+z) Face.
BC_Faces_of_Cells = Reflective_BC
end where
Phi_BC_Faces_of_Cells = zero
! Set Phi and X to an initial guess to avoid division by zero error in
! LAMG in the stopping test.
Phi = one
Phi_MV = Phi
X_BNV = one
X_MV = X_BNV
! Define Steady State Diffusion Equation.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, BC_Faces_of_Cells, &
Phi_BC_Faces_of_Cells, Phi_MV, 'Cells', Mesh, &
'Diffusion Term', Extrapolation, Equation, NEquations, &
status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, Variable, &
Equation, NEquations, status)
! Output Terms.
NRows = NCells_PE(Mesh)
!call Output ( Constant_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output ( Linear_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
!call Output (Ortho_Diff_Term, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of Equation Terms.
VERIFY(Valid_State( Constant_Term),0)
VERIFY(Valid_State( Linear_Term),0)
VERIFY(Valid_State(Ortho_Diff_Term),0)
! Populate Matrix Equation.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
! ---------------------------------------------------------------------------
! This is a difficult problem, and it turned out that setting the diagonals
! to unity on zero rows, as is done for the other problems, resulted in a
! fairly good answer, but the gradients didn't pass the current tests. So,
! instead of taking that approach, we enter the same equation multiple times.
! So we repeat the steps above for the other equations. This works well.
do eq = 1, NEquations
if (eq /= Equation) then
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Define Steady State Diffusion Equation for Equation = eq.
call Initialize (Ortho_Diff_Term, Diffusion_Coefficient, &
BC_Faces_of_Cells, Phi_BC_Faces_of_Cells, Phi_MV, &
'Cells', Mesh, 'Diffusion Term', Extrapolation, eq, &
NEquations, status)
Exponent = one
call Initialize (Linear_Term, Sigma_a, Exponent, Phi_MV, 'Cells', Mesh, &
'Absorption Term', Derivative, Delta_t, Phi_MV, &
Variable, eq, NEquations, status)
Exponent = zero
call Initialize (Constant_Term, Source, Exponent, Phi_MV, 'Cells', Mesh, &
'Source Term', Derivative, Delta_t, Phi_MV, &
Variable, eq, NEquations, status)
! Populate Matrix Equation for Equation = eq.
call Add_to_Matrix_Equation (Ortho_Diff_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Linear_Term, Matrix_ELLM, RHS_MV)
call Add_to_Matrix_Equation (Constant_Term, Matrix_ELLM, RHS_MV)
end if
end do
! End of repeated equation.
! -------------------------
! Solve Matrix Equation.
VERIFY(Valid_State(Matrix_ELLM),0)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Initialize (Solver, 'LAMG')
call Set (Solver, 'Epsilon', 1.d-10)
call Solve (Solver, Matrix_ELLM, X_MV, RHS_MV)
call Finalize (Solver)
! Output.
VERIFY(Valid_State(X_MV),0)
!call Output (X_MV)
!call Output (Matrix_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
VERIFY(Valid_State(RHS_MV),0)
!call Output (RHS_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check answers.
! Note: I tested problem X.7 and it converged by 2nd Order when
! refined in the X-direction.
X_BNV = X_MV
Phi = X_BNV(Equation::NEquations)
Phi_MV = Phi
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (X < 0.5) then
Exact(c) = Q0 / (twelve*D2) * (one + eight * D2 + &
(D2 - D1) / (sixteen*D1) - (D2/D1) * X**4)
else
Exact(c) = Q0 / (twelve*D2) * (one + eight * D2 - X**4)
end if
! write (6,*) c, Coordinates_Cells(1,c), Phi(c), Exact(c), Exact(c) - Phi(c)
end do
Error_MV = Phi - Exact
call Output (Error_MV, 0, 0)
Exact_MV = Exact
Relative_L2_Error_Norm = Norm(Error_MV) / Norm(Exact_MV)
if (this_is_IO_PE) &
write (6,102) ' Relative L2 error norm = ', Relative_L2_Error_Norm
! Plotting output (turn on if desired).
if (.false.) then
call Dump_GMV ('GMVout.gmv', Mesh, Phi_MV, Exact_MV, Error_MV)
end if
! Check gradient.
call Evaluate_Gradient_Cells (Ortho_Diff_Term, Phi_MV, Gradient_Cells)
do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (X < 0.5) then
if (ABS(Gradient_Cells(1,c) - &
(-Q0 / (three*D1) * X**3) ) > 1.1d0 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 1.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 1.d-8) then
write (6,*) 'Error (Problem X.7): Grad = ', Gradient_Cells(:,c)
write (6,*) 'Error (Problem X.7): Real = ', &
(-Q0 / (three*D1) * X**3), X
write (6,*) 'Error (Problem X.7): Diff = ', &
ABS(Gradient_Cells(1,c) - (-Q0 / (three*D1) * X**3) )
end if
else
if (ABS(Gradient_Cells(1,c) - &
(-Q0 / (three*D2) * X**3) ) > 1.1d0 .or. &
ABS(Gradient_Cells(2,c) - (zero) ) > 1.d-8 .or. &
ABS(Gradient_Cells(3,c) - (zero) ) > 1.d-8) then
write (6,*) 'Error (Problem X.7): Grad = ', Gradient_Cells(:,c)
write (6,*) 'Error (Problem X.7): Real = ', &
(-Q0 / (three*D2) * X**3), X
write (6,*) 'Error (Problem X.7): Diff = ', &
ABS(Gradient_Cells(1,c) - (-Q0 / (three*D2) * X**3) )
end if
end if
end do
! Finalize terms.
call Finalize (Ortho_Diff_Term)
call Finalize (Linear_Term)
call Finalize (Constant_Term)
! Reset entire Matrix Equation.
!
! Note that we must finalize X_MV before Matrix_ELLM because the call to
! Solver may use a MatVec to calculate a residual, which would create an
! OV in X_MV based on the Data_Index in Matrix_ELLM. Tricky.
call Finalize (X_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Initialize (Matrix_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
'Coefficients', status)
call Initialize (RHS_MV, Row_Structure, 'Right-Hand Side', &
status)
call Initialize (X_MV, Column_Structure, 'Solution', status)
!~~~~~~~~~~~End of test problems.~~~~~~~~~~~
! Finalize variables.
call Finalize (BC_Faces_of_Cells)
call Finalize (Flag_Faces_of_Cells)
call Finalize (Phi_BC_Faces_of_Cells)
call Finalize (Coordinates_Faces_of_Cells)
call Finalize (MaxNonzeros)
call Finalize (Phi_MV)
call Finalize (Error_MV)
call Finalize (RHS_MV)
call Finalize (Matrix_ELLM)
call Finalize (Phi)
call Finalize (Source)
call Finalize (Sigma_a)
call Finalize (Diffusion_Coefficient)
call Finalize (Mesh)
call Finalize (NDimensions)
call Finalize (Lengths)
call Finalize (NCells_X)
call Finalize (NCells_Y)
call Finalize (NCells_Z)
call Finalize (Comm)
! Format statements.
100 format (/,a,:,a,:,a,:,a)
102 format (a,1pe22.15)
end