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 domain set-up.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

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

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
    LIQUID_PRESSURE dirichlet
  /
  TEMPERATURE 1.5d0 C
  LIQUID_PRESSURE 101325 Pa
END

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

FLOW_CONDITION right_end
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
  /
  TEMPERATURE 2.d0 C
  LIQUID_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))
  nxyz_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)

    # load data from hdf5
    output_str = '/1D_steady_thermal_BC_1st_kind.h5'
    index_str = 'Time:  1.00000E+02 y/Temperature [C]'
    T_pflotran[:] = read_pflotran_output_1D(path+output_str,index_str,remove)
    ierr = check(T_pflotran)
    
    max_percent_error = calc_relative_error(T_soln[1:nx+1],T_pflotran,ierr)
    record_error(error_analysis,nxyz_record,dxyz_record,max_percent_error,
                 nxyz,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
  plot_1D_steady(path,x_soln,T_soln,x_pflotran,T_pflotran,'Distance [m]',
                 'Temperature [C]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,nxyz_record,dxyz_record,path,try_count,1)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  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 domain set-up.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

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

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
  /
  TEMPERATURE 1.0d0 C
  LIQUID_PRESSURE 101325 Pa
END

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

FLOW_CONDITION right_face
  TYPE
    ENERGY_FLUX neumann
    LIQUID_PRESSURE dirichlet
  /
  ENERGY_FLUX +1.5d0 W/m^2
  LIQUID_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))
  nxyz_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)

    # load data from hdf5
    index_string = 'Time:  1.00000E+02 y/Temperature [C]'
    T_pflotran[:] = read_pflotran_output_1D(path+
                    '/1D_steady_thermal_BC_1st_2nd_kind.h5',index_string,remove)
    ierr = check(T_pflotran)

    max_percent_error = calc_relative_error(T_soln[1:nx+1],T_pflotran,ierr)
    record_error(error_analysis,nxyz_record,dxyz_record,max_percent_error,
                 nxyz,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
  plot_1D_steady(path,x_soln,T_soln,x_pflotran,T_pflotran,'Distance [m]',
                 'Temperature [C]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,nxyz_record,dxyz_record,path,try_count,1)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  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 domain set-up.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

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

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
  /
  TEMPERATURE 1.0d0 C
  LIQUID_PRESSURE 101325 Pa
END

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

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

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

FLOW_CONDITION south_face
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
  /
  TEMPERATURE 0.0d0 C
  LIQUID_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.bytes_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
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.bytes_('Y')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
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))
  nxyz_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)

    # load data from hdf5
    index_string = 'Time:  1.00000E+02 y/Temperature [C]'
    T_pflotran[:,:] = read_pflotran_output_2D(path+
                      '/2D_steady_thermal_BC_1st_kind.h5',index_string,remove)
    ierr = check(T_pflotran)

    max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
    record_error(error_analysis,nxyz_record,dxyz_record,max_percent_error,
                 nxyz,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
  plot_2D_steady(path,x_soln,y_soln,T_soln,T_pflotran,'Distance [m]',
                 'Temperature [C]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,nxyz_record,dxyz_record,path,try_count,2)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  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 domain set-up.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

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

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
  /
  TEMPERATURE 1.0d0 C
  LIQUID_PRESSURE 101325 Pa
END

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

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

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

FLOW_CONDITION south_face
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
  /
  TEMPERATURE DATASET temperature_bc_south
  LIQUID_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.bytes_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [2.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
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.bytes_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [2.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
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.bytes_('Y')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [1.]
# Location of origin
h5grp.attrs['Origin'] = [0.]
# Load the dataset values
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))
  nxyz_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)

    # load data from hdf5
    index_string = 'Time:  1.00000E+02 y/Temperature [C]'
    T_pflotran[:,:] = read_pflotran_output_2D(path+
                      '/2D_steady_thermal_BC_1st_2nd_kind.h5',index_string,remove)
    ierr = check(T_pflotran)

    max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
    record_error(error_analysis,nxyz_record,dxyz_record,max_percent_error,
                 nxyz,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
  plot_2D_steady(path,x_soln,y_soln,T_soln,T_pflotran,'Distance [m]',
                 'Temperature [C]',"{0:.2f}".format(max_percent_error))
  # Plot error analysis
  plot_error(error_analysis,nxyz_record,dxyz_record,path,try_count,2)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  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 domain set-up.

Comparison of the PFLOTRAN vs. analytical solution for GENERAL mode.

Comparison of the PFLOTRAN vs. analytical solution for TH mode.

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
  CHARACTERISTIC_CURVES cc1
  PERMEABILITY
    PERM_X 1.d-14
    PERM_Y 1.d-14
    PERM_Z 1.d-14
  /
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

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
  /
  TEMPERATURE 1.0d0 C
  LIQUID_PRESSURE 101325 Pa
END

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

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

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

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

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

FLOW_CONDITION top_face
  TYPE
    TEMPERATURE dirichlet
    LIQUID_PRESSURE dirichlet
  /
  TEMPERATURE DATASET temperature_bc_top
  LIQUID_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.bytes_('YZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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.bytes_('YZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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.bytes_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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.bytes_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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.bytes_('XZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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.bytes_('XZ')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [0.1,0.1]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
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))
  nxyz_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)
    
    # load data from hdf5
    index_string = 'Time:  1.00000E+01 y/Temperature [C]'
    T_pflotran[:,:,:] = read_pflotran_output_3D(path+
                        '/3D_steady_thermal_BC_1st_kind.h5',
                        index_string,remove)
    ierr = check(T_pflotran)

    max_percent_error = calc_relative_error(T_soln,T_pflotran,ierr)
    record_error(error_analysis,nxyz_record,dxyz_record,max_percent_error,
                 nxyz,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]
  plot_3D_steady(path,x_soln,y_soln,z_levels,T_soln,T_pflotran,'Distance [m]',
                 'Temperature [C]',"{0:.2f}".format(max_percent_error))
  
  # Plot error analysis
  plot_error(error_analysis,nxyz_record,dxyz_record,path,try_count,3)
  
  # Add test result to report card
  add_to_report(path,test_pass,max_percent_error,ierr)

  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.