H.2.11 Ortho_Diffusion Class Unit Test Program

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) :: 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 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))

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

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
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 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))

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

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
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 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))

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

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
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 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))

! 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

do c = 1, NCells_PE(Mesh)
X = Coordinates_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
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 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))

! 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

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
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 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))

! 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

do c = 1, NCells_PE(Mesh)
X = Coordinates_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): 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.

! ---------------------------------------------------------------------------
! 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.

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))

! 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

do c = 1, NCells_PE(Mesh)
X = Coordinates_Cells(1,c)
if (X < 0.5) then
(-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): 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
(-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): 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
```

Michael L. Hall