# 1D Steady Thermal Conduction, BCs of 1st Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The Python Script

## The Problem Description¶

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.1.1, pg.14, “A 1D Steady-State Temperature Distribution, Boundary Conditions of 1st Kind.”

The domain is a 100x10x10 meter rectangular beam extending along the positive x-axis and is made up of 10x1x1 cubic grid cells with dimensions 10x10x10 meters. The domain material is assigned the following properties: thermal conductivity K = 1.0 W/(m-C); specific heat capacity Cp = 0.001 J/(m-C); density rho = 2,800 kg/m^3.

The temperature is initially uniform at T(t=0) = 1.5C. The boundary temperatures are T(x=0) = 1.0C and T(x=L) = 2.0C, where L = 100 m. The simulation is run until the steady-state temperature distribution develops.

The LaPlace equation governs the steady-state temperature distribution,

${{\partial^{2} T} \over {\partial x^{2}}} = 0$

The solution is given by,

$T(x) = ax + b$

When the boundary conditions T(x=0) = 1.0C and T(x=L) = 2.0C are applied, the solution becomes,

$T(x) = T(x=0){x \over L}$

## The PFLOTRAN Input File (GENERAL Mode)¶

The General Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE GENERAL
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 8 1 1
DXYZ
12.5d0
10.0d0
10.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 10.0d0 10.0d0
/
END

REGION right_end
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 10.0d0 10.0d0
/
END

MATERIAL_PROPERTY beam
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

CHARACTERISTIC_CURVES cc1
SATURATION_FUNCTION VAN_GENUCHTEN
LIQUID_RESIDUAL_SATURATION 0.5d-1
M 0.75
ALPHA 1.d-3
/
PERMEABILITY_FUNCTION MUALEM
LIQUID_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.1d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END

STRATA
REGION all
MATERIAL beam
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 2.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 100
NO_PRINT_INITIAL
FORMAT HDF5
/
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.5d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.D0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 2.D0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
REGION left_end
FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
FLOW_CONDITION right_end
REGION right_end
END

END_SUBSURFACE


## The PFLOTRAN Input File (TH Mode)¶

The TH Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE TH
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 8 1 1
DXYZ
12.5d0
10.0d0
10.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 10.0d0 10.0d0
/
END

REGION right_end
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 10.0d0 10.0d0
/
END

MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
SPECIFIC_HEAT 1.d-3 J/kg-C
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

SATURATION_FUNCTION default
SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
RESIDUAL_SATURATION 0.5d-1
LAMBDA 0.75
ALPHA 1.d-3
END

STRATA
REGION all
MATERIAL beam
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 2.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.5d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 2.d0 C
PRESSURE 101325 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
REGION left_end
FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
FLOW_CONDITION right_end
REGION right_end
END

END_SUBSURFACE


## The Python Script¶

def thermal_steady_1D_BC1stkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.1.1, pg.14
# "A 1D Steady-State Temperature Distribution, Boundary Conditions of 1st Kind"
#
# Author: Jenn Frederick
# Date: 06/27/2016
# ******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
test_pass = False
try_count = 0

# initial discretization values
lxyz[0] = 100.    # [m] lx
lxyz[1] = 10.     # [m] ly
lxyz[2] = 10.     # [m] lz
nxyz[0] = 8       # [m] nx
dxyz = lxyz/nxyz  # [m]
T0 = 1.           # [C]
ierr = 0

while (not test_pass) and (try_count < num_tries):
print_discretization(lxyz,nxyz,dxyz)
nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
try_count = try_count + 1

x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)    # [m]
x_pflotran = x_soln                               # [m]
x_soln = np.concatenate(([0.],x_soln),axis=0)     # [m]
x_soln = np.concatenate((x_soln,[Lx]),axis=0)     # [m]
T_soln = np.zeros(nx+2)                           # [C]
T_pflotran = np.zeros(nx)                         # [C]

# create the analytical solution
T_soln = np.array(x_soln/Lx + T0)                 # [C]

# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_str = 'Time:  1.00000E+02 y/Temperature [C]'
ierr = check(T_pflotran)

max_percent_error = calc_relative_error(T_soln[1:nx+1],T_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz[0] = nxyz[0]*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions
'Temperature [C]',"{0:.2f}".format(max_percent_error))

# Plot error analysis
plot_error(error_analysis,dxyz_record,path,try_count)

# Add test result to report card

return;


Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

# 1D Steady Thermal Conduction, BCs of 1st & 2nd Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The Python Script

## The Problem Description¶

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.1.2, pg.15, “A 1D Steady-State Temperature Distribution, Boundary Conditions of 1st and 2nd Kind.”

The domain is a 100x10x10 meter rectangular beam extending along the positive x-axis and is made up of 10x1x1 cubic grid cells with dimensions 10x10x10 meters. The domain is composed of two materials and is assigned the following properties: thermal conductivity K1(x<2L/5) = 100 W/(m-C); thermal conductivity K2(x>2L/5) = 300 W/(m-C); specific heat capacity Cp = 0.001 J/(m-C); density rho = 2,800 kg/m^3.

The temperature is initially uniform at T(t=0) = 1.0C. The boundary temperature at the left end is T(x=0) = 1.0C and a heat flux q(x=L) = -1.5 W/m^2 is applied at the right end, where L = 100 m. The simulation is run until the steady-state temperature distribution develops.

The LaPlace equation governs the steady-state temperature distribution,

${{\partial^{2} T} \over {\partial x^{2}}} = 0$

The solution is given by,

\begin{align}\begin{aligned}T(x<2L/5) = a_1 x + b_1\\T(x>2L/5) = a_2 x + b_2\end{aligned}\end{align}

When the boundary conditions T(x=0) = 1.0C and q(x=L) = -1.5 W/m^2 are applied, the solution becomes,

\begin{align}\begin{aligned}T(x<2L/5) = -{{qx} \over {K1}} + T(x=0)\\T(x>2L/5) = -{{qx} \over {K2}} + T(x=0) + {{2qL} \over {5}}\left({1 \over K2}-{1 \over K1}\right)\end{aligned}\end{align}

## The PFLOTRAN Input File (GENERAL Mode)¶

The General Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE GENERAL
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 10 1 1
DXYZ
10.0d0
10.0d0
10.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left
COORDINATES
0.d0   0.d0   0.d0
40.0d0 10.0d0 10.0d0
/
END

REGION right
COORDINATES
40.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 10.0d0 10.0d0
/
END

REGION right_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 10.0d0 10.0d0
/
END

MATERIAL_PROPERTY beam_left
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

MATERIAL_PROPERTY beam_right
ID 2
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 300 W/m-C
THERMAL_CONDUCTIVITY_WET 300 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

CHARACTERISTIC_CURVES cc1
SATURATION_FUNCTION VAN_GENUCHTEN
LIQUID_RESIDUAL_SATURATION 0.5d-1
M 0.75
ALPHA 1.d-3
/
PERMEABILITY_FUNCTION MUALEM
LIQUID_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.1d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END

STRATA
REGION left
MATERIAL beam_left
END

STRATA
REGION right
MATERIAL beam_right
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 2.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.0d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION left_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 1.d0 C
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION right_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX +1.5d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial_left
REGION left
FLOW_CONDITION initial
END
INITIAL_CONDITION initial_right
REGION right
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_face
REGION left_face
FLOW_CONDITION left_face
END

BOUNDARY_CONDITION right_face
REGION right_face
FLOW_CONDITION right_face
END

END_SUBSURFACE


## The PFLOTRAN Input File (TH Mode)¶

The TH Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE TH
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 10 1 1
DXYZ
10.0d0
10.0d0
10.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left
COORDINATES
0.d0   0.d0   0.d0
40.0d0 10.0d0 10.0d0
/
END

REGION right
COORDINATES
40.0d0 0.0d0 0.0d0
100.0d0 10.0d0 10.0d0
/
END

REGION left_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 10.0d0 10.0d0
/
END

REGION right_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 10.0d0 10.0d0
/
END

MATERIAL_PROPERTY beam_left
ID 1
POROSITY 1.d-5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
SPECIFIC_HEAT 1.d-3
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

MATERIAL_PROPERTY beam_right
ID 2
POROSITY 1.d-5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
SPECIFIC_HEAT 1.d-3
THERMAL_CONDUCTIVITY_DRY 300 W/m-C
THERMAL_CONDUCTIVITY_WET 300 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

SATURATION_FUNCTION default
SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
RESIDUAL_SATURATION 0.5d-1
LAMBDA 0.75
ALPHA 1.d-3
END

STRATA
REGION left
MATERIAL beam_left
END

STRATA
REGION right
MATERIAL beam_right
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 2.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.0d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION left_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION right_face
TYPE
ENERGY_FLUX neumann
PRESSURE dirichlet
/
ENERGY_FLUX +1.5d0 W/m^2
PRESSURE 101325 Pa
END

INITIAL_CONDITION initial_left
REGION left
FLOW_CONDITION initial
END
INITIAL_CONDITION initial_right
REGION right
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_face
REGION left_face
FLOW_CONDITION left_face
END

BOUNDARY_CONDITION right_face
REGION right_face
FLOW_CONDITION right_face
END

END_SUBSURFACE


## The Python Script¶

def thermal_steady_1D_BC1st2ndkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.1.2, pg.15
# "A 1D Steady-State Temperature Distribution, Boundary Conditions of 1st and
# 2nd Kind"
#
# Author: Jenn Frederick
# Date: 06/27/2016
# *****************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
lxyz[0] = 100.    # [m] lx
lxyz[1] = 10.     # [m] ly
lxyz[2] = 10.     # [m] lz
test_pass = False
try_count = 0

# initial discretization values
nxyz[0] = 10      # [m] nx
dxyz = lxyz/nxyz  # [m]
K1 = 100.         # [W/m-C]
K2 = 300.         # [W/m-C]
q = -1.5          # [W/m^2]
T0 = 1.           # [C]
ierr = 0

while (not test_pass) and (try_count < num_tries):
print_discretization(lxyz,nxyz,dxyz)
nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
try_count = try_count + 1

x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)    # [m]
x_pflotran = x_soln                               # [m]
x_soln = np.concatenate(([0.],x_soln),axis=0)     # [m]
T_soln = np.zeros(nx+1)                           # [C]
T_pflotran = np.zeros(nx)                         # [C]

# create the analytical solution
k = -1
for j in x_soln:
k = k + 1
if j <= (2.*Lx/5.):
T_soln[k] = np.array( -((q/K1)*j) + T0 )
else:
T_soln[k] = np.array( -((q/K2)*j) + T0 + (q*(2.*Lx/5.)*((1/K2)-(1/K1))) )

# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+02 y/Temperature [C]'
ierr = check(T_pflotran)

max_percent_error = calc_relative_error(T_soln[1:nx+1],T_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz[0] = nxyz[0]*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions
'Temperature [C]',"{0:.2f}".format(max_percent_error))

# Plot error analysis
plot_error(error_analysis,dxyz_record,path,try_count)

# Add test result to report card

return;


Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

# 2D Steady Thermal Conduction, BCs of 1st Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The Dataset

The Python Script

## The Problem Description¶

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.1.3, pg.16, “A 2D Steady-State Temperature Distribution, Boundary Conditions of 1st Kind.”

The domain is a 1x1x1 meter slab extending along the positive x-axis and y-axis and is made up of 10x10x1 hexahedral grid cells with dimensions 0.1x0.1x1 meters. The domain material is assigned the following properties: thermal conductivity K = 1.0 W/(m-C); specific heat capacity Cp = 0.001 J/(m-C); density rho = 2,800 kg/m^3.

The temperature is initially uniform at T(t=0) = 1.0 C. The boundary temperatures are:

\begin{align}\begin{aligned}T(0,y) = 0 \hspace{0.25in} x=0 \hspace{0.15in} face\\T(L,y) = T0 {y \over L} \hspace{0.25in} x=L \hspace{0.15in} face\\T(x,0) = 0 \hspace{0.25in} y=0 \hspace{0.15in} face\\T(x,L) = T0 {x \over L} \hspace{0.25in} y=L \hspace{0.15in} face\end{aligned}\end{align}

where L = 1 m and T0 = 1.0 C. The simulation is run until the steady-state temperature distribution develops.

The LaPlace equation governs the steady-state temperature distribution,

${{\partial^{2} T} \over {\partial x^{2}}} + {{\partial^{2} T} \over {\partial y^{2}}}= 0$

The solution is given by,

$T(x,y) = T_{t=0} {x \over L}{y \over L}$

## The PFLOTRAN Input File (GENERAL Mode)¶

The General Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE GENERAL
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 10 10 1
DXYZ
0.1d0
0.1d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
1.0d0 0.d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0  0.d0  0.d0
1.0d0 0.d0 1.0d0
/
END

MATERIAL_PROPERTY plate
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

CHARACTERISTIC_CURVES cc1
SATURATION_FUNCTION VAN_GENUCHTEN
LIQUID_RESIDUAL_SATURATION 0.5d-1
M 0.75
ALPHA 1.d-3
/
PERMEABILITY_FUNCTION MUALEM
LIQUID_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.1d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END

STRATA
REGION all
MATERIAL plate
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 5.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 100
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_x
HDF5_DATASET_NAME x_line_node_centered
FILENAME ../dataset.h5
END

DATASET temperature_bc_y
HDF5_DATASET_NAME y_line_node_centered
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.0d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 0.0d0 C
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE DATASET temperature_bc_y
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE DATASET temperature_bc_x
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 0.0d0 C
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

END_SUBSURFACE


## The PFLOTRAN Input File (TH Mode)¶

The TH Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE TH
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 10 10 1
DXYZ
0.1d0
0.1d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
1.0d0 0.d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0  0.d0  0.d0
1.0d0 0.d0 1.0d0
/
END

MATERIAL_PROPERTY plate
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
SPECIFIC_HEAT 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END

SATURATION_FUNCTION default
SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
RESIDUAL_SATURATION 0.5d-1
LAMBDA 0.75
ALPHA 1.d-3
END

STRATA
REGION all
MATERIAL plate
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 5.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 100
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_x
HDF5_DATASET_NAME x_line_node_centered
FILENAME ../dataset.h5
END

DATASET temperature_bc_y
HDF5_DATASET_NAME y_line_node_centered
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.0d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 0.0d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_y
PRESSURE 101325 Pa
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_x
PRESSURE 101325 Pa
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 0.0d0 C
PRESSURE 101325 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

END_SUBSURFACE


## The Dataset¶

The hdf5 dataset required to define the initial/boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

filename = 'dataset.h5'
h5file = File(filename,mode='w')

# 1d line in x
# Temperature boundary condition
h5grp = h5file.create_group('x_line_node_centered')
h5grp.attrs['Dimension'] = np.string_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
nx = 2
rarray = np.zeros(nx,'=f8')
for i in range(nx):
rarray[i] = float(i)
h5dset = h5grp.create_dataset('Data', data=rarray)

# 1d line in y
# Temperature boundary condition
h5grp = h5file.create_group('y_line_node_centered')
h5grp.attrs['Dimension'] = np.string_('Y')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
ny = 2
rarray = np.zeros(ny,'=f8')
for i in range(ny):
rarray[i] = float(i)
h5dset = h5grp.create_dataset('Data', data=rarray)


## The Python Script¶

def thermal_steady_2D_BC1stkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.1.3, pg.16
# "A 2D Steady-State Temperature Distribution, Boundary Conditions of 1st Kind"
#
# Author: Jenn Frederick
# Date: 06/28/2016
# ******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
lxyz[0] = 1.      # [m] lx
lxyz[1] = 1.      # [m] ly
lxyz[2] = 1.      # [m] lz
test_pass = False
try_count = 0

# initial discretization values
nxyz[0] = 10      # [m] nx
nxyz[1] = 10      # [m] ny
dxyz = lxyz/nxyz  # [m]
T0 = 1.           # [C]
ierr = 0

while (not test_pass) and (try_count < num_tries):
print_discretization(lxyz,nxyz,dxyz)
nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
try_count = try_count + 1

x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)     # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny)     # [m]
T_pflotran = np.zeros((nx,ny))                     # [C]
T_soln = np.zeros((nx,ny))                         # [C]

# create the analytical solution
i = -1
for x in x_soln:
i = i + 1
j = -1
for y in y_soln:
j = j + 1
T_soln[i,j] = T0*(x/Lx)*(y/Ly)

# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+02 y/Temperature [C]'
ierr = check(T_pflotran)

max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz[0] = nxyz[0]*2.
nxyz[1] = nxyz[1]*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions
'Temperature [C]',"{0:.2f}".format(max_percent_error))

# Plot error analysis
plot_error(error_analysis,dxyz_record,path,try_count)

# Add test result to report card

return;


Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

# 2D Steady Thermal Conduction, BCs of 1st & 2nd Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The Dataset

The Python Script

## The Problem Description¶

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.1.4, pg.17, “A 2D Steady-State Temperature Distribution, Boundary Conditions of 1st and 2nd Kind.”

The domain is a 2x1x1 meter slab extending along the positive x-axis and y-axis and is made up of 20x10x1 hexahedral grid cells with dimensions 0.1x0.1x1 meters. The domain is composed of a single material and is assigned the following properties: thermal conductivity K = 1 W/(m-C); specific heat capacity Cp = 0.001 J/(m-C); density rho = 2,800 kg/m^3.

The temperature is initially uniform at T(t=0) = 1.0 C. The boundary temperatures are:

\begin{align}\begin{aligned}T(2L,y) = {T0 \over L} (2L+2y) \hspace{0.25in} x=2L \hspace{0.15in} face\\T(x,L) = {T0 \over L} (2L+x) \hspace{0.25in} y=L \hspace{0.15in} face\\T(x,0) = {T0 \over L} x \hspace{0.25in} y=0 \hspace{0.15in} face\\{{\partial T} \over {\partial x}}(0,y) = {T0 \over L} \hspace{0.25in} x=0 \hspace{0.15in} face\end{aligned}\end{align}

or

$q(0,y) = -K {T0 \over L} = -1 {W \over {m^2}} \hspace{0.25in} x=0 \hspace{0.15in} face$

where L = 1 m and T0 = 1.0 C. The simulation is run until the steady-state temperature distribution develops.

The LaPlace equation governs the steady-state temperature distribution,

${{\partial^{2} T} \over {\partial x^{2}}} + {{\partial^{2} T} \over {\partial y^{2}}} = 0$

The solution is given by,

$T(x,y) = {T0 \over L} (x+2y)$

## The PFLOTRAN Input File (GENERAL Mode)¶

The General Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE GENERAL
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 20 10 1
DXYZ
0.1d0
0.1d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
2.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
2.0d0 0.d0 0.d0
2.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
2.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0  0.d0  0.d0
2.0d0 0.d0 1.0d0
/
END

MATERIAL_PROPERTY plate
ID 1
POROSITY 1.d-8
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END

CHARACTERISTIC_CURVES cc1
SATURATION_FUNCTION VAN_GENUCHTEN
LIQUID_RESIDUAL_SATURATION 0.5d-1
M 0.75
ALPHA 1.d-3
/
PERMEABILITY_FUNCTION MUALEM
LIQUID_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.1d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END

STRATA
REGION all
MATERIAL plate
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 5.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 100
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_north
HDF5_DATASET_NAME x_line_node_centered_north
FILENAME ../dataset.h5
END

DATASET temperature_bc_south
HDF5_DATASET_NAME x_line_node_centered_south
FILENAME ../dataset.h5
END

DATASET temperature_bc_east
HDF5_DATASET_NAME y_line_node_centered_east
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.0d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END

FLOW_CONDITION west_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX -1.0d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_east
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_north
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_south
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

END_SUBSURFACE


## The PFLOTRAN Input File (TH Mode)¶

The TH Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE TH
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 20 10 1
DXYZ
0.1d0
0.1d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
2.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
2.0d0 0.d0 0.d0
2.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
2.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0  0.d0  0.d0
2.0d0 0.d0 1.0d0
/
END

MATERIAL_PROPERTY plate
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
SPECIFIC_HEAT 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END

SATURATION_FUNCTION default
SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
RESIDUAL_SATURATION 0.5d-1
LAMBDA 0.75
ALPHA 1.d-3
END

STRATA
REGION all
MATERIAL plate
END

TIME
FINAL_TIME 100 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 5.d0 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 100
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_north
HDF5_DATASET_NAME x_line_node_centered_north
FILENAME ../dataset.h5
END

DATASET temperature_bc_south
HDF5_DATASET_NAME x_line_node_centered_south
FILENAME ../dataset.h5
END

DATASET temperature_bc_east
HDF5_DATASET_NAME y_line_node_centered_east
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.0d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION west_face
TYPE
ENERGY_FLUX neumann
PRESSURE dirichlet
/
ENERGY_FLUX -1.0d0 W/m^2
PRESSURE 101325 Pa
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_east
PRESSURE 101325 Pa
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_north
PRESSURE 101325 Pa
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_south
PRESSURE 101325 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

END_SUBSURFACE


## The Dataset¶

The hdf5 dataset required to define the initial/boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

filename = 'dataset.h5'
h5file = File(filename,mode='w')

L = 1.

# 1d line in x
# Temperature boundary condition; SOUTH
h5grp = h5file.create_group('x_line_node_centered_south')
h5grp.attrs['Dimension'] = np.string_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [2.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
nx = 2
rarray = np.zeros(nx,'=f8')
for i in range(nx):
# T = x
x = (float(i)*2*L)/(nx-1)
rarray[i] = x
h5dset = h5grp.create_dataset('Data', data=rarray)

# 1d line in x
# Temperature boundary condition; NORTH
h5grp = h5file.create_group('x_line_node_centered_north')
h5grp.attrs['Dimension'] = np.string_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [2.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
nx = 2
rarray = np.zeros(nx,'=f8')
for i in range(nx):
# T = 2 + x
x = (float(i)*2*L)/(nx-1)
rarray[i] = 2 + x
h5dset = h5grp.create_dataset('Data', data=rarray)

# 1d line in y
# Temperature boundary condition; EAST
h5grp = h5file.create_group('y_line_node_centered_east')
h5grp.attrs['Dimension'] = np.string_('Y')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
ny = 2
rarray = np.zeros(ny,'=f8')
for j in range(ny):
# T = 2 + 2*y
y = (float(j)*L)/(ny-1)
rarray[j] = 2 + 2*y
h5dset = h5grp.create_dataset('Data', data=rarray)


## The Python Script¶

def thermal_steady_2D_BC1st2ndkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.1.4, pg.17
# "A 2D Steady-State Temperature Distribution, Boundary Conditions of 1st and
# 2nd Kind"
#
# Author: Jenn Frederick
# Date: 06/29/2016
# ******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
lxyz[0] = 2.      # [m] lx
lxyz[1] = 1.      # [m] ly
lxyz[2] = 1.      # [m] lz
test_pass = False
try_count = 0

# initial discretization values
nxyz[0] = 20      # [m] nx
nxyz[1] = 10      # [m] ny
dxyz = lxyz/nxyz  # [m]
T0 = 1.           # [C]
ierr = 0

while (not test_pass) and (try_count < num_tries):
print_discretization(lxyz,nxyz,dxyz)
nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
try_count = try_count + 1

x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)  # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny)  # [m]
T_pflotran = np.zeros((nx,ny))                  # [C]
T_soln = np.zeros((nx,ny))                      # [C]

# create the analytical solution
i = -1
for x in x_soln:
i = i + 1
j = -1
for y in y_soln:
j = j + 1
T_soln[i,j] = (T0/Ly)*(x+(2*y))

# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+02 y/Temperature [C]'
ierr = check(T_pflotran)

max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz[0] = nxyz[0]*2.
nxyz[1] = nxyz[1]*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions
'Temperature [C]',"{0:.2f}".format(max_percent_error))
# Plot error analysis
plot_error(error_analysis,dxyz_record,path,try_count)

# Add test result to report card

return;


Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.

# 3D Steady Thermal Conduction, BCs of 1st Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The PFLOTRAN Input File (TH Mode)

The Dataset

The Python Script

## The Problem Description¶

This problem is adapted from Kolditz, et al. (2015), Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions, Springer International Publishing, Switzerland. Section 2.1.5, pg.18, “A 3D Steady-State Temperature Distribution, Boundary Conditions of 1st Kind.”

The domain is a 1x1x1 meter cube extending along the positive x-axis, y-axis, and z-axis and is made up of 10x10x10 cubic grid cells with dimensions 0.1x0.1x0.1 meters. The domain material is assigned the following properties: thermal conductivity K = 1.0 W/(m-C); specific heat capacity Cp = 0.001 J/(m-C); density rho = 2,800 kg/m^3.

The temperature is initially uniform at T(t=0) = 1.0 C. The boundary temperatures are:

\begin{align}\begin{aligned}T(0,y,z) = T0 \left( {{0}+{y \over L}+{z \over L}} \right) \hspace{0.25in} x=0 \hspace{0.15in} face\\T(x,0,z) = T0 \left( {{x \over L}+{0}+{z \over L}} \right) \hspace{0.25in} y=0 \hspace{0.15in} face\\ T(x,y,0) = T0 \left( {{x \over L}+{y \over L}+{0}} \right) \hspace{0.25in} z=0 \hspace{0.15in} face\\ T(L,y,z) = T0 \left( {{L}+{y \over L}+{z \over L}} \right) \hspace{0.25in} x=L \hspace{0.15in} face\\ T(x,L,z) = T0 \left( {{x \over L}+{L}+{z \over L}} \right) \hspace{0.25in} y=L \hspace{0.15in} face\\ T(x,y,L) = T0 \left( {{x \over L}+{y \over L}+{L}} \right) \hspace{0.25in} z=L \hspace{0.15in} face\end{aligned}\end{align}

where L = 1 m and T0 = 1.0 C. The simulation is run until the steady-state temperature distribution develops.

The LaPlace equation governs the steady-state temperature distribution,

${{\partial^{2} T} \over {\partial x^{2}}} + {{\partial^{2} T} \over {\partial y^{2}}} + {{\partial^{2} T} \over {\partial z^{2}}} = 0$

The solution is given by,

$T(x,y,z) = T_{t=0} \left( {x \over L}+{y \over L}+{z \over L} \right)$

## The PFLOTRAN Input File (GENERAL Mode)¶

The General Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE GENERAL
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 8 8 8
DXYZ
0.125d0
0.125d0
0.125d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
1.0d0 0.d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0 0.d0 0.d0
1.0d0 0.d0 1.0d0
/
END

REGION bottom_face
FACE BOTTOM
COORDINATES
0.d0 0.d0 0.d0
1.0d0 1.0d0 0.d0
/
END

REGION top_face
FACE TOP
COORDINATES
0.d0 0.d0 1.0d0
1.0d0 1.0d0 1.0d0
/
END

MATERIAL_PROPERTY cube
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-20
PERM_Y 1.d-20
PERM_Z 1.d-20
/
END

CHARACTERISTIC_CURVES cc1
SATURATION_FUNCTION VAN_GENUCHTEN
LIQUID_RESIDUAL_SATURATION 0.5d-1
M 0.75
ALPHA 1.d-3
/
PERMEABILITY_FUNCTION MUALEM
LIQUID_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.1d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END

STRATA
REGION all
MATERIAL cube
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 10 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.10 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 10
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_west
HDF5_DATASET_NAME node_centered_surf_west
FILENAME ../dataset.h5
END

DATASET temperature_bc_east
HDF5_DATASET_NAME node_centered_surf_east
FILENAME ../dataset.h5
END

DATASET temperature_bc_north
HDF5_DATASET_NAME node_centered_surf_north
FILENAME ../dataset.h5
END

DATASET temperature_bc_south
HDF5_DATASET_NAME node_centered_surf_south
FILENAME ../dataset.h5
END

DATASET temperature_bc_bottom
HDF5_DATASET_NAME node_centered_surf_bottom
FILENAME ../dataset.h5
END

DATASET temperature_bc_top
HDF5_DATASET_NAME node_centered_surf_top
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 1.0d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_west
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_east
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_north
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_south
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION bottom_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_bottom
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

FLOW_CONDITION top_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_bc_top
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

BOUNDARY_CONDITION bottom_face
REGION bottom_face
FLOW_CONDITION bottom_face
END

BOUNDARY_CONDITION top_face
FLOW_CONDITION top_face
REGION top_face
END

END_SUBSURFACE


## The PFLOTRAN Input File (TH Mode)¶

The TH Mode PFLOTRAN input file can be downloaded here.


SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE TH
/
/
END

SUBSURFACE

GRID
TYPE structured
NXYZ 8 8 8
DXYZ
0.125d0
0.125d0
0.125d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 1.0d0
/
END

REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END

REGION east_face
FACE EAST
COORDINATES
1.0d0 0.d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION north_face
FACE NORTH
COORDINATES
0.d0 1.0d0 0.d0
1.0d0 1.0d0 1.0d0
/
END

REGION south_face
FACE SOUTH
COORDINATES
0.d0 0.d0 0.d0
1.0d0 0.d0 1.0d0
/
END

REGION bottom_face
FACE BOTTOM
COORDINATES
0.d0 0.d0 0.d0
1.0d0 1.0d0 0.d0
/
END

REGION top_face
FACE TOP
COORDINATES
0.d0 0.d0 1.0d0
1.0d0 1.0d0 1.0d0
/
END

MATERIAL_PROPERTY cube
ID 1
POROSITY 0.01d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
SPECIFIC_HEAT 1.d-3
THERMAL_CONDUCTIVITY_DRY 1 W/m-C
THERMAL_CONDUCTIVITY_WET 1 W/m-C
SATURATION_FUNCTION default
PERMEABILITY
PERM_X 1.d-14
PERM_Y 1.d-14
PERM_Z 1.d-14
/
END

SATURATION_FUNCTION default
SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
RESIDUAL_SATURATION 0.5d-1
LAMBDA 0.75
ALPHA 1.d-3
END

STRATA
REGION all
MATERIAL cube
END

TIME
FINAL_TIME 10 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.10 yr at 0.d0 yr
END

OUTPUT
SNAPSHOT_FILE
TIMES yr 10
NO_PRINT_INITIAL
FORMAT HDF5
/
END

DATASET temperature_bc_west
HDF5_DATASET_NAME node_centered_surf_west
FILENAME ../dataset.h5
END

DATASET temperature_bc_east
HDF5_DATASET_NAME node_centered_surf_east
FILENAME ../dataset.h5
END

DATASET temperature_bc_north
HDF5_DATASET_NAME node_centered_surf_north
FILENAME ../dataset.h5
END

DATASET temperature_bc_south
HDF5_DATASET_NAME node_centered_surf_south
FILENAME ../dataset.h5
END

DATASET temperature_bc_bottom
HDF5_DATASET_NAME node_centered_surf_bottom
FILENAME ../dataset.h5
END

DATASET temperature_bc_top
HDF5_DATASET_NAME node_centered_surf_top
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE 1.0d0 C
PRESSURE 101325 Pa
END

FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_west
PRESSURE 101325 Pa
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_east
PRESSURE 101325 Pa
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_north
PRESSURE 101325 Pa
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_south
PRESSURE 101325 Pa
END

FLOW_CONDITION bottom_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_bottom
PRESSURE 101325 Pa
END

FLOW_CONDITION top_face
TYPE
TEMPERATURE dirichlet
PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_bc_top
PRESSURE 101325 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END

BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END

BOUNDARY_CONDITION north_face
REGION north_face
FLOW_CONDITION north_face
END

BOUNDARY_CONDITION south_face
FLOW_CONDITION south_face
REGION south_face
END

BOUNDARY_CONDITION bottom_face
REGION bottom_face
FLOW_CONDITION bottom_face
END

BOUNDARY_CONDITION top_face
FLOW_CONDITION top_face
REGION top_face
END

END_SUBSURFACE


## The Dataset¶

The hdf5 dataset required to define the initial/boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

filename = 'dataset.h5'
h5file = File(filename,mode='w')

T0 = 1.  # [C]
L = 1.   # [m]

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition x=0 face
h5grp = h5file.create_group('node_centered_surf_west')
h5grp.attrs['Dimension'] = np.string_('YZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
ny = 11
nz = 11
rarray = np.zeros((ny,nz),'=f8')
for j in range(ny):
for k in range(nz):
y = (float(j)*L)/(ny-1)
z = (float(k)*L)/(nz-1)
rarray[j][k] = T0*(0+(y/L)+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition x=L face
h5grp = h5file.create_group('node_centered_surf_east')
h5grp.attrs['Dimension'] = np.string_('YZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
ny = 11
nz = 11
rarray = np.zeros((ny,nz),'=f8')
for j in range(ny):
for k in range(nz):
y = (float(j)*L)/(ny-1)
z = (float(k)*L)/(nz-1)
rarray[j][k] = T0*(1+(y/L)+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition z=0 face
h5grp = h5file.create_group('node_centered_surf_bottom')
h5grp.attrs['Dimension'] = np.string_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
nx = 11
ny = 11
rarray = np.zeros((nx,ny),'=f8')
for i in range(nx):
for j in range(ny):
x = (float(i)*L)/(nx-1)
y = (float(j)*L)/(ny-1)
rarray[i][j] = T0*((x/L)+(y/L)+0)
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition z=L face
h5grp = h5file.create_group('node_centered_surf_top')
h5grp.attrs['Dimension'] = np.string_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
nx = 11
ny = 11
rarray = np.zeros((nx,ny),'=f8')
for i in range(nx):
for j in range(ny):
x = (float(i)*L)/(nx-1)
y = (float(j)*L)/(ny-1)
rarray[i][j] = T0*((x/L)+(y/L)+1)
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition y=0 face
h5grp = h5file.create_group('node_centered_surf_south')
h5grp.attrs['Dimension'] = np.string_('XZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
nx = 11
nz = 11
rarray = np.zeros((nx,nz),'=f8')
for i in range(nx):
for k in range(nz):
x = (float(i)*L)/(nx-1)
z = (float(k)*L)/(nz-1)
rarray[i][k] = T0*((x/L)+0+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature boundary condition y=L face
h5grp = h5file.create_group('node_centered_surf_north')
h5grp.attrs['Dimension'] = np.string_('XZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
nx = 11
nz = 11
rarray = np.zeros((nx,nz),'=f8')
for i in range(nx):
for k in range(nz):
x = (float(i)*L)/(nx-1)
z = (float(k)*L)/(nz-1)
rarray[i][k] = T0*((x/L)+1+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)


## The Python Script¶

def thermal_steady_3D_BC1stkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
# Based on:
# Kolditz, et al. (2015) Thermo-Hydro-Mechanical-Chemical Processes in
# Fractured Porous Media: Modelling and Benchmarking, Closed Form Solutions,
# Springer International Publishing, Switzerland.
# Section 2.1.5, pg.18
# "A 3D Steady-State Temperature Distribution"
#
# Author: Jenn Frederick
# Date: 06/28/2016
# ******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
lxyz[0] = 1.      # [m] lx
lxyz[1] = 1.      # [m] ly
lxyz[2] = 1.      # [m] lz
test_pass = False
try_count = 0

# initial discretization values
nxyz[0] = 8       # [m] nx
nxyz[1] = 8       # [m] ny
nxyz[2] = 8       # [m] nz
dxyz = lxyz/nxyz  # [m]
T0 = 1.     # [C]
ierr = 0

while (not test_pass) and (try_count < num_tries):
print_discretization(lxyz,nxyz,dxyz)
nx = int(nxyz[0]); ny = int(nxyz[1]); nz = int(nxyz[2])
dx = dxyz[0]; dy = dxyz[1]; dz = dxyz[2]
Lx = lxyz[0]; Ly = lxyz[1]; Lz = lxyz[2]
try_count = try_count + 1

T_soln = np.zeros((nx,ny,nz))                   # [C]
T_pflotran = np.zeros((nx,ny,nz))               # [C]
x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx)  # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dx/2.),ny)  # [m]
z_soln = np.linspace(0.+(dz/2.),Lz-(dx/2.),nz)  # [m]

# create the analytical solution
i = -1
for x in x_soln:
i = i + 1
j = -1
for y in y_soln:
j = j + 1
k = -1
for z in z_soln:
k = k + 1
T_soln[i,j,k] = T0*((x/Lx)+(y/Ly)+(z/Lz))

# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+01 y/Temperature [C]'
index_string,remove)
ierr = check(T_pflotran)

max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz[0] = nxyz[0]*2.
nxyz[1] = nxyz[1]*2.
nxyz[2] = nxyz[2]*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions
z_levels = [0,math.floor(nx/3),math.floor(nx/1.5),nx]

Refer to section Python Helper Functions for documentation on the qa_tests_helper module, which defines the helper functions used in the Python script above.